Desalination, 28 (1979) 1-11
© Elsevier Scientific Publishing Company, Amsterdam - Printed in The Netherlands
DYNAMIC MODELLING AND PARTIAL SIMULATION OF A PILOT SCALE COLUMN CRYSTALLIZER *
R . AKHTAR, L . McGRATH AND P . D . ROBERTS The City University, London (United Kingdom) (Received January 10, 1979)
SUMMARY
A dynamic model of a continuous crystallizer has been attempted as an extension of our work on the continuous freezing of ice from saline solutions and subsequent melting of the ice to give a potable water . In pilot plant equipment three distinct zones are designated as freezing, purification and melting, and our ultimate aim is to achieve a unified model applicable to the total system . The scope of this paper covers the freezing section in which basic equations for heat and mass transfer were applied covering the different components .
SYMBOLS
Aw CB CR Dp E Ex f G gb HF Hp A L A
Nucleation rate equation empirical constant Area of freezing wall chamber, m 2 Specific heat of liquid brine, J kg 1 K-1 Specific heat of refrigerant, J kg-1 Ki Particle diameter, m Empirical constant in Eq . 10 Ratio of actual thermal driving force to apparent driving force Mass transfer from liquid to solid, kg s i Crystal growth rate constant Constant determining entrained liquid Height of freezing chamber, m Planck's constant, J s-1 Latent heat of solid (ice), J kg F' K 1 Liquid brine feed, kg s 1
* Presented at the Sixth International Symposium on Fresh water from the Sea, September 17-22, 1978, Las Palmas, Gran Canaria .
2 LD Ls LF ML MR
MS PS Q
r RG S S,,, S„ Tb TFi T1 TLi Ta, U V WL
WLi Wq
N. AKHTAR, L. McGRATH AND P. D. ROBERTS
- Adhered liquid leaving freezing section, kg s 1
- Concentrated brine take-off, kg s 1 - Flow rate of refrigerant, kg s 1 - Mass of liquid brine hold up in freezing section, kg - Mass of refrigerant surrounding freezing chamber, kg - Mass of solid (ice) hold-up in freezing chamber, kg - Density of solid (ice), kg m3 - Reflux mass flow rate, kgs1 - Average crystal radius, m - Ideal gas constant, J K-1 mole- ' - Crystal mass rate out of freezing chamber, kg s 1 - Mass of a single crystal, kg - Nucleation rate, No_ of particles s-1 - Bulk temperature of liquid brine, K - Freezing temperature of brine used, K - Inlet temperature of refrigerant, K - Ice/liquid interface temperature, K - Inlet temperature of brine, K - Temperature of freezing chamber wall, K - Heat transfer co-efficient, Jm2 K-1 S-1 - Vertical velocity of screw, m s-1 - Mass fraction of salt in liquid brine - Mass fraction of salt in brine feed - Mass fraction of salt in reflux
INTRODUCTION
Early work on column crystallization was carried out by Arnold [1] whose equipment was modified and improved by Schildknecht [2] for the separation of organic mixtures, the equipment herein described being of similar form . In the process of desalination by freezing, an ice slurry is transported to a melter, some of the melt water being used to wash away adhered brine solution from the ice crystals . The use of a column crystallizer for the purpose of desalination was described by Hobson and McGrath [3] and further developments given later [4] . Modelling of the system in a column crystallizer was attempted by Albertins and Powers [5] and later by Gates and Powers [6] ; Henry and Powers [7] presented a steady state model based on the analogy to the continuous distillation process . A further contribution to the steady state modelling was made by Gladwin [8) . The dynamic modelling and simulation of a continuous column crystallizer here presented takes account of the earlier steady state models [5-8] and methods of Franks [9] and Rose [10] in the analysis of continuous
PILOT SCALE COLUMN CRYSTALLIZER
3
systems_ The process of nucleation is considered to follow the Eyring rate equation for nucleation from a melt which was used by Umano and Kawasaki [111 for the formation of ice from sea water . Considerations of the work by Wey and Estrin [121 on diffusion and mass transfer effects based on the modelling of the ice-brine batch crystallization process and of Bolsaitis [131 and Devyatykh [141 for continuous countercurrent crystalliz~ation have also been taken into account . PROCESS DESCRIPTION
A photograph of the column showing the various components is presented in Fig. 1 . As explained earlier, the column consists of three sections : freezing, purification and melting sections, but is in fact joined together into a single piece of equipment. The feed point is close to the freezing section and the inlet feed goes almost directly into the freezing section . The initial separation of the components takes place here by means of fractional solidification, a refrigerant being circulated around the freezing section of the column . Adherence of mother liquor to the ice crystal and occlusion within the lattice gives rise to overall impurity . By means of an Archimedean spiral conveyor, the crystals are transported out of the freezing section, into the purification section and finally into the melting section . Reflux is generated in the melting section due to a heat input ; the major part of the melt is taken off as product and the remainder is returned as a countercurrent wash stream . This pure reflux serves to wash adhered liquid
Fig. 1 . 4 in . diameter pilot plant crystallizer.
N.AKHTAR, L. McGRATH AND P . D . ROBERTS
4
from around the crystals within the purification section . The crystals melt in the melting section due to the pressure of accumulated crystals and application of heat by means of hot water circulation . The purification section of the column is also jacketted to avoid loss of heat and a temperature gradient is established inside the column between the freezing and melting sections . DYNAMIC MODEL OF FREEZING SECTION
A schematic diagram of the freezing section is shown in Fig . 2. Ref Jut
Adhereo
Brine feed L
i~ ILD
Refngerant Lr
liquid
I
Liqui d Solid phase , phase S
S
Ice crystal
F Liquid
take off
• LE
Fig. 2 . Schematic of the freezing section . Brine of composition WL ; and temperature TL ; is fed into the section at a rate L mass per unit time . Liquid reflux of brine composition W o enters from the purification section at a rate Q mass per unit time and concentrated brine is taken off at a rate L E . Ice crystals leave the freezing section and enter the purification section at a rate S together with adhered liquid at a rate LD _ Brine is also entrained with the ice in the amount gb mass per unit mass of ice . The dynamic mathematical model is developed by applying conservation laws of mass and heat to the various phases and components of the freezing section . However, in order to simplify the solution of the equations it is necessary to make the following assumptions concerning the nature of the process : (a) The contents of the freezing section are considered to be completely mixed . Thus the compositions and temperatures of the liquid streams leaving the section are considered to be the same as the composition W L and temperature T b of the bulk liquid in the freezing section . (b) The reflux is considered to enter the freezing section at bulk temperature Tb . (c) The heat capacity of the walls is assumed to be negligible in comparison to those of the refrigerant and bulk liquid. (d) Variations in density, specific heat, latent heat and heat transfer coefficient, are considered to be negligible over the temperature range of interest. (e) Crystal size may be defined by a single dimension and, for the purpose of determining crystal mass, crystals are assumed to be spherical . (f) Nucleation rate can be approximated by adapting the Eyring rate theory to nucleation from a melt [15] .
PILOT SCALE COLUMN CRYSTALLIZER
5
An overall mass balance on the liquid phase determines the mass M L of brine- in the freezing section according to the first-order differential equation dML dt
L-L E - LD - I - Sgb
(1)
+ Q
where f is the rate of mass transfer from liquid to solid during the formation of the ice crystals . The mass MS of crystals is determined from a mass balance on the solid phase which gives d Ms _ f -S
(2)
dt
A salt balance gives the composition W L of brine in the bulk liquid in terms of the salt compositions W Lj and W Q of the brine feed and reflux according to the equation dWL _ 1 dt
ML
{L(WLI -
WL)
+ Q(WQ -
WL)
+ fWL }
(3)
The bulk temperature Tb of the liquid depends upon the temperature TLj of the brine feed and the temperature T R of the refrigerant surrounding the section wall according to a heat balance on the liquid phase giving dTb 1 { LCB (TLI - Tb) dt ML CB
- fX
-UA W (Tb
-
TW )}
(4)
Similarly, the temperature Tjv is determined by applying a heat balance to the refrigerant to obtain dTA.
1
dt
MR cR
{L F C R (TFj -TW ) + UA W (Tb -TW )}
(5)
Crystal growth is calculated from a thermal driving force relationship [16] obtained from a dynamic heat balance on the solid phase . This determines the rate of change of diameter Dp of a single crystal with respect to the difference between the temperature Tj of the ice/liquid interface and temperature Tb of the bulk liquid, according to dDp _ a (T, TO dt Dp
(6)
6
N. AKHTAR, L . McGRATH AND P . D. ROBERTS
where G is a rate constant defined in terms of a driving force ratio E,, as -8
G = 0 .6129 10
Ex
(7)
and Ex is dependent on salt concentration W L of the bulk liquid according to Ex
= 1/ {1 + 77WL /(1 - WL )}
(8)
The ice/liquid interface temperature T; is defined in terms of the liquid bulk temperature and the driving force ratio using
T, = (Tf - Tb )Ex +
(9)
Tb
where Tf is the freezing temperature of brine used (271 .7K). In order to complete the mathematical model, expressions are required for rate of adhered liquid LD , the mass transfer rate f and the crystal mass rate S to the purification section . The adhered liquid leaving the freezing section with the crystals depends upon crystal site and the crystal mass rate and is assumed to be of the form
LD = E S f(r)
(10)
where F = 0 .5 Dp is the average crystal radius and E is a constant of proportionality . For the purpose of this preliminary investigation the functional expression f(r) is assumed to be linear. The rate of mass transferred from liquid to solid depends upon the nucleation rate S„ and the mass of a single crystal S m according to the relationship f = Sn Sm
(11)
where S„ is the nucleation rate given as follows,
S„ =
R` H
Tb
A(Tf z exp ) Tb ) 2 l, Tb < Tf Tb (Tf J
[16) (12)
, Tb > Tf
0
R e, is the ideal gas constant, Hp is Planck's constant and A is an empirical constant which may be calculated by assuming that S n = 1 at nucleation . Assuming that the crystals are spherical, the mass of a single crystal is obtained from
Sm = 47rr 3 ps /3 where ps is the density of ice .
(13)
Finally, the rate at which crystals are transported to the purification section depends upon the rotational speed of the screw and is given by
MS S = H
V
(14) F
PILOT SCALE COLUMN CRYSTALLIZER
7
where v is the vertical velocity of the screw and Hp is the height of the freezing section . Computer simulation
Eqs . (1) to (14) define a set of differential and algebraic relationships which may be employed for simulating the freezing section of the column crystallizer. Digital computer simulation has been chosen using a modular subroutine structure with standard FORTRAN . A master routine is initially called which defines the storage area required, reads in the model data and initial conditions, and calls a library numerical integration routine based on Hamming's modified predictor corrector method [171 . This is a fourth order technique using four preceding points to compute a new output state at each instance of simulated time . At each simulated time instance, the integration routine is supplied with the right-hand side of each differential equation and returns with the appropriate integrated values . This is accomplished through calls to a user supplied subroutine containing all the algebraic relationships and the right-hand sides of the first-order differential equations which comprise the mathematical model and whose purpose is to supply the numerical integration routine with the values of all the derivatives at each simulated time instance. The numerical integration routine adjusts the integration step size during the whole computation to satisfy a given error criterion . When accurate values of the integrated states are obtained the numerical integration routine calls another user supplied subroutine which controls the printing of results . The complete time response of the simulated system is obtained, entirely under the control of the numerical integration routine, by incrementing the simulated time parameter from zero to a chosen final time value . DISCUSSION AND RESULTS
Experiments have been performed using the programmed computer simulation of the mathematical model with the aim of determining the dynamic response of the freezing section during start-up . Data for the mathematical model have been obtained from the constructional dimensions of the pilot plant crystallizer and physical properties of brine and ice . These data are given in Table I. Feed, and take-off rates and conditions, constitute boundary inputs to the model and these are listed in Table II . Table III gives the initial conditions for the start of a simulation experiment. The freezing section is considered to contain an initial quantity of liquid at feed temperature, void of ice, and surrounded by an empty refrigerant chamber at room temperature . Refrigerant is then fed through the chamber in order to start the freezing process and subsequent ice crystal formation . Fig . 3 shows the variation of the temperature of the bulk liquid as the liquid is cooled to freezing conditions . The initial small rise in temperature occurs because the walls of the freezing chamber are at room temperature,
8
N . AKHTAR, L. McGRATH AND P. D. ROBERTS
TABLE I DATA A Aw CB CR
E gb N_ F
Hp MR RG Tf U
x Ps
Nucleation rate equation empirical constant Area of freezing chamber wall Specific heat of liquid brine Specific heat of refrigerant Empirical constant in adhering liquid equation Constant determining entrained liquid Height of freezing chamber Planck's constant Mass of refrigerant around freezing chamber Ideal gas constant Freezing point of brine used Heat transfer coefficient Latent heat of ice Density of ice
-8.88 0.0405 m2 4017Jkg 1 K1 3138 J kg1 K-1 0.06 0 .5 x 105 0.33 m 6.6256 x 10-34 J S71 0.45 kg 8.314 J K 1 mole1 271 .7 K 142Jm 2 K-1 s-1 3.33 x 10 5 J kg' 920 kg m-3
TABLE II BOUNDARY INPUTS L LE LF Q TLi
TFi V IV Li PVq
Liquid brine feed Concentrated brine take-off Flow rate of refrigerant Reflux mass flow rate Inlet temperature of brine Inlet temperature of refrigerant Vertical velocity of screw Mass fraction of salt in brine feed Mass fraction of salt in reflux
3.45 X 10 4 1 .53 X 10-4 0.02 kg s 1 3.45 X 10-' 275 K 258 K 3.70 X10 4 0 .0" 5 0.0035
TABLE III INITIAL CONDITIONS Dp ML MS
Tb
fly IVL
Crystal mean diameter Mass of liquid brine Mass of ice crystals Bulk temperature of liquid brine Temperature of freezing chamber wall Mass fraction of salt in liquid brine
0.00 m 0 .3795 kg 0.00 kg 275 K 292 K 0.035
-1
kg S kg S-1 kg
s1
ms 1
9
PILOT SCALE COLUMN CRYSTALLIZER
275
a E r 265 0
500
2000
1000 1500 Time (s)
Fig. 3. Variation of bulk temperature with time .
103
102
10 1 I I 1500 1000 0 500 Time (s)
Fig. 4. Variation of nucleation rate with time. 0.10
0.10
t
n 7
N
1 9 O
s 0.05
v 0 N
r 1
I I I
5000
10 000 Time (s)
15 000
1 I 1 1 1 i 1 40 20 60 80 Time (s)
I
0
Fig. 5. Variation of mass of solids hold-up with time . Fig. 6 . Variation of transfer of ice crystals to purification section with time .
10
N. AKHTAR, L . McGRATH AND P. D. ROBERTS
which is higher than the brine liquid feed temperature, at the beginning of the response. Nucleation is arranged to take place only when a finite minimum undercooling occurs (3 .5° C) and Fig. 4 shows a rapid rise in nucleation rate after 300 seconds followed by a gradual fall . The rate of nucleation is governed by Eq. (12) which implies that the rate can change by several orders of magnitude within a small range of bulk temperature . The consequent sensitivity of the model to change in bulk temperature can give rise to numerical difficulties within the computer integration routine requiring control over the allowable range of the exponential in Eq . (12) . As soon as nucleation occurs, crystals start to form as shown by Figs . 5 and 6 which illustrate the accumulation of the solid mass and transfer of ice crystals to the purification section as functions of time . CONCLUSIONS
The behaviour of the column may be predicted from the results obtained for the freezing section so far as nucleation and growth of solids present are concerned . Similarly, liquid behaviour - its temperature with respect to the kind and rate of refrigerant circulated around the freezing section - is predicted. An optimum temperature range is essential to run the column efficiently . Above this temperature level there is very little nucleation and below it there is over-production of crystals which, practically, results in a blockage of the column . In the case of a numerical analysis this over-protection is given in terms of a numerical overflow and the program is thrown out . The graphs of heat and mass balance, show similar behaviour to the pilot plant experiments of Gladwin [8] although they were not discussed extensively by him . The optimum working temperatures of the column obtained experimentally were confirmed by using the dynamic model . REFERENCES 1. 2. 3. 4. 5. 6. 7. 8. 9.
P . M . Arnold, U .S. Patent 2, 540, 977, 1951 . H. Schildknecht, J . Breiter and K. Maas, Separation Sci ., 5 (1970) 99-112 . M. D. Hobson and L . McGrath, Proc. Fourth Intern . Symp- on Fresh Water from the Sea, Heidelberg, Sept . 1973, 3 (1973) 357-369 . C . Bates, R. P. Gladwin and L . McGrath, Proc . Fifth Intern. Symp. on Fresh Water from the Sea, Alghero, May 1976, 3 (1976) 201-212 . R. Albertins and J. E. Powers, Am. Inst. Chem . Eng . J ., 15 (1969) 554-560 . W. C. Gates and J . E. Powers, Am . Inst. Chem . Eng. J., 16 (1970) 648-658 . J . D. Henry, Jr. and J. E. Powers, Am . Inst. Chem . Eng. J., 16 (1970) 1055-1063 . R. P . Gladwin, Ph .D . Thesis, City University, London, 1974R . G. E . Franks, Modelling and Simulation in Chemical Engineering, Wiley, New York, N.Y., 1961 .
PILOT SCALE COLUMN CRYSTALLIZER
11
10 . L. M. Rose, The Application of Mathematical Modelling to Process Development and Design, Applied Science Publishers, London, 197611 . S. Umano and S . Kawasaki, Tokyo Kogyo, Shikensho Hekoku, 54 (1959) 233-245 . 12. J. S. Wey and J. Estrin, Ind . Eng. Chem. Process Des . De•7elop., 12 (1973) 236-245. 13. P. Bolsaitis, Chem . Eng. Sci., 24 (1969) 1813-1825 . 14. G. G. Devyatykh, Russ. J. Phys. Chem ., 41 (1967) 507-50915. A . Vanhook, Crystallization, Reinhold, New York, N.Y., 1961 . 16. J. C. Orcutt, Fractional Solidification, Vol. 1, Chap . 17, pp 441-459, Arnold, London, 1967 . 17. A. Raison and H. S. Wilkf, Mathematical Methods for Digital Computers, Wiley, New York, N.Y., 1960, pp. 95-109 .