Markov Chain Modeling of Fluidized Bed ... - ACS Publications

ARTICLE pubs.acs.org/IECR

Markov Chain Modeling of Fluidized Bed Granulation Incorporating Simultaneous Aggregation and Breakage Muammer Catak,*,† Kevin Cronin,‡ and Dario Medina-Tellez‡ † ‡

Department of Electronics and Communication Engineering, I_zmir University, Turkey Department of Process and Chemical Engineering, University College Cork, Ireland ABSTRACT: Fluidized bed granulation to achieve particle size enlargement is a significant process in chemical engineering. A key aim of the modeling of such a system is prediction of the evolution of the particle size distribution with time. The traditional tool to achieve this is through the use of population balances. This Article proposes an alternative modeling approach of Markov chain, which is a well-known tool in stochastic theory. Specifically, in this Article, it is employed to model bottom spray Wurster type fluidized bed granulation and a one-dimensional chain where the particle size is the state variable was used. Three important elements of Markov chains are the initial state vector a(t), the transition time step τ, and the transition matrix P. The transition matrix was assembled from consideration of breakage and aggregation dynamics. Experiments were conducted to validate the approach using glass beads as the elementary particles and polyethylene glycol as the liquid binder. From the results, it is shown that the Markov chain method is an efficient and accurate tool to analyze the particle size evolution of a size enlargement process.

1. INTRODUCTION Fluidized bed processing has been used in the pharmaceutical industry for more than 50 years. This technology was originally designed for the rapid drying of pharmaceutical granulations.2 Over the years, fluid bed processing has become an important application for other processes such as granulation, coating, and layering. Fluidized bed granulation (FBG) is a particle size enlargement technique especially prevalent in industries concerned with pharmaceuticals, detergents, fertilizers, and food. Common reasons for granulation are to obtain a narrow band particle size distribution, to improve flow properties, and to avoid cake formation during storage. In FBG, a liquid binder solution is injected into the bed usually by a nozzle. There are three different types of FBG according to the spray method employed; they are top spray, bottom spray, and tangential spray. The purpose of the binder is to produce wet and sticky particle surfaces and hence promote agglomeration when interparticle collisions occur. Whether a collision successfully results in granulation is determined by the granulation mechanism, which in turn depends on particle and bed properties. Ennis and Litster3 presented a general approach for wet granulation that describes the whole granulation process using three subprocesses: (1) wetting and nucleation, (2) consolidation and growth, and (3) attrition and brakeage. Although these three processes provide an overview of granulation, specific granulation operations may require more detailed models. Tan et al.4 deduced that the initial nucleation stage can be very short and may not necessarily require examination from knowledge of the fact “the particle collisions in a dense enough fluidized bed happen in very short time about microseconds” as reported by Goldschmit et al.5 A detailed view of Tan et al.4 on some of the phenomena is displayed as a sequence of events in Figure 1. The approach of this Article on the nature of the granulation process in bottom spray FBG is slightly different from that of r 2011 American Chemical Society

Tan et al.4 Before introducing the approach, the terminology used in this Article will be outlined. The starting point is a population of individual particles, which can have different sizes. They can combine to form multiparticle bodies that we call granules. Assuming the elementary particles cannot break, the number of particles is conserved in the granulation process, although the number of granules is dynamically changing. According to our view, the first process is that of body
droplet collisions, which result in wet and sticky bodies. The second process is body
body collisions, which may produce bigger bodies provided at least one of the bodies is wet and sticky enough to contribute the liquid bridge to hold the bodies together. The third process is the liquid bridge rupture that we call unsuccessful aggregation. The fourth process is binder solidification that is denoted successful aggregation. The fifth and final process is solid bridge rupture, which means breakage of the bodies. These processes are illustrated in Figure 2. Population balance modeling (PBM), which is a technique primarily based on the law of mass conservation and represented by integro-differential equation, is widely used to model the operations of particulate systems. There is an enormous number of papers on the application of PBM on granulation processes in the literature.6
9 Alternatively, discrete element modeling (DEM), which incorporates numerical techniques to compute the motion of discrete particles that interact with neighboring particles, has been used by some researchers to model granulation.10,11 For any kind of model, the validity (i.e., consistency between the model output and the experimental results) and the complexity (ease of understanding, ease of implementation, ease of computation, etc.) are very important issues. A continuous number density function in PBEs makes the PBM method quite powerful in analyzing the Received: December 20, 2010 Accepted: August 17, 2011 Revised: August 15, 2011 Published: August 17, 2011 10811

dx.doi.org/10.1021/ie102513v | Ind. Eng. Chem. Res. 2011, 50, 10811–10823

Industrial & Engineering Chemistry Research

ARTICLE

Figure 1. Sequential events in fluidized bed granulation.

In this Article, Markov theory will be employed to model and analyze particle size enlargement in fluidized bed granulation where aggregation and breakage occur simultaneously.

2. THEORY

Figure 2. Notional rate processes in fluidized bed granulation.

dynamics of a process. However, the structure of these equations is complex due to the intrinsic partial integro-differential equations, and hence analytical solutions may not often be possible except for rather simple cases. DEM is also a very complex technique for modeling, and it may be inadequate particularly when the number of parameters of the model increases. The tracking of a very large number of discrete particles, each undergoing many simultaneous and interacting physical processes, is hugely challenging to capture accurately. In addition to this, current computer technology is not capable of dealing with DEM simulations of production scale units.12 A stochastic model based on Markov chain theory is a very valuable tool because of its simplicity and flexibility.

2.1. Markov Chains. The principle of the Markov process is one of the most applied and well-known mathematical modeling methods in stochastic theory. The name comes from the Russian mathematician Andrei Andreyevich Markov (1856
1922) who released the first work13 on this topic. The basic property of a Markov process is that it has no memory. That is, the past and the future of the process are mutually independent, and only the present can influence the future. 14 Where the process is a finite set of states, then it is referred to as Markov chain. If the total process time is divided into the finite intervals of duration τ, then a representation of the process, which is discrete in time, is obtained. This Markov process is called discrete time Markov chains (DTMC). 15 In this study, the term Markov chains refers to DTMC. To constitute a Markov chain, three important notations are used: a state vector a(t), a transition time step τ, and a transition matrix P. The vector a(t) with components a 1 (t), a 2 (t), ..., a n(t) is called the state probability distribution vector of the system at time t. The transition matrix P has entries p ij, which represents transition probabilities from state i to state j at each time step. If the probability for being in the state j at time t is denoted by a j (t), then the state probability distribution of state i for the next time step, a i(t + 1), is given by the sum of product of all probabilities. This is formulated as follows:

