- Email: [email protected]
Research into the kinetics of the mass crystallization process and the design of MSMPR-type crystallizers
Chemical Engineering and Processing, 31
(I 992) 57-61
57
Technical Note Research into the kinetics of the mass crystallization process and the design of MSMPR-type crystallizers Karol Machej Institute of Chemical Engineering nnd Apparatus Construction, Silesim Technical University, ul. Kuczewskiego 7, 44- 100 Gliwice (Poland) (Received
October
29, 1991)
Abstract The idea of utilizing the modified Larson-Garside kinetic equation in the elaboration of measurement data is presented. Application of the data to a design process is also discussed. A computer program verifying the method and simulating the mass crystallization process is proposed.
assumption is undoubtedly true, in practice it may lead to uncertainty in the calculation results. Many expressions in the measurement processing data equations or in the design calculations are a function of L, [l-3]. Unfortunately, a determination of this parameter experimentally is very difficult and sometimes charged with a large error of the order of L,. For this reason, there is always some uncertainty that the calculations may be charged with considerable errors. The proposal, of Larson and Garside, presented in their classic work [4], avoids this inconvenience. The authors, in their considerations, imply a size of zero for the seed. This assumption, though not strict from a theoretical point of view, considerably simplifies the calculations and does not cause large errors. These authors derive equations for the MSMPR-type crystallizer [4], based on the population balance and using the well-known relations L* = k, AC”
Many articles devoted to the mass crystallization process appeared in technical periodicals at the turn of the 1960s. Most attention was paid to the theory of the process kinetics. The creation of the theoretical bases of the mass crystallization process, especially in relation to MSMPR-type crystallizers?, was a result of the work initiated by two schools: in Europe it was started, among others, by NLvlt, and in North America by Larson and Garside. Both of these scientific schools paid a lot of attention to finding practical applications of the theoretical solutions. In this approach, the McCabe rule is used for design purposes. Distinct differences exist between the two above-mentioned schools. The European school applies the concept of mass to the description of basic quantities, such as, for example, crystal growth rate and nucleation rate, and uses the classical size analysis to describe particle size distributions. On the other hand, the American school describes these quantities by employing the concept of crystal number. Of the two roughly characterized trends, the first is less popular. One reason for this is probably the assignment of a definite size L, > 0 to crystal seeds. Although, from the theoretical point of view, this
+hixing
suspension
0255-2701/92/$5.00
mixed product
removal
(1)
Ac”
B, = k,W
(2)
and assuming the independence of L* from the crystal size (the McCabe rule), which, after simple algebraic transformations, can be rewritten in the form L*
(2f;;I,;-J(‘-‘)
=
(3)
where i =mjn
and
k, = k,/k,’
(4)
For a crystallizer (Fig. 1) with a partial linear classification the crystal growth rate is described as
I/(i1) ’-i L*=2LD4Wrp exp( 1 27M
3&y ‘1-W
(5)
where y’ = 5, /T;
(6)
and 5,
=
v,/v:
(7)
7; = (V, + V*)/(V
+ V:)
(8)
The above relationships, contrary to expectation, are not widely applied in designing. This is due to frequently occurring discrepancies between the obtained and predicted results. There are various reasons for these discrepancies. The difficulties connected with the
~0 1992 ~ Elsevier Sequoia.
All rights reserved.
58
which is slightly different from eqn. (8). Accordingly, the definition of y’ (eqn. (6)) was changed to Y = t1lr2 The following formations:
L* = L,/(3r,)
(a)
(b) Fig. 1. Schematic diagrams of the MSMPR-type vacuum crystallizer; (b) evaporative crystallizer.
crystallizers:
(a)
determination of the kinetic constants k,, k,, n, m, and j, are undoubtedly one of the reasons. The sophistication of the problem is confirmed by the existence of many measuring methods to determine the kinetic constants [5]. This situation is made worse by the fact that the measurements are devoted mainly to the theoretical aspects and are carried out under idealized conditions which are rather different from those existing in industrial crystallizers. An outline of the measurement of the crystallization kinetic constants, according to design requirements, is presented in ref. 6, taking into account the above facts. The calculations are based principally on the transformed equations (3) and (5). Introducing a parameter p, where
P = (Wvpkr)P it is possible in the form
(9) to rearrange
l/Lo4=pMi+i(L*)i~1
the above-mentioned
relations
(10)
and
(11) where zz=
22 is defined V,/(Vf
+ VT)
as (12)
(13) relation
was also employed
in the trans-
(14)
Equation ( 14) exactly satisfies eqns. (3) and ( lo), and is a good approximation of the relationships (5) and ( 11). Using equalization methods for the experimental data, it is relatively easy to determine the constants p, i, andj simultaneously. However, this method causes eqn. ( 10) to become semi-empirical, and the same applies to eqn. (11). This means that the constants p, i, and j lose their theoretical interpretation (e.g. a calculated value of j does not fulfil the condition 0
59
opinion, the factor L,, has a typically conventional character. Also taking into account that the estimation of the value of L,, is practically impossible by direct measurement, due to too fine graining of the solid phase carried by the stream V,*, it is obvious that the obligatory sedimentation equation w = f(L)
provide a stochastic model of the process and to examine the history of each individual crystal separately [ 121. The example presented below describes the procedure. Using the measurement data dealing with the crystallization process of CuSO,. 5H,O from the acidified solution, carried out in a vacuum crystallizer with outer circulation and a volume of 36 1 (together with the tubes), the following equations were derived:
(15)
should be employed. The simple relation proposed by N$vlt, based on the statement that the oversize fraction equal to 64.7% in the crystal product just corresponds to the size L,, can be used to determine the value of L,
L* = 3.117 x 1O-7 AcO.‘*~
(18)
L,e4
(19)
= 1.2014 x 1023M3.92’(L*)Z.5’4
ill].
which gave
important parameters, describing the quality of the crystal product, is the mean crystal size L,. Moreover, when McCabe’s rule is fulfilled, the following dependence is also valid:
B, = 3.991 x 10-3M492’ Ac2-h02
One of the more
L* = L,/(42,)
Taking
(16)
into consideration
eqn. (14),
4Lo = 3L,
(17)
may be obtained. Knowing L, and calculating z, from (7), the last value, L*, can be found from (14). The method of performing the measurement and calculation operations leading to the determination of the value L* makes an adequate estimation of the level of accuracy of this calculation difficult. Therefore, an’ attempt at an approximate estimation of the accuracy A method of the calculation of L* was undertaken. using a computer simulation of the crystallization process was employed. Equations ( 1) and (2) and the population, mass and heat dependences were utilized to
volume of vacuum crystallizer volumetric feed rate volumetric vapour rate feed temperature feed concentration of CuS04.5H20
500 I
3
iT
I
I
7
1I
I
-
I
-LOO
20-
7L
3
0
% II
2 15-
-300
lo-
-200
I -100
5-
i
I
0 0
I I I
i I
250
I I I 500
z
,
I I I
z
I
B and M.
0 1000
750
time [II
Fig. 2. Effect cm time on the parameters
for the
2.7 m3 3.3 x 10e4 m3/s 1.5 m3js 80 “C 0.45
The behaviour of the simulated crystallizer was characterized by the curves in Fig. 2. The curves, on which M and B were placed, showed the fluctuations of the suspension composition indirectly. Measurements of some averaged parameters (presented in Table l), for the ranges marked in Fig. 2, were carried out. Furthermore, the particle size distributions were determined in the particular ranges. The mean particle size distribution of all the ranges (I-IV) in the Lz coordinate system is presented in Fig. 3. The following well-known
25 I
(20)
Employing eqns. (18) and (20), a simulation following parameters was carried out:
60
TABLE
1. General
experimental
Parameter
results
Range
* (s)
t (“C) P (N/m*) L,, x 10’ (m)
(as in Fig. 2)
I
II
III
IV
I-IV
8662 25.92 3055.5 0.6613
8665 25.93 3056.9 0.6169
8668 25.89 3049.6 0.6212
8668 25.95 3059.9 0.60 14
8665 25.92 3055.5 0.6252
Z
6-
1 .‘ ‘*
L Fig.
3. Characteristics
relationship z is valid:
of the product
between
g=l-(1+z+z2/2+z3/6)e-’
particle
the measured
size distribution.
mass fraction
g and (21)
TABLE 2. Comparison and (16) Range
where z = L/(,5*7,)
mm
of the L* values
According to the theoretical assumptions, the points should lie along a straight line. Starting from the origin of the coordinates, the line should pass through the point (L,, 4). This point and the theoretical line were plotted on the diagram. It is clear that the measurement data do not fit the straight line, although the investigated system fulfilled McCabe’s rule (eqn. ( 18)). The existing discrepancies are partly caused by the theoretical grain size, the maximum value of which was assumed to be infinite, which can never be obtained in practice. It was mentioned earlier that the value of L* should be determined on the basis of (14) or (16). The values of L* calculated from eqn. (16) are collected together in
from
eqns.
( 16)
L* x 108 Eqn.
(22)
calculated
(16)
Eqn.
( 18)
I II III IV
I .9086 1.7799 1.7916 1.7345
1 .X066 1.7868 1.7470 1.7470
I-IV
1JO38
1.7736
Table 2. For comparision, the values determined from eqn. (18) are also presented in Table 2. It can be seen that the agreement is entirely satisfactory. Equations ( 14), ( 16) and (22) show that the characteristic dimensions L, and L, correspond to values z = 3 and z = 4. Using this information, the parameters LD and L, were calculated from the measurement data by the linear interpolation method and collected in
61 TABLE
3. Characteristic
m, n, P
sizes of the crystals
Range
Lo x 102
L, x 101
LDIL,
L, x IO’ (measured)
I
II III IV
0.5193 0.4669 0.4871 0.4899
0.6675 0.6375 0.6221 0.6337
0.7779 0.7325 0.7830 0.7731
0.6613 0.6169 0.6212 0.6014
I-l\’
0.49 14
0.6400
0.7679
0.6252
Table 3. The calculated values of L, are close to those measured directly. On the other hand, the ratio L,/L, in particular ranges differs slightly from the number 0.75 resulting from eqn. (17). To summarize, it can be said that the results of the simulated experiments presented here confirm the possibility of employing the proposed method in industrial practice. It is important to mention that, in order to utilize this method for experimental purposes, the measurements carried out should cover a period of time long enough to determine the required average values.
Nomenclature number of crystals in unit volume nucleation rate concentration (mass fraction) kinetic constants kinetic constants linear crystal growth rate dominant size ‘cut’ crystal size average crystal size seed size mass of crystals per unit volume of suspension
V, VT v: ‘i
P 71.2
*,
kinetic constants working volume of crystallizer volumetric product rate (suspension) volumetric rate of circulating suspension mean retention time ratio crystal density mean retention time volumetric shape factor
References 1 J. Nivlt, In&sfriol Cr~stullisation, Verlag Chemie, Weinheim New York, 1978. 2 J. Nivlt and M. Broul. Generalized design relations for continuous mixed crystallizer, Collect. Czech. Chem. Commun., 41 ( 1976) 536-539. 3 H. KoEovB, Z. Blechta, M. cern$ and J. NLvlt, Kontinualni laboratorni krystalirator pro stanoveni kinetiky krystalizace, Chem. Prtim., 27152 (9) (1977) 443-446. 4 M. A. Larson and J. Garslde, Crystallizer design techniques using the population balance. Chem. Eng. (London), 6 (1973) 318-328. 5 J. Njvlt. SzybkoSi: wzrostu krysztaibw. Cz@ I i II, In!. Apar. Chrm., 29/21 (2-3) (1990) 3-10. 6 K. Machej, Koncepcja badali modelowych na potrzeby projektowania krystalizator6w pr6iniowych i wyparnych. Int. C!!enl. Prows., 4 ( 1986) 627-636. 7 Z. Rojkowski, Teoretyczne podstawy procesu krystalizacji masowej, Int. Apar. Chem., 23/15 (4) (1984) 3-10. 8 E. V. Khamski, Kristnllixtsiycr u Zhimicheskoi Promyshlennosti, Khimya. Moscow, 1979. 9 R. H. Hedrick, Checking crystallizer performance, C/tern. Eng., 7 (1988) 116-l 18. 10 E. J. de Jone. Entwicklune van Kristallisatoren. Chem.-ln~.Tech., 54 (3;i 1982) 193-202. 11 K. Machej and J. Wbjcik. Opracowanie danych fizykochemicznych i kinetycznych dla potrzeb projektowania krystalizator6w przy pomocy EMC, In!. Apar. Chew., 26/18 (4) (1987) 13-15. 12 K. Machej, Symulacja procesu krystalizacji masowej na elektronicznej maszynie cyfrowej, Inf. Apar. Chem., .?9/21 (2-3) (1990) 19-22.
(I 992) 57-61
57
Technical Note Research into the kinetics of the mass crystallization process and the design of MSMPR-type crystallizers Karol Machej Institute of Chemical Engineering nnd Apparatus Construction, Silesim Technical University, ul. Kuczewskiego 7, 44- 100 Gliwice (Poland) (Received
October
29, 1991)
Abstract The idea of utilizing the modified Larson-Garside kinetic equation in the elaboration of measurement data is presented. Application of the data to a design process is also discussed. A computer program verifying the method and simulating the mass crystallization process is proposed.
assumption is undoubtedly true, in practice it may lead to uncertainty in the calculation results. Many expressions in the measurement processing data equations or in the design calculations are a function of L, [l-3]. Unfortunately, a determination of this parameter experimentally is very difficult and sometimes charged with a large error of the order of L,. For this reason, there is always some uncertainty that the calculations may be charged with considerable errors. The proposal, of Larson and Garside, presented in their classic work [4], avoids this inconvenience. The authors, in their considerations, imply a size of zero for the seed. This assumption, though not strict from a theoretical point of view, considerably simplifies the calculations and does not cause large errors. These authors derive equations for the MSMPR-type crystallizer [4], based on the population balance and using the well-known relations L* = k, AC”
Many articles devoted to the mass crystallization process appeared in technical periodicals at the turn of the 1960s. Most attention was paid to the theory of the process kinetics. The creation of the theoretical bases of the mass crystallization process, especially in relation to MSMPR-type crystallizers?, was a result of the work initiated by two schools: in Europe it was started, among others, by NLvlt, and in North America by Larson and Garside. Both of these scientific schools paid a lot of attention to finding practical applications of the theoretical solutions. In this approach, the McCabe rule is used for design purposes. Distinct differences exist between the two above-mentioned schools. The European school applies the concept of mass to the description of basic quantities, such as, for example, crystal growth rate and nucleation rate, and uses the classical size analysis to describe particle size distributions. On the other hand, the American school describes these quantities by employing the concept of crystal number. Of the two roughly characterized trends, the first is less popular. One reason for this is probably the assignment of a definite size L, > 0 to crystal seeds. Although, from the theoretical point of view, this
+hixing
suspension
0255-2701/92/$5.00
mixed product
removal
(1)
Ac”
B, = k,W
(2)
and assuming the independence of L* from the crystal size (the McCabe rule), which, after simple algebraic transformations, can be rewritten in the form L*
(2f;;I,;-J(‘-‘)
=
(3)
where i =mjn
and
k, = k,/k,’
(4)
For a crystallizer (Fig. 1) with a partial linear classification the crystal growth rate is described as
I/(i1) ’-i L*=2LD4Wrp exp( 1 27M
3&y ‘1-W
(5)
where y’ = 5, /T;
(6)
and 5,
=
v,/v:
(7)
7; = (V, + V*)/(V
+ V:)
(8)
The above relationships, contrary to expectation, are not widely applied in designing. This is due to frequently occurring discrepancies between the obtained and predicted results. There are various reasons for these discrepancies. The difficulties connected with the
~0 1992 ~ Elsevier Sequoia.
All rights reserved.
58
which is slightly different from eqn. (8). Accordingly, the definition of y’ (eqn. (6)) was changed to Y = t1lr2 The following formations:
L* = L,/(3r,)
(a)
(b) Fig. 1. Schematic diagrams of the MSMPR-type vacuum crystallizer; (b) evaporative crystallizer.
crystallizers:
(a)
determination of the kinetic constants k,, k,, n, m, and j, are undoubtedly one of the reasons. The sophistication of the problem is confirmed by the existence of many measuring methods to determine the kinetic constants [5]. This situation is made worse by the fact that the measurements are devoted mainly to the theoretical aspects and are carried out under idealized conditions which are rather different from those existing in industrial crystallizers. An outline of the measurement of the crystallization kinetic constants, according to design requirements, is presented in ref. 6, taking into account the above facts. The calculations are based principally on the transformed equations (3) and (5). Introducing a parameter p, where
P = (Wvpkr)P it is possible in the form
(9) to rearrange
l/Lo4=pMi+i(L*)i~1
the above-mentioned
relations
(10)
and
(11) where zz=
22 is defined V,/(Vf
+ VT)
as (12)
(13) relation
was also employed
in the trans-
(14)
Equation ( 14) exactly satisfies eqns. (3) and ( lo), and is a good approximation of the relationships (5) and ( 11). Using equalization methods for the experimental data, it is relatively easy to determine the constants p, i, andj simultaneously. However, this method causes eqn. ( 10) to become semi-empirical, and the same applies to eqn. (11). This means that the constants p, i, and j lose their theoretical interpretation (e.g. a calculated value of j does not fulfil the condition 0
59
opinion, the factor L,, has a typically conventional character. Also taking into account that the estimation of the value of L,, is practically impossible by direct measurement, due to too fine graining of the solid phase carried by the stream V,*, it is obvious that the obligatory sedimentation equation w = f(L)
provide a stochastic model of the process and to examine the history of each individual crystal separately [ 121. The example presented below describes the procedure. Using the measurement data dealing with the crystallization process of CuSO,. 5H,O from the acidified solution, carried out in a vacuum crystallizer with outer circulation and a volume of 36 1 (together with the tubes), the following equations were derived:
(15)
should be employed. The simple relation proposed by N$vlt, based on the statement that the oversize fraction equal to 64.7% in the crystal product just corresponds to the size L,, can be used to determine the value of L,
L* = 3.117 x 1O-7 AcO.‘*~
(18)
L,e4
(19)
= 1.2014 x 1023M3.92’(L*)Z.5’4
ill].
which gave
important parameters, describing the quality of the crystal product, is the mean crystal size L,. Moreover, when McCabe’s rule is fulfilled, the following dependence is also valid:
B, = 3.991 x 10-3M492’ Ac2-h02
One of the more
L* = L,/(42,)
Taking
(16)
into consideration
eqn. (14),
4Lo = 3L,
(17)
may be obtained. Knowing L, and calculating z, from (7), the last value, L*, can be found from (14). The method of performing the measurement and calculation operations leading to the determination of the value L* makes an adequate estimation of the level of accuracy of this calculation difficult. Therefore, an’ attempt at an approximate estimation of the accuracy A method of the calculation of L* was undertaken. using a computer simulation of the crystallization process was employed. Equations ( 1) and (2) and the population, mass and heat dependences were utilized to
volume of vacuum crystallizer volumetric feed rate volumetric vapour rate feed temperature feed concentration of CuS04.5H20
500 I
3
iT
I
I
7
1I
I
-
I
-LOO
20-
7L
3
0
% II
2 15-
-300
lo-
-200
I -100
5-
i
I
0 0
I I I
i I
250
I I I 500
z
,
I I I
z
I
B and M.
0 1000
750
time [II
Fig. 2. Effect cm time on the parameters
for the
2.7 m3 3.3 x 10e4 m3/s 1.5 m3js 80 “C 0.45
The behaviour of the simulated crystallizer was characterized by the curves in Fig. 2. The curves, on which M and B were placed, showed the fluctuations of the suspension composition indirectly. Measurements of some averaged parameters (presented in Table l), for the ranges marked in Fig. 2, were carried out. Furthermore, the particle size distributions were determined in the particular ranges. The mean particle size distribution of all the ranges (I-IV) in the Lz coordinate system is presented in Fig. 3. The following well-known
25 I
(20)
Employing eqns. (18) and (20), a simulation following parameters was carried out:
60
TABLE
1. General
experimental
Parameter
results
Range
* (s)
t (“C) P (N/m*) L,, x 10’ (m)
(as in Fig. 2)
I
II
III
IV
I-IV
8662 25.92 3055.5 0.6613
8665 25.93 3056.9 0.6169
8668 25.89 3049.6 0.6212
8668 25.95 3059.9 0.60 14
8665 25.92 3055.5 0.6252
Z
6-
1 .‘ ‘*
L Fig.
3. Characteristics
relationship z is valid:
of the product
between
g=l-(1+z+z2/2+z3/6)e-’
particle
the measured
size distribution.
mass fraction
g and (21)
TABLE 2. Comparison and (16) Range
where z = L/(,5*7,)
mm
of the L* values
According to the theoretical assumptions, the points should lie along a straight line. Starting from the origin of the coordinates, the line should pass through the point (L,, 4). This point and the theoretical line were plotted on the diagram. It is clear that the measurement data do not fit the straight line, although the investigated system fulfilled McCabe’s rule (eqn. ( 18)). The existing discrepancies are partly caused by the theoretical grain size, the maximum value of which was assumed to be infinite, which can never be obtained in practice. It was mentioned earlier that the value of L* should be determined on the basis of (14) or (16). The values of L* calculated from eqn. (16) are collected together in
from
eqns.
( 16)
L* x 108 Eqn.
(22)
calculated
(16)
Eqn.
( 18)
I II III IV
I .9086 1.7799 1.7916 1.7345
1 .X066 1.7868 1.7470 1.7470
I-IV
1JO38
1.7736
Table 2. For comparision, the values determined from eqn. (18) are also presented in Table 2. It can be seen that the agreement is entirely satisfactory. Equations ( 14), ( 16) and (22) show that the characteristic dimensions L, and L, correspond to values z = 3 and z = 4. Using this information, the parameters LD and L, were calculated from the measurement data by the linear interpolation method and collected in
61 TABLE
3. Characteristic
m, n, P
sizes of the crystals
Range
Lo x 102
L, x 101
LDIL,
L, x IO’ (measured)
I
II III IV
0.5193 0.4669 0.4871 0.4899
0.6675 0.6375 0.6221 0.6337
0.7779 0.7325 0.7830 0.7731
0.6613 0.6169 0.6212 0.6014
I-l\’
0.49 14
0.6400
0.7679
0.6252
Table 3. The calculated values of L, are close to those measured directly. On the other hand, the ratio L,/L, in particular ranges differs slightly from the number 0.75 resulting from eqn. (17). To summarize, it can be said that the results of the simulated experiments presented here confirm the possibility of employing the proposed method in industrial practice. It is important to mention that, in order to utilize this method for experimental purposes, the measurements carried out should cover a period of time long enough to determine the required average values.
Nomenclature number of crystals in unit volume nucleation rate concentration (mass fraction) kinetic constants kinetic constants linear crystal growth rate dominant size ‘cut’ crystal size average crystal size seed size mass of crystals per unit volume of suspension
V, VT v: ‘i
P 71.2
*,
kinetic constants working volume of crystallizer volumetric product rate (suspension) volumetric rate of circulating suspension mean retention time ratio crystal density mean retention time volumetric shape factor
References 1 J. Nivlt, In&sfriol Cr~stullisation, Verlag Chemie, Weinheim New York, 1978. 2 J. Nivlt and M. Broul. Generalized design relations for continuous mixed crystallizer, Collect. Czech. Chem. Commun., 41 ( 1976) 536-539. 3 H. KoEovB, Z. Blechta, M. cern$ and J. NLvlt, Kontinualni laboratorni krystalirator pro stanoveni kinetiky krystalizace, Chem. Prtim., 27152 (9) (1977) 443-446. 4 M. A. Larson and J. Garslde, Crystallizer design techniques using the population balance. Chem. Eng. (London), 6 (1973) 318-328. 5 J. Njvlt. SzybkoSi: wzrostu krysztaibw. Cz@ I i II, In!. Apar. Chrm., 29/21 (2-3) (1990) 3-10. 6 K. Machej, Koncepcja badali modelowych na potrzeby projektowania krystalizator6w pr6iniowych i wyparnych. Int. C!!enl. Prows., 4 ( 1986) 627-636. 7 Z. Rojkowski, Teoretyczne podstawy procesu krystalizacji masowej, Int. Apar. Chem., 23/15 (4) (1984) 3-10. 8 E. V. Khamski, Kristnllixtsiycr u Zhimicheskoi Promyshlennosti, Khimya. Moscow, 1979. 9 R. H. Hedrick, Checking crystallizer performance, C/tern. Eng., 7 (1988) 116-l 18. 10 E. J. de Jone. Entwicklune van Kristallisatoren. Chem.-ln~.Tech., 54 (3;i 1982) 193-202. 11 K. Machej and J. Wbjcik. Opracowanie danych fizykochemicznych i kinetycznych dla potrzeb projektowania krystalizator6w przy pomocy EMC, In!. Apar. Chew., 26/18 (4) (1987) 13-15. 12 K. Machej, Symulacja procesu krystalizacji masowej na elektronicznej maszynie cyfrowej, Inf. Apar. Chem., .?9/21 (2-3) (1990) 19-22.