Residence time distributions of different size particles in the spray zone ofa Wurster fluid bed studied using DEM-CFD
Liang Li1, 2, Johan Remmelgas2, *, Berend G.M. van Wachem3, Christian von Corswant2, Mats Johansson2, Staffan Folestad2, Anders Rasmuson1
1Department of Chemical and Biological Engineering, Chalmers University of Technology, SE-412 96, Gothenburg, Sweden
2Pharmaceutical Development, AstraZeneca R&D, Mölndal, Sweden
3Division of Thermofluids, Department of Mechanical Engineering, Imperial College London, Exhibition Road, London SW7 2AZ, United Kingdom
*Correspondence should be addressed to Johan Remmelgas at
Tel +46-(0)-31-7065838 / Fax +46-(0)31-7763729
.
Abstract
Particle cycle and residence time distributions in different regions, particularly in the spray zone, play an important role in fluid bed coating. In this study, a DEM-CFD (discrete element method, computational fluid dynamics) model is employed to determine particle cycle and residence time distributions in a laboratory-scale Wurster fluid bed coater. The calculations show good agreement with data obtained using the positron emission particle tracking (PEPT) technique. The DEM-CFDsimulations of different sizeparticles show that large particles spend a longer time in the spray zone and in the Wurster tube than small particles. In addition, large particles are found on average to move closer to the spray nozzle than small particles, which implies that the large particles could shield small particles from the spray droplets. Both of these effects suggest that large particles receive a greater amount of coating solution per unit area per cyclethan small particles. However, the simulations in combination with the PEPT experiments show that this is partly compensatedfor by a longer cycle time for large particles. Large particles thus receive more coating per unit area per pass through the spray zone, but also travel through the spray zone less frequently than small particles.
Keywords: fluid bed, coating, spray zone, residence time distribution, discrete element method
1.Introduction
The Wurster process[1] has been widely used to coat particles in industry for a number of different purposes[2].In the pharmaceutical industryit is used to coat tablets or pellets with drug substances or functional films that, e.g. control the release of the drug substance [3-6]. The process typically takes place in a fluid bed, as illustrated in Figure 1,in which particles are kept near minimum fluidization by a fluidization air flow that is supplied through a distributor plate at the bottom of the fluid bed. One or more two-fluid nozzleslocated at the bottom of the bed supply an atomization air flow and a liquid suspension in the form of droplets. The liquid droplets are deposited onto the particlesthat pass through the spray zone. The particles are then dried by the air flow as they move upwards and as they settle back towards the bottomof the fluid bed.
Figure 1Schematic of the Wurster process and the different regions: 1) the spray zone, 2) the Wurster tube, 3) the fountain region, 4) the downbed region, and 5) the horizontal transport region.
The fluid bed can bedivided intodifferent regions as shownin Figure 1[7-9]. The particle coating cycle, i.e.the sequence of coating and subsequent drying,is also illustrated in Figure 1. It is here definedas the trajectory in which the particle travels through the spray zone, the Wurster tube, the fountain regionand the downbed region before it reappears in the spray zone to begin the next coating cycle. The total amount of coating deposited onto a particle depends on the number of times the particle passes the spray zone and on the amount of coating solution that the particle receives per pass[10-13]. Thus it is clear that the cycle time distribution (CTD) and the residence time distribution (RTD) in the spray zone of particles are important factors in determining the quality of the product such as the film thickness and the film thickness variability. Furthermore, the RTDs in different regions determine the drying rate of the particles in different regions, which plays a key role in the film formation process.
There have been a number of experimental studies on CTDs and RTDs in previous work. Mann and Crosby[14] andShelukar et al.[12]measured the CTDs using radioactive or magnetic particles in spouted beds and Wurster beds respectively, and San José et al.[15] investigated the CTDs by monitoring colored particles in spouted bed dryer.While the CTD certainly affects the amount of coating solution deposited onto the particles, the RTD in the spray zone and the detailed motion of particles in the spray zone also play an important role.Moreover, it has been proposed that particles of one size receive a smaller portion of the coating solution due to shielding by particles of a different size. For example, Cheng and Turton[16] reported that the major cause of coating variation can be attributed to particles receiving different amounts of coating when passing through the spray zone. Recently, Li et al.[17]presented an experimental study of particle cycle and residence time distributions in different regions using the positron emission particle tracking (PEPT) technique. The results from the PEPT technique confirmed clear differences in the RTD for different size particles, also for mixtures of particles. Unfortunately, because of limitations with the PEPT technique the particle motion in the spray zone could only partially be characterized. Additional studies weretherefore suggested by Li et al.[17] to provide further insight into the detailed particle motion in the spray zoneand its contribution to particle coating process.
In order to shed light on the detailed particle motion in the spray zone, DEM-CFD (discrete element method, computational fluid dynamics) modeling may be employed. In this modeling approach, every particle is modeled (as a discrete element) while the gas phase is treated as a continuum. When the number of particles is large, as is typical in particle coating processes, the computational time can become very long. For many particle coating systems used in the pharmaceutical industry, however, batch sizes of a few hundred grams are often employed during development. Depending on the particle size, these systems thus contain between tens of thousands and a few million particles. Simulations involving tens of thousands of particles are fairly straightforward, while simulations involving a few million particles are challenging but still feasible (see e.g. [18]). In the pharmaceutical industry, CFD-DEM modeling thus offers a unique opportunity to study process systems that are commonly used for small-scale production.
The discrete element method was proposed for granular assemblies by Cundall and Strack[19] and was first introduced into research in fluid beds by Tsuji et al. [20].The DEM-CFD approach is increasingly applied in modeling of particle/powder processes [21-29] and has become popular [30] because interparticle interactions are accounted for in a straightforward manner. For example, Link et al. [22] assessed the capability of DEM to reproduce several important flow regimes observed in a 3D spout-fluid bed and van Buijtenenet al. [31, 32] investigated the effect of elevating the spout on the dynamics of a (multiple) spout-fluidized bed.In the work by Frieset al.[26], an example of the RTD in the spray zone of a Wurster coater was reported by means of DEM-CFDsimulations. In addition, Yanget al. [29, 33, 34] presented the hydrodynamics in a 3D spouted bed and found that the solid residence time is shortest in the spout region.
The aim of the present study is to establish and validate a DEM-CFD model for simulations of particle motion in a Wurster fluid bed. This is performed as an important step in developing a simulation tool for the entire particle coating process. By following particle trajectories, the CTDs and RTDs in different regions of the Wurster fluid bed process are determined and compared to data from recent PEPT experiments [17]. In these recent experiments, it was pointed out that it is difficult to detect the tracer particle in the spray region towards the bottom of the bed because the -rays tend to be absorbed by the bulky metal parts in this region. Since it was not possible to employ the PEPT technique to measure the particle residence time in the spray zone, special attention is in the present study paid to evaluate the predicted particle motion in the spray zone. The validated DEM-CFD model is thus employed to simulate particle motion in the spray zone to gain an improved understanding of the coating process. These latter simulations are used to predict the residence time in the spray zone and to examine whether particles of one size may preferentially shield particles of a different size from the spray droplets. Lastly, the results are discussed in the context of coating and the effect on the overall coating process.
2.Mathematical model
2.1.Equations of particle motion
The motion of an individual particle is described using Newton’s second law
/ (1)where is the mass of the particle,is the particle velocity, is the gas velocity,is the gravitationalacceleration,is the particle volume,is the gas pressure gradient,is the particle volume fraction, is the contact force during particle-particle and particle-wall collisions and is the interphase momentum transfer coefficient.
The particle volume fraction is calculatedusing
/ (2)where is the volume of the computational cell. The coupling between the discrete and continuous phase is handled based upon a multigrid technique [35], as described by e.g.Bruchmülleret al. [36]. The particles are tracked on a so-called particle mesh, which is a Cartesian mesh with the size larger than the particle. All the coupling terms, including the volume fraction and the drag forces, are determined on the length-scale of this particle mesh. In case the local fluid cell is smaller than a particle, the coupling between that fluid cell and a particle occupying that fluid cell will be length-scale weighted and thus remain physical. The procedure has been described and tested by Mallouppas and van Wachem [37], and has been further worked out in the recent work of Capecelatro and Desjardins [38].
The interphase momentum transfer coefficient is based on the Ergun equation [39]in the dense regime and onthe Wen and Yu drag correlation in the dilute regime [40],
/ (3)where is the gas volume fraction, is the dynamic viscosity of air, and is the particle diameter. In equation (3), is the drag coefficient, which is written as
/ (4)where is the particle Reynolds number
/ (5)In equation (5) is the air density.
The angular momentum of the particle is calculated by
/ (6)where is the moment of inertia of the particle, is the rotational velocity and is the total torque acting on the particle.
2.2.Soft sphere model
In general, particle-particle and particle-wall collisions can be modeled using the hard sphere or soft sphere model. In the hard sphere model, the interaction forces are assumed to be impulsive and all other forces are negligible during collisions [41]. In the soft sphere model, the deformation of particles in contact is modeled by relating the local linear deformation in the normal and tangential directions to a force in these directions, using Hertz-Mindlin theory [42]. The hard sphere model is easier to use but applicable only to binary collisions. The soft sphere model directly implements the physics of collisions but is computationally costly. In this study, since there are dense regions where many particles can be in contact for a long time, the soft sphere model is used.
In the soft sphere model, deformation of particles is replaced with overlap of two particles and the energy loss during collisions is taken into account through a spring-dashpot system. This system is characterized using the spring stiffness, , the damping coefficient, , and the friction coefficient, . The former two quantities are calculated using Hertz contact theory [42, 43], as explained below.
The normal and tangential contact forces acting on the particle, and , are given by
/ (7)/ (8)
where and are the normal and tangential displacements of the particle, and are unit vectors normal and tangential to the contact plane, respectively, is the relative velocity between the and particles, and is the slip velocity of the contact point.
The normal and tangential stiffnesses, and , are expressed by
/ (9)/ (10)
where is the Young’s modulus of the particle, is the shear modulus, is the Poisson’s ratio and is the radius. The suffixes and denote the components in the normal and tangential directions, while the suffixes and denote the particle and the particle respectively.
According to Tsuji et al. [20], the damping coefficients and represent the visco-elastic dissipation of kinetic energy in the normal and tangential directions, respectively. These are given as
/ (11)/ (12)
where denotes a constant related to the coefficient of restitution [44], and represents the effective particle mass and is calculated by
/ (13)More detailed description of the model can be found in literature [19, 20, 30, 41].
2.3.Equations of gas flow
A characteristic time scale of the flow may be calculated based on the atomizer flow rate and the nozzle diameter. The Reynolds number thus calculated is approximately 104 and single-phase flow of air may therefore be expected to be turbulent.However, since the Stokes number of the particles is quite large, turbulence is not expected to have a direct influence on the particle trajectories. In addition, for dense gas-solid flows, such as the one in this study, the particle stress is expected to be much greater than the stress due to turbulence [45]. A turbulence model is therefore not used in the DEM-CFD model. The continuity and momentum equations for gas flow are then writtenas follows
/ (14)/ (15)
where is the viscousstress tensor for incompressible flow,
/ (16)3.Experimental
3.1.The Wurster bed configuration
The model of the fluid bed used in this study is based on the geometry of the Strea-1 fluid bed. The fluid bed is 380 mm high, with top and bottom diameters of 250 mm and 114 mm respectively. Detailed dimensions are given in the experimental paper [17].
The atomization air flow was supplied through a nozzle with a diameter of 5 mm located 3 mm above the bottom of the fluid bed. This nozzle was only a big circular inlet with no opening for any liquid suspensions. Even in this case, the nozzle air velocity is high, e.g. 50 m/s. This nozzle was used in the experiments in anticipation of future computational work so as to avoid the numerical difficulties associated with simulating sonic air flow while at the same time maintaining a high velocity jet.
The fluidization air flow passes through the bowl-shaped distributor, which consists of a central base and an outer annular region. While the central base of the distributor is fully open, the spray nozzle at the center implies that the region is also annular. The outer annulus region of the distributor has a number of orifices with a diameter of 3 mm, as illustrated in Figure 2. A wire mesh screen was put over the distributor to prevent particles from falling down below the distributor.
Figure 2A sketch of the bowl-shaped distributor: 1) nozzle, 2) solid, 3) fully-opened central base, and 4) outer annulus
3.2.PEPT data
In this work, calculations using the DEM-CFD model are compared to recent PEPT experimental data[17]; no additional experiments are performed. In the PEPT measurement system, a radioisotope is incorporated into a tracer particle. As the tracer particle moves around in the vessel the radioisotope decays in beta-decay and releases back-to-back γ-rays. These γ-rays are detected via two large position sensitive detectors. After a sufficient number of γ-rays are collected the location of the tracer particle can be found through three dimensional triangulation. More details about the technique and the algorithms used to determine the location have been reported by Parker et al. [46-48].
3.3.Material
As a common material in the pharmaceutical industry, microcrystalline cellulose (MCC) spheres were employedin the PEPT experiments, as well as a model particulate material in the DEM-CFD simulations. The pellets that typically are used in practice have a diameterthat ranges from 200 µm to 1 mm. However, smaller or larger particles, such as tablets, are also coated using the Wurster process. In order to limit the number of particles, the particle size of 1749 and 2665 µm, which is still relevant and applicable to product development, was selected. In the PEPT experiments, the particle size was measured using a QICPIC Particle Size Analysis (Sympatec GmbH), and the sphericity ranged from 0.85 to 0.95 (see [17] for a photograph of samples).
4.Numerical
4.1.Meshing
To solve the equations of particle motion and gas flow, a fully coupled multi-phase flow solver, MultiFlow ( is employed. The model of thefluid bed is the same as in the PEPT experiments. Figure 3shows the global mesh and close-ups of the surface mesh at the bowl-shaped distributor and nozzle.Approximately 65000 hexahedral cells are required based on a grid independence check via simulations of single phase gas flow. The computational mesh near the distributor and inside the Wurster tube is locally refinedin order to properlycapture the flow characteristics in this region.
Figure 3(a) The surface mesh of the fluid bed and (b) the close-up of the bowl-shaped distributor (from top view)
4.2.Boundary and initialconditions
In the recent PEPT experiments,the flow ratesof the atomization and fluidization air were measured separately using rotameters[17]. These measurements provide the global flow rates, but they do not give any information about the flow distribution at the distributor.In order to specify the flow distribution at the distributor, a single-phase CFD model including the air supply chamber, the distributor, the wire mesh screen, and the fluid bed was developed in Ansys Fluent. In this single-phase CFD model, the distributor was modeled with all its orifices, but the wire mesh screen was included as a region with a specified pressure loss (see e.g. [49]). The results obtained using this single-phase CFD model (not shown) show that theair flow passing through the central base of the distributor is in between 53% and 64% of the total fluidization air flow.