ai ðt þ 1Þ ¼

n

∑ aj ðtÞpji j¼1

ð1Þ

If eq 1 is written for all states, then it is possible to show transitions between the states using a square matrix, which is called a Markov or transition matrix. On this basis, the 10812

dx.doi.org/10.1021/ie102513v |Ind. Eng. Chem. Res. 2011, 50, 10811–10823

Industrial & Engineering Chemistry Research

ARTICLE

Figure 3. Classification of Markov chains according to the transition matrix.

transition matrix P is given by n  n matrix as follows: 16 2 3 p11 p12 ::: p1n 6 6 p21 p22 ::: p2n 7 7 7 P¼6 6 ::: ::: ::: ::: 7 4 5 pn1 pn2 ::: pnn nn

The classification of Markov chains according to the dependency of the transition matrix on the elapsed time and the state distribution vector is shown in Figure 3, which is adapted from Berthiaux and Mizonov.14 The general form of eq 1 for both homogeneous and nonhomogeneous transition matrices at various time steps can be given by the following equations, respectively: aðt þ kτÞ ¼ aðtÞPk

ð2aÞ

aðt þ τÞ ¼ aðtÞPðtÞ

ð2bÞ

A Markov chain is called an ergodic Markov chain if it is possible to go from every state to every other state (not necessarily in a single transition). This implies that the transition matrix is a regular ergodic matrix. Where the granulation process corresponds to an ergodic Markov chain, granule growth continues until the asymptotic state vector point, a(η), is reached, where a(η + 1) = a(η)P. In actuality, the asymptotic limit is a dynamical balance point. In other words, aggregation and breakage processes occur, but there is no net size growth or size reduction. The asymptotic behavior of the process can be calculated using matrix algebra. The limit state probability distribution can be written as: aðηÞ ¼ að0ÞPη

ð3Þ

a(η) is an eigenvector of the transpose matrix of P (i.e., P0 ), and the corresponding eigenvalue of P0 equals 1. With a sufficiently large value of η, the converging state vector can be calculated numerically with an acceptable error using: aðη þ 1Þ ¼ aðηÞP ¼ aðηÞ þ εr

ð4Þ

where εr is the error component. A Markov chain simulation can, in principle, have any number of dimensions where each dimension corresponds with a property of the system under study. For this work, a one-dimensional chain will be analyzed with particle size (in fact equivalent diameter) as the state property. For nonspherical particles, equivalent diameter, d, is determined from  1=3 6V ð5Þ d ¼ π

A key implicit assumption underlying this approach is that instantaneous particle size is the most important particle property in determining future evolution of the particle size. The approach could be extended so that if particle porosity or particle sphericity was a significant determinant of the evolution of a particle size (for instance, through the mechanism of porosity influencing granule strength and hence breakage behavior), a two-dimensional chain could be constructed. Finally, selection of time step is important for the Markov model. Some process systems have a natural periodicity such as rotating mixers (period of rotation), vibratory conveyors (period of excitation) for which the transition time step can be chosen from basic physics. For systems that have no obvious period, the time step τ should be chosen taking into consideration process kinetics and total processing time. Also, the state property under analysis (particle size in this work) must be divided into a sufficient number of intervals, n. Ultimately, a suitable pair of (τ, n) values should be selected to represent the time and particle property in discrete form efficiently. Markov theory has been used in many fields such as astronomy, biology, computer science, communications, forecast, game theory, and radio engineering.14 There also are several applications of Markov chains in chemical engineering area such as (i) particle flow in fluidized beds,17
19 (ii) particle mixing,20
22 and (iii) particle breakage.23,24 However, it can be stated that Markov theory is not well-known in chemical engineering. There are some processes for which Markov chains have almost never been applied such as crystallization, agglomeration, coating, drying, and membrane separation.14 Readers can find more detailed information about Markov chains and their application in the following works of Kemeny and Snell,25 Turner,26 and Tamir.27 2.2. Application of Markov Chains to Fluidized Bed Granulation. The granulation process displays very random behavior because the outcomes of many of the subprocesses are uncertain. Probability density functions are commonly used to represent processes, which have a random nature. In this Article, such functions will be used to model aggregation and breakage as a means to capturing the stochasticity in fluidized bed granulation. Aggregation and breakage are the two important subprocesses to explain the fluidized bed granulation process. The term aggregation is used for size enlargement processes whereby small particles are gathered into larger, relatively permanent masses in which the original particles can still be distinguished.28 The term breakage can be defined as separation of fragments from a whole particle. Because aggregation and breakage are the two main subprocesses in fluidized bed granulation, to achieve a successful model, it is important to take into account both processes. Also, because particle size is the Markovian state property, the sensitivity of aggregation and breakage to particle size must be determined. For this work, all of the aggregation and breakage events are binary, meaning that aggregation results from the simultaneous collision of at most two particles, and breakage only gives two daughter particles. 2.2.1. Aggregation Modeling. Sastry29 suggested that the aggregation kernel can be divided into two parts, one depending on the particle size, and the other determined by all other process and material properties: βðε, v
εÞ ¼ β0 gðε, v
εÞ

ð6Þ

β0 is the size-independent part of the kernel, and it usually represents all of the process parameters. g is a function of size 10813

dx.doi.org/10.1021/ie102513v |Ind. Eng. Chem. Res. 2011, 50, 10811–10823

Industrial & Engineering Chemistry Research

ARTICLE

Figure 4. Graphical representation of eq 7.

