403
Ind. Eng. Chem. Process Des. Dev. 1904, 23, 463-460
Modeling the Batch Crystallization Process via Shifted Legendre Polynomials Rong-Yeu Chang and Maw-Llng Wang' DepaHment of Chemicei Engineering, Natbnai Tsing Hue University, Hslnchu, Taiwan, ROC
The population balance equation, which is used to describe the crystal size distribution within the batch crystallization process, is solved. The technique of the shifted Legendre polynomials function is used to solve the particle conservation equation. The key method is that the population density function is expressed by a series of shifted Legendre functlons. The partlai differential equation of population balance is transformed into a series of ordinary differential equations of expansion coefficients. A more generalized sizedependent growth rate expression and a complicated moment form of nucleation rate expression are considered in the solution of population balance density function. One illustrative example of batch ice crystaUitation from brine solution, which involves population balance, material, and energy balances, Is given to demonstrate the validity of the presently proposed method.
Introduction The well-known population balance equation has been extensively applied to many particulate systems in chemical engineering (Himmelblau and Bischoff, 1968; Randolph and Larson, 1971). The transient and steady-state cyrstal size distributions were studied (Randolph and Larson, 19621, from which a general population balance was used to derive equations relating transient population density to ita size. The 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. The moments method (Hulburt and Katz, 1964) was then employed to approach the population balance equation describing either a batch or a continuous system of the crystallization processes. As indicated, the solution of population density function is to employ the so-called moments which relate to the crystal number, length, area, and volume. By use of such a method, a partial differential equation of population balance can be transformed into a series of ordinary differential equations of moments with physical significance. The j t h moments of population balance are calculated from those ordinary differntial equations. The population balance density function is thus obtained by inverting the moments with the aid of Laguerre function expansion in terms of the coefficients that were a linear combination of moments. However, as pointed out by Randolph and Larson (1971), the main disadvantage of using the moments method is that the higher order corrections of m+ ments involving the Laguerre polynomial terms produce oscillations about the r weighting function and the Laguerre series expansion, for the higher order term is unstable. The moment equations are normally not closed. Therefore, the other disadvantage of using the moments method to solve the population balance equation is that the calculation of the lower order moment will frequently require the information of the higher order moments when the expression of size-dependent growth rate is complicated. In order to overcome this difficulty, an appropriate expression for the higher order moments in terms of the leading moments (say po to p 3 ) associated with the expression of Laguerre polynomial functions was proposed (Hulburt and Katz, 1964). However, as pointed out previously, the expansion of Laguerre polynomial functions is unstable. Recently, a powerful method to analyze a class of problems of crystallization processes, which are described
by the population balance equation, was proposed (Chang and Wang, 1982). The population balance, either described by functional differential equation or ordinary differential equations, is solved. The computational results are satisfactory, compared with that of the exact solution and other numerical values given in the literature. Nevertheless, the conventional technique of the shifted Legendre functions (Chang and Wang, 1982) did not provide satisfactory computational results for such a complicated system as the ice-brine system. In modeling the batch crystallization process, one has to consider material and energy balance equations as well as population balance equation. Additionally, the growth rate shows infinite behavior a t zero size and the nucleation rate expression involves complicated moments form. Therefore, those conservation equations, which represent the crystallization process, need further modification before the technique of shifted Legendre polynomial functions is applied. In the following section, the model of batch crystallization process is briefly discussed. Population Balance Equation in a Batch Crystallization Process Crystallization from solution is characterized by the production of a spectrum of different size crystals. Combining population balance with the mass and energy balances and the relevant rate equations, we can predict the performance of a batch crystallizer. For a batch crystallizer with no net inflow or outflow of crystals, the population balance can be written as (Wey and Estrin, 1973)
where n is the population density of crystals at size L, and G is the growth rate of crystal. The boundary condition, n(0,t)is the nuclei population density, no,and is related to the nucleation rate, Bo by Bo = lim Gn (nucleation) (2) L-0
If the crystallizer is initially seeded, the initial population density, n(L,O),for a batch suspension crystallizer is usually expressed by the population density of seeds n(L,O) = no(L) (seeding) (3) In addition to the population balance equation, suitable kinetic equations are needed for the evaluation of the
0196-4305/84/1~23-0463$01.50/0 0 1984 Amerlcan Chemical Society
484
Ind. Eng. Chem. Process Des. Dev., Vol. 23,No. 3, 1984
nucleation rate, B O , and the growth rate, G. Empirical power law forms are frequently used to correlate G and Bo (Moyers and Randolph, 1973). For the growth rate, G G = KG(AC)d
(4)
where AC is the supersaturation. The exponent d is an empirical constant. The growth rate constant, K G , which may be influenced by temperature and degree of agitation, is a function of crystal size when the crystal growth is affected by bulk diffusion (Wang, 1974) or the GibbsThomson effect (Wey, 1977). In general, the way to fit the growth rate constant, KG, into a polynomial form of crystal size is described by J
G = (AC)d
u~L'
1=-1
(5)
where aldepends upon the physical properties, operating parameters, and growth mechanism. J is a positive integer. Equation 5 shows an infiiite growth rate at zero size which will make the population balance equation difficult to be solved (Wey and Estrin, 1973). The nucleation rate, B O , may be represented by a semiempirical power law model (Wey and Estrin, 1973).
Bo = K,(AC)"(P/V)@(l "0 U n a)'
G=
bo
+ BlX +
+ ...
and bo
a0 -.a-1 dl = -L*2' a, = a,; 63 = a,L* L*'
In eq 10, G is assumed to be expressed by an infiiite series of shifted Legendre functions of x . In general, the expansion coefficients, gI, can be obtained from Gaussian integration formula even when the expression of G(x,t) is of a very complicated form. By using such a transformation, as shown in eq 9 and 10, the singular point for G of infinite growth rate at zero size is thus prevented. The population balance density function, n(x,t), is assumed to be expressed by an infinite series of the shifted Legendre polynomials, Pj(x). 0)
n(x,t) = XC(2j j=O
+ l)LYj(t)Pj(X)
(12)
Multiplying each term of eq 8 by P,(x) and integrating the equation with respect to x from 0 to 1 gives
(6)
where P / V is the power dissipation of the crystallization process. The nucleation rate constant, KN,is likely to depend on the temperature and degree of agitation. a,P, and y are parameters of nucleation rate and are determined by experiments. Equations 5 and 6 provide a general model for expressing growth rate and nucleation rate of crystal within a batch crystallization process. The nucleation rate, which is shown in eq 6, is related to the jth moment of population density function, i.e., Som nu' dL. The upper integration limit, which is set at infinite, is convenient for mathematical manipulation using the method of moments. Physically, the upper integration limit of population balance density function should be a value of L*, which is the maximum crystal size for a real crystallization process within the crystallizer. Therefore, the population balance density function will vanish when the sizes of those crystals are greater than L*. Hence, eq 6 becomes
Equation 13 is integrated to an appropriate form
(14) One of the properties of shited Legendre function is its recurrence (i + I)P~+~(X) = (2i + 1)(2x - l)Pi(x) - iPi-l(x)
or =
XPi(X)
1 2(2i
+ 1)((i + l)Pi+l(x)+ (2i + l)Pi(x) + iPi-l(x)j (15)
Substituting eq 15 and 12 into the first term of eq 14 yields d d " 1 - l l n P i ( x ) dx= dt ;=n C(2j + l)ajJ xPi(x)P,(x)dw] = dt o
-[
In the following section, the shifted Legendre polynomial function is employed to solve the population balance equation subject to eq 5 and 7 and the initial condition. The Solution of Shifted Legendre Function Series In order to apply the shifted Legendre functions to solve the population balance equation, the independent variable of L is transformed to x by letting L = L*x. Thus, eq 1 reduces to
Employing eq 9 and 10, and 2, the second term of eq 14 becomes
The above equation is formulated based on the fact that
Pi(l)= 1 and Pi(0)=(-l)i. From eq 9 where L* is the maximum crystal within the crystallizer. The population density function will tend to zero for crystal size beyond L*. The corresponding growth rate expression of crystal, as shown in eq 4 or 5, will be
-G(l't) - - G(l,t)(AC)d= ( g g l ) ( A C ) d L*
l=O
Thus eq 17 becomes
_ - %(AC)d L* where
x
where n(1,t) is directly obtained from eq 12, i.e.
(18)
Ind. Eng. Chem. Process Des. Dev., Vol. 23, No. 3, 1984 465 m
n(1,t)= E ( 2 j + 1)aj j=O
Furthermore, substituting eq 12 and 10 into the third term of eq 14 yields
g(2j j=O
+ 1)(AC)daj5 JIP[(x)pj(x)&(x)dx = l=O 0 m
m
Equations 16, 19, and 21 are substituted into eq 11 dai dai-1 dai+l (2i 1)i= 2(2i 1) x (i 1)dt dt dt m BO(t) (-l)iC201j(2i + l ) ( Z j + l)[Cgl(Bij,l-l)](AC)d L* j=o 1=0 (22) where Bij,l is defined as
+ +
+
+
-
+
(23)
A detailed derivation for evaluating Bi.,lin order to save computer time, is given in the Appendix. In order to write eq 22 into a matrix form, a tridiagonal matrix, T ,is defined (Shih, 1971)
T=
[:
2
6
o... o... 3 ...
0
0
0 2
1 3
.........
0
...
0 0 0
0 0 0
N-1
2h'-1
I
Define the following matrix and column vectors A = [Aij] m
Aij = 2(2i + 1)(2j + l)[Cgl(Bf,l- l ) ] 1=0
D = [Di] Dj = 2(-1)'(2i =
[a09
a17
where
t,
is the volume fraction of solution, defined as
+
Bij,l= JIP[(x)Pj(x)Pl(x) dx
1 1
However, as pointed out previously, the supersaturation of solution may be affected by the temperature and the concentration of solute, which exists in the solution during the batch run. Therefore, the combination of heat and mass balance equations with population balance is essential in modeling the batch crystallization process. The crystallization of ice from brine solution in a Couette flow crystallizer is studied. In the crystallization of ice from brine solution, the salt is not incorporated into the crystal, and the conservation of salt can be written as d --(Vt,C,) = 0 dt
+ 1)-L*1 **.?
ffN-1IT
Thus, eq 22 can be written as
T* = Aa[AC(t)ld + DBo(t)
e,(t) = 1 - SmKJ3n(L,t) 0 dL
If the selection of L* value is large enough or the population density function vanishes as L > L*, eq 33 becomes e&)
The inversion of tridiagonal matrix can be obtained by the well-known numerical analysis (Shih, 1971). The initial conditions of a are obtained by making transformation of no(L)to a(0) from eq 3. Therefore, the simulation of batch crystallizationprocess is obtained by solving eq 31, with material and energy equations. In the following section, a typical example of crystallization of ice from brine solution in a Couette flow crystallizer is studied (Wey and Estrin, 1974). Illustrative Example In general, eq 1through eq 6 constitutes the entire description of crystal-size distribution in a batch crystallizer.
(34)
(35)
The energy balance equation for the crystallization of ice from brine solution in a Couette flow crystallizer may be expressed as (Wang, 1974)
The initial condition for Tm(t)is Tm(0)= Ti
(37)
Generally, the energy balance can be obtained by applying the first law of thermodynamics to the batch crystallizer. Equation 36 was set up by taking into consideration the sensible heat of solution, the heat of crystallization, the heat removed by the coolant through the crystallizer, the heat lost to the surroundings, and the power dissiptation due to the agitation of the rotor. For the batch crystallizing system used in this study, Ti and Ci are related by the freezing point curve for the ice-brine system. The supersaturation is thus zero at zero time and develops quickly since Tc < Ti. Applying the technique of Legendre polynomials, eq 34 becomes 1
ev(0) = 1- K&*4$ 0 x3n0 dx
or
dt
= 1 - xL'KvL3n(L,t)dL
The initial condition for Cm(t)is CJO) = ci
dt
-d a-- TIAa[AC(t)ld+ T I D B o ( t )
(33)
K,L*4 cv(t) = 1 - -[a4+ 7a3 20az 70 The energy equation reduces to
+
U84
-(Tar V
-
+ 28al + 14a0]
(38) (39)
P ( K > * 4 ) [ 2- + 7 - + Tm)+ - + PLXR 70 V
daz
20-
dt
+ 28-da1 + 14- dt dt
(40)
Equations 31,39, and 40 are solved simultaneously to give the expansion coefficients, ak. Once 'Yk are obtained, the
466
Ind. Eng. Chem. Process Des. Dev., Vol. 23, No. 3, 1984
A
/ I
0 0 4 7 2 0 0/ 1 7
6000~
N L == 2n 1. 1 1 x
L*=o.::
4800
ti \
io0
L * = 0.11
NTl t
/
10
L/L*
Figure 1. Population balance density function vs. normalized length of crystal.
10
Figure 2. Total number of crystals vs. time.
population balance is obtained from eq 12. The following physical data and operating conditions are adapted from the system (Wey and Estrin, 1973).
Bo = 1.8 X 106(AC)2p2(az= 2; y
30
t. (secl
5000
-
t=30 sec N=26
= 1; j = 2)
N =21 N
=11
U A / V = 0.08 cal/(OC cm3)
U,A,/V = 0.0125 cal/("C cm3) d = l L* = 0.11 cm
G=
+
(0.168 X 10-2L-1 0.192 - 0.57L
+ 0.926L2) (AC) cm/s
C,(O) = 0.03 g/cm3 T,(O) = -1.78 "C I
K , = 1.4 P / V = 2.83
X
cal/(cm3
0
I
I
0.2
I
_.....
OTTY-?.O
"
L/L*
8)
C, = 1.0 cal/(g "C)
I
0.4
Figure 3. Population balance density function vs. normalized length of crystal at t = 30 s.
C, = 0.5 cal/(g "C)
-
N = 2 l
n&)
--_-- N = 1 6
= 1.01852 X 10SL2e-509.26L2
...........
p c = 0.917 g/cm3
N =I1
r
pL = 1.0 g/cm3
XR = 80 cal/(g "C)
T,,= 0.0 "C C , = -0.00826 - 0.0261T- - 0.00278Tm2- 0.000194Tm3 The computational results are shown in Figures 1-4. Results and Discussion In evaluating aifrom eq 31, the leading number of equations to be calculated is determined by i. However, the upper limit of the summation, j , should also be equal to the value of the number of equations. Otherwise, there is great difficulty in getting the solution, i.e., i = j = N 1. Thus,the value of N is the leading number of equations to be calculated. A larger number of leading equations involved in simulation will make the results convergent rapidly or more closely to the exact solution. Figure 1 reveals the relationship between population density
-
6
.
2
-
io
20
30
40
50
t, (sec)
Figure 4. Temperature trajectory of solution vs. time.
function n(x,t) and characteristic size x. The population density function is obtained by solving 21 ordinary differential equations (N= 21) of expansion coefficients; i.e., the first 21 terms of expansion coefficients, ai (ai,i = 0, 1, ..., 20) are calculated. The portion of the curve of Figure 1which includes the main peak indicates that the population density number
Ind. Eng. Chem. Process Des. Dev., Vol. 23,No. 3, 1984 407
is generated due to nucleation. The tailing part a t larger crystal size is due to the growth of seeded crystal. This special phenomenon can be explained by the solution of the population balance equation using the method of characteristics. The solution of the first-order partial differential equation of population balance is closely related to the seeding and nucleation rate expression. In general, the seeds within the crystallizer will grow independently and their size distribution will retain their shape as time proceeds if a simple model is calculated to obtain an exact solution (Chang and Wang, 1982, 1983). Nevertheless, the present model is very difficult to solve due to the complicated forms of growth rate and nucleation rate expressions. In addition, due to epecial oscillating characteristics of higher-order shifted Legendre polynomials, the "tailing part" of the population density function which is orginally contributed by the growth of seeds will also shown oscillatory behavior. NeveFtheless, this oscillation of the tailing part will die out if the number of leading terms of expansion coefficients to be calculated is chosen sufficiently large and there is no effect of round-off error of the computer. The total number of particles NT = JT;'n dL can be obtained with different values of N. Figure 2 shows the computational results for NT vs. t with N = 11 and N = 21. As shown in the figure, eleven ordinary differential equations of expansion coefficient will be enough to describe the system. However, the integration time interval will be chosen to a smaller value when a larger N is chosen. Otherwise it will lead to the divergence of computation. As indicated in eq 6, the original nucleation rate expression involves the integration of the population density function from 0 to m , i.e., JGnL2 dL. However, the upper limit of integration set at infinity is convenient for calculation using the method of moment (Wey and Estrin, 1973). Physically, it is more sujtable to integrate the population density function from zero size to maximum size, L*,within the crystallizer, i.e., JT;'nL2dL. Meanwhile, L* is just used to be a transformation parameter. Using such a transformation as shown in eq 7, the growth rate expression which originally involved a singular point a t zero size now will be nonsingular. Therefore, eq 7 can be solved without difficulty by using the shifted Legendre approximation. In simulating the crystallization process with nucleation rate in terms of moment forms, the L* value is chosen to be equal to the maximum crystal size if the largest crystal within the crystallizer is known. However, if the size of the largest crystql is not possible to determine, an appropriate increase of the value of L* is necessary until there is no significant change of the calculated population balance function. Increasing the value of L* lea& to the rapid coqvergence of the computational results. In the present example, a choice of L* greater than 0.11, which corresponds to the normalized crystal size x = 1,will not affect the calculated value n any further. Therefore, L* is set a t 0.11 for studying the ice cfystal size distribution from brine solution at transient state. As indicated previously, one of the main purposes of the present paper is to provide an effective tool to solve the most general population balance equation. As shown in Figures 1-4,it reveals that the presently proposed method is an effective tool in solving the population balance equation involving complicated moments forms of nucleation rate and infinite growth rate a t zero size. Figure 3 shows the computational results of population balance vs. crystal size a t t = 30 s for N = 11,21, and 26, respectively. As shown in ~e figure, the curve will oscillate
when the smaller values of N are chosen (for example: N = 11). The curve of n will go smoothly for N = 26 (see Figure 4). This indicates that the shifted Legendre polynomial function is adequately convergent. The corresponding temperature trajectory of solution in the crystallizer is also shown in Figure 4. From the figure, it is seen that the solution temperature dropped as time increased up to 20 s. At time, t = 20 s, the solution temperature is -2.3 "C, which is enough to induce nucleation of ice. The design work of the crystallization process is concerned with the nucleation rate and growth rate expression. In general, the crystal growth rate, which is usually subject to heat and mass transfer, is rather well-known. N e v p theless, the theory of nucleation is still in ita infancy. Thus, nucleation rate is usually assumed to be expressed by a simple power law model with operating condition parameters. Therefore, the presently effective approximate method can be used directly to deal with the experimental crystal size distribution vs. size d@a. At the beginning, the parameters of nucleation rates are assumed. Then, the crystallization process is simulated by the population balance equation using the present method. The system parameters are adjusted and are obtained to make the minimum difference between the calculated and experimental population density function.
Conclusion The population balance and the material and energy balance equations, which represent the batch crystallization process, are simulated. The solutions are obtained by employing the shifted Legendre polynomial functions. A generalized size-dependent growth rate expression and moment form of nucleation rate expression are considered in this paper. The proposed method leads one to results which are indeed reasonable as judged by the nature of the recovered distributions. Therefore, the technique of the shifted Legendre function can serve as a powerfql tool in designing the crystallization process. Appendix Properties of Legendre Polynomial Functions. The shifted Legendre function, P~(T), is related to the wellknown Legendre function, Pi(x), by changing the independent variable to x = 27 - 1. Thus, the binomial expansion of the shifted Legendre polynomial function can be written as k = ~k ! [ ( i-
= --[d(~ 1 di i! dTi
k)!I2
(A-1)
The products of the shifted Legendre polynomial function can be expressed in terms of its polynomials (Gradshteyn and Ryzhik, 1977)
where dij,i-j+21 is defined as
and al are a0 = 1
21
al+l
+1
= (-)a1
(A-4a) (A-4b)
Applying the orthogonality property of the shifted Le-
488
Ind. Eng. Chem. Process Des. Dev., Vol. 23,No. 3, 1984
gendre polynomial functions, the integration of a triple products of shifted Legendre polynomial functions is thus obtained; i.e., for j Ii , k # i - j + 21, 1 = 0, 1, ...,j (A-5)
for j Ii , k = i - j
+ 21, 1 = 0, 1, ...,j (A-6)
Define
Cij,k
1
Cij,k
J f j i ( 7 ) f j j ( 7 ) P k ( 7d7)
=
(A-7)
Therefore, the following identities hold C k j , k = 0; k # i - j 21 (I = 0, 1, 2 ,
(A-8)
C i j h = dij,k; k = i
(A-9)
+
-j
...,j ) + 21 ( i = 0, 1, 2 , ..., j )
Using eq A-5 to A-9, the integration of the triple product of the shifted Legendre polynomial functions can, therefore, be obtained. The other property of shifted Legendre polynomial functions is that the derivative of its function can also be expressed in terms of its polynomials (Farrell and Ross, 1963) P:(7) = P’+J7) + 2(2i - 1)Pi-1(7) (A-10) Thus, the other form of the integration of the triple products of shifted Legendre functions is obtained by employing the above equation, i.e. Bf,k = Bi-l,k 2(2i - l ) C i - l j , k (A-11)
+
where Bij,kis defined as
The following simple identities are directly obtained from eq A-11 Boi,k = B’j,k
0; i = 0
=
(A-13) (A-14)
2cOj,k
Bij,o= J1P[(7)Pj(7)d7 = 0 = 1-
(j L i ) (-l)i+j
(j < i) (A-15)
Using eq A-11 to A-15, the integration of the triple multiplication of Legendre polynomial functions in which one of them is the derivative can be calculated. Nomenclature A = heat transfer area from coolant to solution A, = heat transfer area from coolant to surroundings B’. = defined in eq A-12 ~ d ” =nucleation rate
Ci = initial salt concentration (brine solution) C, = equilibrium salt concentration, C,(T,) C, = salt concentration C, = specific heat of ice crystal C,, = specific heat of salt solution CiJ,k = defined in eq A-7 d = parameter of growth rate of crystal dij,k= defined in eq A-3 G = growth rate of crystal G(0,t) = zero size growth rate G(1,t) = growth rate of crystal at x = 1 G = defined in eq 9 KG = growth rate constant K, = surface shape factor of crystals K, = volume shape factor of crystals
L = characteristic length of crystals L* = maximum characteristic length of crystals n = population density function no(L)= population density function of seeding pi(x) = shifted Legendre polynomial function t = time T , = coolant temperature T,, = surroundings temperature T , = bulk solution temperature Ti = initial bulk solution temperature T = tridiagonal matrix U = overall heat transfer coefficient from coolant to solution U, = overall heat transfer coefficient from coolant to surroundings Greek Symbols r = independent variable XR = heat of crystallization pL = density of solution p, = density of crystal t, = volume fraction of solution a,&y = parameters of nucleation rate aj = expansion coefficient of Legendre polynomial AC = supersaturation of solution, C, - C, Literature Cited Chang, R. Y.; Wang, M. L. Symposium on Transport Phenomena and AppUcations, Natlonal Talwan Unlversity; Taiwan 1982, 111. Chang, R. Y.; Wang, M. L. ChICM J . 1983, 14, 25. Farrell, 0. J.; Ross, B. “Solved Problems: Gamma 8 Beta Functions, Legendre Polynomials, Bessel Functions”: Macmlllan: New York, 1963; pp 118-234. . . - -.. .
Gradshteyn, 1. S.; Ryzhik, I. M. “Tables of Integrals Series, and Products”; Academic Press: New York. 1977: D 1026. Hlmmelbiau, D. M.; Bischoff, K. B. “Proc‘ess Analysis and Simulation”; Wiley: New York, 1968. Hulburt, H. M.; Katz, S. Chem. Eng. Sci. 1964, 19, 554. Moyers, C. G., Jr.; Randolph, A. D. AIChE J. 1973, 19, 1089. Randolph, A. D.; Larson, M. A. “Theory of Partlculate Processes”; Academlc Press: New York, 1971. Randolph, A. D.; Larson, M. A. A I C M J. 1962, 8, 639. Shih, Y. P. “Numerlcal Analysis”; National Cheng Kung University: Taiwan, 1971; p 40. Wang, M. L. Ph.D. Dissertation, Clarkson C o l l w of Technology, Potsdam, NY, 1974. Wang, M. L.; Chang, R. Y. Chem. Eng. Commun. 1983, 22, 115. Wey, J. S.; Estrln. J. Ind. Eng. Chem. F’rocess Des. D e v . 1973, 12, 236. Wey, J. S.; Estrin, J. Desalination 1974, 14, 103. Wey. J. S. Phofogr. Sci. Eng. 1977, 21, 248.
Received for review June 8, 1982 Revised manuscript received August 15, 1983 Accepted September 12, 1983
