(5) Deposition of metal in the feed line and condensation salt components in the vacuum lines must be preof vented for long-term operation of a molten salt still. Nomenclature
liquid yi
=
average mole fraction of compound i in still-pot
mole fraction of component i in vapor effective relative volatility of component i with respect to LiF, defined by eq 1 =
Hightower, J. R., LTcNeese, L. E., J . Chem. Eng. Data, 17 (3), 342 (1972). Hi htower, J. R., hlcNeese, L. E., U . 8. At. Energy Comm. Rep., 8RNL-4434 (1971). Hightower, J. R., McNeese, L. E., Hannaford, B. A., Cochran, H. D., U . S. At. Energy Comm. Rep., ORNL-4577 (1971). Lindauer, R. B., U . S . At. Energy Comm. Rep., ORNL-TM-2578 (1969). Smith, F. J., Ferris, L. M., Thompson, C. T., U . S . At. Energy Comm. Rep., ORNL-4415 (1969).
CY~-L,F =
literature Cited
Carter, W. L., Lindauer, R. B., McNeese, L. E., C. S . A t . Energy Comm. Rep., ORNL-TM-2213, 1968.
RECEIVED for review January 20, 1972 ACCEPTED January 19, 1973 Research sponsored by the U. S. Atomic Energy Commission under contract with Union Carbide Corporation.
Modeling the Batch Crystallization Process. The Ice-Brine System Jong-Shinn Wey and Joseph Estrin* Chemical Engineering Department, Clarkson College of Technology,Potsdam, New York 13676
The cooling-down suspension batch crystallization process was modeled by considering the population balance and material and energy balance equations. A general size-dependent growth model based upon diffusion-controlled growth mechanisms and a simple power law for the nucleation rate involving the supersaturation and total surface area were included. These were applied to the ice-brine system. The Hulburt and Katz numerical technique involving the moments of the distribution was used to solve the batch conservation equations. The importance of the maximum in the distribution in interpreting size distributions from batch experiments i s discussed. The results reveal that the heat transfer properties of the crystallizer, the level of agitation, the coolant temperature, and the number of seeds sensitively determine the final distribution. Data from batch crystallization experiments are also discussed.
T h e techniques for interpreting crystal size distributions from crystallization processes to ascertain the nature of the controlling mechanisms developed rapidly over the period from 1962 to 1964 (Randolph and Larson, 1962; Hulburt and Katz, 1964). Many experiments have since made use of these techniques using steady mixed suspension crystallization methods. Steady-state experiments provide size distribution data which are readily interpreted to give meaningful results concerning the growth and nucleation mechanisms. Relatively little has been written concerning the analysis of size distributions from batch crystallizer operations (Hulburt and Katz, 1964) which do offer advantages of simpler experimental methods. Further, many industrial suspension crystallizations are carried out batch-wise. This work describes results obtained from mathematically modeling batch crystallizations. Emphasis has been placed upon the system consisting of a suspension of ice crystals in brine; some experimental size distributions for this system are available. The results of the modeling illustrate which experimental parameters are probably important in determining the final size distribution, and could apply reasonably generally to any batch process for which diffusion-controlled growth of the dispersed phase is a reasonably accurate growth model. The 236
Ind. Eng. Chem. Process Des. Develop., Vol. 12,
No. 3, 1973
mode of the distribution is suggested to be a particularly meaningful characteristic in what follows. Other growth mechanisms are implied in the discussion which follows and pertain to the work of Sauter (1972) and of Baliga and Larson (1970). A. Description of the Crystallization Process
Shroff (1969) has demonstrated that the use of the following operating procedure for carrying out experimental batch crystallizations employing the ice-brine system is practical. That is, they may be carried out without ice deposition occurring on the heat transfer surface. Ordinarily this deposition phenomenon prevents the use of a solid heat transfer surface for suspension crystallization studies involving this system. The procedure is based upon the following schedule. ,4 known quantity of salt solution is cooled by transferring heat through malls of a Couette flow crystallizer until nucleation of ice crystals occurs. The apparatus and contents are permitted to warm by stopping coolant flow until the solution is above the freezing point and the crystals begin to melt. This procedure removes any ice deposit adhering to the walls. The coolant flow is initiated and the solution and suspended solids decrease in temperature because of heat transfer to the
refrigerant. When the saturation temperature is attained the run is considered to have been initiated, Le., in the model, this moment is taken to be time zero. The seeds present apparently serve as a preferential surface on which the ice deposits as the level of undercooling increases and deposition on heat transfer surfaces is avoided. Secondary nucleation and ice crystal growth spontaneously occur. Equilibrium is attained when the slurry contents are a t constant temperature and a constant, final distribution of crystal sizes exists. B. Nucleation and Growth Models
A size-dependent growth mechanism is considered here. Small ice crystals are considered to grow due to the convective diffusion process described by Harriott (1967); larger sizes are assumed to grow according to the empirical expression offered by Calderbank and Moo-Young (1961) to describe the heat and mass transfer to a dispersed phase in stirred vessels. These growth descriptions are in reasonable agreement with the work of Brian, et al. (1969), who studied heat and mass transfer from dissolving ice spheres in a stirred tank. The absence of inherent growth resistance effects (at least for the faster growing edges of the disk-shaped crystals) is reasonably well justified (Fernandez and Barduhn, 1967; Wey and Estrin, 1972). The growth model is described in greater detail in the Appendix and is illustrated in Figure 1. The secondary nucleation rate is assumed here to be given by a power law
Bo = KN(AC*) LIp2
(1)
K N and a are the important parameters to be determined as a consequence of the modeling procedure and experiments. Success in showing consistency in the behavior of K N and a over a number of experiments of varying operating parameter value would give credence to the functional form of the expression. Inadequacies of this oversimplified mechanism are suggested later; it is used here as a reasonable first estimate for illustrative reasons. The value of K N used here is estimated from the work of Margolis, Brian, and Sarofim (1971) who measured nucleation rates of ice produced in a well-stirred continuous crystallizer. From their experimental data, one knows Bo,4C*, and ~2 so that the value of K N can be estimated from
where ( K N ) is ~ the value of K N when a = 1. The value of ( K N )estimated ~ from their work is 1.8 X lo2g-l. The nucleation model used here is consistent with some evidence described in the literature (Sadek, 1966; Clontz and McCabe, 1971). We visualize the nucleation mechanism for the ice-brine suspension system as a microbreakage mechanism. I t is further assumed that these breakage events do not disturb the growth rates which occur by a heat and mass transfer process. We assume that some nuclei generated by the fracture process are near to the critical size but their growth rates are still not appreciably altered from those implied in Figure 1. This point has been theoretically demonstrated elsewhere when interfacial kinetics are not significant (Wey and Estrin, 1972). The estimate of a zero initial size associated with the nuclei is permissible due to the rapid growth rates assumed for small size crystals (hence their time at any size near to zero is very short as it very rapidly grows to larger sizes relative to the growth rates of larger crystals). The effect of surface energy upon growth for this system is referred to later when describing some available experimental data.
I\
E
s8
I
-
.06
.04
-
.02
-
00
4 ,02
,04
CRYSTAL
.06
.08
.IO
.I2 4
DIAMETER (cm)
Figure 1 . Size dependent growth rate coefficient as a function of crystal size (3% NaCl solution)
K N is assumed to depend upon the level of agitation or specific power dissipation P / V . The dependence of K N on P/V may be determined experimentally by changing the rotational velocity of the rotor and/or the annular gap of the Couette flow crystallizer. The dependence of K N upon P / V is not discussed further here. K N would display temperature dependence; the range of temperatures permitted by the brine-ice system is limited and is not considered to be a major independent experimental variable. C. Characteristier of the Maximum of the Number Distribution
The maximum in the number density distribution, should it exist, is a very prominent feature and is useful in revealing what mechanisms are involved in determining the distribution. The location of the peak is related to growth in the following way. Recognizing that n is a function of the two variables L and t, the chain rule may be written in the following form dn
dn dt
dt
dndL bL dt
The derivative dL/dt implies a relationship between L and t such that this derivative exists. Now the population balance holds in general dn - + - =dGn o at bL On combining the two equations, the following is obtained
dn _ _ - an(dL _ _ - G) - dG (3) dt dL dt We now identify the derivative dL/dt with the velocity of the maximum along the L coordinate and note that the function, n, a t the maximum has the property ~
Thus, the magnitude of according to
bn/bL = 0 (4) a t the maximum changes with time
TI.
where Lp is the value of L a t which the maximum is located. Ind. Eng. Chem. Process Des. Develop., Vol. 12, No. 3, 1 9 7 3
237
To obtain the motion of the peak as dL,/dt, we carry out an argument similar to that above and consider the derivative bn/bL to be a function of two variables d ( an) - b2ndL -% bL2 dt
dt
+
”(”)
bt bL
The order of differentiation of the right-hand term may be interchanged and the time derivative of bn/bL, on the left, is zero at the peak; on replacing bn/bt by -(dGn/dL) the following is ultimately obtained
Expression 5 reveals how the magnitude of n a t L, changes with time. This rate of magnitude change depends upon the grom-th mechanism and particularly upon the nature of the dependence ot G upon L. For the ice-brine system, the growth function assumed here is a nonincreasing one as L increases so that (5) implies n a t L, will not decrease with time. For size-independent growth the peak will remain unchanged in value. The same kinematic properties are true for the minima of the distribution. Thus, if experimentally one can follow the maximum and minimum values as the crystallization occurs, growth characteristics of the dispersed phase may be inferred. Expression 6 reveals how the location of the peak changes with time. For the growth function shown in Figure 1, the second term on the right is negative and dL,/dt is less than G a t L,. Depending upon the value of b2n/bL2a t L, one peak may tend to coalesce with another or one may conceivably travel in the direction of smaller sizes. For size-independent growth, the peak remains distinct and moves toward the larger sizes a t a velocity equal to the growth rate, G. For the growth function employed here, minima move to the right more rapidly than maxima (this is demonstrated by the skewness of the data of Wey (1972) in Figure 5 which shows the distribution of sizes of ice crystallizing from brine). We can continue the above discussion by considering conditions under which a peak is generated due to nucleation a t size Lo. Generation of a maximum (or minimum) refers to conditions a t which bn/bL = 0 and dL,/dt > 0, as L + LO from the right. The population balance is multiplied through by G, the derivatives of products are expanded or reformed, and the resulting equation is solved for bn/bL to yield
A further requirement would be that dL,/dt > 0 (11) but this can not be usefully introduced into this presentation in any straightforward way. On examination of eq 10 it is seen that for the growth mechanism depicted here, in Figure 1, for the ice-brine system, the term
is negative and is proportional to (AC*)@so that for the type of operation considered here involving a build up of supersaturation from a zero value this term is initially zero. On the other hand, when AC* + 0, the second term [(CY- P)/AC*](dAC*/ dt) is large, assuming (CY- p) > 0; in general the first term d In fiz/dt is nonnegative as is evident from physical considerations (that is, dissolution of crystals does not occur). An equation for the second moment is obtained directly from the population balance equation by multiplying through by L2 and integrating from LOto infinity. If LOis considered small, the moment equation becomes approximately
Now as the maximum value of G is approached as L LnG dL
< (AC*)@K(Lo)pl
bn I
dBO dt
G2-
dLl~+~o
-
Bo and the above equa-
Bo dG1 G bt !L.-L~
+ Bo bL
L+Lo
]
(8)
.A necessary condition for a maximum (or minimum) in n just
to the right of L = Lo is that the right-hand side of eq 8 be equal t o zero. When the growth and nucleation mechanisms are G
K(L)(AC*)@ /3 2 1
(9)
Bo = K ~ p z ( A c * ) ~
(1)
=
then eq 8, as the necessary condition, becomes d In PZ dt
___ + ( a - PI d IndtAC* + (Ac*1@dK(L)l dL IL+L~
=
0
~
(10) 238
Ind. Eng. Chem. Process Des. Develop., Vol. 12, No. 3, 1 9 7 3
Lo (13)
Furthermore, if the growth of a crystal of size L + LOmay be estimated from simple diffusion considerations to be (AC*)@K(L)
(k/L)(AC*)@
then
The left-hand side of eq 10 is approximated by
+
2(hC*)@K(Lo)!!! L(2
- a)
(CY
d In AC* dt ~
-
and for distributions which are not extremely spiked 2(Pl/P2)
Now as L -+ Lo from the right, Gn tion becomes
-
0. This conclusion holds only for the type of operation considered here. When the operation of the crystallizer is modified so that AC* is a t its maximum value a t t = 0 and subsequently decreases, monotonically, after being seeded, the above discussion implies that no maximum would be generated so long as the diffusion growth model is as assumed, no interfacial kinetics and no surface energy effects are involved, CY 2 p, and the simple nucleation model provided by eq 1 is applicable. In fact this kind of experiment has been made and the resulting distribution displayed the same peaked characteristic as the standard operation displayed. The most serious shortcoming among the listed assumptions above is probably the simplicity of the nucleation model (eq 1). This shortcoming
CRYSTAL
SIZE.
L (cm)
Figure 2. Particle size distribution data from a typical batch run (aqueous MgSOd.7HzO system (data of Weed, 1972)) CRYSTAL
is evident from the appearance of the experimental data of Schneider (1970) which is discussed later. Another example of experimental data illustrates the above discussion concerning the presence or absence of the mode of the distribution. This relates to the distribution of sizes from continuous MgSO4.7H20 crystallizations and from isothermal batch experiments. The continuous operation yielded distributions which were decidedly concave upward when plotted as In n us. L (Sauter, 1972). This type of distribution would result from growth rates which increase with size; this was believed to be the case for this system; the smaller, less blemished crystals were subject to strong influence due to interfacial kinetics which retard growth while the larger crystals, with highly blemished surfaces and consequently many growth centers, were less subject to this influence (Sauter, 1972). Furthermore, analysis of the data from both continuous and batch runs indicated that quite closely a = p. Consideration of this information, viz. a = p
dK(L)/dL
>0
together with the eq 10 suggests that no peak should be observed in the distributions resulting from the batch experiments. This appeared to be the case as is seen in a typical distribution of sizes provided in Figure 2 (Weed, 1972). This kind of analysis may also be applied to the distributions of sizes of potassium sulfate crystals as obtained by Baliga and Larson (1970). From their data and discussion the factors (a - p) and d In AC*/dt are both negative and [ d K ( L ) / d L ] ( L + ~isopositive so that all terms in eq 10 are positive and the equality can not be satisfied. The implication is that their distributions should permit no mode a t L + Lo as was indeed so. D. Effect of the Nucleation Model Parameters
Figure 3 shows the effects of the parameters ( K N )and ~ a upon the final size distributions. As ( K N )increases, ~ the peak shifts to the small size region and has a higher value. This is because ( K N ) ,determines the magnitude of nucleation rate and enters into the competition between nucleation and growth processes for the solute available in solution. A higher ~ that nucleation predominates, ie., more value of ( K N )means
DIAMETER (cm)
Figure 3. Particle size distribution curves, effect of the secondary nucleation power law model: Cb(0) = 3 g/100 ml, po(0) = 100 no./ml, Tb(O) = -1.78", T, = -5.5", UA = 64 cal/sec "C,P / V = erg/cma sec, ( K N ) l = 1.8 X 1 O2 9 - l in Figures 4,7,9,and 1 1
lo4
crystals exist in the small size region for given total precipitated solids. The parameter a shows the relative dependence of nucleation and growth processes on supersaturation. When nucleation has a greater subcooling dependence, Le., a > 1, the peak shifts to the large size region and fewer crystals exist in the small size region. The driving force AC* is very small g/ml), so that a larger value of a yields a lower nucleation rate for given AC*, especially as Ac* -+ 0, as occurs near the end of the batch operation.
E. Size Distribution
as a Function of l i m e
Figure 4 shows the behavior of the particle size distribution variation with time. The behavior of the curves strongly depends upon the value of a, which substantially determines the outcome of the competition of nucleation and growth processes for available supersaturation. For a = 1 a second peak occurs a t early time (t = 1 min). Later, the peak a t the larger size disappears and the peak a t the smaller size remains in the smaller size region and nIL,continues to increase in value. For a = 2 only one peak occurs which persists in moving to the larger size region. The distributions generated over a period of one run were obtained using the size-dependent growth function illustrated in Figure 1. For a equal to 2 the peak is seen to move to the right and persistently to larger values until equilibrium is attained. The increase in n a t L, is small for the larger sizes (later times) as is expected, since ?)G/BL+ 0 for these sizes and dG/bL directly influences dn/dt a t L,. In general these peak characteristics are in accord with the equations for dLp/dt and dn/dt given as eq 5 and 6. The a = 1 distributions display a different behavior. The magnitude of n a t the maximum a t 1 min is less in value than the corresponding seeds value. This is counter to what is expected according to the discussion above and can be attributed only to the inadequacy of the computational method. After 1 min the distribution maximum increases in value monotoniInd. Eng. Chem. Process Des. Develop., Vol. 12, No. 3, 1973
239
m
2 MINUTE 0 5 MINUTE
t I O MINUTE 0 STEADY STATE I> I5 YII1.1
C
CRYSTAL
DIAMETER (cm)
Figure 4. Particle size distribution curves, time dependence (conditionsof Figure 3)
0 8lEADY STATE
5
Figure 6. Particle size distribution data from a typical batch run (ice-brine system (data of Wey, 1972))
Figures 5 and 6 show some short-time distributions obtained from the Couette flow crystallizer for the ice-brine system (Wey, 1972). The suspended crystals were photographed as they passed through an external cell. The Zeiss Particle Counter was employed to determine distributions from the photographs. The apparatus and procedures used will be described in a future paper. It is apparent from Figures 5 and 6 that the mode does move in the direction of large size and takes on greater values as time increases. Qualitatively, the appearance of mode behavior suggests an CY value in the vicinity of 2, say, and certainly greater than unity. The leveling off of the modes a t about L = 0.06 cm, which is especially apparent in Figure 6, is consistent with the growth model described in Figure 1. F. Effect of Refrigerant Temperature on Size Distribution
CRYSTAL DIAMETER (cm)
Figure 5. Particle size distribution data from a typical batch run (ice-brine system (data of Wey, 1972))
cally with time. The appearance of the second peak (left-hand peak) suggests the importance of nucleation in determining the distribution. The right-hand peak originated from the seeds and ultimately is to be insignificant. Evidently the numerical method employed is not applicable to a bimodal distribution with good accuracy as is evidenced by the decreasing first peak value with increasing time, which is not consistent with the growth model. However, since the distribution of seeds and short-time distribution behavior do not influence the long-time distribution, the comparison of equilibrium distributions is permissible. Apparently the mode does not shift to the right for CY = 1 as it does for CY = 2. 240 Ind.
Eng. Chem. Process Des. Develop., Vol.
12, No. 3, 1973
Figure 7 shows the effect of refrigerant temperature on the final equilibrium distribution for fixed initial brine concentration. The size distribution curves shift upward when a lower refrigerant temperature is used. The change in peak value due to a change in T,is not sensitive when CY is 1 ; it is sensitive when a is 2 and presumably larger. It appears then that the refrigerant temperature may be an important experimental parameter for use in discriminating CY values. In the study of secondary nucleation phenomena the refrigerant temperature should be maintained high enough (or AC* small enough) so that the contribution of the homogeneous nucleation process can be neglected. In general, homogeneous nucleation effects are not significant for the ice-brine system if the subcooling of the solution is less than 3' (Sherwood, et al., 1966). The deposition of ice on the heat transfer surface limits the extent to which T,may be lowered below the saturation temperature of the feed brine. 0. Effect of Seeds on Size Distribution
The size distribution of seeds can be transformed into moments which serve as initial conditions for the batch conservation equations in moment form. Figure 8 shows several size distributions of seeds. Seeds C have the same mass as seeds A,
iI
I
I I
I
2x104 0
1
I I I 1
'
,04
I
.OS
.I2
I
.I6
.20
.24
CRYSTAL DIAMETER (cm)
Figure 7. Particle size distribution curves, effect of coolant temperature (conditions of Figure 3)
CRYSTAL
w
B C D
1000 35 35 1.25
PI
20
0.7 2.5 0.09
DIAMETER
(cml
Figure 9. Particle size distribution curves, effect of seeds (conditionsof Figure 3)
DIAMETER (cm)
Figure 8. Particle size distribution of seeds
A
CRYSTAL
P2
0.5
0.018 0.2 0.0067
Pa
0.0154 0.00054 0.0154 0.00054
but the former have a smaller number and larger particles. Seeds B have the same number of particles as seeds C, but the former have a smaller particle size distribution. Seeds D have the same mass as seeds B, but the former has smaller number and larger particles. Figures 9 and 10 show the effects of seeds on the final (equilibrium) particle size distributions for different values of the constants in the secondary nucleation expression. In the case of low ( K N ) the ~ , final particle size distributions are not sensitive to the number or size distribution of the
CRYSTAL
DIAMETER
(em)
Figure 10. Particle size distribution curves, effect of seeds (conditionsof Figure 3)
seeds so long as the number of seeds is not too great (say, less than 50/ml as in seeds types B, C, and D). When the number of seeds is large, the final particle size distribution becomes very different. These results are exhibited in Figure 9. An explanation is that the seeds grow very rapidly during the initial period of time because little or no secondary nucleation competes for the supersaturation which is large. The rapid growth of the seeds increases the total exposed crystal surface which subsequently gives rise to an active secondary nucleation during this early period and it is these nuclei which serve to determine the final particle size distributions. When the number of seeds is great, some of the supersaturation is consumed by the fast growing seeds during the initial period Ind. Eng. Chem. Process Der. Develop., Vol. 12, No. 3, 1973
241
in Figure 10. This comes about because, in the case of large (&)I, the nucleation rate remains relatively high throughout the run and the influence of newly generated crystals, rather than the seeds, predominates in determining the distribution. H. Effect of Heat Transfer Characteristics on Size Distribution
CRYSTAL
DIAMETER (cm)
Figure 1 1 . Particle size distribution curves, effect of UA (conditions of Figure 3)
Figure 11 shows the effect of heat transfer characteristics through the wall of the Couette flow crystallizer ( U A ) on the final particle size distribution. For a = 1, nucleation has the same supersaturation dependence as growth and there is no effect of U A as can be shown directly from the equations. However, when nucleation has a greater supersaturation dependence, Le., a > l, the heat transfer characteristics ( U A ) do bring about changes in final size distribution. A smaller U A value yields a flatter final size distribution. Distributions for CY = 3 as well as for CY = 2 are displayed to show the sensitivity of the peak magnitude to the value of ( U A ) . Experimentally one could determine the parameter CY through modeling procedures by changing the total heat transfer area or the overall heat transfer coefficient U.A reduction in U could be obtained by partially inhibiting the heat transfer by the addition of a thermally nonconducting coating on the heat transfer surface. 1. consideration of Experimental Data
The experimental batch crystallization data of Schneider (1970) are shown in Figure 12. The crystallization experiments involved the ice-brine system and employed isobutylene as a primary, direct contact refrigerant. These experiments involved a typical stirred tank configuration with the liquid refrigerant introduced below the stirrer. Crystal size analyses were carried out using a modified Royco Liquidborne particle size analyzer. -4description of the apparatus is given in the reports by Schneider, et al. (1968; 1971). These data are very interesting in that they include sizes well below 100 p and display a bimodal distribution a t later times. Schneider attributes the smaller size distributions to an agglomeration mechanism. An alternative description is suggested here for the nature of these distributions. When we include death and nucleation terms into the crystal population balance, this becomes dn - + - dGn at b~
=B-D
and the appropriate expansion of n (t,L) according to the chain rule is
CRYSTAL SIZE, L (microns)
Figure 12. Particle size distributions determined from batch crystallization experiments (ice-brine system (data of Schneider, 1970))
so that the nucleation rate is relatively small in the case of ~ . one compares the equilibrium size distribution low ( K N ) As from seeds A with those from B, C, and D, it is obvious that virtually no nucleation occurs in the seeds A experiment. This is especially true for the CY = 2 distribution. For CY = 1, the growth still appears to predominate but there is some contribution of nucleation though it is still too small to change the nature of the final size distribution very significantly. In the case of large ( K N ) , the , final size distributions are not determined by the nature of seeds. These results are given 242 Ind. Eng. Chem. Process Des. Develop., Vol. 12, NO. 3, 1973
The death term describes the sudden disappearance of crystals of size L; the nucleation term describes the sudden appearance of crystals of size L. We may identify dL/dt with dL,/dt but we need not restrict eq 16 to apply only to size L = L,; that is dn/dt is applied to the entire distribution to the left of, say 150 p . If we accept the right-hand term as
-(-
dn dL, - G ) bL dt it may certainly be argued that it is small over the range 0 < L 5 100 p in which dn/dL is zero a t several locations; the three remaining terms then determine the uniform decrease in the value of n which occur with time over this size range. The third term on the right-hand side may be accepted as contributing to the decrease in value of n when bG/dL > 0
This may well be the case when particles are actually dissolving rather than growing. Certainly the growth behavior due to diffusion alone (with no interfacial energy influences) is not an acceptable mechanism here as then dG/bL < 0 and is large in magnitude in this size range. This will be argued later. It is suggested here that the death and growth (i.e., positive or negative growth) are dependent upon the supersaturation and the thickness of the ice disk-like crystals. The characteristic L is not ordinarily considered to be descriptive of this dimension, yet it is the thickness times the disk diameter product which is the surface area a t which linear growth velocity occurs most rapidly. Also from the arguments of Fernandez and Barduhn (1967) one anticipates that this surface is the one which will display surface energy effects due to curvature when the thickness is sufficiently small, and further, when the thickness is small enough one anticipates that the crystal dissolves rather than grows. Therefore, in this range of small sizes, it might well be that the usual size L and the effective thickness are not correlated strongly so that even a relatively “large” crystal, say 100 p , may show rapid dissolution characteristics because of a very thin platelet geometry. That is, the very thin platelets dissolve very rapidly and may be considered to disappear suddenly. These very thin ones are described in the above balance equation by the death term, D. When the thickness and the usual size dimension, L , are correlated only weakly, the apparent death of a larger crystal (say, L = 100 p ) could occur as frequently as a smaller crystal (say, L = 30 p ) and the death function would appear to be independent of L. When D is the dominant term on the right of eq 16, the uniform decrease in value for n(L), 0 5 L < -100 p , observed from Schneider’s data, is expected. This assumes implicitly that B is also approximately uniform over 0 to 100 p . We may note that the distributions over the sizes larger than 150 p or so (Figure 12) appear to display the momentary occurrence of a zero second derivative a t about 250 I.( a t times between those corresponding to curves 3 and 4.This momentary event then apparently gives rise to the minimum in the distribution (curve 4) appearing between the modes. If we assume that the birth function, B, and the death function, D, are not applicable in these regions of larger L, then eq 5 and 6 hold and describe the behavior of both the maximum and minimum in the distribution. The right-hand peak retains its value over a period of time since it was generated a t the second derivative which occurred a t about the same value for n. From (5), size independent growth, dG/dL = 0, should give rise to this behavior, Le., (dn/dt)(L, = 0. Size-independent growth (or nearly so) is anticipated in this size range (Figure 1). A similar examination of the minimum shows that dn/dt < 0 and therefore dG/dL > 0. This is not the growth behavior displayed in Figure 1, which ignores surface energy effects. Mass transport (growth or dissolution) which increases with size is representative of growth rate-size behavior of growing spheres subject to growth retarding surface energy influences (Wey and Estrin, 1972). This suggests that the minimum lies in a size region in which the crystals either dissolve or grow very slowly due to capillarity effects. Consistent with this, then, even smaller crystals do dissolve which was the claim made earlier when describing those in the smaller size region. Equation 6 may be rewritten in the form it had before rearrangement
From this it is seen that a necessary condition for a second derivative to be equal to zero momentarily a t L, is b2G/bL2 = 0 Thus, the growth characteristics of the crystals in the vicinity of the location a t which (b2n/bL2)/L-,L, = 0 are described by
dG/bL
>0
bG/dL
-
0
as concluded from earlier discussion and bZG/bL2 = 0 Such behavior could be exhibited by a particle in a region of sizes in which surface energy influences and possibly turbulence influences upon mass transport are important. From eq 6, if G dominates the right-hand side, the above characteristics are consistent with the observation that (dLp/dt)max> (d&/dt),,, and especially when G passes through zero, the minimum and maximum are in fact distinct (which is the case in Schneider’s data). Conclusions
The preceding discusses the nature of the distribution of sizes which may be obtained from the batch operation of a laboratory crystallizer in which ice crystallizes from a brine solution. Here the nucleation expression is assumed to be adequately represented by
Bo = K N ( A C * ) ~ ~ ~ Results of preliminary experiments however suggest this description may be over simplified. Growth is based on a diffusion-controlling mechanism so that AC* is involved to the power unity. This growth mechanism appears to be substantially verified by experiments with single ice crystals which have indicated the absence of significant growth-retarding mechanisms (at least in directions lying in the basal plane (Fernandez and Barduhn, 1967)). The location and magnitude of the mode or peak in the distribution is a clearly observed characteristic of the distribution and as such its behavior, when the operating parameters are varied, may provide very useful information concerning the growth and nucleation mechanisms. The value of the mode in the final size distribution is sensitive to the value of K N in the nucleation expression. As K N is considered to be dependent upon the power being dissipated in the agitated slurry, the peak in the distribution potentially provides a measure of the dependence of secondary nucleation upon the specific power dissipation, P / V . The motion of the peak, ie., Lp(t),is sensitive to the growth function. This should be revealed clearly when number distribution vs. time data are obtained experimentally. The data of Schneider (1970) are discussed from this point of view. These data are of interest in that they involve small sizes for the ice-brine system. The number of seeds may influence the final size distribution markedly; the distribution of sizes of the seed crystals has virtually no influence. However when the seeds number is small, this influence of number does not exist. This complex dependence of the distribution upon initial seeds suggests that batch crystallization nucleation experiments should be carried out a t low concentrations of initial seeds to avoid final distributions which may be difficult to interpret. Acknowledgment
This work was financially supported by the Office of Saline Water (Grant No. 14-30-2573). The authors are indebted to Ind. Eng. Chem. Process Der. Develop., Vol. 12, No. 3, 1973
243
Dr. G. R. Schneider for making available to them a computer printout of the reduced data from his automatic particle size analyzer. Appendix
The equations used for modeling the batch crystallizer are briefly described in this Appendix. A more detailed presentation can be found in reference (Wey, 1970).
A. Batch Conservation Equations. (a) Population Balance. bn
+
b
Equation A.1 is derived under the following assumptions: (1) the system is perfectly mixed, (2) particle and solution densities are constant and nearly equal, and (3)all new particles are generated a t a small characteristic length, Lo, which is assumed to be zero compared to practically measurable sizes.
cb =
ci/(1
- Ks
La
(A.2)
(c) Energy Balance. (PL
- Cb)eCp, f (1 - @)PcCp,l(dTb/dt)+ - Cb)hR(de/dt) - &(dCb/dt) = ( U A / V ) ( T C- Tb) 64-31
where e =
1
- J,
KsLan(L, t ) dL
(A.4)
Equation A.3 is derived by applying the first law of thermodynamics to the closed batch system. (d) Boundary Conditions. Equation A.l,the population balance, requires two conditions: (1) an initial condition which states the distribution of sizes a t time t = 0, n(L,O) = ni(L); and (2) the value of the distribution function for zero size as a function of time, n(0,f) = no(t). ni refers to the seeds introduced to initiate the process or to the “initially bred” secondary nuclei which derive from the seeds. From the definitions of population density and growth rate one can relate no((t)to the nucleation rate Bo no = Bo/Go
64.5)
Size “zero” refers to the size of the nuclei which are generated by some mechanism involving an existing crystal surface and which are capable of further growth. Its true size is assumed to be small relative to the bulk crystals in the slurry and in the mathematics is considered to be zero. The initial conditions associated with the material and energy balances are Cb(0)
=
=
1-
La
KsLani dL
Ind. Eng. Chem. Process Des. Develop., Vol. 12,
KAC*
(A.6)
=
Cs(Tb)
-
(AS)
cb
The mass and heat transfer coefficients, in general, vary with particle size. For small particles (NR, < 1) K , and h are described by (Harriott, 1967)
+ 0 . 5 2 N ~ , ’ / ~ ( N por,
Nse)l13
(A.9)
A correlation for heat and mass transfer coefficients in which transfer is affected by turbulence in the surrounding fluid is given by Calderbank and Moo-Young (1961) based on the experimental data of many investigators. The expressions of Calderbank and Moo-Young are
K,
O.l3[(P/V)p/p~~]”~(Ns,)-*/~ (A.10)
=
h = 0.13[(P/V)p/pLZ]1’4CpLpL(NP,)
-2’8
(All)
C. Moment Equations. The conservation equations are difficult to solve since they contain simultaneous partial differential equations and nonlinear boundary conditions. One solution method involves the introduction of moments for size distribution, n. The ith moment of n(L,t) is designated pf and is defined pt =
la
Lin(L, t) dL
(A.12)
Since the growth rate of ice crystals as described here is size dependent, one can fit the growth coefficient, K , into the polynomial form w K = AnLn (A.13)
2
-P
I n this expression the A , depend upon the physical properties and operating parameters and .M and P are positive integers. The function K(L) which involves equations A.7, A.9, A.10, and A . l l may be combined to give a nonincreasing function of L and is illustrated in Figure 1. Inserting the expression for G, the population balance can be written in terms of moments of n(L,t) .M
dpi/dt
=
iAC*
dpo/dt
=
-P
Anpn+t-i i # 0
Bo = KN(AC*)=W
(A.14)
(h.15)
eters are to be determined by the modeling process. The material and energy balances can also be written in moment form
where Ti and Ci are the known initial temperature and concentration, respectively, and in the type of experimental run conducted by Shroff (1969) these are equilibrium values. 244
=
Equation A.15 contains the boundary condition a t L = 0. It is the secondary nucleation rate expression whose param-
ci
Tb(0) = Ti e(Oj
Ac*
N N or~ N g h = 2
L8n dL)
Equation A.2 is obtained from the conservation of salt species in the ice-brine system in which the salt is not incorporated into the ice crystals. [(PL
G where
(nG) = 0
(b) Material Balance.
B. Growth Mechanisms. Basic to the growth of a n ice particle suspended in brine solution are the simultaneous processes diffusion of salt away from the ice-brine interface and diffusion of heat away from the crystal interface. The resistance to growth associated with the intrinsic growth process is so small as to be unimportant (Fernandez and Barduhn, 1967; Wey and Estrin, 1972) and is neglected here. The growth expression which described the two basic steps is given as
No. 3, 1973
cb =
dTb dt
J-
=
Ci/(I - Kspaj
dp3 XRPLK~ dt
UA +(To - Tb) V
(-1.16)
(h.17)
a
size distributions are obtained when seven moments are used. This does not occur when only four moments are used, and strongly suggests that the Laguerre series expansion for n is unstable and a judiciously chosen truncation point is necessary. For this reason, when comparing experiment with theory, i t is apparently more accurate to compare moments rather than size distributions, per se. This characteristic of the Laguerre series expansion method is discussed elsewhere (Randolph and Larson, 1971).
z 0
Nomenclature
c
Ea
z Y
5a
A Bo
g
C
= heat transfer surface area, cm2 = nucleation rate at size LO,no./sec ml
salt concentration, g/ml specific heat of the crystals, cal/"C g CPL - specific heat of solution, cal/"C g CPW = specific heat of water, cal/"C g (dC/dT),, = freezing point depression constant, g/"C ml D = death rate, no./sec ml D , = mass diffusivity, cm2/sec G = growth rate, cm/sec h = heat transfer coefficient, cal/sec cm2 "C K = growth coefficient given in eq A.6, cm4/g sec K , = mass transfer coefficient, cm/sec K N = proportionality constant in the secondary nucleation expression Ks = shape factor, dimensionless KT = thermal conductivity, cal/sec cm "C L = characteristic length of a crystal, cm Lo = characteristic length at which nucleation occurs, cm N N ~= Nusselt number, ZhR/KT Npr = Prandtl number, C,L~/KT N R , = Reynolds number, 2Rup~/p N s c = Schmidt number, p / p ~ D ~ N s ~= Sherwood number, 2K,R/D, n = population density distribution, no./cm ml P = power dissipation, erg/sec R = particle radius, cm T = temperature, "C T , = coolant temperature, "C t = time, sec C = overall heat transfer coefficient, cal/sec em2 "C V = working volume of the crystallizer, ml u = terminal velocity of falling particles, cm/sec
Cpc
CRYSTAL
DIAMETER
(cm)
Figure 13. Particle size distribution curves, effect of the number of moments
where J
= (PL
-
Cb)Cp,
+ Ks!-Q[PcC~, -
(PL
- C,)Cp,I
(X.18)
The initial condition moments are derived directly from the seeds' distribution, nl. D. Description of Numerical Work. The moment equations associated with the conservation equations used for modeling the batch crystallization process can not be closed off a t any value of i when a general size-dependent growth expression is used. Hulburt and Katz (1964) have suggested a procedure for estimation of the unknown distribution from a few of its leading moments by a suitably truncated expansion of the distribution in terms of the r distribution and associated Laguerre functions. The coefficients of the Laguerre function expansion are related to the moments of the distribution in a simple manner. Thus, the moments, which appear in the moment equations but can not be closed off, can be expressed in terms of the leading moments of n and the solution of moment equations would then yield values of as many leading moments as may be of interest. I n this study only the first four moments ( p o - p 3 ) are used t o solve the sequence of moment equations since they are the "basic" ones which are directly involved in the moment equations. The numerical results show that the first four moments obtained from the moment equations with truncation a t pa are insignificantly different from those with truncation a t p6. These moments represent end-of-the-run distributions and therefore the differences display maximum "error" accumulation. Presumably, then, the first four moments give reasonable results when used in the above procedure. This computational method is only a n approximate one and weakness is displayed in some earlier discussions. T o qualitatively illustrate the results obtained here the distributions have been graphed using the method of reconstructing the distribution from its moments as described by Hulburt and Katz (1964). The reconstruction of the distribution introduces error especially when the distribution differs significantly from the applied weighting function. Figure 13 shows the effects of the number of moments applied and thus of the number of terms in the expansion used in the calculations upon the final computed particle size distributions. Wavey
=
=
GREEKLETTERS cx = exponent constant associated with supersaturation in secondary nucleation expression p = exponent constant associated with supersaturation in growth expression XR = heac of solution, cal/g pc = density of crystals, g/ml PL = density of solution, g/ml E = fractional volume occupied by solution in the crystallizer, dimensionless p = viscosity, P p t = i t h moment of distribution function AC* = c, - C b , g/ml SUBSCRIPTS b = bulk condition i = initial condition s = thermodynamic equilibrium condition Literature Cited
Baliga, J. B., Larson, 31. A., Paper presented at 63rd AIChE Annual Meeting, Chicago, Ill., Nov 1970. Brian, P. L. T., Hales, H. B., AIChE J., 15, 419 (1963). Brian. P. L. T.. Hales, H. B., Sherwood, T. K., AZChE J., 15, 727'(1969).
'
Calderbank, P. H., Moo-Young, M. B., Chem. Eng. Sci., 16, 39 (1961).
Clontz, N. A., McCabe, W. L., Chem. Eng. Progr. S y m p . Ser., 67, No. 110 (1971).
Fernandez, R., Barduhn, A. J., Desalination, 3, 330 (1967). Harriott, P., AZChE J.,13, 7 5 5 (1967). Hulburt, H. M., Katz, S., Chem. Eng. Sci., 19, 555 (1964). Ind. Eng. Chem. Process Des. Develop,, Vol. 12, No. 3, 1973
245
hfargolis, G., Brian, P. L. T., Sarofim, A. F., Ind. Eng. Chem., Fundam., 10, 439 (1971). Randolph, A. D., Larson, M. A., AIChE J.,8, 639 (1962). Randolph, A. D., Larson, M. A., "Theory of Particulate Processes," Academic Press, Kew York, N. Y., 1971, pp 58-61. Sadek, S. E., Sc.D. Thesis, Massachusetts Institute of Technology, 1966. Sauter, W. A., Ph.D. Thesis, Clarkson College of Technology, .?,"O
IYiZ.
Schneider, G. R., Newton, P. R., Sheehan, D. F., U . S . Off. Saline Water, Res. Develop. Prog. Rep., NO.R-7717 (1968). Schneider, G. R U.S . Off. Saline Water, Res. Develop. Prog. Rep., No. R-7i79-6 (1970).
Schneider, G. R., U . S. O j . Saline Water, Res. Develop. Prog. Rep., No. R-7979-10 (1971). Sherwood, T. K., Brian, P. L. T., Sarofim, A. F., Smith, K. A., M.I.T. Desalination Research Laboratory, Report No. 295-9 (1Qfifi) ~-~--,. Shroff, H. R., M.S. Thesis, Clarkson College of Technology, 1969. Weed, D. R., M.S. Thesis, Clarkson College of Technology, 1972. Wey, J. S., M.S. Thesis, Clarkson College of Technology, 1970. Wey, J. S.,Estrin, J., AIChE Symp. Ser., 68, No. 121 (1972). Wey, J. S.,Ph.D. Thesis, Clarkson College of Technology, 1972. RECEIVED for review March 16, 1972 ACCEPTEDApril 2, 1973
A Least-Squares Solution to Mass Balance around a Chemical Reactor A. K. S. Murthy Allied Chemical Corporation, Morristown, N . J. 07960
A simple procedure to adjust the measured flow rates of the various species entering and leaving a chemical reactor so as to satisfy the element conservation equations is developed by using Lagrange multipliers and linear algebra. Since this procedure requires the solution of only 1 simultaneous linear equations, where 1 is at most the number of chemical elements in the system, it may even b e used for manual process calculations. A numerical example is given to illustrate the computational procedure.
I n experimental studies of chemical reactors at steady state, due to experimental errors, often the measured flow rates of the various species entering and leaving the reactor do not satisfy the law of conservation of elements. It is a common practice in manual process computations to attribute all the error to a few suspected measurements and adjust only these flows. When care is taken to eliminate all systematic errors, a more desirable procedure is to assume that all measurements are subject to random experimental errors and adjust the raw data in such a may that the required adjustments are minimum in some sense. This can be achieved in several ways; the methods differ in the objective function minimized and the technique used to minimize the objective function. In discussing computer control of chemical processes, Kuehn and Davidson (1961) have shown how the method of least squares can be employed to adjust raw data to satisfy material and energy balances. Vaclavek (1969) has considered the problem of material balance around a chemical reactor specifically, but his solution requires the inversion of a large matrix. Since element balances are linear constraints on the flow rates, linear programming can also be used to achieve material balance if statistical validity of the objective function is not important. The purpose of this communication is to show how a simple procedure to obtain mass balance around a chemical reactor can be developed by using the weighted least-squares procedure.
contained in a molecule of the ith species. If xi and yi are, respectively, the moles of the i t h species entering and leaving the system in unit time, then the fact that the amount of the j t h element must be conserved corresponds to N
N
etjzi = i=l
C eijvt i=l
or A'
- yi)
eij(zr i=l
=
0
j = 1,2,.. . , L
(1)
If 5 , and g i are the measured flow rates of the ith species entering and leaving the system, respectively, and ut and vi are the weights associated with the measurements, the problem is to minimize the objective function x S(z*,vd = %(Et - l d 2 Vi(!/$ - ad2 (2)
+
i=l
subject to the constraints given by eq 1. According to the method of Lagrange multipliers this is equivalent to minimizing the following unconstrained objective function .v S * ( ~ i , y i , ~= j) (Ui(zi - ~ i ) * V i b i i=l
+
+
Mathematical Formulation
Consider a chemical reactor involving hr species and L elements. Let ei3 be the number of atoms of the j t h element 246 Ind.
Eng. Chem. Process Des. Develop., Vol. 12, NO. 3, 1973
The constants z j are the Lagrange multipliers. A necessary condition for the function S* to be a miriimum