such as v, ε. An aggregation kernel based on the equi-partition of kinetic energy (known as the EKE model), which was derived by Hounslow,30 is used in this Article. The EKE kernel assumes that the interparticle collisions occur as a result of the random component of the particle velocities. Furthermore, the kinetic energy associated with the random motion of the particles is shared equally across the particle size range. The validity of the EKE kernel depends on whether particle velocity conforms to the characteristics predicted by kinetic theory. Specifically, Hounslow states that the EKE kernel is correct for collision statistics if the underlying (random) velocity distribution in each direction is Gaussian and the velocities are isotropically distributed. This means the overall random component of velocity is described by the Maxwell speed distribution. More detailed information can be obtained from Tan et al.31 and Boerefijn and Hounslow.32 As a check in this work, the velocity of five particles in the system was examined by image analysis. Particle displacements versus time histories were converted into velocity versus time histories. The systematic or deterministic component of the velocity was filtered out to leave the residual random component. The distribution in velocity was plotted in frequency histogram form, and the normality of the data was tested with the software Minitab 14 (Minitab Inc., U.S.) using the Kolmogorov
Smirnov test with 0.01 as the level of significance. Three of the five particles had a Gaussian distributed velocity, while the remaining two just lay outside the confidence limits; this was regarded as solely reflecting the limited amount of sampled data. The following discrete aggregation kernel, which is explained in two functions, the product of a size-independent term β0 and a

term that is a function of the agglomerating particle sizes li and lj, respectively, is given by Tan et al.:4 sffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 1 2 βij ¼ β0 ðli þ lj Þ ð7Þ 3 þ 3 li lj Specifically, li and lj are the representative particle diameters in the size intervals i and j, respectively. A physical interpretation of the dependency on size of eq 7 is that the more the two colliding particles differ in size, the more likely they are to aggregate; conversely, equal sized particles are less likely to agglomerate. A visual implementation of eq 7 can be seen in Figure 4. The size-independent term β0 is the aggregation frequency or the rate at which particles aggregate. It, in turn, will be the product of the size-independent collision rate, the geometric success factor, and the physical success factor.33 The collision rate will primarily reflect the average collision velocity, which itself depends on the magnitude of the variance in the random component of particle velocity (as measured by the granular temperature of the system). The geometric success factor is a measure of the fraction of the total surface area of the particles in the system that is wetted by the binder liquid, which is a prerequisite for agglomeration. Surface area wetness will be a function of the binder liquid addition rate less the drying rate. The physical success factor will depend on the Stokes numbers (for Type I and Type II coalescence) and their levels with respect to the corresponding critical Stokes numbers. According to Type I coalescence, elastic properties of the particles are dominant. The physics of particle
particle contact are as follows: Two 10814

dx.doi.org/10.1021/ie102513v |Ind. Eng. Chem. Res. 2011, 50, 10811–10823

Industrial & Engineering Chemistry Research

ARTICLE

colliding particles have certain kinetic energies. It is assumed that these particles are coated with a well-distributed layer of liquid binder film. A consequent granule formation occurs if the initial kinetic energy of the particles has been dissipated; otherwise, particles will rebound.34 Type II coalescence models assume that the plastic properties are dominant over the elastic properties of the particles. Therefore, direct particle
particle contact happens at least for a finite time period as a consequence of collisions. During the contact, a bond whose strength is dependent on initial amount of plastic deformation and contact time develops. Coalescence happens if this bond is strong enough to avoid being broken by subsequent collisions and agitation forces.35 While generalizations can be difficult, in terms of direct operational variables, β0 should have a positive dependency on binder liquid addition rate, binder liquid viscosity, and fluidizing air turbulence intensity and an inverse dependency on air temperature as this last factor promotes drying and decreases effective viscosity. Equation 7 can be written to obtain a discrete probability distribution as: sffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 1 2 ðli þ lj Þ 3 þ 3 li lj sffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð8Þ βij ¼ β0 n 1 1 2 ðli þ lj Þ þ 3 l3i lj j¼1

observed. It was seen that the smaller was the granule, the longer was the time between breakage events, meaning smaller granules have a lower breakage frequency as indicated by eq 10. The constant coefficient k in eq 11 will depend on the operational variables such as air flow rate, binder viscosity, temperature, etc. For instance, a higher binder viscosity, meaning stronger intragranule bonds, would give a lower value for the coefficient. With further analysis, a specific functional relationship could be established between k and the input operational variables. The fragment particle size distribution should be chosen bearing in mind process parameters. Berthiaux et al.24 used the following formulas for particle breakage in which the particle can transit only into neighboring intervals as: bjm ¼

∑

ð9Þ 2.2.2. Breakage Modeling. There are two breakage functions that constitute the breakage kernel: the breakage frequency and the fragment particle size distribution. The breakage frequency gives the rate particles of size x break per unit time. In the literature, it is common to assume that the breakage frequency is proportional to the volume of the particle as:36,37 SðxÞ ¼ kx

ð10Þ

where k is an arbitrary constant and x is the particle volume. Accordingly, the discrete breakage frequency can be written as: Si ¼ kxi  kl3i

ð11Þ

where xi and li are representative volume size and diameter of interval i, respectively. The breakage frequency function of eq 11 implies that the larger is the granule, the weaker it is and so it is more likely to break. This assumption conforms to the deformation Stokes number criterion, which states that larger particles have a greater deformation Stokes number and hence are more likely to break in impact.38 Furthermore, this assumption was supported by the image analysis of the experiments of intergranule collisions in the fluidized bed. Ten individual granules were arbitrarily selected and were tracked over a time period. For each granule, its behavior in successive impacts with other granules was

ð12Þ

If the restriction on only permitting passage to adjacent states is removed, then the distribution function in eq 12 can be generalized as: bij ¼

∑

The probability that particles in the ith interval then may transit to interval j is written as: 8 sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi > 1 1 2 > > ðli þ lj
i Þ > 3 þ 3 > l l > i j
i > < sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi if i e n 2n gij ¼ 1 1 2 > ðli þ lj
i Þ þ 3 > > l3i lj
i > j¼i þ 1 > > > :0 else

xRm xRj
1 þ xRm

¼

xRj xR1 þ xR2 þ ::: þ xRj þ ::: þ xRi
1 xRj i
1

∑

k¼1

for all i

ð13Þ

xRk

