1064
Ind. Eng. Chem. Res. 2006, 45, 1064-1073
Simulation of Semibatch Precipitation of Sodium Perborate Tetrahydrate in an Industrial Crystallizer D. Teslic*,† and C. Pohar‡ Belinka Perkemija Chemical Company, ZasaVska cesta 95, 1001 Ljubljana, SloVenia, and Faculty of Chemistry and Chem. Technology, UniVersity of Ljubljana, AskerceVa 5, 1000 Ljubljana, SloVenia
A model of an industrial semibatch sodium perborate tetrahydrate (SPT) precipitation process, based on mass, population, and energy balances, is presented. The model parameters of growth, nucleation, and agglomeration were obtained by nonlinear parametric optimization from the experimental data for the sodium perborate tetrahydrate concentration in the mother liquid, the mass of the solid phase, the crystal size distribution, and the suspension temperature as a function of time. A combined integral-differential strategy was used to determine the kinetic parameters. The original set of five kinetic parameters was reduced to two subproblems, each of which was solved consecutively, and in this way, the parameters were calculated much faster. The model was verified by comparison with measurements obtained in precipitation experiments, which were carried out under process conditions different from those used in the reference run. A rather good agreement between the model predictions and the experimental data was found. 1. Introduction
2. Theory
Sodium perborate tetrahydrate (SPT) and its anhydrous form sodium perborate monohydrate have found many applications on account of their stability as oxidizing agents, especially as bleaching agents. The process of preparing sodium perborate tetrahydrate involves the chemical reaction between hydrogen peroxide and sodium metaborate (SMB) in mother liquor saturated with sodium metaborate and sodium perborate, with appropriate cooling of the exothermic reaction:
Sodium Perborate Tetrahydrate Precipitation Model. The mathematical model for semibatch precipitation is based on the population balance approach. The kinetic processes of growth, nucleation, and agglomeration depend on supersaturation and temperature. Thus, for a complete description of the crystallization process, mass and energy balances should be incorporated into the model. (A) Mass Balance. The volume change of solution for a given system is given by an equation which assumes a constant density of solution
NaBO2 + H2O2 + 3H2O f NaBO2H2O23H2O
(1)
The reaction is supposed to be instantaneous, whereas the crystallization rate is slower. SPT has a tendency to agglomeration. The solubility of SPT is a function of the temperature and the amount of metaborate present in the solution.1 The outstanding characteristic of the mother liquor is its tendency to supersaturate with respect to sodium perborate; the relevant factors, which modify this tendency, are the compositions of the reactants, the composition and temperature of the mother liquor, and the efficiency of mixing. An involuntary change of any of these parameters can be responsible for a severe loss in product quality. Because of the complex influences of all the above-mentioned process parameters, a mathematical model of the industrial crystallization process would be desirable. In recent years, considerable progress has been achieved in understanding the analysis, design and performance of all kinds of precipitation systems. However, precipitation models were applied mostly to laboratory scale crystallizers,2-5 while for real industrial scale systems successfully applied and validated models are rare. The purpose of the present work was to apply process simulation analysis using the principles of mass, population, and energy balances to investigate the semibatch nonisothermal reactive precipitation system of sodium perborate tetrahydrate, carried out in a 7000 L industrial crystallizer. * To whom correspondence should be addressed. E-mail:
[email protected]. Fax number: +38615886303. † Belinka Perkemija Chemical Company. ‡ University of Ljubljana.
FSMB FML dV FH2O2 fH2O2 + fSMB + f + ) dt Fsusp Fsusp Fsusp ML dms 1 1 (2) dt Fs Fsusp
[
]
whereas the mass change of solid phase is given by
dms ) 3FskvV dt
∫LL 1
N+1
ψL2G dL + RkVFSL13VB0n
(3)
where the terms on the right-hand side represent the increase of the solid phase due to crystal growth and nucleation. The nucleation term is nil except for class n ) 1. In the case of no excess of hydrogen peroxide, the mass balance of dissolved sodium perborate present can be stated as
MSPB dms dml ) cH2O2fH2O2 dt MH2O2 dt
(4)
The excess mass of sodium metaborate present in solution is calculated as
MSMB dmSMB ) fSMBcSMB - fH2O2cH2O2 dt MH2O2
(5)
The solubility of sodium perborate tetrahydrate depends on temperature T and the excess of sodium metaborate wSMB present in solution. Using the experimental data obtained by Van
10.1021/ie049039t CCC: $33.50 © 2006 American Chemical Society Published on Web 01/07/2006
Ind. Eng. Chem. Res., Vol. 45, No. 3, 2006 1065
Gelder1, it was correlated as follows (T in °C):
Nuc )
c* ) 0.7171wSMB2 e-0.06261T + wSMB(0.0009938T2 + 0.004752T - 0.1946) + 5.895 e-0.08564T wSMB + 0.000914T + 0.005685T + 0.817 2
0.825Re0.56Pr1/3
mcp
dt
)
∑i
η ηw
0.14
dms dml fiFiHi ∆Hk ∆Hr - φmcp(Tout - Tin) dt dt (7)
(8)
The fourth term on the right-hand side of eq 7 represents the heat transferred to the cooling liquid:
φmcp(Tout - Tin) ) UA
(T - Tin) - (T - Tout) T - Tin ln T - Tout
(9)
where A is the area available for heat transfer and U is the overall heat transfer coefficient,
()
do ln(do/di) 1 1 1 do ) + + U hi di λ ho
(10)
The heat transfer coefficients on the inner and outer sides of the cooling coil, needed for a numerical evaluation of eq 9, were estimated by standard correlation equations. The inner side heat transfer coefficient can be expressed by the Dittus-Boelter equation
Nu ) 0.023Re0,8Pr0,4
(11)
The Reynolds number, Re, is defined as
Re )
FVdi η
Nuc ) 1.05Re0.62Pr1/3
cpη λ
hidi λ
ni
(np)0.15
Ci
0
do D
nil
(13)
(14)
The outer side heat transfer coefficient can be estimated by the standard correlation6
-0.3
(15)
( ) () ( ) η ηw
Re )
0.14
d D
-0.25
nib D
0.15
()
(np)0.15
Dc D (16)
Nrd2Fsusp η
(17)
and the Nusselt number, Nu, is given by
Nu )
hoD λ
(18)
One should note that the characteristic dimension in the Nusselt number expression is the diameter of the crystallizer D, while in the Reynolds number expression it is related to the impeller speed Nr and diameter d. (C) Population Balance. Modeling of the precipitation process followed a procedure suggested by David et al.,2-5 which is based on the population balance equation introduced by Randolph and Larson.7 The computational approach used the method of classes proposed by Marchal et al.8 The population balance for this case may be written as
d(VNn) + V[G(Ln+1)ψ(Ln+1) - G(Ln)ψ(Ln)] ) dt VB0n + VRA,n; n ) 1, ..., nc (19) The first term on the left-hand side of eq 19 represents the accumulation of crystals in the nth size class. The second term represents the rate of change of the number of particles in the differential size range L to L + dL to the rates of growth in to and out of that range. The terms on the right-hand side are the source and sink terms due to primary nucleation and agglomeration of particles. The nucleation term is nill except for the case of the class n ) 1. The elementary crystallization processes are described as follows. The overall nucleation rate is the sum of the primary and secondary nucleation rates: 0 ) Bn)1
and the Nusselt number, Nu, is given by
Nu )
( ∑ )( )
0.15
which is valid when the impeller is located underneath the coil array. In the latter case, the Reynolds number in the crystallizer, Re, is defined by
(12)
The Prandtl number, Pr, is given by
Pr )
ni b D
which is used when the impeller is located inside the coil array, or by the correlation6
where Hi is the enthalpy of the ith reactant flow of mother liquor, sodium metaborate, and hydrogen peroxide:
Hi ) cpi(Ti - Tr)
-0.25
d D
(6)
(B) Energy Balance of the Crystallizer. In a larger, industrial crystallizer, the ratio between the cooling coil area and the suspension volume is lower than that in a pilot or laboratory scale crystallizer. That is why the temperature is difficult or even impossible to keep constant during precipitation. Thus, beside the mass and population balances, energy balances accounting for the temperature inside the crystallizer and the temperature inside the cooling coil are required to complete the model. For semibatch precipitation of sodium perborate, the change of the temperature T inside the vessel is given by
dT
( ) () ( )
R [kn(c - c*)n + kns(c - c*)mS] kvFsL13
(20)
Crystal growth results from a combination of solute integration into the crystal and transfer from the bulk solution to the crystal surface. A quantitative measure of the degree of diffusion or surface integration control is made through the concept of the effectiveness factor ηr:9
G)
ks dL ) k η (c - c*)g dt 3kvFs g r
where ηr is calculated from
(21)
1066
Ind. Eng. Chem. Res., Vol. 45, No. 3, 2006
[
]
kg (c - c*)g-1 ηr + ηr1-g - 1 ) 0 kd
(22)
Newton’s method was used to find the roots of eq 22. The mass transfer coefficient kd is a function of the hydrodynamics and of the dimension of the crystal. According to Mersmann et al.,10 it can be calculated by
[ ( )( ) ]
DAB jL4 kd(L) ) 0.8 3 L νl
1/5
νl DAB
1/3
+2
(24)
The impact on class n of agglomeration between particles of classes i and j g i is represented by stoichiometric coefficients νn,i,j in analogy with a chemical reaction system. The corresponding stoichiometric coefficients that stand for the impact of agglomeration (i, j) on particles of class n are12
δn,i+1 - δn,i 2
νn,j,iK(i, j)
(27)
(28)
The expression for K(i, j) for large crystals (Sj > ηk) is given by David et al.2 for the chemical regime:
Si + Sj λe
2
+2
(30)
( ) 60νl
1/2
(31)
where u′ and can be evaluated from velocity and power dissipation measurements in a discharge flow14
u′ ) 0.3πNrd, ) 10P
(32)
The agglomeration rate of small crystals (Sj < ηk) was calculated by Higashitani’s15 expression, developed for the case of spherical crystals of classes i and j.16
Ra2,n )
[]
i
∑ ∑ i)1 j)1
νn,j,ikalNiNj(Si + Sj)3
ν
1/2
ηa
(33)
where ηa is the efficiency of the agglomeration defined in terms of characteristic times16
ηa )
tcr )
i
[ ( )]
1/3
1 tcr 1+ td
(34)
(26)
The term δn,i is an element of the Kronecker matrix (δn,i ) 0 if n * i and δn,i ) 1 if n ) i). These coefficients are calculated in order to balance the solid volume, which is equivalent to the conservation of the 3rd moment of the distribution, and to remove a single particle for each agglomeration (except for agglomeration (i, i) where only 1/2 of a particle disappears due to symmetry) to permit compliance with the 0th moment equation of the particle size distribution (PSD). The resulting agglomeration rate for class n is12
K(i, j) ) kaSi(1 + F)2f(F)Nrd 1 -
νl DAB
Crystallization time may be expressed as a function of the growth rate G of the crystal
Vj Vi δ δ νn,j,i ) Vi - Vj n,i Vi - Vj n,j
∑ ∑ i)1 j)1
λg ) u′
(25)
δn,i+1 δn,i - δn,i-1 νn,i-1,i ) 2 2
nc
1/5
The Lagrangian microscale λl is assumed to have the value of the Taylor microscale, which can be calculated by the approximation reported by Costes and Couderc:13
nc
RA,n ) Ra1,n + Ra2,n
Ra1,n )
[ ( )( ) ]
DAB jSj4 0.8 3 kd νl
(23)
The symbols j and νl in eq 23 denote the specific energy dissipation rate within the crystallizer and the kinematic viscosity of the liquid phase, respectively. The density, viscosity, and diffusivity of the sodium perborate tetrahydrate solution as a function of concentration have been reported by Frances et al.11 The overall agglomeration rate was represented as the sum of the agglomeration rates of class n for small and large crystals:
νn,i,i )
δb )
where F represents the crystal size ratio Sj/Si and f(F) is the relative shape function of both crystals. For spheres, f(F) always ranges between 10.8 and 12.2 H is the Heaviside step function (H(x) ) 1 for x g 0, H(x) ) 0 for x < 0). The boundary layer thickness δb can be estimated from standard correlations, for instance, that of Mersmann et al.:10
(35)
The disruption time is expressed as17
td )
σ*Lc [AlFsuspSi]
(36)
The term σ*Lc/Al was fixed at 1 in accordance with the litterature,17 although it is appreciated that σ* and Lc can vary from polymorph to polymorph. According to the literature, ka2 should be equal to 0.16.18 Small crystals are those smaller than the Kolmogoroff microscale size, which can be estimated from
ηk ) [(η/Fl)3/j]1/4
(37)
In the case of slow processes, crystals are rapidly carried around the recirculation path within the tank and are, thus, exposed to the average power dissipation. The mean energy dissipation rate can be calculated using the power number Np, and the tank and impeller size:
×
H(Si - δb)H(λl - Si - Sj)NiNj(c - c*)g (29)
Sjf(F) G
j )
4NpNr3d5 πD2l
(38)
For present case, the Kolmogoroff microscale size and mean energy dissipation rate were in the range ηk ) 20-45 µm and j ) 0.2-0.85 m2 s-3. The additional assumptions of the model are the following: (i) The system is well mixed with negligible hydrodynamic effects. There was little modification observed in the concentra-
Ind. Eng. Chem. Res., Vol. 45, No. 3, 2006 1067
tion curve and crystal size distribution (CSD) as a result of changing the reactant feed points from that on the surface to that close to the impeller. Therefore, according to the discussion of Marcant and David19 about the importance of good micromixing in the reactor, we assume that micromixing effects can be neglected. (ii) The redox reaction of hydrogen peroxide with sodium metaborate is supposed to be instantaneous. (iii) The temperature change in the observed time interval is so small that constant values of the kinetic parameters could be assumed. (iv) The crystal shape is spherical. (v) No breakage occurs. 3. Parameter Identification Once the structure of the model is known, the next step in the modeling process is identification of the five unknown kinetic parameters included in the model using the method of nonlinear parametric optimization. An objective function is selected to calculate the parameters of the model using the method of nonlinear programming. The function is a measure of how well the model describes the experimental data. The objective function can be calculated in two ways. The first uses the so-called integral approach, where for a certain set of parameters the equations of the model are integrated to obtain the time course of the dependent variables. Then taking into account the experimental values of the dependent variables, the values of the objective function is calculated. The difference between the predicted and the real course directs the selection of values in the following step. In this way, a relatively accurate parameter calculation is possible. However, the integration of the model equation is very time consuming. The special advantage of the integral method is that it allows kinetic parameters to be calculated using incomplete data on the particle size distribution because the value of the objective function can be calculated at any stage by simply using the information on concentrations and large particle classes. Another way of determining the kinetic parameters is the application of the differential approach introduced by Bramley et al.20 The advantage of the differential method is that instead of repeatedly integrating the model equations for every new set of kinetic parameters, only the derivatives are calculated. This method is much faster but less accurate since there is no feedback information about model predicted values of concentrations and crystal size distributions, as provided by the integral method. The original task can be divided into three subjobs, each of which can be solved consecutively, and in this way, the parameters can be calculated much faster.21 The problem when using the differential method occurs in determining the kinetics of nucleation since in actual practice it is difficult to get reliable experimental data on the number of particles in the small size classes. Given that we have relatively reliable information about the concentration of the solute and to some extent less reliable data about particle size distribution, we decided to use the following combined integral-differential strategy to determine the kinetic parameters.22 Calculation of the Particle Agglomeration Constant and Nucleation Parameters using the Integral Method. The agglomeration and nucleation parameters can be calculated theoretically by using the faster differential method. However, the integral approach was used to determine these parameters, because the influences of the agglomeration and nucleation processes are especially noticeable in the smallest particle classes, as well as in the early stage of the process for which
Figure 1. Schematic diagram of the crystallizer.
we have no reliable experimental data on the number of particles. The objective function to be minimized was defined as: Ns
Φc,m ) wfc
∑
N Ns
(cexp µ
∑ ∑(mµ,iexp - mµ,i)2 i)1 µ)1
- cµ) + wfm 2
µ)1
(39)
On completion of this step, we have an estimate of the set of kinetic parameters, as well as of the entire course of the crystal mass in each particle size class. Determination of Growth Parameters kg and g using the Differential Method. Because the approximate course of the crystal mass has been reconstructed in the earlier step (this being determinative for the course of concentration), new growth rate parameter values can now be calculated using a concentration criterion containing experimentally determined (c˘ exp µ ) and model calculated (c˘ µ) concentration derivatives:
Φdc )
2 ∑µ)1N [c˘ exp µ - c˘ µ] s
(40)
This criterion contains information about the mass of solute in solution, i.e., the mass of solid phase at a certain time.
c˘ )
dc d(ml/Vl) 1 dml ml dVl ) ) - 2 dt dt Vl dt V dt
(41)
l
The two kinetic processes contributing to enlarged crystal mass during crystallization are crystal growth and nucleation. Taking into account the fact that the influence of nucelation on the total crystal mass is negligible in comparison with its influence on crystal growth and can, therefore, be disregarded, the incorporation of eqs 2-4 into the above equation leads to the next correlation for the time derivative of the solute concentration in solution:
(
MSPT dc 1 ≈ cH2O2fH2O2 - 3FskVV dt Vl M H2 O2 ml
(
Vl Fl 2
∫LL
N+1
)
ψL2G dL -
1
∑ Fifi - 3FskVV∫L
LN+1
ψL2G dL) (42)
1
Index i denotes sodium metaborate, hydrogen peroxide, and mother liquor. Values Vl, V, and ml can easily be derived by using the mass balance equations and experimental data. The experimental data
1068
Ind. Eng. Chem. Res., Vol. 45, No. 3, 2006
Figure 2. Concentration of sodium perborate solution at different reactant flow rates as a function of time. Table 1. Expressions of Nucleation and Growth Rates according to the Literature Nucleation ref no.
nucleation rate 10-9
24
RN ≈
23 22 25 5
RN ) (3.44 × 10-10)∆w4.71 RN ) (2.768 × 10-10)∆c4.741 RN ) (2 × 10-7)∆c3.93 rn ) (89.5 × 107) × exp(-0.434/ ln(cn/cn*)2) RN ) (2.999 × 10-11)∆c4.785
this work
units
exp method
RN (kg m-3 s-1)
18
0.01
1 × 10-9
kg m-3 s-1 kg m-3 s-1 kg m-3 s-1 m-3 s-1
measmnt of the mass precipitated in a solution kept at const supersat metastable range width measmnt semibatch metastable range width measmnt MSMPR
20-30 15 8.8-13.5 23
S ∼ 1.4-1.6 0-0.065 0-0.010 0-0.070
kg m-3 s-1
semibatch
26-28
0-0.065
3.19 × 10-4 1.41 × 10-6 2.36 × 10-4 7.82 × 10-12 for L1 ) 1 µm 1.48 × 10-7
kg
m-3 s-1
supersat domain ∆w (kgsalt kgwater-1)
T (°C)
Growth
27 28 24 26 25 22 29 5 this work
growth rate G or specific growth rate RG
units
exp method
T (°C)
supersat domain ∆w (kgsalt kgwater-1)
RG ) (1.34 × 10-3)∆c RG ) (1.28 × 10-3)∆c1.132 G ) 4.44 × 10-5(cn/cn* - 1)2 G ) (8.53 × 10-6)∆w RG ) (1.49 × 10-6)∆c1.5 RG ) ηr(4.812 × 10-6)∆c1.086 RG ) (4.15 × 10-3)∆w G ) 5.13 × 10-11(cn - cn*)2 RG ) ηr(6.404 × 10-6)∆c1.0784
g cm-2 s-1 g cm-2 s-1 cm s-1 m s-1 kg m-2 s-1 kg m-2 s-1 kg m-2 s-1 m s-1 kg m-2 s-1
rotating disc fluidized bed face growth measmnt batch face growth measmnt semibatch fluidized bed MSMPR semibatch
23 25 19 20 21 15 20 23 26-28
0-0.015 0-0.015 0.003-0.008 0-0.015 0-0.01 0-0.065 0.0-0.014 0-0.070 0-0.065
we have on particle size distribution in the lower class sizes are incomplete. To calculate the growth rate, we therefore used an estimated particle size distribution calculated in the previous step by the integral method. The local approximate values of experimentally determined derivatives of concentration were obtained by numerically differentiating the experimental data exp exp cµ+1 ∂cexp - cµ-1 µ = ∂t tµ+1 - tµ-1
(43)
exp exp where cµ+1 and cµ-1 represents points next to the value cexp µ , measured at times tµ+1 and tµ-1.
RG (kg m-2 s-1) 7.99 × 10-5 3.88 × 10-5 1.64 × 10-5 4.48 × 10-5 1.05 × 10-5 3.40 × 10-5 ηr ) 1 2.49 × 10-5 6.36 × 10-5 4.46 × 10-5 ηr ) 1
The above two steps were repeated until the parameter values converged. It should be pointed out that the experimental errors could have influence on the accuracy of the determined kinetic parameters. On the basis of repeated experimental runs, the relative experimental errors of about 5% for the concentration and about 10-15% for crystal size distribution were estimated. These errors certainly affect the accuracy of the determined kinetic parameters, but the precise influence has not been determined yet. This problem was theoretically tackled by Livk et al.22 They showed that the experimental error of 5% does not significantly influence the values of the determined parameters by the differential method.
Ind. Eng. Chem. Res., Vol. 45, No. 3, 2006 1069
homogeneity, the stirring speed was increased 45 s before the sample was taken. No significant differences were noticed in the CSD of samples taken at different locations in the crystallizer. The sample was filtered, and the crystals were washed with a 2-propanol solution. The crystal size distribution (CSD) of the dry solid product was determined by sieving. A set of sieves of 100, 150, 200, 250, 425, and 850 µm was used according to ISO Standard 3118. The concentrations of hydrogen peroxide, SMB, and SPT were determined by titration of filtered samples using standard procedures. The solution was diluted to avoid any additional crystallization. Other quantities such as overall mass, solid phase mass, suspension volume, and excess of SMB, as well as supersaturation, were easily derived from the mass balance. 5. Results and Discusion Figure 3. Mass of the product at different flow rates as a function of time.
4. Experimental Section Crystallization experiments were carried out in a flat-bottomed 7000 L industrial crystallizer. The stainless-steel cylindrical vessel was equipped with a cooling coil and a 4 pitched-blade turbine agitator. The mechanical stirrer (d ) 1.17 m, Np ) 1.36) was operated in the range from 25 to 50 rpm. For temperature measurement, a Pt100 resistance thermometer was used. It was attempted to keep the temperature in the crystallizer at the desired value by control of the flow of cooling water. The outlet from the crystallizer was connected to a centrifuge, which served as the separation unit. All measured variables were displayed and stored in a process computer. The process computer was also used to monitor and control the temperature and reactant flow rates. The complete experimental setup is schematically presented in Figure 1. A semibatch experiment was conducted as follows. SMB and hydrogen peroxide solutions were fed into the initial volume of the SPT solution. The addition of hydrogen peroxide started several minutes after the addition of SMB to preserve the prescribed molar ratio of the reagents. The feed rates were constant during each experiment. At predetermined time intervals, suspension samples were taken by means of a vacuum pump. To achieve suspension
Figure 4. Evaluation of the particle size distribution as a function of time.
With the information available about the process, it is almost impossible to define the secondary nucleation constant. Within the selected parameter values, the secondary nucleation constant has only a limited effect on the course of concentration. We have no data, however, on the number of particles in lower classes where its influence is the largest. The secondary nucleation constant was, therefore, left out of the parameters we wished to reconstruct on the basis of experimental data, i.e., it is supposed to have a zero value. The model parameters for sodium perborate tetrahydrate precipitation were already determined for a laboratory scale system and a slightly different agglomeration model by Livk.22 The following values were obtained:
kn ) 2.768 × 10-10 kg-1 m3n-3 s-1
n ) 4.741
kg ) 4.812 × 10-6 kg1-g m3g-2 s-1 ka ) 0.595 × 10
-15
6
m kg
g ) 1.086 -1 -1
s
These values were chosen as an initial guess for our optimization procedure, although we were aware that a different agglomeration model was being used. The constrained optimization problem was solved by the sequential quadratic programming method, which is one of the components of the MATLAB
1070
Ind. Eng. Chem. Res., Vol. 45, No. 3, 2006
Figure 5. Concentration of sodium perborate solution as a function of time in the presence of seeds.
software. The kinetic parameters obtained by the optimization procedure are given as follows:
estimated with the same value of supersaturation, ∆w ) 0.006 kgsalt/kgwater. This is the same value as that used in the comparison made by David et al.,5 for reasons of simplicity. It could be seen that the orders of magnitude of the growth rate all fall in the same range. There are several factors, which could be responsible for the different orders of magnitude of the nucleation rates reported in the literature. Nucleation in industrial crystallizers is in most cases heterogeneous, while most authors examined homogeneous nucleation (Chianese et al.,23 Dugua & Simon,24 and Livk et al.25). It has also been shown24 in the case of sodium perborate tetrahydrate that the stirring rate has a strong influence upon the nucleation rate of the solution. It should also be pointed out that secondary nucleation could play a major role in the second part of the process, when the crystals have reached a sufficient number and range of sizes; nevertheless, we believe that primary nucleation is crucial for the course of the process. However, the influence of secondary nucleation could be included in the primary nucleation constants. Chianese et al.23 also reported a strong decrease of the SPT nucleation rate due to the presence of SMB in aqueous solution.
kn ) 2.999 × 10-11 kg-1 m3n-3 s-1
n ) 4.785
6. Comparison of the Model and the Experimental Results
kg ) 6.404 × 10-6 kg1-g m3g-2 s-1
g ) 1.078
In this section, simulation results for the crystallizer model are presented. By comparison of the simulated behavior with measured data, the model for the industrial scale crystallizer may be validated. The following experimental parameters were varied: (i) the flow of reactants; (ii) the initial volume of the mother liquid; (iii) seeds/no seeds. Concentration differences of the reactants in various runs are mainly due to technological factors in the real industrial system. However, all concentrations were kept constant during each experiment. From the estimated kinetic parameters and the accepted precipitation model, the time variation of the CSD, mass, concentration, and temperature in the reactor was calculated. The experimental and model predicted values for the mass, concentration, particle size distribution, and temperature courses are presented in Figures 2-8. From Figures 2 and 3, a good
ka ) 1.655 × 10-12 kg-g m3g+1 It is important to note that the evaluated kinetic parameters are accurate only for a given temperature range determined by the experimental conditions. The exact temperature dependence of the kinetic parameters is the subject of our further investigations. Comparison with the Literature Data. Different growth and nucleation laws collected from the literature are summarized in Table 1. Generally, different authors, using various measurement techniques, proposed various forms of the growth and nucleation equations for sodium perborate tetrahydrate. To compare these proposed growth and nucleation rates with those predicted in the present work, all rates have been calculated or
Figure 6. Evaluation of the particle size distribution as a function of time in the presence of seeds.
Ind. Eng. Chem. Res., Vol. 45, No. 3, 2006 1071
Figure 7. (a-c) Differential particle size distributions at the end of crystallization at different reactant flow rates and (d) in the presence of seeds for feed flow conditions f.
Figure 8. Suspension temperature as a function of time at different reactant flow rates.
agreement between the experimental and the model predicted concentration and mass curves for different reactant flow rates is evident. Here, the flow of reactants was chosen as a process variable. The following values for the flow of reactants were taken: 0.75f (Run 1), 1.25f (Run 3), and nominal flow, 1.00f (Run 2), where f represents the nominal flows of 10 L/min of H2O2 and 50 L/min of SMB. The experimental and calculated evaluation of the crystal size distribution during the precipitation process is presented in Figure 4. One can observe a rather good model prediction of
the crystal size distribution during the precipitation process. Particles pass to higher classes due to growth and mainly to agglomeration. It is obvious from the results presented in Figures 2-4 that the accuracy of the prediction of the course of the particle size distribution is much more sensitive to errors in estimating the kinetic parameters than the prediction of the mass or concentration course. This can be attributed to the neglible effects of nucleation and agglomeration on the total mass of the crystals.
1072
Ind. Eng. Chem. Res., Vol. 45, No. 3, 2006
Measurement of the particle size distribution is still the major problem in monitoring the crystallization process. The smallest particles form and disappear in the early stages of the precipitation process. In this early stage of nucleation, it is practically impossible to obtain the amount of material required to conduct a reliable measurement of the particle size distribution. Also, as mentioned before, the secondary nucleation constant has been omitted from the parameters we wish to reconstruct. These facts can lead to an inaccurate determination of nucleation and agglomeration parameters, which can be reflected in disagreement between model and experimental values of the crystal size distribution. Nucleation already begins before t ) 13 min, but it is not detected before that time. On the basis of the rather good agreement between model predicted and experimental values for crystal size distribution at later times, the correct prediction of CSD at the beginning of the precipitation process is assumed. In Figure 7a-c, the predicted and experimental final differential CSD as a function of the flow of reactants are compared. The capability of the time evaluation and final CSD prediction of the model in the case of seeded and nonseeded precipitation is presented in Figures 6 and 7d. As it is seen from Figure 7, a very good agreement with experimental data was obtained. It can be seen that different flow rates of reactants do not change the final crystal size distribution. A good prediction of temperature during the precipitation is of paramount importance for proper process control. In the nucleation period, the fastest increase of the solid phase concentration occurs and, also, the highest amount of crystallization heat is released, while the heat flow due to chemical reaction is constant all the time chemical reaction takes place. Thus, at the beginning of the process the suspension temperature rises because of insufficient cooling due to the low suspension level. Since the process is semibatch, the suspension level and, therefore, the exposed surface of the cooling coil rises as the process continues. The reason for a slight temperature decrease in the second half of the process is the increase of the cooling area. The experimental and predicted temperature profiles during precipitation are presented in Figure 8. A reasonable agreement between the experimental and model predicted temperature course is obvious. 7. Conclusions A kinetic model based on growth, nucleation, and agglomeration has been applied to the semibatch precipitation of sodium perborate tetrahydrate on an industrial scale. Kinetic parameters were estimated from the supersaturation curves and crystal size distribution. A combined integral-differential strategy was used to determine the kinetic parameters. The original set of five kinetic parameters was reduced to two subjobs, each of which was solved consecutively, and in this way, the parameters were calculated much faster. It is important to note that the evaluated kinetic parameters are accurate only for the given temperature range determined by the experimental conditions. This should be further investigated. The energy balance had to be incorporated into the model because of the variation of temperature in the industrial crystallizer. On the basis of the information available about the process, it is almost impossible to define the secondary nucleation constant. The secondary nucleation constant has, therefore, been left out of the parameters we wished to reconstruct on the basis of the experimental data, i.e., it is supposed to have a zero value. The simulation results were compared to measurements of the temperature, supersaturation, mass of crystals, and crystal
size distribution. These comparisons lead to a rather good agreement between simulations and experiments. Acknowledgment The authors thank Dr. Iztok Livk for helpful discussions. This work was supported by Belinka Perkemija, Chemical Company, and the Slovenian Ministry of Education, Science and Sport under Grant No. P1-0201 (Physical Chemistry). Nomenclature ∆Hk ) heat of crystallization, J kg-1 ∆Hr ) heat of reaction, J kg-1 A ) total area of heat transfer, m2 Al ) laminar efficiency constant b ) impeller width, m B0i ) nucleation rate, m-3 s-1 c ) concentration, kg m-3 c* ) solubility, kg m-3 Ci ) ith impeller height, m cp ) specific heat, J kg-1 K-1 d ) impeller diameter, m di ) inside diameter of the cooling coil, m do ) outside diameter of the cooling coil, m D ) vessel diameter, m Dc ) loop diameter of the cooling coil, m f ) volume flow rate, m3 s-1 g ) growth rate parameter G ) growth rate, m s-1 ∆H ) enthalpy, J kg-1 H ) Heaviside function hi ) heat transfer coefficient on inner side of cooling coil, W m-2 K-1 ho) heat transfer coefficient on outer side of the cooling coil, W m-2 K-1 K(i, j) ) elementary agglomeration rate, m3 s-1 ka1 ) agglomeration rate parameter for small crystals ka ) agglomeration rate parameter, m3g+1 kg-g kg ) growth rate constant, kg1-g m3g-2 s-1 kn ) primary nucleation parameter, kg1-n m3n-23 s-1 kns ) secondary nucleation parameter, kg-2 s-1 kv ) volume shape factor l ) height of liquid in the reactor, m L ) crystal size, m Lc ) contact length between particles, m M ) molar weight, kg kmol-1 ML ) mother liquor ml ) mass of dissolved sodium perborate, kg ms ) mass of solids, kg mµ,i ) model calculated mass of crystals of granulometric class i in experimental point µ exp mµ,i ) experimentally determined mass of crystals of granulometric class i in experimental point µ n ) primary nucleation parameter N ) crystal concentration, m-3 Nr ) stirring speed, s-1 Ns ) number of experimental points ni ) number of impellers np ) number of paddles nc ) total class size number Np ) power number P ) power input from the stirrer per unit mass of fluid, m2 s-3 RA ) overall agglomeration rate, m-3 s-1
Ind. Eng. Chem. Res., Vol. 45, No. 3, 2006 1073
Ra1 ) agglomeration rate of large particles, m-3 s-1 Ra2 ) agglomeration rate of small particles, m-3 s-1 RG ) growth rate, kg m-2 s-1 rN ) nucleation rate, m-3 s-1 S ) average class size, m SMB ) sodium metaborate SPT ) sodium perborate tetrahydrate t ) time, s tcr ) crystallization time, s td ) disruption time, s T ) temperature, K Tin ) input temperature of cooling water, K Tout ) output temperature of cooling water, K Tr ) reference temperature, K u′ ) terminal velocity, m s-1 U ) overall heat transfer coefficient, W m-2 K-1 V ) suspension volume, m3 V ) velocity of cooling water, m s-1 Vl ) volume of the liguid phase, m3 wf ) objective function weighting factor w ) mass fraction, % Greek Letters R ) volume part of liquid phase in suspension δb ) thickness of the boundary layer surrounding the crystal, m δn,i ) element of the Kronecker matrix ) energy dissipation rate, m2 s-3 j ) mean energy dissipation rate, m2 s-3 η ) dynamic viscosity, kg m-1 s-1 ηa ) agglomeration efficiency ηk ) Kolmogoroffov microscale, m ηr ) effectiveness factor θ ) paddle angle, deg λ ) thermal conductivity, W m-1 K-1 λe ) Lagrangian microscale, m λg ) Taylor microscale, m ν ) stoichiometric coefficient F ) cooling water density, kg m-3 Fi ) density of the ith reactant flow, kg m-3 Fl ) density of the liquid phase, kg m-3 Fs ) solid density, kg m-3 Fsusp ) suspension density, kg m-3 σ* ) yield stress, Pa Φ ) objective function ψ ) crystal population density function function, m-3 ωi ) weighting factor for the ith measurement Subscripts c ) concentration exp ) experimental i, j, n ) class index m ) mass w ) wall µ ) time index Literature Cited (1) Van Gelder, D. W. Equilibria in the ternary system sodium metaborate - sodium perborate - water. Recl. TraV. Chim. Pays-Bas 1958, 75, 117. (2) David, R.; Marchal, P.; Klein, J. P.; Villermaux J. Crystallization and precipitation engineering - IV. Kinetic model of adipic acid crystallization. Chem. Eng. Sci. 1991, 46, 1129.
(3) David, R.; Marchal, P.; Klein, J. P.; Villermaux J. Crystallization and precipitation engineering - III. A Discrete formulation of the agglomeration in industrial crystallization from solution. Chem. Eng. Sci. 1991, 46, 205. (4) David, R.; Marchal, P.; Marcant, B. The modeling of agglomeration in industrial crystallization from solution. Chem. Eng. Technol. 1995, 18, 302. (5) David, R.; Bossoutrot, J. M. Crystallization and precipitation engineering - VII, The modeling of sodium perborate tetrahydrate crystallization from solution. Chem. Eng. Sci. 1996, 51, 4939. (6) Nagata, S. Mixing; Principle and Applications; Wiley: New York, 1975. (7) Randolph, A. D.; Larson, M. A. Theory of Particulate Processes; Academic Press: New York, 1971. (8) Marchal, P.; David, R.; Klein, J. P.; Villermaux, J. Crystallization and precipitation engineering - I. An efficient method for solving population balance in crystallization with agglomeration. Chem. Eng. Sci. 1988, 43, 59. (9) Garside, J. The concept of effectiveness factors in crystal growth. Chem. Eng. Sci. 1971, 26, 142. (10) Mersmann, A.; Rennie, F. W.; Design of crystallizers and crystallization processes. In Crystallization Technology Handbook, 1st ed.; M. Dekker: New York, 1995; pp 215-325. (11) Frances, C.; Biscans, B.; Laguerie, C. Some physicochemical data on tetrahydrate sodium perborate in aqueous solutions. J. Chem. Eng. Data 1990, 35, 423. (12) David, R.; Paulaime, A.; Espitalier, F.; Rouleau, L. Modelling of multiple-mechanism agglomeration in a crystallization process. Powder Technol. 2003, 130, 228. (13) Costes, J.; Couderc, J. P. Study by laser Doppler anemometry of turbulent flow - II. Spectral analysis and scales of turbulence. Chem. Eng. Sci. 1988, 43, 2765. (14) Mahouast, M.; David, R.; Cognet, G. Caracterization des champs hydrodinamique et de concentration dans une cuve standard alimentee en continu. Entropie 1987, 133, 7. (15) Higashitani, K.; Yamauchi, K.; Matsuno, Y.; Hosokawa, G. Turbulent coagulation of particles dispersed in a viscous fluid. J. Chem. Eng. Jpn. 1983, 16, 299. (16) David, R.; Espitalier, F.; Cameirao, A. Developments in the understanding and modeling of the agglomeration of suspended crystals in crystallization from solutions. Kona 2003, 21, 40. (17) Liew, T. L.; Barrick, J. P.; Hounslow, M. J. A micro-mechanical model for the rate of aggregation during precipitation from solution. Chem. Eng. Technol. 2003, 26, 282. (18) Mersmann, A.; Braun, B. Agglomeration. In Crystallization Technology Handbook, 2nd ed.; M. Dekker: New York, 2001; pp 135-137. (19) Marcant, B.; David, R. Experimental evidence for a prediction of micromoxing effects in precipitation. AIChE J. 1991, 37, 1698. (20) Bramley, A. S.; Hounslow, M. J.; Ryall, L. R. Aggregation during precipitation from solution: A method for extracting rates from experimental data. J. Colloid Interface Sci. 1996, 183, 155. (21) Livk, I.; Pohar, C.; Ilievski, D. Estimation of batch precipitation kinetics by a simplified differential method. AIChE J. 1999, 45, 1593. (22) Livk, I. Dynamic Modelling of the process of Crystallization (in Slovene). Ph.D. Thesis, University of Ljubljana, Ljubljana, Slovenia, 1997. (23) Chianese, A.; Contaldi, A.; Mazzarota B. Primary nucleation of sodium perborate in aqueous solutions. J. Cryst. Growth 1986, 78, 279. (24) Dugua, J.; Simon, B. Crystallization of sodium perborate from aqueous solutions. J. Cryst. Growth 1978, 44, 265. (25) Livk, I.; Smodisˇ, M.; Golob, J.; Pohar, C. Crystallization kinetics of sodium perborate. Cryst. Res. Technol. 1995, 30 (7), 911. (26) Sohnel, O.; Bravi, M.; Chianese, A.; Mazzarotta, B. Growth kinetics of sodium perborate from batch crystallization. J. Cryst. Growth 1996, 160, 355. (27) Chianese, A. Growth and dissolution of sodium perborate in aqueous solutions by using the RDC technique. J. Cryst. Growth 1988, 91, 39. (28) Chianese, A.; Condo, A. The growth and dissolution of sodium perborate crystals in a fluidized bed crystallizer. J. Cryst. Growth 1989, 97, 375. (29) Frances, C.; Biscans, B.; Laguerie, C. Modelling of a continuous fluidized bed crystallizer. Effects of mixing and segregation on crystal size distribution during the crystallization of tetrahydrate sodium perborate. Chem. Eng. Sci. 1994, 49, 3269.
ReceiVed for reView October 4, 2004 ReVised manuscript receiVed November 8, 2005 Accepted November 30, 2005 IE049039T
