Article pubs.acs.org/IECR
Discretized n‑Dimensional Population Balance for Agglomeration S. M. H. Hashemi Amrei and Asghar Molaei Dehkordi* Department of Chemical and Petroleum Engineering, Sharif University of Technology, P.O. Box 11155-9465, Tehran, Iran S Supporting Information *
ABSTRACT: In this article, a new method for solving multidimensional population balance (PB) equations for agglomeration is presented. First, discretized 2-dimensional PB (DPB) is derived and by applying probability approach, this method is extended to n-dimensional PB. The discretization of size intervals is in the form of v(j) = 2j‑i v(i) and the volume based PB is used. Because of conservative nature of volume, all internal coordinates are converted to some independent volumes. A simple and clear procedure for solving multidimensional PB is proposed using probability concept. In addition, various types of birth and death mechanisms involved in PB are analyzed separately. Moreover, results of obtained DPB in 2 and 3-dimensional are compared with time variations of zeroth moment of continuous 1D PB and excellent agreement is found. Finally, the developed method is applied for the simulation of granulation processes and the predicted results are presented and discussed.
1. INTRODUCTION Population balance is the continuity equation for particulate systems to describe the changes in the number of particles during processes. Mass balance, thermal energy equation, and momentum balance equations are governing equations in process modeling and; population balance equation (PBE) is required to model systems dealing with particles (i.e., particulate systems). PBE is used in modeling various processes such as polymerization, crystallization, granulation, etc. Particulate nature of a particle can be characterized in terms of various specifications such as size, shape, liquid saturation, porosity, and its age. These specifications are known as internal coordinates while spatial coordinates such as rectangular coordinates, cylindrical coordinates, and spherical coordinates that can be used to characterize the location of a particle in a system are denoted as external coordinates. The most important internal coordinate of a particle includes its size1 and various processes can be modeled and described by this single characterization. However, for granulation processes more coordinates need to be taken into account to model the process well; for example, size, liquid content, and porosity are all important particle’s parameters in the granulation processes. Generally, PBE is in the form of nonlinear integro-partialdifferential equation, which is quite difficult to solve analytically. A useful numerical technique for solving PBE is to discretize the particle size domain into various intervals and supposing that the particle size distribution (PSD) is constant within each interval.2 This technique yields a discretized PBE (DPBE), which can be solved using various numerical schemes. A number of methods for solving 1D PBE have been reported in the literature to date. Hounslow et al.3 derived a discretized 1D PBE with the advantage of being simple to solve and ensuring accurate prediction of total number and the volume of particles. An extension to Hounslow’s discretization has been presented by Litster et al.2 Litster et al.2 upgraded Hounslow’s domain discretization (i.e. v(i + 1)/v(i) = 2) to an adjustable discretization (i.e. v(i + 1)/v(i) = 21/q). Kumar and Ramkrishna4 developed a discretization procedure with more flexible method using small and large intervals in various size © 2013 American Chemical Society
ranges. Nevertheless, there are scarce efficient works dealing with solving multidimensional PBE. Iveson5 in a comprehensive work studied the limitations of 1D population balance models for predicting the wet granulation processes. Verkoeijan et al.6 proposed the use of truly mutually independent particle specifications as the internal coordinates; thus, in a 3D formulation, they proposed the use of volumes of solid, liquid, and gas as the internal variables rather than the particle total volume, binder content, and porosity that are not mutually independent of each other. Hounslow et al.7 applied reduction/ reformulation techniques to transform a 2-D PB into two 1-D PBs by exploiting a correlation between the two internal coordinates. The weak point of this method is correlations required to correlate the internal coordinates together that can be difficult to construct. Henson8 and Henson et al.9 applied Monte Carlo techniques to simulate multidimensional PB. Although, Monte Carlo techniques are successful for the simulation of n-D PB, this method is time-consuming and, moreover, some of situations that have small probability to occur can be lost. Immanuel and Doyle10 developed an efficient solution technique to solve 3-D PB for granulation processes and, later, they developed an extension to their method to solve higher dimensional systems.11 In their method, the PB is first discretized into subpopulations and the PB is formulated for each of these semilumped subpopulations. The major objective of the present work was to propose a new, simple, and clear method for solving multidimensional PB. In this regard, using the concept of probability, a discretized n-D population balance (DPB) for agglomeration was developed. The obtained DPB can accurately predict the time variations of total number of particles compared to those obtained by the analytical solution of zeroth moment of the 1-D continuous PB. This article is organized as follows: in Section 2, we present a method for solving 2-D PB and a Received: Revised: Accepted: Published: 17487
April 23, 2013 October 18, 2013 October 30, 2013 October 30, 2013 dx.doi.org/10.1021/ie401287w | Ind. Eng. Chem. Res. 2013, 52, 17487−17500
Industrial & Engineering Chemistry Research
Article
In Figures 3−7 used for derivations of the discretized 2D PB in the following subsections, primary and secondary cells are shown with solid and shaded patterns, respectively. 2.1. Birth Mechanisms. Four mechanisms can be considered for possible births in the Eulerian cell (i1, i2) where i1 and i2 denote interval counters of v1 and v2, respectively. Possible birth mechanisms can be expressed as follows (see Table 1): • Mechanism 1: Agglomerations that can lead to produce particles in the class (i1, i2) or in smaller classes. • Mechhanism 2: Agglomerations that can lead to produce particles in the class (i1, i2) and it would be possible to produce particles larger than this class with respect to v1. • Mechanism 3: Agglomerations that can lead to produce particles in the class (i1, i2) and it would be possible to produce particles larger than this class with respect to v2. • Mechanism 4: Agglomerations that can lead to produce particles in the class (i1, i2) or in larger classes with respect to both the internal coordinates. In the following subsections, discretized relations for each mechanism are derived. 2.1.1. Mechanism 1. Mechanism 1 can be taken place in five types and each type should be analyzed separately (Figure 2). 2.1.1.1. Mechanism 1.1 ((i1 − 1,i2 −1) + (1 ≤ j1 ≤ i1 − 2,1 ≤ j2 ≤ i2 − 2)). Consider a particle with volumes a1 < v1 < a1 + da1 and a2 < v2 < a2 + da2 in class (j1, j2). This particle could be aggregated with a particle with volumes 2i1 − a1 < v1 < 2i1 and 2i2 − a2 < v2 < 2i2 from class (i1 − 1, i2 − 1) to produce a new particle in the class (i1, i2). Therefore, the rate of successful agglomeration can be expressed as
mathematical derivation has been developed to obtain a DPB; in Section 3, we discuss the concept of using probability in the derivation of DPB; in Section 4, the developed method for 2-D PB is extended to n-D PB using probability approach; in Section 5, a case study of spray granulation is simulated and the model predictions are presented and discussed; and finally, the summary of the present work and the most important conclusions and future works are presented in Sections 6 and 7, respectively.
2. DISCRETIZED 2D PB In the present method, it is assumed that all the internal coordinates dealing with the system can be transformed into some independent volumes or masses.6 For example; size, liquid content, and the porosity of a particle can be replaced by solid volume, liquid volume, and gas volume. Note that because of the conservative nature of mass and volume, this assumption can be so useful in the formulations. Another important concept of this assumption is that the volumes should be independent and this means that in the agglomeration of two particles, their internal coordinates (i.e., volumes) can be added separately to yield the corresponding volumes of the produced particle. Now consider a 2D PB with internal coordinates v1 and v2, which are both volume and independent. Discretization in terms of domains of the internal coordinates is just like that Hounslow et al.3 used in their work for 1-D PB (i.e., v(j) = 2j−iv(i)) and numbering of the intervals is started with 1. Moreover, v(1) shows the minimum value of the corresponding internal coordinate over the total particles of the system. However, this discretization can be easily extended to adjustable discretization (i.e., v(j) = 2(j−i)/q v(i)). The discretization in 2D space is illustrated in Figure 1. As it is shown in this figure,
dR i1.1 = β × N (i1 − 1, i2 − 1) × N (j1 , j2 ) 1− 1, i 2 − 1, j , j 1 2
×
da1 da 2 a1 a 2 2 j1 2 j2 2i1− 1 2i2 − 1
(1)
Now, to evaluate the rate of successful agglomeration between classes (j1, j2) and (i1 − 1, i2 − 1), eq 1 should be integrated over the entire class (j1, j2) R i1.1 = β × N (i1 − 1, i2 − 1) × N (j1 , j2 ) 1− 1, i 2 − 1, j , j 1 2
× (32 × 2 j1 + j2 − i1− i2)
(2)
Finally, taking the summation over all classes (j1, j2) of mechanism 1.1, the rate expression of this mechanism can be obtained as i1− 2 i 2 − 2
Figure 1. 2-dimensional discretization of the size domains.
R i1.1 = 1, i2
∑ ∑ β × N(i1 − 1, i2 − 1) × N(j1 , j2 ) j1 = 1 j2 = 1
dividing internal coordinate values by 1/2 × v(1), results in dimensionless internal coordinate values. Accordingly, interval “i” contains particles with dimensionless volumes ranging from 2i to 2i+1 (v(i) to v(i + 1)). First, we introduce some definitions that may be useful to provide a framework for better understanding of the formulation. • Secondary cells: In the derivation procedure of the present work, cells that the elements are taken in them and, also, integrations are performed over their boundaries are termed as “secondary” classes. • Primary cells: cells considered for seeking a particle in them for a successful agglomeration with particles in the secondary classes are termed as “primary” classes.
× (32 × 2 j1 + j2 − i1− i2)
(3)
2.1.1.2. Mechanism 1.2 ((i1 − 1, i2 − 1) + (1 ≤ j1 ≤ i1 − 2, j2 = i2 − 1)). If a particle with volumes of a1 < v1 < a1 + da1 and a2 < v2 < a2 + da2 in class (j1,j2), agglomerates with a particle of volumes 2i1 − a1 < v1 < 2i1 and any v2 of class (i1 − 1, i2 − 1), the resulting particle would be classified in the class (i1, i2). Therefore, the rate of successful agglomeration can be evaluated by dR i1.2 = β × N (i1 − 1, i2 − 1) × N (j1 , i2 − 1) 1− 1, i 2 − 1, j , j 1 2
× 17488
da1 da 2 a1 2 j1 2 j2 2i1− 1
(4)
dx.doi.org/10.1021/ie401287w | Ind. Eng. Chem. Res. 2013, 52, 17487−17500
Industrial & Engineering Chemistry Research
Article
Table 1. Birth and Death Mechanisms of Agglomeration for 2D PB mechanisms
patterns
Birth Mechanism 1
•1.1 •1.2 •1.3 •1.4 •1.5
(i1 (i1 (i1 (i1 (i1
− − − − −
1, 1, 1, 1, 1,
i2 − 1) + (1 ≤ j1 ≤ i1 − 2, 1 ≤ j2 ≤ i2 − 2) i2 − 1) + (1 ≤ j1 ≤ i1 − 2, j2 = i2 − 1) i2 − 1) + (j1 = i1 − 1, 1 ≤ j2 ≤ i2 − 2) i2 − 1) + (j1 = i1 − 1, j2 = i2 − 1) 1 ≤ j2 ≤ i2 − 2) + (1 ≤ j1 ≤ i1 − 2, j2 = i2 − 1)
Mechanism 2
•2.1 (i1, i2 − 1) + (1 ≤ j1 ≤ i1 − 1, 1 ≤ j2 ≤ i2 − 2) •2.2 (i1, i2 − 1) + (1 ≤ j1 ≤ i1 − 1, j2 = i2 − 1) •2.3 (i1, 1 ≤ j2 ≤ i2 − 2) + (1 ≤ j1 ≤ i1 − 1, i2 − 1)
Mechanism 3
•3.1 (i1 − 1, i2) + (1 ≤ j1 ≤ i1 − 2, 1 ≤ j2 ≤ i2 − 1) •3.2 (i1 − 1, i2) + (j1 = i1 − 1, 1 ≤ j2 ≤ i2 − 1) •3.3 (1 ≤ j1 ≤ i1 − 2, i2) + (i1 − 1, 1 ≤ j2 ≤ i2 − 1)
Mechanism 4
•4.1 (i1, 1 ≤ j2 ≤ i2 − 1) + (1 ≤ j1 ≤ i1 − 1, i2)
Death Mechanism 5
•5.1 (i1, i2) + (1 ≤ j1 ≤ i1 − 1, 1 ≤ j2 ≤ i2 − 1) •6.1 (i1, i2) + (i1 ≤ j1 < ∞, i2 ≤ j2 < ∞) •6.2 (i1, i2) + (i1 ≤ j1 < ∞, 1 ≤ j2 ≤ i2 − 1) •6.3 (i1, i2) + (1 ≤ j1 ≤ i1 − 1, i2 ≤ j2 < ∞)
Mechanism 6
where the coefficient 1/2 is introduced to avoid counting the rate twice. 2.1.1.5. Mechanism 1.5 ((i1 − 1,1 ≤ j2 ≤ i2 − 2) + (1 ≤ j1 ≤ i1 − 2, i2 −1)). As it may be obvious from Figure 2, mechanism 1.5 can be occurred in two situations that are equal and there is no need to consider both of them. Now consider the agglomeration of a particle with volumes a1 < v1 < a1 + da1 and a2 < v2 < a2 + da2 from class (j1, i2 −1) with particles from the class (i1 − 1, j2); in this case, depending on the range of distance 2i2−1 to 2i2 that a2 is located, three different situations for the successful agglomeration can be considered as follows: 1. Mechanism 1.5.1: If 2i2−1 < a2 < 2i2 − 2j2+1, then a successful agglomeration cannot be occurred and that is because the largest v2 of the class (i1 − 1, j2) is 2j2+1. 2. Mechanism 1.5.2: If 2i2 − 2j2+1 < a2 < 2i2 − 2j2, then this particle should be aggregated with particles with volumes 2i1 − a1 < v1 < 2i1 and 2i2 − a2 < v2 < 2i2+1 from class (i1 − 1, j2) to produce a particle in the class (i1, i2). 3. Mechanism 1.5.3: If 2i2 − 2j2 < a2 < 2i2, with respect to v2 all the agglomerations of this particle with the particles classified in class (i1 − 1, j2) would be successful to produce a new particle in the class (i1, i2). But about v1, this particle should be aggregated with a particle of volume 2i1 − a1 < v1 < 2i1 from class (i1 − 1, j2) to lead a successful birth. Therefore, considering three mechanisms mentioned above, the rate expression is given by
Then, integration over class (j1, j2) and taking the summation over the classes (j1, j2) of the mechanism 1.2, the rate expression of this mechanism can be expressed as i1− 2
R i1.2 = 1, i2
∑ β × N(i1 − 1, i2 − 1) × N(j1 , i2 − 1) j1 = 1
× (3 × 2 j1 − i1)
(5)
2.1.1.3. Mechanism 1.3 ((i1 − 1, i2 − 1) + (j1 = i1 − 1,1 ≤ j2 ≤ i2 − 2)). This mechanism is just the same as the mechanism 1.2 and the rate expression is as follows: i2 − 2
R i1.3 = 1, i2
∑ β × N(i1 − 1, i2 − 1) × N(i1 − 1, j2 ) j2 = 1
× (3 × 2 j2 − i2)
(6)
2.1.1.4. Mechanism 1.4 ((i1 − 1, i2 − 1) + (j1 = i1 − 1,j2 = i2 − 1)). Consider a particle of volumes a1 < v1 < a1 + da1 and a2 < v2 < a2 + da2 from class (i1 − 1, i2 − 1). Agglomeration of this particle with all particles in class (i1 − 1, i2 − 1) results in producing a new particle in the class (i1, i2) dR i1.4 = β × (N (i1 − 1, i2 − 1))2 × 1− 1, i 2 − 1, j , j 1 2
da1 da 2 2i1− 1 2i2 − 1 (7)
Then, integrating over the class (i1 − 1, i2 − 1) we get 1 R i1.4 = β × (N (i1 − 1, i2 − 1))2 1, i2 2
(8)
⎧ 0 2i2 − 1 < a 2 < 2i2 − 2 j2 + 1 ⎪ ⎪ ⎛ da da a 2 j2 + 1 − 2i2 + a ⎞ 2 1 2 ⎪⎜ 1 ⎟ 2i2 − 2 j2 + 1 < a 2 < 2i2 − 2 j2 1.5 dR i1− 1, j , j , i2 − 1 = β × N (j1 , i2 − 1) × N (i1 − 1, j2 ) × ⎨ ⎝ 2 j1 2i2 − 1 2i1− 1 2 j2 ⎠ 2 1 ⎪ ⎪ ⎛ da1 da 2 a1 ⎞ ⎜ j i −1 i −1 ⎟ 2i2 − 2 j2 < a 2 < 2i2 ⎪ ⎝ 2 1 2 2 21 ⎠ ⎩
17489
(9)
dx.doi.org/10.1021/ie401287w | Ind. Eng. Chem. Res. 2013, 52, 17487−17500
Industrial & Engineering Chemistry Research
Article
Figure 2. Patterns of the mechanism 1 (Birth). dR i2.1 = β × N (i1, i2 − 1) × N (j1 , j2 ) 1, i 2 − 1, j , j 1 2
Integrating over the corresponding ranges we get R i1.5 1− 1, j2 , j1 , i 2 − 1
×
= β × N (j1 , i2 − 1) × N (i1 − 1, j2 ) 2
× (3 × 2
j1 + j2 − i1− i2
)
(10)
i1− 1 i2 − 2
R i2.1 = 1, i2
i1− 2 i2 − 2
∑ ∑ β × N(j1 , i2 − 1) × N(i1 − 1, j2 ) × (3 × 2
∑ ∑ β × N(i1, i2 − 1) × N(j1 , j2 ) j1 = 1 j2 = 1
j1 = 1 j2 = 1 2
× (3 × 2 j2 − i2 − 32 × 2 j1 + j2 − i1− i2 − 1) j1 + j2 − i1− i2
)
(12)
Integrating over class (j1, j2) and taking the summation over all classes (j1, j2), the rate equation of mechanism 2.1 is
Finally, taking the summation, the rate expression of the mechanism 1.5 can be expressed as R i1.5 = 1, i2
da1 da 2 2i1+ 1 − 2i1 − a1 a 2 2 j1 2 j2 2i1 2i2 − 1
(11)
(13)
2.1.2.2. Mechanism 2.2 ((i1, i2 − 1) + (1 ≤ j1 ≤ i1 − 1,j2 = i2 −1). Consider a particle with volumes a1 < v1 < a1 + da1 and a2 < v2 < a2 + da2 from class (j1, j2). In terms of v2, agglomeration of this particle with all particles present in class (i1, i2 − 1) leads to successful birth in the class (i1, i2). But with v1, the mentioned particle should be aggregated to a particle with volume 2i1 < v1 < 2i1+1 − a1 to produce a particle located in the class (i1, i2). Therefore, the birth rate in the class (i1, i2) can be evaluated by
2.1.2. Mechanism 2. Mechanism 2 can be occurred in three types and each should be analyzed separately (Figure 3). 2.1.2.1. Mechanism 2.1 ((i1, i2 − 1) + (1 ≤ j1 ≤ i1 − 1,1 ≤ j2 ≤ i2 − 2)). A particle with volumes a1 < v1 < a1 + da1 and a2 < v2 < a2 + da2 from class (j1, j2) should be aggregated with a particle of volumes 2i1 < v1 < 2i1+1 − a1 and 2i2 − a2 < 2 < 2i2 from class (i1, i2 −1) to produce a particle in the class (i1, i2). Therefore, the rate of successful agglomeration becomes 17490
dx.doi.org/10.1021/ie401287w | Ind. Eng. Chem. Res. 2013, 52, 17487−17500
Industrial & Engineering Chemistry Research
Article
Figure 3. Patterns of the mechanism 2 (Birth).
1. Mechanism 2.3.1: If 2i2 < a2 < 2i2 − 2j2+1, then there is no chance for successful agglomeration of this particle with particles in class (i1, j2). 2. Mechanism 2.3.2: If 2i2 − 2j2+1 < a2 < 2i2 − 2j2, then this particle should be aggregated with a particle of volumes 2i1 < v1 < 2i1+1 − a1 and 2i2 − a2 < v2 < 2j2+1 from class (i1, j2) to produce a particle in the class (i1, i2). 3. Mechanism 2.3.3: If 2i2 − 2j2 < a2 < 2i2, then this particle should be aggregated to a particle of volume 2i1 < v1 < 2i1+1 − a1 and any v2 from class (i1, j2) to yield a new particle in the class (i1, i2).
dR i2.2 = β × N (i1, i2 − 1) × N (j1 , i2 − 1) 1, i 2 − 1, j , j 1 2
×
da1 da 2 2i1+ 1 − 2i1 − a1 2 j1 2i2 − 1 2i1
(14)
Integrating over the class (j1, j2) and the summation of these classes results in i1− 1
R i2.2 = 1, i2
∑ β × N(i1, i2 − 1) × N(j1 , i2 − 1) j1 = 1
× (1 − 3 × 2 j1 − i1− 1)
(15)
2.1.2.3. Mechanism 2.3 ((i1, 1 ≤ j2 ≤ i2 − 2) + (1 ≤ j1 ≤ i1 − 1, i2 − 1)). Consider a particle with volumes a1 < v1 < a1 + da1 and a2 < v2 < a2 + da2 from class (j1, i2 − 1), depending on the range of distance 2i2−1 to 2i2 which a2 located in, different situations can be considered:
Considering above statements, the rate of the agglomeration can be evaluated by
⎧ 0 2i2 − 1 < a 2 < 2i2 − 2 j2 + 1 ⎪ ⎪ ⎛ da da 2i1+ 1 − 2i1 − a 2 j2 + 1 − 2i2 + a ⎞ 1 2 ⎪ ⎟ 2i2 − 2 j2 + 1 < a 2 < 2i2 − 2 j2 ⎪⎜ j1 2 2.3 2i1 2 j2 dR i1, j , j , i2 − 1 = β × N (j1 , i 2 − 1) × N (i1, j2 ) × ⎨ ⎝ 2 1 2i2 − 1 ⎠ 2 1 ⎪ + i i 1 ⎛ da da 2 1 − 2 1 − a ⎞ ⎪ 1 ⎜ j 1 i −21 ⎟ 2i2 − 2 j2 < a 2 < 2i2 ⎪ ⎪ 2i1 ⎝ 2 1 22 ⎠ ⎩
Again, like the other mechanisms, integrating and taking summation yields
2.1.3. Mechanism 3. This mechanism is just the same as the mechanism 2 (see Figure 4) and final obtained equations are just presented here. 2.1.3.1. Mechanism 3.1 ((i1 − 1, i2) + (1 ≤ j1 ≤ i1 − 2,1 ≤ j2 ≤ i2 − 1)).
i1− 1 i2 − 2
R i2.3 = 1, i2
∑ ∑ β × N(j1 , i2 − 1) × N(i1, j2 ) j1 = 1 j2 = 1
× (3 × 2 j2 − i2 − 32 × 2 j1 + j2 − i1− i2 − 1)
(16)
(17) 17491
dx.doi.org/10.1021/ie401287w | Ind. Eng. Chem. Res. 2013, 52, 17487−17500
Industrial & Engineering Chemistry Research
Article
Figure 4. Patterns of the mechanism 3 (Birth). i1− 2 i2 − 1
R i3.1 = 1, i2
one type, for keeping the continuity of derivations, this one type is named as mechanism 4.1 (Figure 5). 2.1.4.1. Mechanism 4.1 ((i1, 1 ≤ j2 ≤ i2 − 1) + (1 ≤ j1 ≤ i1 − 1, i2)). Consider a particle with volumes a1 < v1