and eq 13 will be subsequently used. The physical interpretation of eq 13 is that a granule can break and produce fragment particles in any size class lower than it, and the closer the size class is to the parent size class, the greater is the proportion of fragment particles that will be found in that class. Granule tracking using image analysis was also used to assist in the selection of the fragment size distribution. Despite some attrition that was observed, the main breakage type was (binary) fragmentation where a granule could break into any two sizes less than it. This breakage mechanism may be present because the granules are quite porous and weak and show no sign of consolidating with time (as measured by the evolution of average granule porosity). Also, as time progresses, granule sphericity continually decreases, which is more suggestive of fragmentation rather than attrition; if attrition were the dominant mechanism, spherical granules could be expected to remain approximately spherical. The value of R is selected as an adjustment parameter of the model. Berthiaux et al.24 reported that in applications to model a real process that R should lie in the range 1 < R < 2. As the magnitude of the parameter R is increased, the greater proportion of fragment particles will lie in size classes closer to the parent size class. The effect of changes to R is illustrated in Figure 5 for a parent granule of size 10 units. This figure shows the fragment size distribution by volume fraction, which is the proportion of fragment particles (for this example, from size 1 unit to size 9 unit) resulting from the breakup of the parent particle. Increasing R accentuates the fragment particle distribution and biases it toward the size classes closer to the parent class, which is more indicative of attrition. For the system under study here, the R value will be taken as the default level of 1. 2.2.3. Formulation of the Markovian Transition Matrix, P. Assuming the initial size range is divided into n state spaces, which are enough to represent the process, the discrete size domain is obtained. If the maximum particle size that can aggregate in the system is represented by the nth interval, then the number of the intervals will be 2n using a uniform discretization procedure. 10815

dx.doi.org/10.1021/ie102513v |Ind. Eng. Chem. Res. 2011, 50, 10811–10823

Industrial & Engineering Chemistry Research

ARTICLE

Figure 5. Illustration of the physical sense of the parameter, R.

On this basis, breakage and aggregation functions must be converted into a discrete form. A particle in any interval can go to a higher interval by aggregating, or can go down to a lower interval as a result of a breakage, or it can stay in the same interval after a transition time step. The transition matrix can thus be divided into three parts regarding the modeling algorithm: • Diagonal entries contain the probability that the particle can stay in the same interval after a transition time step. If βi and Si are the aggregation and breakage frequencies of the ith interval, then the probability that the particle stays in the same interval is: pii ¼ ð1
βi
Si Þ

ð14Þ

• Upper diagonal entries represent the aggregation process. The probability that the particle will go to size interval j from size interval i can be written as: pij ¼ βi gij if j > i

ð15Þ

• Entries below the diagonal comprise the breakage process. The probability that the particle goes to interval j from the interval i due to a breakage can be written as pij ¼ Si bij if j < i

ð16Þ

The transition matrix Pnn combines the entries that are described in eqs 14
16. 8 > < pii ¼ ð1
βi
Si Þ P ¼ pij ¼ βi gij if j > i ðaggregationÞ > : pij ¼ Si bij if j < i ðbreakageÞ

ð17Þ

2.2.4. Markov Chain Homogeneity and Linearity. The issues of chain homogeneity and linearity are central to correct implementation of the Markov approach. In a homogeneous Markov chain, the entries of the transition matrix (giving the probability a granule has of transiting from one size class to another) have no dependence on elapsed time (i.e., number of

time steps, k) and solely depend on the current state (size) of the granule. In this sense, the process history that any particular granule has experienced does not affect the transition probability. As an example, consider two particles each of an arbitrary size 50. One may have been generated by the aggregation of a size 40 and size 10 particle, while the other by the aggregation of a size 30 and size 20 particle. From the Markovian perspective, this history is irrelevant, and both particles have the same future possibility of aggregating or breaking. Considering the process to be homogeneous imposes certain limits on the phenomena that can be represented by the model. For example, the liquid bridges of the binder that hold a granule together tend to solidify with time, meaning the granule becomes more difficult to break. Hence, the probability of transiting to a lower size (by breakage) depends on the elapsed time since the formation of that granule. There are a number of possible techniques to handle such effects. If the characteristic time of such a subprocess is very small as compared to the time step chosen for the simulation, the effect may not be important. Alternatively, instead of using a single transition matrix to simulate the entire process, a number of different transition matrices may be used to simulate sections of the total process duration with the introduction of the changed matrices reflecting some fundamental alteration in system dynamics. Finally (with some restrictions), a continuous time Markov chain approach may be used to deal with time effects where transition rates replace transition probabilities. For a linear Markov chain, the entries of the transition matrix do not depend on the current value of the state vector, that is, on the instantaneous distribution in the property that is being analyzed. For a granulation process with granule size as the property being studied, this means the transition probabilities from one size class to another cannot have sensitivity to the fractional number of particles in any size class. Nonlinear effects are not explicitly present in any of the equations used to formulate the Markov chain in this Article. That is not to say nonlinear effects are not present in the system; rather that if present they are not included in the model. One potential nonlinear phenomenon is that the size-independent aggregation constant, β0, will be a function of the fraction of wet particles that are present in the system. This wetness fraction will depend on the absolute number of particles in the system (which is taken into account) and may depend on the proportion of particles in each size class depending on the model of particle wetting that is adopted. Alternatively, if the transition matrix entries depended on the initial state vector, this would clearly be a nonlinear effect. Technically, this is the case at the start of the process because if the initial size distribution is monodispersed (and all in state 1), then given binary collisions, the only allowable transition is to state 2. Once state 2 is populated with particles, the transitions to states 3 (2 +1) and 4 (2 + 2) become possible, and in addition a transition from state 2 back to state 1 is allowed (breakage). Very quickly though, as the particle size distribution evolves, all possible transitions become permissible and the transition matrix has validity. This does mean though that for early times in the simulation, model output is not accurate. Because this Article has the objective of introducing the concept of Markov modeling to granulation, hence a simplified approach is adopted, and the process is assumed linear. For nonlinear effects, the transition matrix would need updating according to the nature of the state vector, similar to the means of accommodating nonhomogeneous effects. 10816

dx.doi.org/10.1021/ie102513v |Ind. Eng. Chem. Res. 2011, 50, 10811–10823

Industrial & Engineering Chemistry Research

ARTICLE

Figure 6. Initial particle size distribution of glass beads.

Figure 8. D10, D50, and D90 evolutions of the granulation process.

