Article pubs.acs.org/crystal
The Application of Two-Dimensional Population Balance Model To Study the Effect of Temperature Profile on the Crystal Size Distribution and Aspect Ratio Mitsuru Shoji*,† and Hiroshi Takiyama‡ †
R&D Center, Mitsubishi Chemical Corporation, 17-1 Towada, Kamisu-shi, Ibaraki 314-0102, Japan Department of Chemical Engineering, Tokyo University of Agriculture and Technology, 2-24-16 Nakacho, Koganei-shi, Tokyo 184-8588, Japan
‡
ABSTRACT: The effect of cooling profiles on the average crystal size and aspect ratio was studied as an example of the application of two-dimensional population balance model. The nonseeded cooling type batch crystallization system was chosen for this study. The simulation successfully provided the information about the crystal size distribution and the average aspect ratio of each crystal product obtained through different cooling profiles. In addition, using the proposed model, it was expected that the modulate cooling in which undersaturation is introduced during the cooling operation, that has already been reported in previous research as an effective method to improve crystal size distribution, is also applicable for the improvement of aspect ratio. However, it was also implied that the results were influenced by various factors such as nucleation rates, growth rates, and dissolution rates of length and width.
■
INTRODUCTION In the field of crystallization, many researchers have worked to improve and control the quality of the product crystals in terms of crystal purity, crystal shape, and crystal size distribution (CSD). For the description of CSD, one-dimensional distribution using a characteristic size is usually employed even for the system of nonspherical crystals (e.g., needle). There are some reports about the method to improve such a one-dimensional CSD,1−4 whether the studies are based on experiments or simulations. Crystallization processes were often modeled using one-dimensional population balances. In industrial crystallization, on the other hand, many of the crystal shapes are not spherical. The problem is that the shape of crystals is often considered as a key factor of the process. For example, it is known that the size of crystals and the ratio of length to width (aspect ratio) affect the efficiency of specific unit operations such as filtration, drying, and packing. For solid−liquid separation, which is usually conducted after crystallization, the crystals are required to be monodisperse and within the range of specific values of aspect ratio. Therefore, considering two or multidimensional crystal size and its distribution should be one of the key issues for the improvement of crystallization technology. Turning to the present state of the crystallization field, however, the research related to multidimensional crystal size and shape is limited. The numerical method for solving the multidimensional PB model is often discussed, and some researchers have proposed an application of various methods to crystallization simulation.5−7 Doherty et al. have discussed a mathematical model to predict the shape of a crystal that grows and dissolves.8 Although this research has provided very basic and important information, it is still difficult to find research in © 2012 American Chemical Society
which multidimensional crystal shape is considered to evaluate and to design the crystallization conditions for the control of crystal shape in industrial crystallizers. The purpose of this research is to study the influence of cooling profiles on the aspect ratio of crystals and the CSD in nonseeded batch type cooling crystallization based on simulation using the two-dimensional PB model. The reason to focus on nonseeded crystallization is that in the pharmaceutical field, where cooling type batch crystallization is widely used, seed crystals are not used because addition of foreign substances can cause contamination. It is expected that different cooling profiles provide products with different two-dimensional CSD, as many researchers did with one-dimensional systems.1−4 It will be very useful in industry when the two-dimensional PB model can be applied to design the crystallization conditions that provide crystal products with required specifications.
■
MATHEMATICAL MODEL Population Balance and Kinetics. Nucleation Rate. The number of total crystals NT [no./m3] is calculated by integrating the population density n(t,LS,LL) [no./m5] for all crystal sizes as eq 1. ∞
NT =
∫0 ∫0
∞
n(t , LS , L L) dLS dL L
(1)
where LL is the length and LS is the width. Then the crystal generation rate is expressed as the time derivertive of NT. Received: May 18, 2012 Revised: September 13, 2012 Published: October 1, 2012 5241
dx.doi.org/10.1021/cg300680f | Cryst. Growth Des. 2012, 12, 5241−5246
Crystal Growth & Design
Article
∞ ∞ dNT d = n(t , LS , L L) dLS dL L dt dt ∞ 0 ∞ 0 = {B(LS , L L) + D(LS , L L)} dLS dL L
Width: Birth and Spread model with size dependence
∫ ∫ ∫0 ∫0
(2)
where jP and jS are the primary and secondary nucleation rates represented using temperature T, solution concentration C, saturated concentration C*, and magma density MC. The size of nuculei formed by both primary and secondary nuculeation was assumed to be critical size, LC. Moreover, when it is assumed that the surface energy of a crystal is constant on each face, LC should have the same value for the direction of LS and LL because LC is now determined only by supersaturation. Therefore, the following expressions were used in this study.
(4)
⎧ k M (C − C*) LS = L L = LC jS = ⎨ b C ⎩0 otherwise
(5)
= 0
AL = A L0 exp(−E /(RT ))
(11)
σ1 = σ0/(RT )
(12)
D = kd(C − C*)
(13)
⎧ ⎛ εL4 ⎞1/5⎛ ν ⎞1/3⎫ ⎪ DAB ⎪ ⎨ kd = 2 + 0.8⎜ 3 ⎟ ⎜ ⎟ ⎬ ⎪ ⎪ L ⎩ ⎝ ν ⎠ ⎝ DAB ⎠ ⎭
(14)
kBT 2πηdm
(15)
In eqs 14 and 15, ε is specific power input, η is dynamic viscosity, ν is kinematic viscosity, ρ is density, kB is the Boltzmann constant, and dm is molecular size of diffuser. Since the calculation using eqs 13−15 requires cyrstal size L as one-dimensional, sphere equivalent size was substituted for eq 14 to obtain D as the dissolution rate of length. The dissolution rate of width was calculated with the assumption that aspect ratio of each crystal will not be changed by dissolution.8,14 Cross Moment. To simplify the evaluation of the results that should include two-dimensional distribution function, cross moments μij described below were also considered.15
(6)
∞
μij =
∫0 ∫0
∞
n(L1 , L 2)L1iL 2j dL1dL 2
(16)
With the definision of L1 as width and L2 as length, the following properties can be calculated from the CSD obtained as the solution of the PB model.
n(0, LS , L L) = 0
n(t , LS , 0)
(10)
DAB =
The boundary and the initial conditions are shown as eq 7.
= 0
⎛A ⎞ G L = ⎜ L ⎟σ 2 tanh(σ1/σ ) ⎝ σ1 ⎠
An estimate for the diffusion coefficient DAB can be obtained using the Stokes−Einstein equation, eq 15.
∂n(t , LS , L L) dL ⎫ ∂ ⎧ ⎨n(t , LS , L L) S ⎬ + ∂t ∂LS ⎩ dt ⎭ dL ⎫ ∂ ⎧ ⎨n(t , LS , L L) L ⎬ + ∂L L ⎩ dt ⎭
n(t , 0, L L)
(9)
where σ0 is a parameter that represents the depencency of the growth rate on temperature.10 Dissolution Rate. The dissolution rate is described as mass transfer using mass transfer coefficient kd estimated by eq 14.11,12
The death rate should be considered related to the breakage of crystals. The model of the breakage of needle-shaped crystals has already been proposed9 and also considered. As a result, it was found that in the condition assumed in this study the effect of the breakage on the CSD and aspect ratio was small enough to be neglected. Therefore, the death rate D(LS,LL) was assumed to be zero. Population Balance. The population balance equation is derived from eqs 1−5 as shown in eq 6.
= jP (T , S , LS , L L) + jS (T , S , LS , L L)
β = β0 /(RT )2
Length: BCF model
(3)
⎧ ⎤ ⎡ jb ⎪ j exp⎢ ⎥ L S = L L = LC 3 2 a jP = ⎨ ⎣ T ln(C /C*) ⎦ ⎪ ⎩0 otherwise
(8)
where α and β0 are parameters that represent the dependency of the growth rates on width13 and temperature,10 respectively.
where B(LS,LL) and D(LS,LL) are the birth and the death rates, respectively. The birth rate can be expressed as the nucleation rate. B(LS , L L) = jP + jS
GS = A S0(1 − e−αLS)σ 5/6 exp(−β /σ )
(7)
The number of crystals: NT [no./m3]
The growth rates should be substituted into dLS/dt and dLL/ dt when the condition is supersaturated (C − C* > 0) or the dissolution rates in undersaturation (C − C* ≤ 0). Growth Rate. The following equations were chosen for the description of the growth rate of width (GS) and length (GL), respectively, with the help of previous research.10−13 σ is the relative supersaturation.
NT = μ00
(17)
Number base average length: L̅L
L̅ L = μ01 /μ00
5242
(18)
dx.doi.org/10.1021/cg300680f | Cryst. Growth Des. 2012, 12, 5241−5246
Crystal Growth & Design
Article
Number base average width: L̅S L̅S = μ10 /μ00
(19) 3
Mass of crystals: MT [kg/m ]
M T = ϕv ρc μ21
(20)
where ϕv is volume shape factor and ρc is crystal density. Mass base average length: L̅L,m L̅ L,m = μ22 /μ21
(21)
Mass base average width: L̅ S,m L̅S,m = μ31 /μ21
(22)
Number base average aspect ratio, ar̅ , and mass base average aspect ratio, ar̅ ,m, are also calculated by eqs 23 and 24. ar̅ =
1 NT
∞
∫0 ∫0
∞
⎛L ⎞ n(L1 , L 2)⎜ 2 ⎟ dL1 dL 2 = μ−11 /μ00 ⎝ L1 ⎠ (23)
ar,m ̅ =
■
1 MT
∞
∫0 ∫0
∞
⎛L ⎞ ϕv ρc n(L1 , L 2)L12L 21⎜ 2 ⎟ dL1 dL 2 ⎝ L1 ⎠
= μ12 /μ21
Figure 1. Temperature profiles. (24)
many nuclei” at the beginning of the operation as the alternative of seed addition, because the target of this study is nonseeded batch cooling crystallization. The profiles were defined to examine the cooling curve of convex upward, convex downward, linear, and modulated. The reason a cubic function (third) was selected is that the cubic function is very typical and had already been known as the cooling profile that improves CSD. Natural cooling was selected as an example of uncontroled cooling, which is known to be an exponential function. The modulate cooling profile used here is just an example, in which undersaturation was introduced during cooling operation. The purpose of applying the modulate cooling profile was to examine and discuss the influence of modulation on the aspect ratio of crystals. According to previous research,34 it is known that modulate cooling provides better CSD (large average size and narrow distribution) because of the dissolution of small crystals. With the consideration by analogy, it is also assumed that modulation and dissolution must influence the aspect ratio because, in general, the values of aspect ratio are different between small crystals and larger ones. In Table 2, variated kinetic parameters are listed as cases 1, 2 and 3. Case 1 is the condition represented by the parameters
CALCULATION METHOD AND CONDITIONS There is much research about calculation methods to solve the multidimensional PB model. One work concluded by some researchers as an appropriate method, in terms of accuracy and short computing time, is “method of characteristic”.16,17 In this study, CIP method18,19 was used. CIP method is a kind of “method of characteristic” and was originally developed to solve the advection equation that includes population balance equation. The parameters used in this study are listed in Table 1. The values were estimated from the experimental data of DTable 1. Kinetic Parameters Used for Case 1 value
equation
ja [no./s] jb [K3] kb [no./s] AS0 [μm/s]
1.53 × 105 5.61 × 106 5.45 × 109 2.53
eq eq eq eq
α [1/μm] β0 [J2/mol2] AL0 [μm/s]
0.03 1.28 × 104 2.03 × 109
eq 8 eq 9 eq 11
E [J/mol] σ0 [J/mol] ε [J/(s·kg)]
35.7 × 103 75.0 × 103 0.2
eq 11 eq 12 eq 14
4 4 5 8
Table 2. Kinetic Parameters Changed in Each Case ja [no./s] kb [no./s] AL0 [μm/s]
mannitol−water system crystallization. D-Mannitol is a typical needle-shaped crystal. The two-dimensional PB model was solved, and the change of CSD and cross moment values in time were calculated in a nonseeded batch cooling type crystallizer. The cooling profiles studied here are shown in Figure 1. The profile “modulate” features the introduction of undersaturation. Total operation time and initial and final temperature were decided based on the existing experiments in the laboratory. The purpose of the first heating is to avoid “too
case 1
case 2
case 3
1.53 × 105 5.45 × 109 2.03 × 109
1.53 × 108 0 2.03 × 109
1.53 × 105 5.45 × 109 2.03 × 1010
listed in Table 1. Case 2 is a hypothecial situation in which small crystals are hardly generated. The values of the parameters related to the nucleation rates are changed. As mentioned above, the dissolution of small crystals influences the CSD. It is easily predicted that the dissolution of small crystals may also have an influence on the aspect ratio. The 5243
dx.doi.org/10.1021/cg300680f | Cryst. Growth Des. 2012, 12, 5241−5246
Crystal Growth & Design
Article
different profiles in cases 1 and 2. However in case 3, natural cooling resulted the smallest value of the aspect ratio. If only case 1 was examined, it can be concluded that modulate cooling is a better profile than three other profiles to obtain crystal products that have large size and small aspect ratio. However, it is not always true as seen above. It means that an appropriate cooling profile to obtain the products with specific properties is not always appropriate for other systems. One of interesting things is when and how the modulate cooling works to make product crystals have large size and small aspect ratio. To look into details, the change of the average sizes and aspect ratio in time during modulate cooling operations were investigated for cases 1, 2, and 3. Figures 2 and 3 show the change in average length and aspect ratio with time, respectively. In Figure 2, it was found that the average length increases when undersaturation was introduced in cases 1 and 3, conditions in which secondary nucleation occurs. On the other hand, in case 2 without secondary nucleation, the average length was decreased by the introduction of undersaturation. In addition, the change of the average length along with cooling was also opposite between the conditions with (cases 1 and 3) and without (case 2) secondary nucleation. In case 2, the average length got larger by cooling where in cases 1 and 3 the average length became smaller. With the consideration of eq 13, which tells that the smaller crystal dissolves faster than the larger one, all of the observations mentioned above imply that the average length is influenced by the number of small crystals generated by secondary nucleation. It can be summarized that the dissolution of small crystals caused by undersaturation is the reason why the modulate cooling works to produce crystals with larger average length. Since there was no secondary nucleus in case 2, the introduction of undersaturation did not work. In case 1, the dissolution of small crystals was considered as a main reason of larger average length, but it does not explain the reason the modulate cooling causes the smallest aspect ratio within three cases. It should be noted again that all of the calculations were conducted with the assumption that “aspect ratio of each crystal is not changed by dissolution” as described in the previous section. It was found that the influence of modulation on the average aspect ratio was different among cases 1, 2, and 3. The introduction of undersaturation causes lower average aspect ratio in case 1, no change in case 2, and larger average aspect ratio in case 3. The reason will be discussed in the rest of this section. As seen in Figures 2 and 3, the average length became smaller with cooling while the average aspect ratio became larger in case 1. This implies that the small crystals have relatively larger aspect ratio than the larger crystals. It should be noted again that the average length is influenced by the number of small crystals. In case 3 on the other hand, the situation is opposite. The average length became smaller with cooling as in case 1, but the average aspect ratio became smaller, unlike in case 1. In addition, when undersaturation was introduced in case 3, the average aspect ratio became larger. All of the observations mean that in case 3, the small crystals have relatively smaller aspect ratio than the larger crystals. It can be predicted whether the average aspect ratio increases or decreases if the average aspect ratio of the crystals is known. When undersaturation was introduced, the average aspect ratio
parameters varied in Case 2 were chosen to examine the influence of the number of small crystals on how the aspect ratio is affected by the modulation. In Case 3, another hypothetical set of parameters was selected in which the growth rate of length is larger than the rate in Case 1. If the dissolution of small crystals influences the aspect ratio, it is simply predicted that the aspect ratio is influenced by the ratio of growth rate of length to that of width. When small crystals have smaller aspect ratio than the crystals with larger size, the modulation would make the average aspect ratio larger. On the other hand, when small crystals have larger aspect ratio, the modulation would make the average aspect ratio smaller. In Case 3, the varied parameters were chosen to examine the effect of the growth rates.
■
RESULTS AND DISCUSSION The average crystal length, width and aspect ratio obtained through the cooling profiles shown in Figure 1 are summarized in Tables 3,4, and 5. Table 3. Average Length and Aspect Ratio of the Product in Case 1 L̅ L [μm] L̅ S [μm] ar̅ L̅ L,m [μm] L̅ S,m [μm] ar̅ ,m
modulate
linear
third order
natural cooling
52 22 2.37 166 63 2.86
54 21 2.47 153 52 3.20
57 24 2.36 159 57 2.93
51 21 2.37 154 54 3.07
Table 4. Average Length and Aspect Ratio of the Product in Case 2 L̅ L [μm] L̅ S [μm] ar̅ L̅ L,m [μm] L̅ S,m [μm] ar̅ ,m
modulate
linear
third order
natural cooling
1238 740 1.67 1251 749 1.67
1302 737 1.77 1305 740 1.76
1248 741 1.69 1252 744 1.68
1294 741 1.75 1298 744 1.75
Table 5. Average Length and Aspect Ratio of the Product in Case 3 L̅ L [μm] L̅ S [μm] ar̅ L̅ L,m [μm] L̅ S,m [μm] ar̅ ,m
modulate
linear
third order
natural cooling
159 20 7.66 512 26 17.4
144 20 7.00 426 24 16.1
191 20 9.33 449 22 19.3
137 20 6.68 437 25 15.9
When the mass base average length and width are compared between different cooling profiles, it was found that “modulate” profile results in the largest crystal size in cases 1 and 3. The results seem reasonable since the effect of modulation to improve the CSD has already reported in previous research.3,4 On the other hand in case 2, modulation did not work to make crystals larger. It is also predictable according to the previous research. Turning to the mass base average aspect ratio, the modulate cooling provided the smallest aspect ratio within the four 5244
dx.doi.org/10.1021/cg300680f | Cryst. Growth Des. 2012, 12, 5241−5246
Crystal Growth & Design
Article
Figure 2. Change in average length with modulate cooling operation (cases 1−3 correspond to those described in Table 2).
Figure 3. Change in average aspect ratio with modulate cooling operation (cases 1−3 correspond to those described in Table 2).
changes as small crystals vanish; if the small crystals have small average aspect ratio, dissolution of them results in larger average aspect ratio, and vice versa. In case 2, ar and ar,m were almost always the same value. This means that small crystals and large crystals have same aspect ratio. Therefore, the introduction of undersaturation, which causes crystal dissolution, has no influence on the average aspect ratio. In addition, the result of case 2 shows that the effect of nuclei formed by primary nucleation is almost negligible. Seeing through the results, it was implied that the number and shape of small crystals have a big influence on the average size and aspect ratio. How to control the number of small crystals must be one of the key issues for the design of crystallization conditions. The modulation cooling, in which undersaturation is introduced, may be one of the methods to
control the number of small crystals. However, the effect needs to be considered carefully because how the modulation would work depends on other factors, for example, the growth rates along each axis as examined in this study. The kinetic study of crystallization will become more and more important for advanced operation design. This is an example of how we can examine the effect of cooling profiles on crystal properties using two-dimentional PB model. This is also an example of how we can find a way to control the properties of product crystals. What was found through the discussion above could provide us with the possible strategy to control the average crystal size and the average aspect ratio of needle-shaped crystals using modulate cooling. 5245
dx.doi.org/10.1021/cg300680f | Cryst. Growth Des. 2012, 12, 5241−5246
Crystal Growth & Design
■
Article
(19) Yabe, T.; Aoki, T. Comput. Phys. Commun. 1991, 66, 219−232.
CONCLUSIONS The effect of cooling profiles on the average crystal size and aspect ratio were studied using two-dimentional population balance model. The nonseeded cooling type batch crystallization was chosen for the simulation. The simulation provided the information about the differences of each crystal produced by different cooling profiles. It was expected that modulate cooling, which has already been reported in previous research as an effective method to improve crystal size distribution, is also able to be applied to improve the aspect ratio of crystals in specific conditions. In addition, it was found that the effects of the modulation depend on various factors such as crystal size distribution itself, nucleation rates, growth rates of length and width, ratio of growth rates, and dissolution rates. The population balance model and its solution based on kinetic data will help engineers predict the results of operations and can become an effective tool to design operation conditions for the improvement of properties such as the crystal size and aspect ratio of product crystals under the conditions of nonseeding and nonadditives. Since all the results in this study are based only on the simulations, the model will need proper validation before it is applied to real processes. Even though the application of the proposed model is limited, authors strongly believe that the basic strategy to discuss and design the crystallization conditions using PB model for the control of the crystal shape and the size distribution is useful in many industrial scenes.
■
AUTHOR INFORMATION
Corresponding Author
*E-mail:
[email protected]. Notes
The authors declare no competing financial interest.
■
REFERENCES
(1) Heffels, S. K.; de Jong, E. J. AIChE Symp. Ser. 1991, 87, 170−181. (2) Moscosa-Santillán, M.; Bals, O.; Fauduet, H.; Porte, C.; Delacroix, A. Chem. Eng. Sci. 2000, 55, 3759−3770. (3) Takiyama, H.; Shindo, K.; Matsuoka, M. J. Chem. Eng. Jpn. 2002, 35, 1072−1077. (4) Shoji, M.; Eto, T.; Takiyama, H. J. Chem. Eng. Jpn. 2011, 44, 191−196. (5) Heineken, W.; Flockerizi, D.; Voigt, A.; Sundmacher, Kai Comput. Chem. Eng. 2011, 35, 50−62. (6) Ma, D. L.; Tafti, K.; Braatz, R. D. Ind. Eng. Chem. Res. 2002, 41, 6217−6223. (7) Briesen, H. Chem. Eng. Sci. 2009, 64, 661−672. (8) Snyder, R. C.; Doherty, M. F. Proc. R. Soc. A 2009, 465, 1145− 1171. (9) Sato, K.; Nagai, H.; Hasegawa, K.; Tomori, K.; Kramer, H. J. M. Chem. Eng. Sci. 2008, 63, 3271−3278. (10) Mersmann, A., Ed. Crystallization Technology Handbook; Marcel Dekker Inc.: New York, 2001. (11) Gahn, C.; Mersmann, A. Chem. Eng. Sci. 1999, 54, 1273−1282. (12) Gahn, C.; Mersmann, A. Chem. Eng. Sci. 1999, 54, 1283−1292. (13) Mydlarz, J. Cryst. Res. Technol. 1996, 31, 541−565. (14) Snyder, R. C.; Doherty, M. F. AIChE J. 2007, 53, 1337−1348. (15) Briesen, H. Chem. Eng. Sci. 2006, 61, 104−112. (16) Mesbah, A.; Kramer, H. J. M.; Huesman, A. E. M.; Van den Hof, P. M. J. Chem. Eng. Sci. 2009, 64, 4262−4277. (17) Hasebe, S.; Hamamura, T.; Naito, K.; Sotowa, K.-I.; Kano, M.; Hashimoto, I; Betsuyaku, H.; Takeda, H. J. Process Control 2000, 10, 441−448. (18) Takewaki, H.; Yabe, T. J. Comput. Phys. 1987, 70, 355−372. 5246
dx.doi.org/10.1021/cg300680f | Cryst. Growth Des. 2012, 12, 5241−5246
