7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Simulation of detailed flow behaviour in a Wurstertype fluidized bed based upon
Discrete Element and TwoFluid Modelling
Berend van Wachem,* Johan Remmelgast and
Ingela NiklassonBjntm
Department of Mechanical Engineering, Imperial College London, Exhibition Road, London SW7 2AZ
t Pharmaceutical Development, AstraZeneca R&D, MOlndal, Sweden
Keywords: Wurstertype fluidized bed, TwoFluid Modelling, Discrete Element Modelling
Abstract
Discrete element or Lagrangian models have become a very useful and versatile tool to study the hydrodynamical
behaviour of pariculate flows. For example, they are very attractive because phenomena such as coating, drying,
agglomeration, and breakup can be incorporated in a straightforward manner. In these models, Newton's equations
of motion are solved for each particle, and a collision model is applied to handle particle encounters. The motion
of the fluid phase is determined from the volume averaged governing equations in an Eulerian framework. Although
recent increases in computational power has significantly increased the applicability of discrete element models, their
usefulness is still limited to systems involving on the order of one or a few million particles. For commercialscale
production processes involving large numbers of particles, other models such as twofluid models must be employed.
In the twofluid approach, where both phases are treated as volume averaged, fully interpenetrating continue, it is
much more difficult to include phenomena such as agglomeration and breakup.
In this study, fluid mechanics calculations using a discrete element model (DEM) are employed to study the flow
behaviour in a laboratoryscale Wurstertype coating process computationally. The results from the DEM calculations
are compared to predictions based upon simulations using a twofluid model in the context of studying the flow
behavior in and around the jets that emanate from the orifices in the distributor plate. For the overall flow behavior,
it is found that the two different models give predictions that are similar. However, for the detailed flow behavior in
the immediate vicinity of the distributor plate and the atomization nozzle, the results differ. The DEM model predicts
that the air jets near the distributor plate and the atomization nozzle persists further into the particle bed than for the
twofluid model. However, the effect of the jets on the behavior of the particles is greater for the twofluid model.
While the effect of the atomization jet on the timeaveraged particle distribution for the twofluid model, it is hardly
discernible for the DEM model. For both models, however, the effect of the jets near the distributor plate is limited to
the region in the immediately downstream of the distributor plate.
Introduction
The Wurster bed coating process is used frequently in
the pharmaceutical industry. For example, it is em
ployed to deposit a drug substance or a film on pellets
for controlledrelease formulations. The development of
such a process, from the laboratory to the commercial
scale, is difficult, and scaleup is always a challenge.
There have been a number of studies dealing with
scaleup of fluidized beds. For example, Glicksman and
coworkers (Glicksman 1984, 1988) analyzed the equa
tions of motion to develop a set of scaling rules. These
scaling rules were later considered by van Ommen et al.
(2006), who evaluated them computationally and found
that similitude between scales could not be obtained us
ing any of the proposed rules. Despite the extensive re
search in this field, scaleup of fluidized bed processing
remains a challenge. Given the fact that it is challenge
for a fluidized bed, it is obviously also a challenge for
a Wurster bed, i.e. a spouted bed with a draft tube, that
involves a droplet spray.
One reason for the difficulty in scaleup of fluidized
beds is obvious: the particle size does or should not
change upon scaleup while the geometrical configura
tion of the processing equipment does. This is obviously
the case also for Wurster beds, though in this case it can
be expected to be more difficult since the flow profile is
not uniform throughout the cross section. That is, for
Wurster beds the distributor plate has been specifically
designed in order to create or at least enhance the circu
lating behavior that is characteristic of spouted beds.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
A nonuniform flow profile throughout the cross
section is also obtained immediately downstream of the
distributor plate due to the gas jets that emerge from the
orifices in the distributor plate. These jets may penetrate
the bed and result in bypass or attrition. The effect of the
jets in determining the flow behavior obviously depends
on the length of the jets compared to the dimensions of
the vessel. In particular, the influence of the jets may
be expected to be greater for small vessels and they are
thus a complication that must be dealt with in order to
understand how scaleup can be expected to affect pro
cess and, more imlli'i .lm.ll \ product performance.
In this work, the flow behavior in a laboratoryscale
Wurster bed is studied computationally using discrete el
ement and twofluid models in order understand how the
predictions between the two models differ. For the lat
ter model, the solids phase is modeled using a constant
particlephase viscosity. The results are employed to un
derstand how the atomization jet and the gas jets formed
by the the distributor plate affect the performance of the
process.
Computational models
The computational model includes twophase flow of air
and particles in a model of a Wurster bed that includes
an inlet section, a distributor plate including a wire mesh
screen, an atomization nozzle, and a product chamber. A
sketch of the cross section of the (axisymmetric) geom
etry is shown in Figure 1. Figure 1 also shows a close
up of the geometry near the spray nozzle. In reality, the
nozzle consist of an annular air inlet and a circular liquid
inlet, as also indicated in Figure 1.
Following Karlsson et al. (2009) the nozzle is mod
eled as an inlet with a velocity of 40 m/s. In order to
study the influence of the nozzle on the global flow be
havior, calculations are carried out for two nozzle diam
eters: 0.4 mm and 4 mm corresponding to volumetric
flow rates of 0.02 and 2 m3/h at at the relevant con
ditions for the atomization flow, i.e. a temperature of
20 C and a pressure of 1 atm.
In previous work, Karlsson et al. (2009), the model
of the distributor plate was simplified by representing it
as three annular regions, as shown in Figure 2, where
each region is homogeneous. A different velocity was
specified in each region in order to create a nonuniform
flow profile at the bottom of the bed. Specifying the
velocity at the bottom of the bed is a good approximation
if the pressure drop across the distributor plate is large
compared to the pressure drop across the particle bed
because then the flow distribution through the plate will
not be affected by the flow behavior in the bed. As an
alternative to specifying the velocity at the bottom of the
bed, one may include an inlet plenum and model flow
H =60mm
r
L,=320mm
GdF 15 mm
D1 =100 mm
(a) Vessel geometry
Dnoe
(b) Geometry near spray nozzle
Figure 1: Sketches of geometry of Wurster bed and ge
ometry near spray nozzle.
(a)
Reflection plane
(b)
Figure 2: Sketch of the distributor plate.
through the three annular regions of the plate using a
different pressure loss coefficient for each annulus. A
case in which there is such a coupling between the inlet
and the product chambers has, in fact, been considered
by Johansson et al. (2006), though only for a distributor
plate the creates a uniform flow profile throughout its
cross section.
In both the aforementioned cases the presence and be
haviour of the orifices is neglected. While this may be
a reasonable approximation in most cases, it is not ob
vious that the effect of the jets on the behavior of the
bed can be neglected for smallscale beds. In this work,
the detailed geometry of the distributor plate is included.
The distributor plate is modeled as sketched in Figure 2,
which shows the smallest repeatunit (a 30degree slice)
from which the distributor plate may be recreated by ro
tational translations. In Figure 2 it may be noted that the
30degree slice has reflection symmetry with respect to
the indicated plane so that it may be possible to model
the flow using as little as 1/24th of the full geometry.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Air Particles
Boundary U[m/s] a [] U [m/s] a []
Jet 44 1.0 0.0 0.0
Inlet 1.4 1.0 0.0 0.0
Table 1: Inlet boundary conditions.
Condition Air Particles
Wall slip No slip Free slip
Outlet volume fraction 1.0 0.0
Outlet pressure 1 atm
Table 2: Outlet and wall boundary conditions.
FluidPhase Model
In the calculations, air is assumed to have properties
corresponding to a temperature of 60 C and a pressure
of 1 atm (for both the atomizing and fluidizing air (see
Karlsson (2007))). Air is thus modeled as an incom
pressible fluid with a density of 1.1 kg/m and a viscosity
of 2.0 x 10 5 Pas (see Table 3). The volumetric flu
idization flow rate of 40 m3/h at a temperature of 20 C
and a pressure of 1 atm corresponds to a a flow rate
of 40 m3/h at a temperature of 60 "C and atmospheric
pressure, which gives an average velocity of 1.6 m/s for
an inlet with an outer diameter of 100 mm diameter and
an inner diameter of 15 mm. The outlet at the top of the
expansion chamber is modeled as a pressure outlet. The
boundary conditions are summarized in Tables 1 and 2.
The interaction of the air with the particles is included
by accounting for the volume that the particles occupy
and through the drag force that the particles exert on the
fluid, which to a first order approximation is linear in the
relative velocity between the fluid, vf, and the particles,
vp. Details of the derivation of the governing equations
may be found in Anderson and Jackson (1967). The con
tinuity equation becomes
d(afpf) dOf(aPfv)
at 9xi
Condition
Density [kg/m3]
Air
1.1
Particles
1.5 x 103
Viscosity [Pa s] 2.0 x 10 0.10
Maximum packing [] 1.0 0.65
Diameter [pm]
N/A
480
Table 3: Properties of air and particles.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
and the momentum equation becomes
O(afpfv) 09 afPf v3 ) 09oxa 3)
at p Oxp Ox(T
Sf OxJ
a _X
phases J f
+TfV3f+S3+ 1 1(f,P) :
where af is the local volume fraction of the fluid phase;
pf is the density of the fluid; p is th pressure; Tf the
source term linear in the velocity field; S} the additional
source terms; and T7 the stress tensor of the fluid, given
by
T + 7) + A f W (3)
a \x3 0x^ ) \ 3 ) ax"
where pf is the shear viscosity and A1 represents the
bulk viscosity of the fluid, which is assumed to be zero.
In equation (2), ,3(,p) is the momentum transfer coeffi
cient, assumed in this work to be in the form due to Wen
and Yu Wen and Yu (1966)
S3 CD apaf vP 265, (4)
4 dp f
where CD is the drag coefficient for an individual parti
cle (assumed to be spherical)
S24(1 + 0.15Re 687) Re < 1000
S0.44, Re > 1000 (5)
The equations are discretized by the finite volume
method and the subsequent set of linearized equations
is solved for the whole domain, in conjunction with the
appropriate equations for the particle phase.
EulerianEulerian ParticlePhase Model
In the EulerianEulerian twophase flow model, the par
ticle and fluid phases are assumed to be continuous and
fully interpenetrating Newtonianlike fluids; the mate
rial parameters are summarized in Table 3. The parti
cle phase is thus modelled by equations similar to Equa
tions 1 and 2. However, in the particle phase momen
tum equation a term involving the gradient of a solids
pressure has also been included to prevent the parti
cle volume fraction from exceeding its maximum value,
ap,max
The particlephase stress is modeled using a constant
particlephase viscosity of p = 0.10 Pas based upon
values reported on in the literature Grace (1970). The
particle phase pressure is modeled using the solids pres
sure model proposed by Gidaspow,
Pp = Goe( ..), (6)
Figure 3: The computational mesh employed for two
phase EulerianEulerian flow predictions.
where Go is the reference modulus of elasticity and c
is the compaction modulus. In this work, a reference
modulus of elasticity of 1 Pa, a compaction modulus of
20, and a maximum particle volume fraction of 0.65 are
employed. The material parameters are summarized in
Table 3. The computational model for twophase flow
does not include flow through the distributor plate. In
stead, inlet boundary conditions are imposed at locations
that correspond to the the orifices in the distributor plate.
Figure 3 shows a typical mesh for this case.
EulerianLagrangian TwoPhase Flow Model
A discrete element model (DEM) proposed by Cundall
and Strack Cundall and Strack (1979) is used to model
the particles. The individual trajectories of the particles
are determined in the Lagrangian framework, where par
ticle collisions are modelled via the softsphere model
proposed by Tsuji et al Tsuji et al. (1991), which ac
counts for some nonreversible deformation.
The trajectories of individual particles are considered
(i.e. described in a Langrangian framework) and New
ton's 2nd law is solved for each particle, accounting for
the fluidparticle, particleparticle, and particlewall in
teractions, and approximating the integral with the Ver
let algorithm.
Newton's 2nd law for the particles is written as
Vp m+ pg
vpVp
N
+ [Fw+ Fpp]
MP = 13(f,p)V v
dv Vp (
m^= ^[v
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
and for the rotational momentum
dwc
I = Tp (8)
where mp is the mass of the particle; Vp is the volume of
the particle; II is the moment of inertia; Tp the torque
acting on the particle; wp is the rotational velocity; vp is
the translational velocity; and 3 is the momentum trans
fer coefficient (see Equation 4). In addition, Fpw and
Fpp are the particlewall and particleparticle interaction
forces, respectively. The right hand side terms of equa
tion (7) are outlined in table 4.
Table 4: Terms of right hand side of equation(7)
Term Force type
/v (vf vp Drag force
mpg Body force due to gravity
VpVPf Force due to the pressure gradient
Fp Particlewall contact force
Fpp Particleparticle contact force
Implementation of particle collisions
The particleparticle and particlewall interactions as
taken into account in this work are assumed to be local;
i.e. no longrange forces are included. If each particle
pair is considered in order to determine this interaction,
the resulting computational effort would scale as N2,
N being the number of particles. In order to develop
a more efficient scheme, a particlemesh algorithm is
adopted in which each of the particles is assigned a cell
in the particle mesh based upon its location. Using this
particlemesh, each particle is only tested for interaction
against particles located in the same or directly neigh
bouring particle mesh cells. Although this method re
quires a more complex algorithm it leads to a scaling of
N logN, making it a lot more favourable for large num
bers of particles. The particle collisions are modelled
by the softsphere model as described by Cundall and
Strack Cundall and Strack (1979). In brief, this model
uses a springdashpotslider arrangement to describe the
particle behaviour before, during and after a collision.
The damping coefficient used, described by Tsuji and
coworkers Tsuji (1994); Tsuji et al. (1991) is, however,
used since it is related to the coefficient of restitution.
The normal and tangential contact forces are given by
the sum of forces due to the spring and dashpot. From
Hertzian contact theory the normal and tangential con
tact forces are,
Fij = (kn,6 .1. G n)n (9)
Ftij = ktet I' Get (10)
where G is the velocity vector of particle i relative to
particle j, Get is the slit velocity vector at the contact
point. The subscripts n and t represent the normal and
tangential components respectively and 6 is the displace
ment, or overlap, of particles i and j during collision.
Three parameters are required by the softsphere model;
stiffness (k), damping coefficient (r) and the coefficient
of friction (p). The parameter p is a well known empiri
cal quantity, stiffness and the damping coefficient can be
estimated using equations
n = a.Mkn (11)
Tpt = aVMkt (12)
where M = ; and m is the mass of the particle.
Note that for wall collisions m,  oc, hence M = p.
The relationship between a and the coefficient of restitu
tion is well defined by Tsuji et al Tsuji et al. (1991). The
spring constants are, based upon elastic deformation,
4 1 a4( 1
56 \ +
( 2 (Ti 2 7j ri + rj 2' 1
kt \ i + G ri (14)
where r is the radius of the particles, a the Poisson's ra
tio; E the Young's modulus, G the shear modulus (given
by G 2(E( ) and 6, the magnitude of the deforma
tion in the normal direction. Note as above, when a par
ticle reaches a wall, r,  oc, hence, =+ 1
Results
Calculations of the EulerianEulerian simulations were
performed with Ansys CFX versions 11.0 and 12.0,
whereas the EulerianLagrangian simulations were per
formed with the inhouse code Multiflow.
Singlephase flow of air
Computations of the singlephase flow of air are em
ployed to determine the details of the flow around the
distributor plate and the atomizer nozzle in order to cal
culate inlet boundary conditions for the multiphase flow
simulations. Figure 4 shows an example of the flow
visualized through the bottom plate. Although Figure
4 shows results for a case with a nozzle diameter of
4.5 mm, the details of the atomization flow do not affect
the flow distribution through the distributor plate. For
example, the same flow distribution is obtained when
the atomization flow is turned off. In any case, the in
let boundary conditions for the cases including the two
phase flow will use the singlephase velocity distribution
Ej rirj
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Velocity
10.0
8.8
7.5
6.2
5.0
3.8
2.5
1.2
0.
[m s^1]
Figure 4: The air velocity
phase flow calculations
as predicted by the single
Orifice Small Medium Big
Velocity 10 9.4 8.0
Table 5: Inlet boundary conditions for twophase flow
calculations.
across the plate, as shown in Table 5. The resulting flow
profile of these boundary conditions is shown in Figure
5 for the EulerianLagrangian simulations.
EulerianEulerian twophase flow predictions
Results are presented for a case for which the atomiza
tion flow enters the vessel through an opening with a
diameter of 0.4 mm and with a velocity of 40 m/s. The
flow in this case is expected to be dominated by the in
fluence from the jets that emanate from the holes in the
distributor plate, though the atomization flow still has
an effect due to its high velocity. Figure 6 shows a snap
shot of the particle distribution at one instant in time that
later may be compared with a similar snapshot for the
DEmodel.
The flow field may also be examined in terms of
the mass flow of particles into and out of the Wurster
tube. Figure 7 shows the instantaneous as well as time
averaged particle mass flow rate into the bottom of the
Wurster tube and out of top of the Wurster tube. Al
though there are considerable fluctuations in the mass
flow rates, the particle mass flow rate into the Wurster
tube at the bottom is clearly greater than the particle
mass flow rate out of the tube at its top. That is, par
Figure 5: The modelled air velocity 0.5 mm above the
distributor plate as modelled in the EulerianLagrangian
simulations.
Figure 6: Snapshot of a particle isovolume correspond
ing to a volume fraction of 0.15 at t 3.79 s.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Figure 7: The particle mass flow rate through the bot
tom of the Wurster tube. The timeaverage mass flow
rates for an averaging period of 2 s are shown with thick
curves.
 : 4s
t =2s
Time Is]
Figure 8: The particle accumulation in the Wurster tube.
tides that enter the Wurster tube through the bottom do
not make it out through the top but fall back down inside
the tube. Since particles enter the Wurster tube but do
not always make it out, there must be an accumulation
of particles inside the tube. Figure 8, which shows the
mass of particles in the Wurster tube is, clearly shows
that particles accumulate in the Wurster tube. The mass
of particles in the Wurster tube may be of interest since
this quantity can possibly be measured using tomogra
phy. Figure 8 also shows the timeaveraged mass of par
ticles in the Wurster tube, calculated using an averaging
period of 2, 3, and 4 s. Figure 8 thus shows that the pro
cess requires at least 6 s to reach a steady state and that
a steady, timeaveraged solution can be calculated with
reasonably precision by forming an average over 2 s.
The fact that particles tend to accumulate in the
Wurster tube may also be seen by examining Figure
Figure 9: Contour plots of the timeaveraged particle
volume fraction in a vertical crosssection through the
centerline of the vessel.
9, which shows the timeaveraged particle volume frac
tion at a vertical plane that passes through the center
line of the vessel. In Figure 9, it can be seen that parti
cles have accumulated in the atomization flow. Figure
10 shows the timeaveraged air velocity for the verti
cal plane shown in Figure 9. The terminal velocity of
a sphere with a density of 1500 kg/m3 and a diameter
of 480 pm can be estimated to 0.5 m/s and Figure 10
indeed shows that the particles in the atomization flow
tend to remain there because there is a balance between
interphase drag and gravity.
The jets that emanate from the distributor plate for
single phase flow of air may be clearly seen in Figure
4. However, Figure 9 shows that when the particles are
present the jets only penetrate a short distance into the
bed, at least on average. Whether this is problem or not
obviously depends on the number or mass of particles in
the vessel. If the amount of particles is small, the jets
can possible penetrate the entire bed which can result in
a fair amount of gas bypass. It must, however, be noted
that Figure 9 shows the timeaveraged particle volume
fraction. At any point in time, bypass may be much
more severe that indicated by Figure 9, as shown in Fig
ure 11, which shows a snapshot of the particle volume
fraction at = 6.9 s.
In Figure 9, it may be seen that the particle volume
fraction is slightly higher next to the wall in the inside
of the Wurster tube. One mechanism by which parti
cles end up near the inside of the Wurster tube may be
seen by examining the sequence of images in Figure 12,
which shows vector plots of the particle velocity super
imposed on contour plots of the particle volume fraction.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
3 ~~:I
(a) t 3.67 s
A
Figure 10: Contour plots of the timeaveraged air veloc
ity fraction in a vertical crosssection through the center
line of the vessel.
Figure 11: Contour plot of the particle volume fraction
at t = 2.1 s for mesh B.
(c) t= 3.75 s (d) t= 3.79 s
Figure 12: Vector plots of the particle velocity superim
posed on contour plots of the particle volume fraction at
t = 3.67, 3.71, 3.75, and 3.79 s.
EulerianLagrangian twophase flow predictions
In the case of the EulerianLagrangian simulations, no
symmetry plane was assumed and the computational
domain considered is the complete geometry. The mesh
of the geometry consists of approximately 400,000 cells
and it is shown in Figure 13. As in the simulations
based upon the EulerianEulerian model, the diameter
of the atomizing nozzle is 0.4 mm and the fluid velocity
through this nozzle is set to 40 m/s.
To have the same number of particles in the Wurster
bed compared to the twofluid simulation case, approxi
mately 1.00 x 106 particles are inserted into the domain.
The same inlet velocity boundary conditions are applied
as in the twofluid model simulation case, see Figure 5.
The simulations were run for a period of 9 s of physical
time, of which the first 4 s are discarded in order to
compute timeaverages. The average fluid velocity
obtained in the EulerianLagrangian simulations is
shown in Figure 14. The fluid velocity is comparable
to the air velocity obtained in the twofluid simulation.
Although the air velocity through the central nozzle
is very high, the mass flow of air through this nozzle
is very small. Two snapshots of the locations of the
particles are shown in Figures 15 and 16 showing the
complete Wurster bed and half of the Wurster bed
respectively. The particles in both figures are coloured
by their velocity. The timeaveraged volume fraction
(b) t 3.71 s
JAM
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Figure 15: A snapshot after t 4 s of the locations
of all the particles in the Wurster bed. The particles are
coloured by their velocity.
Figure 13: Part of the surface mesh employed for two
phase EulerianLagrangian flow predictions.
VelocityPhase1 average
34.52126
110
V h o
1
0.1
0.01
0.005015
Figure 14: A snapshot after t 4 s of the velocity mag
nitude of the air in the center of the Wurster bed. The
snapshot is taken from the twophase flow calculations.
Figure 16: A snapshot after t 4 s of the locations of
all the particles in one side of the Wurster bed, showing
the location and velocities of the particles in the center
region. Particles are coloured by their velocity.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
VFrac _average
1
0.975
0.95
0.925
0.9
Figure 17: The averaged fluid volume fraction shown in
a plane going through the center of the Wurster bed
one side of the Wurster bed, showing the location and
shown in a cross section of the Wurster bed is shown
in Figure 17. The figure is fairly symmetrical, which
implies that the statistics of the linear average is satisfac
tory. This figure shows the particles are moved towards
the Wurster tube in the center of the bed. In Figure 18
the particle mass flow entering the Wurster tube, leaving
the Wurster tube and the difference between these two
fluxes is shown, analogously as in Figure 7 and 8 for the
EulerianEulerian simulations.
Discussion
The comparison between the DE and twofluid model
simulation shows that the two models give approxi
mately the same overall flow behavior, at least on aver
age. For example, both models predict an accumulation
of particle inside the Wurster tube. However, the DE
model predicts that the air jets from the distributor plate
and the atomizer penetrate further into the particle bed
than does the twofluid model. Despite this difference,
the twofluid model predicts that the effect on the jets
on the particles is greater. These differences most likely
depend on the choice of constitutive model for the parti
cle phase, for example through the value of the viscosity
employed for the particle phase.
The fact that particles that enter the Wurster tube
through the bottom but do not exit through the top may
possibly create a problems. For example, if the parti
cle concentration inside the Wurster tube increases there
is an increased risk for agglomeration. Having particles
that enter the Wurster tube but do not make it out may
also create a wider particle size distribution that causes
a corresponding increase in the variability of the final
drug product, such as its content of drug substance or
dissolution time.
In the context of optimizing process parameters, it
Figure 18: The particle mass flow rate in the bottom
and the outlet of the Wurster tube as predicted by the
EulerianLagrangian simulations. Also, the difference
between these is shown; representing the local accumu
lation of particle mass.
is thus of interest to calculate a measure of this back
mixing. If particles that enter through the bottom of the
Wurster tube also make it out through the top, albeit at a
later time, there should be a correlation between the par
ticle mass flow rate into the bottom of the Wurster tube
at some point in time t and the particle mass flow rate
out of the top of the Wurster tube at some later point in
time t + At. It may thus be of interest to calculate the
correlation function
CBT(t) lim 
Too T t Jo
TmB(s) Th(t + s)ds
(15)
where Trbottom is the mass flow rate through the bot
tom of the Wurster tube and rbtop is the mass flow rate
through the top of the Wurster tube. For example, Fig
ure 19, which shows this correlation function normal
ized by the variance of the flow through the bottom of
the Wurster tube for the two cases considered in the
EulerianEulerian simulations, shows that an increase
in the atomization flow decreases the amount of back
mixing inside the draft tube.
Conclusions
Fluid mechanics calculations using a discrete element
model (DEM) were employed to study the flow be
haviour in a laboratoryscale Wurstertype coating pro
cess computationally. The results from calculations us
ing a DEmodel were compared to predictions based
upon simulations using a twofluid model in the context
of the flow behavior in and around the jets that emanate
D=0.4 mm
SD=4.5 mm
0. 10
005 
0.00
*0 200 400 6
Time [ms]
Figure 19: The correlation function for the mass flow of
particles through the bottom and top of the Wurster tube.
from the orifices in the distributor plate. For the overall
flow behavior, it was found that the two different models
gave predictions similar predictions. However, for the
detailed flow behavior in the immediate vicinity of the
distributor plate and the atomization nozzle, the results
differ. It was found that the air jets near the distributor
plate and the atomization nozzle persisted further into
the particle bed for the DEmodel than for the twofluid
model. However, further away from the jet inlets, the ef
fect of the jets on the behavior of the particles was found
to be greater for the twofluid model than for the DE
model. This difference may depend on the choice of
constitute model for the particle phase.
References
T.B. Anderson and R. Jackson. A fluid mechanical de
scription of fluidized beds. Ind. Engng. Chem. Fundam.,
6(4):527539, 1967.
P.A. Cundall and O.D.L. Strack. A discrete numerical
model for granular assemblies. Geotechnique, 29:47
65, 1979.
Leon R. Glicksman. Scaling relationships for fluidized
beds. Chemical Engineering Science, 39:13731379,
1984.
Leon R. Glicksman. Scaling relationships for fluidized
beds. Chemical Engineering Science, 43:14191421,
1988.
J.R. Grace. The viscosity of fluidized beds. The Cana
dian Journal of Chemical Engineering, 48:3033, 1970.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
K. Johansson, B.G.M. van Wachem, and A.E. Almstedt.
Experimental validation of CFD models for fluidized
beds: Influence of particle stress models, gas phase com
pressibility and air inflow models. Chemical Engineer
ing Science, 61:17051717, 2006.
Stina Karlsson. Particle Coating in a Wurster Type Bed.
PhD thesis, Chalmers University of Technology, 2007.
Stina Karlsson, Anders Rasmuson, Berend van Wachem,
and Ingela Niklasson BjOrn. CFD Modeling of the
Wurster Bed Coater. AIChE Journal, 55:25782590,
2009.
Y. Tsuji. Discrete particle simulation of gassolid flows.
0 KONA Powder and Particle, 11:5768, 1994.
Y. Tsuji, T. Tanaka, and T. Ishida. Lagrangian numeri
cal simulation of plug flow of cohesionless particles in a
horizontal pipe. Powder T 1.. .. ,, 71:239250, 1991.
J. Ruud van Ommen, Marnix Teuling, John Nijenhuis,
and Berend G.M. van Wachem. Computational valida
tion of the scaling rules for fluidized beds. Powder Tech
i,. 1... ,, 163:3240, 2006.
C. Y. Wen and Y. H. Yu. Mechanics of fluidization.
Chemical Engineering Progress Symposium Series, 62:
100111, 1966.