Table 1. States Numbers and Corresponding Equivalent Diameters state number

Figure 7. Laboratory scale Wurster fluidized bed granulator: simple schematic diagram (left), active (right).

5

10

15

20

25

30

35

40

equivalent diameter [μm] 208 330 523 829 1315 2084 3303 5244

3. MATERIALS AND METHODS 3.1. Experimental Methods. For this study, the size enlargement process of glass beads is examined. These consist of spherical particles whose mean diameter is 268 μm as illustrated in Figure 6. The average density of the beads is 2500 kg/m3. The binder liquid was an aqueous dissolution of PEG1500 (polyethylene glycol 1500 Da) at 60% w/w, and 10 g was used for a 200 g batch of beads. During the experiments, the overall air flow rate was 0.65 m3/min. The experiments reported here were carried out using the Mini-Airpro fluidized bed granulator (Pro-C-epT, Belgium) with the Wurster configuration. The experimental setup is given in Figure 7. The upper half of the column is denoted the spray zone where most particle and binder droplet interaction occurs. The equipment Camsizer (Retsch Technology, Germany) was used to measure the particle size distribution at every gram of liquid added to the process, that is, approximately every 2.5 min. Particle vertical velocity was determined with dynamic displacement analysis software (ProAnalyst, Xcitex, U.S.) by tracking the particles captured at 500 frame per second (fps), using a high speed camera (X-Motion Pro, AOS Technologies, Switzerland), during their motion from the upper edge of the Wurster column, reaching a maximum height, and falling into the region between the column and the chamber wall. The analysis of these videos provided data about granule breakage, which was accompanied by observations and sphericity measurements of the granules by means of individual images obtained with the equipment PharmaVision 830 (Malvern Instruments, UK).

Figure 9. A typical glass bead PEG granule.

3.2. Determination of Aggregation and Breakage Parameters. The magnitudes of two parameters, the constant or size-

independent aggregation frequency, β0, and the constant or sizeindependent breakage frequency rate, k, must be obtained to run the simulation. While ideally their values should be obtained from physical considerations, practical considerations mean they are better chosen from analysis of the experimental data. Specifically, these constants are found by minimizing the error between predicted and experimental mean granule size versus time. Once the initial estimates of β0 and k that are needed to start the process are reasonable, the system quickly converges to the most accurate values. Using the maximum likelihood approach, which aims to achieve the minimum squared error 10817

dx.doi.org/10.1021/ie102513v |Ind. Eng. Chem. Res. 2011, 50, 10811–10823

Industrial & Engineering Chemistry Research

ARTICLE

between the model and the experimental measurements, the optimum model distribution was obtained as: f ind the optimum values of fβ0 , kg for minf

3.3. Markov Model Implementation. For the implementation of the model, Matlab software package was used for the calculations. The whole size range from 138 to 5500 μm was divided into 40 intervals with a consecutive diameter ratio of 1.1. Next, the representative size of each interval is calculated using geometric means of end of the intervals such as: pffiffiffiffiffiffiffiffiffiffiffi ð18Þ li ¼ di di
1

n

ðEi
Mi Þ2 g ∑ i¼1

Specifically, the magnitudes for the constants β0 and k were found by fitting the model output to the experimentally measured evolution of the D10, D50, and D90 granule diameters versus time. Note D10 is the diameter for which the cumulative distribution function of the particle size dispersion reaches 0.1 or 10%, D50 is the median diameter (50% of granules are larger and 50% smaller), and D90 is the 90th percentile diameter. The experimental histories and fitted curves are given in Figure 8.1 The constant aggregation frequency was calculated as β0 = 5.13  10
3 s
1, while the size-dependent but time invariant breakage frequency rate is found as k = 6.75  10
12 m
3 s
1.

where di and di-
are upper and lower representative size ends of ith interval. To aid understanding, Table 1 gives the equivalent granule diameter for every fifth state. The Markov chain method requires an initial state vector a(0)140, a transition time step τ to be selected, and a transition matrix P(0)4040. The initial vector comes from the discretization of the probability density function for size that is illustrated in Figure 6 as:

að0Þ ¼ ½0:001, 0:003, 0:006, 0:016, 0:042, 0:108, 0:246, 0:317, 0:192, 0:045, 0:013, 0:006, 0:005, 0, 0:::140 Because of the distribution in glass bead diameter, the initial state vector contains entries from size class 1 (144 μm) up to class 13 (435 μm) with a mean bead size of 268 μm in size class 8. Determination of the magnitude of the time step, τ, is more arbitrary. Because the aggregation frequency β0 is time-independent, then the model is relatively insensitive to the time step size (once it is sufficiently small). There are many physical subprocesses that occur sequentially in the overall granulation process such as particle circulation through the spray region, particle surface wetting, surface layer binder liquid drying, interparticle 2 6 6 6 6 6 6 6 6 6 6 6 6 6 P¼6 6 6 6 6 6 6 6 6 6 6 6 4

collisions, etc. Each of these subprocesses has its own unique time constants, and furthermore they vary with particle size. For this work, a number of different values of the time step, τ, in the range of [0.1 s, 60 s] were investigated. In the end, τ was selected as 1 s for the illustration of the process aiming to obtain optimum simulation cost and model accuracy. The entries in the transition matrix are found from the method explained in section 2. It is not possible to show here the 40 by 40 square matrix with good quality, but part of the transition matrix for the first nine states has the appearance of:

9:9810
1

1:6010
6

1:6510
6

1:7210
6

1:8210
6

1:9510
6

2:1110
6

2:3110
6

2:5510
6

:::

2:8310
5

9:9810
1

2:2210
6

2:2610
6

2:3310
6

2:4310
6

2:5710
6

2:7510
6

2:9810
6

:::

1:6010
5

2:1410
5

9:9810
1

3:1310
6

3:1310
6

3:1810
6

3:2810
6

3:4210
6

3:6210
6

9:9810
1

4:4010
6

4:3610
6

4:3710
6

4:4410
6

4:5710
6


5


5


5

2:3610


5

9:9810


1

6:2310


6

6:0910
6

6:0310
6

6:0410
6

2:0910


5

2:7710


5

9:9810


1


6


6

8:3610
6

