Comparison of simulated crystal size distribution control systems

Comparison of simulated crystal size distribution control systems based on nuclei density and supersaturation. Ronald W. Rousseau, and Thomas R. Howel...
0 downloads 0 Views 665KB Size
Ind. Eng. Chem. Process Des. Dev. 1962, 27, 606-610

606

Subscripts j = tray k = tray k,j = influence factor of the specification in the tray k for the variable (liquid or vapor) in the tray j

Hlrose. Y.; Kawase, Y.; Funda, 1. Ind. Eng. Chem. Process D e s . Dev. 1980, 19. 505. Holland, C . D. “Muitlcomponent Dlstlllation”; Prentlce-Hall: Englewocd Cliffs, NJ, 1983; Chapter 3. Wang, J. C.; Henke, G. E. Hydrocarbon Process. lW8, 4 5 , 155.

Received for review March 9, Revised manuscript received April 12, Accepted April 26,

Literature Cited Himmeblau, D.M.; Blschoff, K. 8. “Process Analysls and Simulation”; Wlley: New York. 1968; Chapter 8.

1981 1982 1982

Comparison of Simulated Crystal Size Distribution Control Systems Based on Nuclei Density and Supersaturation Ronald W. Rousseau” and Thomas R. Howell Department of Chemical Engineering, North Carolina State University, Ralelgh, North Carolina 27650

A numerical crystal size distribution (CSD) simulator was developed to test candidate schemes for use in controlling CSD produced in continuous crystallizers. The system model uses a solute balance coupled with a population balance to simulate systems with class I or class I1 behavior. Kinetic data from the literature were used for the simulation of a KCI crystallizer;the simulated behavior of the KCI crystallizer corresponded closely to that observed experimentally; cycling of the product CSD was noted at conditions similar to those found in earlier studies. A control system uslng supersatuniton as the measured variable and flow thrwgh a fines trap as the coniroUed variable was found to be effective In reducing the amplltude of CSD limit cycles in the unstable KCI crystallizer. Potential advantages of this control system are contrasted to difficulties associated with a control system based on measurement of nuclei population density.

Introduction Stability of crystalsize distributions (CSD) in continuous crystallizers with complex flows is not guaranteed, even for steady-state material and energy inputs. For over 20 years, long-term CSD oscillations have been observed in industral ammonium sulfate, ammonium nitrate, and potassium chloride crystallizers (Randolph, 1980; Randolph and Larson, 1962; Saeman, 1956). Control of unstable crystallizen, however, has remained an intractable problem as no CSD control system has yet been successfully operated on a laboratory or industrial crystallizer. Many CSD control systems have been proposed and shown to be effective by population balance analysis. The primary difficulty with all of these control schemes has been that they require an on-line measurement of some property of the crystal size distribution. These measurements have proved to be surprisingly hard to make (Parks, 1979; Rovang and Randolph, 1980). In order to make a control system for a real crystallizer possible, the goal of this study is to devise and evaluate alternative control systems not based on CSD measurements. The R - 2 Crystallizer The R-Z model is a useful representation of industrial crystallizers that expands the MSMPR concept to include process configurations designed to improve CSD. Three additional functions are included in the R-Z model, as shown in Figure 1: clear liquor advance, fines destruction, and classified product removal. The clear liquor advance stream is taken from a quiescent portion of the crystallizer called the fines trap. This gives different residence times for the liquor and solute and results in a higher solids content in the slurry, thereby reducing the volume of liquor that must be handled by a solids recovery unit. A fines destruction loop takes a portion of the clear liquor advance stream, heats it to dissolve the small 0196-4305/82/ 1 121-0606$01.25/0

crystals which are removed with the clear liquor, and returns it to the crystallizer. By selectively removing crystal nuclei the crystallizable solute is forced to grow on fewer crystals, giving a larger average crystal size. In the product classifier, large crystals are allowed to settle out and are removed with the otherwise well-mixed product. A product classifier acts in conjunction with a fines destruction system and clear liquor advance to reduce the variance of the CSD by reducing the number of smaller and larger than average crystals. A dynamic population balance for an R-2 crystallizer is given by (Randolph and Larson, 1971) dn/dt + G(dn/dL) + h(L)n/7 = 0 (1) where R =