1:1910 9:8710
6 8:7310
6 8:0110
6 7:5410
6 7:1910
6 :::

1:5910 1:3210
5 1:1710
5 1:0710
5 1:0110
5 9:6210
6 :::

2:1310 1:7710
5 1:5710
5 1:4410
5 1:3510
5 1:2910
5 :::

8:8510

8:5510

1:9110
5

2:5410
5

3:3910
5

9:9810
1

1:2610
5

1:2010
5

1:8010
5

2:3910
5

3:1910
5

4:2610
5

9:9810
1

7:5410
5

1:7210
5

2:2810
5

3:0510
5

4:0710
5

5:4610
5

9:9410
1

The entries of the matrix represent the transition probabilities between the states in the time step of 1 s. For instance, p1,1 shows that 0.998 (99.8%) of the particles in the size class 1 will remain in the same size class, while p1,2 shows only 1.6  10
6 of them will move to size class 2 by aggregation during the transition time step. The overall pattern of the matrix entries reflects the aggregation and breakage submodels. In each row, entries to the left of the diagonal give the breakage probabilities, and those to the right give the aggregation probabilities. Moving left away from the diagonal entry, along a row, the transition probabilities decrease as a result of the selected breakage kernel. Furthermore, as the breakage rate is taken to be proportional to size, the lower is the row, the greater are the entries left of the diagonal. Moving right from the diagonal entry, the transition probabilities increase as aggregations between different sized granules are more favored

:::

:::

:::

:::

:::

:::

3

7 7 7 7 7 ::: 7 7 7 ::: 7 7 7 ::: 7 7 7 ::: 7 7 7 ::: 7 7 ::: 7 7 7 ::: 7 7 7 5

than for similar sized particles. Each row of the matrix is diagonally dominant (in a given time step, the most likely outcome is that a granule remains in the same size class), although the degree of dominance lessens progressively at lower rows, reflecting the fact that the larger is the granule, the more likely it is to either aggregate or break rather than remain of invariant size.

4. RESULTS AND DISCUSSION 4.1. Results. A typical glass bead-PEG granule is displayed in Figure 9 where the magnification is 5. In the figure, the individual glass beads can also be distinguished. The image analysis software provided with the equipment Camsizer gave a granule diameter of 1.1 mm (placing it in size class 23) and 0.73 for sphericity. The amount of individual particles or beads in the granule cannot be 10818

dx.doi.org/10.1021/ie102513v |Ind. Eng. Chem. Res. 2011, 50, 10811–10823

Industrial & Engineering Chemistry Research obtained directly from the image, but considering an equivalent diameter of 1.1 mm, the number of individual particles would be 64.

Figure 10. The predicted mean size and the experimental mean size evolutions of the granulation process.

ARTICLE

On the basis of average results for the whole batch, the bindersolids mass fraction would be expected to be 2% and the granule porosity 1%, showing most of the space between the beads is filled with PEG. As expected in a granulation process, the average particle size increases with time. The comparison between the predicted mean size and the experimental mean size is displayed in Figure 10.1 Average granule size rises from 270 to 650 over the 25 min process duration. Although the predicted mean size is slightly higher than the experimental (in particular, at later stages of the process at 17.5, 20, and 25 min), agreement is very good in general as measured by the Kolmogorov
Smirnov test with 0.05 confidence interval. The particle size distributions predicted from the model and obtained by experiment are shown in Figure 11 at 2.5 min time intervals. Initially, the particle size has unimodal distribution (see Figure 6). As granulation progresses, a bimodal size distribution is formed with a maximum close to the initial particle size and a second maximum that increases from 500 to above 700 μm. These peaks are more clearly visible especially after 10 min of the process. The first peak is time-invariant and corresponds to the individual glass beads. It is present because some of these beads do not partake in agglomeration, because granule breakage preferentially creates single beads, and because the granulation

Figure 11. Particle size distributions obtained from the model and experiments at various times. 10819

dx.doi.org/10.1021/ie102513v |Ind. Eng. Chem. Res. 2011, 50, 10811–10823

Industrial & Engineering Chemistry Research

ARTICLE

Figure 12. Standard deviation in granule size versus time as predicted by the model and found by experiment. Figure 14. The predicted asymptotic distribution in granule size.

One of the advantages of the Markov chain approach is that it permits the statistics of the size distribution (the higher moments in addition to the mean) to be quantified from the state vectors. Standard deviation in granule size versus time as predicted by the model and found by experiment is plotted in Figure 12. The standard deviation in particle size increases at the initial stages of the experiments as the glass beads begin to aggregate with each other and the range in size increases accordingly. However, after approximately 7.5 min, it tends to decrease slowly as the unimodal size distribution begins to form. This implies that the coefficient of variation in granule size reduces over granulation time. The skewness of the distribution is given (both theoretically and as experimentally measured) in Figure 13. The skewness quantifies the symmetry with increasing positive values, meaning a greater right tail. According to the experimental results and model predictions, the skewness falls from an initial high value at the start of the process and stabilizes out. Hence, the asymmetry in particle size distribution lessens over time until it reaches an asymptotic value, probably because significant numbers of small and large particles appear in the system as a result of simultaneous aggregation and breakage. 4.2. Discussion. From Figure 10, it can be seen that the mean size is increasing to an asymptotic limit although the granulation process is terminated before it is reached (to avoid caking). Using the Markov model, the limit state vector was calculated to be:

Figure 13. The skewness of the distribution (both theoretically and as experimentally measured).

process is terminated relatively early. The second peak reflects the evolution of the mean of the population. Figure 11 shows that the model captures the shape of the particle size distribution and matches the experimental results. Qualitatively, it can be noticed that the agreement between the model and the experiment is better for large and small particles than for intermediate sized particles. In addition, statistical analysis of the results was carried out. According to the Kolmogorov
Smirnov test, the agreement between the model and the experimental results is not less than 90% over the process. 2

0:0000 6 6 0:0193 aðηÞ ¼ 6 6 0:0771 4 0:0063

0:0015 0:0251 0:0707 0:0043

0:0020 0:0321 0:0619 0:0029

0:0026 0:0403 0:0519 0:0020

0:0035 0:0493 0:0418 0:0013

Note that it is in fact a row vector, but due to the visual difficulties it is displayed as a four-row matrix. The asymptotic mean and standard deviation were calculated to be 892 and 342 μm after almost 180 min of granulation processing time. The standard deviation is almost unchanged from its value at 25 min, and this implies that the distribution in granule size cannot be diminished with longer processing times. The predicted asymptotic distribution in granule size is shown in Figure 14. The results are fitted by a two-parameter log-normal distribution with 0.05 confidence interval. Magnitudes of shape and location parameters were 0.49 and 6.77, respectively.

0:0046 0:0587 0:0325 0:0009

0:0062 0:0675 0:0245 0:0005

0:0083 0:0747 0:0180 0:0003

0:0111 0:0791 0:0129 0:0002

3 0:0147 7 0:0800 7 7 0:0091 7 5 0:0001

The Markov chain simulation can also be used to follow the size versus time profile of individual particles by Monte Carlo sampling from the transition probabilities. Four particles, all starting in size class 8 (mean size of the glass beads), had their dynamic size evolution simulated and are given in Figure 15. The successive aggregation and breakage events are clearly seen; after each breakage event, the larger of the two daughter particles became the tracked granule. These simulated realizations predict that granule growth is not steady or continuous with time; early on in the process, large granules can be produced but then can 10820

dx.doi.org/10.1021/ie102513v |Ind. Eng. Chem. Res. 2011, 50, 10811–10823

Industrial & Engineering Chemistry Research

ARTICLE

Figure 15. Dynamic size evolution of four particles, all starting in size class 8 (mean size of the glass beads).

viscosity. If the binder viscosity is very high, the granule will have great strength, which means that the breakage of the granules becomes unimportant.39 Such a system could be modeled by aggregation alone. In our study, the binder viscosity was 36 mPa s, which is relatively low. Figure 16 shows the predicted granule size distribution after 20 min as predicted by the model using aggregation alone as compared to aggregation and breakage. For the former, the agreement with the experimental data is poor with the proportion of smaller sized granules being seriously underestimated. Consequently, both breakage and aggregation functions must be simultaneously included in the model for our case. Finally, the sensitivity of the model output to the magnitudes of the aggregation rate constant β0 and breakage rate constant k was investigated. At early times in the process, the model is more sensitive to the β0 than the k value, while for the later stages, the sensitivity to the k value increases. The main reason for this is because the breakage rate, being positively correlated with granule size, increases with time. Figure 16. The predicted granule size distribution after 20 min as predicted by the model using aggregation alone as compared to aggregation and breakage.

quickly break back to smaller bodies. The importance of breakage in determining size evolution is evident. Fluidised bed granulation is a net size enlargement process where combination of aggregation and breakage takes place competitively. For some systems, the aggregation process may be dominant over breakage. For instance, one of the parameters that affects the granule growth mechanism is the liquid binder

5. CONCLUSION In this Article, the granulation process of glass beads using a bottom spray Wurster-type fluidized bed granulator involving simultaneous aggregation and breakage has been examined. Modeling of this type of particulate process requires the use of stochastic theory to capture the inherent uncertainty in the phenomena. Markov chain theory was the tool employed in this Article. According to the results, Markov chain modeling is an efficient tool to model fluidized bed granulation processes. The transition probabilities, which are the components of the basic operator (transition matrix) of the model, are calculated using 10821

dx.doi.org/10.1021/ie102513v |Ind. Eng. Chem. Res. 2011, 50, 10811–10823

Industrial & Engineering Chemistry Research certain probability density functions that are obtained from the literature. For future work, a Markov model that is more based on physical parameters and dynamics of the granulation process will be built. Specifically, in the current model, the size kernel mainly reflects the physics of collisions. The criteria for whether collisions lead to coalescence may also have a size dependency and inherent randomness; this will require a more detailed treatment of the aggregation rate constant.

’ AUTHOR INFORMATION Corresponding Author

*Tel.: +90 232 246 4949-424. E-mail: muammer.catak@gmail. com.

’ ACKNOWLEDGMENT This publication has emanated from research conducted with the financial support of Science Foundation Ireland. ’ NOMENCLATURE P = Markovian transition matrix a(t) = state probability vector τ = transition time step [s] x = particle volume [mm3] S(x) = breakage frequency (breakage selection function) [s
1] β(x,ε) = aggregation kernel [s
1] β0 = aggregation frequency [s
1] g(x,ε) = size-dependent part of aggregation kernel [mm0.5] gi,j = discrete aggregation probability distribution xi = representative volume of interval i [mm3] li = representative diameter of interval i [mm] d = equivalent diameter [mm] bi,j = discrete fragment particle size distribution εr = error component of the model η = time at which the process reaches its asymptotic limit [s] ’ REFERENCES (1) Catak, M.; Bas, N.; Medina, D.; Cronin, K.; Byrne, E. P.; Fitzpatrick, J. J. Markov chain modelling of fluidised bed granulation. Chem. Eng. J. 2010, 64, 403–409. (2) Wurster, D. E. Means for applying coating to tablets or like. JAPhA 1950, 48, 451. (3) Ennis, B. J.; Litster, J. D. Particle Size Enlargement; Perry’s Chemical Engineers’ Handbook: New York, 1997. (4) Tan, H. S.; Salman, A. D.; Hounslow, M. J. Kinetics of fluidised bed melt granulations I: The effect of process variables. Chem. Eng. Sci. 2006, 61, 1585–1601. (5) Goldschmidt, M. J. V.; Beetstra, R.; Kuipers, J. A. M. Hydrodynamic modelling of dense gas-fluidised beds: comparison of the kinetic theory of granular flow with 3D hard-sphere discrete particle simulations. Chem. Eng. Sci. 2002, 2059–2075. (6) Adetayo, A.; Ennis, B. A new approach to modelling granulation processes for simulation and control purposes. Powder Technol. 2000, 108, 202–209. (7) Lui, L.; Litster, J. Population balance modelling of granulation with a physically based coalescence kernel. Chem. Eng. Sci. 2002, 57, 2183–2191. (8) Peglow, M.; Kumar, J.; Heinrich, S.; Warnecke, G.; Tsotsas, E.; M€orl, L.; Wolf, B. A generic population balance model for simultaneous agglomeration and drying in fluidized beds. Chem. Eng. Sci. 2007, 67, 513–532. (9) Poon, J. M. H.; Immanuel, C. D.; Doyle, I.; Litster, J. D. A three dimensional population balance model of granulation with a mechanistic