+ Qo)/Qo, 2 = Qz/Qo, and h(L) = R (for 0 < L 5 L,) h(L) = 1 (for Lf < L 5 L,)

(QR

h(L) = 2 (for L, < L ) The effect of the product classifier and fines trap is to give three different residence times in the crystallizer. For crystals below the fines cut size the residence time is r/R. For crystals above the product cut size the residence time is T/Z. The steady-state CSD for an R-2 crystallizer is shown in Figure 2. The analytical expression for this CSD is given by the set of equations n = no exp(-RL/GT) (for 0 < L 5 LF) n= no exp(-(R - ~]LF/GT) exp(-L/Gr) (for LF < L 5 L,) (2) n= no exp(-{R - ~ ) L F / Gexp((2 ~ ) -1)L,/G7) exp(-ZL/G.r) (for L, < L ) 0 1982

Amerlcan Chemical Society

Ind. Eng. Chem. Process Des. Dev., Vol. 21, No. 4, 1982

_ _ - - - -

-Fe,

Fines Destruction

I

System

7

‘R

Qcla

I

,

Clear Liquor Advance

v 01

Product

1

c’ m3cla

2 Fines

Trap

-r Solute Balance Boundary

Crystallizer

- - - - - - Figure 1. Schematic diagram of complex crystallizer,

\

J

-*

I

Product and

~ 7 c r\y s t a l l i z e r

I

I

I

I I

I 0

I

Lf

Crystal

LP

Size

Figure 2. Steady-state CSD from an R-2 crystallizer.

The steady state CSD predicted by this model has been used to describe experimentally determined distributions in a variety of cases (Helt and Larson, 1977; Juzaszek and Larson, 1977; Kraljevich and Randolph, 1978; Parks, 1979). Crystal Size Distribution Instability Instability in continuous crystallizers should be distinguished from CSD transients. Crystallizers exhibit a transient CSD in response to an upset in one or more process variables, such as feed rates, temperatures, or concentrations. Some crystallizers, however, have been observed to develop long-term stable limit cycles of the CSD even with constant operating conditions (Saeman, 1956; Randolph and Larson, 1962; Randolph, 1980). This phenomenon is due to an inherent instability of the CSD. Early mathematical analyses of the stabilty of continuous crystallizers were based on population balances and predicted instability by one of two methods. The first method (Randolph and Larson, 1971, pp 91-99) involved writing the dynamic population balance as a set of coupled ordinary differential equations in terms of the moments of the CSD. The moments of the distribution are given by mi = $L’n(L) dL (3) The coupled set of differential moment equations describes the transient response only of the moments and not of the CSD itself. Linearizing and taking Laplace

607

transforms gives a set of linear, homogeneous differential equations. Stability of this set of equations is determined by application of the Routh test. A second method of stability analysis was presented by Lei et al. (1971a) and was later used by Randolph et al. (1973). In their work, the equations given by the original distributed population balance, a boundary condition from nucleation kinetics, and a solute balance were linearized about their steady states. Stability of this set of linearized equations is then determined by spectral (eigenvalue) analysis. The stability criterion for the system is that the roots of a characteristic equation derived from the linearized balance equations must have negative real parts. Randolph and Larson (1965) showed that for an MSMPR crystallizer with nucleation kinetics of the form Bo = KNG’ (4) a value of i greater than 21 is required for the CSD to become unstable. A system exhibiting this extremely high value for the kinetic exponent has not been observed experimentally, Therefore, in an MSMPR crystallizer, the only way that CSD instability will be observed is if there is a discontinuity in the nucleation process (Randolph, 1980). Lei et al. (1971b) showed that the operation of a crystallizer near a metastable concentration, below which nucleation ceases but growth continues, can lead to CSD instability. Randolph et al. (1973) referred to this form of instability as high order cycling, and they showed that it can also be caused by another kind of nucleation discontinuity. If a crystallizer is operated with a high level of supersaturation, as may occur in systems with a very short residence time, homogeneous nucleation may periodically cause a shower of nuclei. As they grow, these new crystals present a large surface area for growth and cause a reduction in supersaturation. When this shower of crystals is removed by natural flow from the system, the supersaturation again begins to increase, so that the cycle is repeated. Although high-order cycling is possible, it is not thought to be a major cause of CSD instability. Industrial crystallizers are frequently operated with fines traps and product classifiers in order to increase the dominant crystal size and narrow the CSD. Sherwin et al. (1969) showed that these associated flows are likely to be the cause of industrial CSD instability. Their analysis was for an MSMPR crystallizer with an ideal product classifier that instantly removed all crystals larger than some cut size. The CSD was shown to be severely destabilized as the cut size was reduced and approached the average particle size. At small product cut sizes they showed that instability could be observed with a nucleation kinetic order of only 2.0. Real crystallizers, of course, do not have ideal product classifiers. Randolph et al. (1973) showed that a crystallizer, with fines destruction and product classification used together, may have an unstable CSD for reasonable values of the nucleation kinetic order. Both fines removal and product classificationwere shown to lower the kinetic order required for instability; when they were used together, the tendency toward destabilization was greatly increased. Fines destruction reduces the number of crystals in the magma and causes an increase in the average crystal size. For a constant product cut size, fines destruction increases the fraction of crystals removed in the product classifier. Stability analyses have, then, pinpointed two modes of CSD instability: high-order cycling which is caused by discontinuous nucleation phenomena, giving an effective nucleation kinetic order greater than 21, and low-order cycling which is primarily due to product classification. In

808

Ind. Eng. Chem. Process Des. Dev., Vol. 21, No. 4, 1982

industrial crystallizers low-order cycling is thought to be the major source of instability. The control problem examined in this study is addressed specifically to stabilizing low-order cycling. Control System Analysis Control systems for unstable crystallizers have been proposed since Saeman (1956) first reported CSD cycling in ammonium sulfate crystallizers. No control system, however, has yet been successfully operated for long periods of time. Almost all proposed control systems require an on-line measurement of some property of the CSD. While the analytical study of CSD control has kept pace with the analytical stability studies, development of on-line CSD measurement techniques has not been so rapid. Han (1969) proposed a control system to reduce CSD transients due to changes in feed concentration to an MSMPR crystallizer; control called for predictive manipulation of the feed rate. Using an analog simulation of the first four moments of the CSD, stability characteristics of the control system were determined without solving for the complete dynamic CSD. It was concluded that this control system reduced the amplitude of transient effects, but it could not stabilize a system exhibiting high-order cycling. On the basis of analog and digital simulation of CSD moments, Gupta and Timm (1971) proposed a complicated scheme for stabilizing high-order cycling in MSMPR crystallizers. Their system initiated fines destruction whenever the number density (as given by the zeroth moment) was above the steady-state value and seeded the crystallizer whenever the number density was below the steady-state value. This control system was shown to eliminate cycling of the moments of the CSD; unfortunately, even 10 years after their study, no instrument is capable of measuring the total number density in a slurry on a real-time basis. Stability analysis by Lei et al. (1971a, 1971b) showed that high-order cycling in a crystallizer with a fines trap can be eliminated by proportional control of the flow to a fines trap. They based their control scheme on a measurement of the area density (second moment, mz) in the fines loop. Instruments that measure surface area directly, such as those monitoring light scattering, may be of great utility where system conditions allow their use. If, however, the surface area must be calculated from a measured population density, the technique would hold little advantage over those requiring measurement of the slurry number density. Beckman and Randolph (1977) have studied the control of low-order cycling in an R-Z crystallizer. As in previous studies, the authors used linearized stability analysis to evaluate their control system. Uniquely, however, Beckman and Randolph also used a digital simulation of the complete time-dependent CSD. Dynamic CSD simulations were not new (Randolph and Larson, 1965; Timm and Larson, 1968),but the work of Beckman and Randolph was different in three important ways. First, a control loop was included in the simulator. With this loop they were able to change process parameters based on current crystallizer conditions. Second, they were able to simulate low-order CSD limit cycles with a frequency close to that observed experimentally. Third, theirs was the only CSD simulator capable of modeling a class I system. All previous simulators (though not all stability analyses) had assumed class I1 behavior in which supersaturation is constrained to a negligible level. Beckman and Randolph confirmed that proportional control of the fines loop flow rate, based on the deviation

of the nuclei density from its steady-state value, can eliminate low-order cycling in an R-Z crystallizer. Subsequent experimental work (Parks and Rousseau, 1979; Rovang and Randolph, 1980) has been directed toward developing an on-line CSD monitor that can be used to measure the nuclei density. This has proven to be a formidable task. The current state of CSD control practice is not very different from what it was when Saeman first discussed the problem in 1956. Holding all material and energy flows constant is the limit of available control technology. Although the system proposed by Beckman and Randolph seems capable of stabilizing an unstable crystallizer, the required instrumentation has not been proven on a commercial installation. Control System Simulation Crystal size distributions are complex and highly nonlinear functions of material and energy flows in a continuous crystallizer. Consequently, CSD control systems are not easily analyzed using the transfer function approach developed for linear systems. In order to be able to observe the effects of control systems on a crystallizer, a computer simulation called GEM was developed. GEM is very similar to the model CYCLER developed by Beckman and Randolph (1977). Both models solve a population balance, a solute balance, and nucleation and growth kinetics equations simultaneously. Both models permit solute to accumulate in the liquor as increased supersaturation (class I behavior). The primary difference between GEM and CYCLER is that GEM uses a finite difference method for the solution of the hyperbolic, partial differential population balance (eq 11, whereas CYCLER uses a linearized form of an analytical solution to the population balance equaton. There are two main reasons that the finite difference method was chosen for GEM. First, no modification of the main program is necessary to accommodate removal functions other than those described by the R-Z model. The removal rate parameter, h(L),is a vector which is set in a subroutine, and any desired removal rate profile can easily be installed. For example, the GEM simulator can incorporate the behavior of a specific classifier with known characteristics. The second reason that the finite difference method was chosen is that this method is not dependent upon the specific form of the nucleation kinetics. It will require no major rewriting of GEM if a crystal system behaving according to Volmer kinetics, for example, is to be modeled. Control by Supersaturation Measurement Using a linearized stability analysis and their model, CYCLER, Beckman and Randolph (1977) showed that a control system based on measurement of nuclei density should be effective in stabilizing low-order cycling. Their control system manipulated the fines loop flow rate in proportion to the deviation of the nuclei density from its steady-state value using the algorithm ( R - R ) / R = K,(nO - i i O ) / i i O (5) where R is the fines removal rate parameter, no is the nuclei density, and K , is the proportional control constant. The overbars designate steady-state values. Beckman and Randolph showed that a control constant of 0.25 will eliminate cycling in an unstable crystallizer. Randolph et al. (1977) have demonstrated low-order cycling in a laboratory KC1 crystallizer. The GEM simulator was tested by running it with the same kinetics and operating parameters as the crystallizer used in their study. GEM clearly predicted CSD cycling as shown in Figure

Ind. Eng. Chem. Process Des. Dev., Vol. 21, No. 4, 1982 609

A

N

r

I

1000

System

T =ZOO

1200

R

KCI/H20

= 28

1400

1600

1800

21

R M S deviation of L a RMS deviation of uncontrolled L D

0.10

1000

1200

1400

1600

1800

2000

Time (minutes)

Time (minutes)

0.05

1

= 28

Z =7.0

min.

Figure 3. Effect of the nuclei density based control scheme on low-order cycling in a KC1 crystallizer.

0.04 0.00

R

System KCI/H,O

0.15

0.20

0.25

0.

Control Constant, K ,

Figure 4. Effect of the control constant on the mean amplitude of oscillations: control by measurement of nuclei density.

3. It showed no tendency to cycle when modeling a stable KN03 crystallizer. The period of cycling in the unstable case is somewhat longer than observed experimentally (300 min vs. an experimental value of 240 min). The controlled response is, however, exactly as observed by Beckman and Randolph (see Figure 3). A control constant of 0.25 reduces cycling to a very low level. For all practical purposes, the damping observed under this control system can not be distinguished from a completely stable system. The success of the tested control scheme was measured quantitatively by a performance index. The performance index was defined as the ratio of the root mean square of the deviation in the dominant crystal size for the controlled crystallizer to the root mean square of the deviation in the dominant size for the uncontrolled crystallizer. The performance index in Figure 4 is not zero when the control constant is 0.25 primarily because of the offset of the controlled dominant size from its steady-state value. Because nuclei density is so difficult to measure, the control system proposed by Beckman and Randolph (1977) has not yet been successfully operated. The nuclei density is, however, only one of many things in a crystallizer that change, and which might be used as the measured variable

Figure 5. Effect of the supersaturation based control system on low-order cycling in a KCl crystallizer.

for a control system. The nuclei density and the crystal growth rate are the two most important parameters of a CSD. It seems reasonable, then, to consider a control system that would stabilize the growth rate rather than the nuclei density. No instrument measures the crystal growth rate directly. The growth rate is, however, generally considered to be linearly proportional to the solution supersaturation. Proportional control based on the deviation of the supersaturation from its steady state is, then, equivalent to control based on the deviation of the growth rate. This concept was tested by using the control algorithm (R - R ) / R = K,(S - S ) / S (6) in the simulator GEM. It was found that this control strategy reduced the amplitude of low-order limit cycles as effectively as had the nuclei density based control scheme (see Figures 5 and 6). It must be pointed out that KC1 follows class I1 behavior and, therefore, MT does not change measurably. Consequently, the effects of this variable on secondary nucleation kinetics are unimportant in this study. The supersaturation in class I1 systems is, by definition, low enough to be ignored in closing a solute balance. Measurement of supersaturation under these conditions is subject to problems with accuracy and sample integrity. Nevertheless, Helt (1976)and Helt and Larson (1977) have demonstrated the potential of using a differential refractometer for this purpose. Measurements with this instrument are not affected by electrical noise or by the presence of crystals in the slurry. In addition, the uncertainty inherent in the determination of no by extrapolation from CSD data is removed. While, under some situations, supersaturation may be easier to measure than nuclei density, a supersaturation control system is more sensitive to measurement error. Figure 6 shows that the supersaturation control system required a control constant of 2.0 to damp CSD oscillations completely. Thii may be compared with the nuclei density control system which required a control constant of only 0.25 for the same degree of damping. Accordingly, an error in the supersaturation measurement will cause a change in the fines loop flow rate that is eight times as large as the same fractional error in a nuclei density measurement. To observe the effect of measurement error, a random noise component was added to the supersaturation mea-

610

Ind. Eng. Chem. Process Des. Dev., Vol. 21, No. 4, 1982 1.o

RMS deviation of LD RMS deviation Of uncontrolled LD

n X

m

u

0.8

-.

- 'F?) Ats - 5,

S(R

0.6 ..

0.4

o.2 V."

..

1 .

0.0

noise = 10% of measured supersaturation no noise

0.5

1 .o

1.5

2.0

Control Constant, K,

Figure 6. E f f e c t o f t h e c o n t r o l constant o n t h e mean amplitude o f oscillations: c o n t r o l by measurement o f supersaturation.

surement in GEM. This noise is a normally distributed random variable with a standard deviation supplied as data to the program. The noise is generated from a random number by the Box-Muller method (Newman and Odell, 1971). The control algorithm including the noise was ( R - R ) / R = KJ(1 varJS- S ) / S (7)

+

where var = a normally distributed random variable = Rand (-2 In (X))0.5cos (2 a Y) and Rand = the width of one standard deviation of the noise variable expressed as a fraction of the measured supersaturation; X,Y = computer supplied random variables uniformly distributed on [0,1]. Figure 6 shows that a noise component with a standard deviation of 10% of the measured supersaturation significantly impairs the effectiveness of the control system. Experiments by Helt and Larson (1977) showed that supersaturation in an MSMPR KN03 crystallizer is hard to measure within 10%. A crystallizer with fines destruction and product classification will, however, operate with higher supersaturation than an MSMPR crystallizer. The supersaturation in an oscillating crystallizer should, then, be easier to measure. In addition, the system used by Helt and Larson had a refractometer connected to the crystallizer by small lines that carried slurry from the crystallizer to the refractometer. They reported considerable difficulty with line blockage and nonrepresentative sampling. If the refractometer measurement cell could be mounted directly inside the crystallizer, these problems would be minimized. GEM indicates that an average measurement of the supersaturation every 5 min will give good control. If measurements are made every few seconds with a refractometer mounted directly in the crystal slurry, and if these measurements are averaged over 5-min periods, then measurements accurate enough for good control should be possible. Conclusions This study leaves the control of an unstable crystallizer undemonstrated. Experimental work in this laboratory and elsewhere leads to the conclusion that control based on measurement of the nuclei densty using an electronic, zone sensing device, such as a Coulter Counter, is not practical. Orifice blocking and electrical noise make these instruments highly unreliable for continuous, automated

CSD measurement in a crystallizing system. Even if the dynamic CSD can be successfully measured, a control system based on determination of nuclei density must overcome severe difficulties. Beckman and Randolph (1977) have shown that if the nuclei density is known to within 50%, then good control is possible. Growth dispersion or size-dependent growth means that even for steady-state conditions the CSD will be curved at small crystal sizes (Janse and DeJong, 1975; White and Wright, 1971). In this case, linear regression may not give a sufficiently accurate estimate of the nuclei density. While a polynomial fit to such data'is possible, the confidence interval for the intercept becomes very wide. Taking the antilog of the intercept gives a value which may not be within 50% of the actual nuclei density. The simulations presented above suggest that an alternative control system based on measurement of the supersaturation is capable of stabilizing low-order cycling. This control system has the advantage of not requiring any distributed measurements of the CSD. If the assumption of a well-mixed crystallizer is valid, the supersaturation can be determined with a single measurement which can be made anywhere in the crystallizer. Although this control system is considerably more sensitive to measurement error, the ease with which the measurements can be made may lead to a control system that is easier to implement than one based on measurement of the nuclei density. These observations presented must be tempered with the recognition that the nucleation and growth kinetics used in the simulations were of a specific type: namely, those observed experimentally by Randolph et al. (1977) for KCl. As pointed out by Shinnar (1981),the accuracy of the model used to select a control algorithm is important. Systems exhibiting forms of nucleation other than that examined in this study may not be as readily stabilized through supersaturation control. Acknowledgment The financial support of the National Science Foundation through Grant ENG-7822656 is gratefully acknowledged. Literature Cited Beckman. J. R.; Randolph, A. D. AIChE J . 1977, 2 3 , 510. Gupta, G.;Ti", D. C. Chem. Eng. Rog. Symp. Ser. 1971, 67(110), 121. Han, C. D. Ind. Eng. Chem. Process Des. Dev. 1969, 8, 150. Helt, J. E.; Larson, M. A. AIChE J . 1977, 2 3 , 822. Helt, J. E. PhD. Thesis, Iowa State University, Ames, IA, 1976. Janse, A. H.; DeJong, E. J. I n "Industrial Crystallization", Mullin, J. W., Ed.; Plenum Press: New York, 1975; p 145. Juzaszek. P.; Larson, M. A. AIChE J . 1977, 2 3 , 460. Randolph, A. D. AIChE J . 1978, 2 4 , 598. KralJevlch, 2. I.; Lei, S.;Shinnar, R.; Katz. S. A I C M J . 1971a, 17, 1459. Lei, S.;Shinnar, R.; Katz, S. Chem. Eng. Rog. Symp. Ser. 1971b, 67(110). 129. Newman. T. G.; Odeli, P. L. "The Generation of Random Variates"; Hafner Publishlng Company: New York, 197 1. Parks, R. M. Ph.D. Thesis, North Carolina State University, Raleigh, NC, 1979. Parks, R. M.; Rousseau, R. W. JACC Roc. 1979, 819. Randolph, A. D.AIChESymp. Ser. 1980, 76(193), 1. Randolph, A. D.; Beckman, J. R.; KraiJevlch,2. I.AIChE J . 1977, 2 3 , 500. Randolph, A. D.; Beer. G. L.; Keener, J. P. A I C M J. 1973, 79, 1140. Randolph, A. D.; Larson, M. A. AIChE J . 1962, 8 , 639. Randolph, A. D.; La",M. A. Chem. Eng. Rcg. Symp. Ser. 1965, 6 1 ( 5 5 ) . 147. Randolph, A. D.; Larson, M. A. "Theory of Particulate Processes"; Acadmeic Press, New York. 1971. Rovang, .R. D.; Randolph, A. D. AIChE Symp. Ser. 1980, 76(193), 18. Saeman, W. C. AIChE J . 1956, 2 , 107. Sherwin, M. 6.; Shlnnar, R.; Katz. S. Chem. Eng. Prog. Symp. Ser. 1969, 65(951. 75. Shinnar, R. Chem. Eng. Commun. 1981, 9 , 73. D. C.; Larson, M. A. AIChE J . 1968, 14, 452. Ti", Whne, E. T.; Wright, P. G. Chem. Eng. Prog. Symp. Ser. 1971, 67(110).81.

Received for review M a y 4 , 1981 Revised manuscript received M a r c h 25, 1982 Accepted April 22, 1982