ARTICLE

representation of the nucleation and aggregation phenomena. Chem. Eng. Sci. 2008, 63, 1315–1329. (10) Muguruma, Y.; Tanaka, T.; Tsuji, Y. Numerical simulation of particulate flow with liquid bridge between particles (simulation of centrifugal tumbling granulator). Powder Technol. 2000, 109, 49–57. (11) Nase, S. T.; Vargas, W. L.; Abatan, A. A.; McCarthy, J. J. Discrete characterization tools for wet granular media. Powder Technol. 2001, 116, 214–223. (12) Salman, A. D.; Hounslow, M. J.; Seville, J. P. K. Handbook of Powder Technology Vol. 11: Granulation; Elsevier: Amsterdam, 2007. (13) Markov, A. Extension of the law of large numbers to dependent events. Bull. Soc. Phys. Mathematics, Kazan, Russia 1906, 15, 155–156. (14) Berthiaux, H.; Mizonov, V. Applications of Markov chains in particulate process engineering: A review. Can. J. Chem. Eng. 2004, 82, 1143–1168. (15) Norris, J. R. Cambridge Series in Statistical and Probabilistic Mathematics; Cambridge University Press: Cambridge, 1997. (16) Farina, L.; Rinaldi, S. Positive Linear Systems: Theory and Applications; Wiley: New York, 2000. (17) Fox, R.; Fan, L. Stochastic Analysis of Axial Solids Mixing in a Fluidised Bed. Proceedings of the 1st World Congress on Particle Technology part 3; Nurnberg, Germany, 1986; pp 581
595. (18) Dehling, H.; Hoffmann, A.; Stuut, H. Stochastic models for transport in a fluidized bed. SIAM J. Appl. Math. 1999, 60, 337–358. (19) Harris, A.; Thorpe, R.; Davidson, J. Stochastic modelling of the particle residence time distribution in circulating fluidised bed risers. Chem. Eng. Sci. 2002, 57, 4779–4796. (20) Aoun-Habbache, M.; Aoun, M.; Berthiaux, H.; Mizonov, V. An experimental method and a Markov chian model to describe axial and radial mixing in a hoop mixer. Powder Technol. 2002, 128, 159–167. (21) Berthiaux, H.; Marikh, K.; Mizonov, V.; Ponomarev, D.; Barantzeva, E. Modeling continuous powder mixing by means of the theory of Markov chains. Part. Sci. Technol. 2004, 22, 379–389. (22) Ponomarev, D.; Mizonov, V.; Gatumel, C.; Berthiaux, H.; Barantzeva, E. Markov chain modelling and experimental investigation of powder-mixing kinetics in static revolving mixers. Chem. Eng. Process. 2009, 48, 828–836. (23) Berthiaux, H. Analysis of grinding processes by Markov chains. Chem. Eng. Sci. 2000, 55, 4117–4127. (24) Berthiaux, H.; Mizonov, V.; Zhukov, V. Application of the theory of Markov chains to model different processes in particle technology. Powder Technol. 2005, 157, 128–137. (25) Kemeny, J. G.; Snell, J. L. Finite Markov Chains; D. Van Nostrand: New York, 1960. (26) Turner, J. C. Modern Applied Mathematics; English University Press Ltd.: London, 1970. (27) Tamir, A. Applications of Markov Chains in Chemical Engineering; Elsevier: Amsterdam, 1998. (28) Snow, R. H.; Allen, T.; Ennis, J. B.; Litster, J. D. Perry’s Chemical Engineers’ Handbook; McGraw-Hill: New York, 1999. (29) Sastry, K. V. S. Similarity size distribution of agglomerates during their growth by coalescence in granulation or green pelletization. Int. J. Miner. Process. 1975, 2, 187–203. (30) Hounslow, M. J. The population balance as a toll for understanding particle rate processes. KONA 1998, 179–193. (31) Tan, H. S.; Goldschmidt, M. J. V.; Boerefijn, R.; Hounslow, M. J.; Salman, A. D.; Kuipers, J. A. M. Building population balance model for fluidized bed melt granulation: lessons from kinetic theory of granular flow. Powder Technol. 2004, 142, 103–109. (32) Boerefijn, R.; Hounslow, M. J. Studies of fluid bed granulation in an industrial R&D context. Chem. Eng. Sci. 2005, 60, 3879–3890. (33) Rajniak, P.; Stepanek, F.; Dhanasekharan, K.; Fan, R.; Mancinelli, C.; Chern, R. T. A combined experimental and computational study of wet granulation in a Wurster fluid bed granulator. Powder Technol. 2009, 189, 190–201. (34) Ennis, B.; Tardos, G.; Pfeffer, R. A microlevel-based characterization of granulation phenomena. Powder Technol. 1991, 65, 257– 272. 10822

dx.doi.org/10.1021/ie102513v |Ind. Eng. Chem. Res. 2011, 50, 10811–10823

Industrial & Engineering Chemistry Research

ARTICLE

(35) Lui, L.; Litster, J.; Iveson, S. M.; Ennis, B. J. Coalescence of deformable granules in wet granulation processes. AIChE J. 2000, 46, 529–539. (36) Braumann, A.; Goodson, M. J.; Kraft, M.; Mort, P. R. Modelling and validation of granulation with heterogeneous binder dispersion and chemical. Chem. Eng. Sci. 2007, 62, 4717–4728. (37) McMillan, J.; Briens, C.; Berruti, F.; Chan, E. High velocity attrition nozzles in fluidized beds. Powder Technol. 2007, 175, 133–141. (38) Tardos, G. I.; Khan, I. M.; Mort, P. R. Critical parameters and limiting conditions in binder granulation of fine powders. Powder Technol. 1997, 94, 245–258. (39) Schæfer, T. Growth mechanisms in melt agglomeration in high shear mixers. Powder Technol. 2001, 117, 68–82.

10823

dx.doi.org/10.1021/ie102513v |Ind. Eng. Chem. Res. 2011, 50, 10811–10823