Ind. Eng. Chem. Process Des. Dev. 1986, 25, 547-551
547
New Algorithm for Calculation of Azeotropes from Equations. of State Shao-Hwa Wang and Wallace B. Whltlng' Department of Chemical Engineering, West Virginia Universiw, Morgantown, West Virginia 26506-6 10 1
By treating vapor-liquid equilibrium for the azeotropic point of a mixture as similar to that for a pure compound, a new algorithm for azeotropic prediction has been developed. Methods of azeotropic prediction for muticomponent systems are also derived from Anderson's vapor-liquid equilibrium calculation method and from Teja's constant-temperature binary azeotrope method. Results from calculations for 36 binary systems and one ternary system show that the new algorithm is the most efficient and robust of all. Azeotropic loci of several mixtures that are weli-described by a simple generalized van der Waals equation of state have been calculated successfully with the new algorithm. This new algorithm is useful for rapid detection of azeotropes and is applicable to any equation of state.
Extensive research has been done on predicting vaporliquid equilibria for industrially important systems. Of particular interest to the chemical engineer, however, is the prediction of azeotropes, a special case of fluid-phase equilibrium, which has not received much attention. General phase-equilibrium algorithms often fail near azeotropes; thus, a specialized algorithm is desirable. Such an efficient and robust algorithm is presented here. There exist various kinds of azeotropic mixtures. However, our attention is restricted here to homoazeotropic systems without chemical reactions. For these systems, two methods previously reported in the literature have been modified for calculation of multicomponent azeotropes. We compared these two methods with our new algorithm for 37 systems. Using the principle of corresponding states along with an experimental equation of state for a reference substance, Teja and Rowlinson (1973) have predicted binary azeotropes at constant temperature. The Powell method (1965) was used to minimize an objective function which is the sum of the squares of the differences of the following thermodynamic properties in the two coexisting phases: pressure, partial derivative of Helmholtz free energy with respect to mole fraction of one component, and molar Gibbs free energy. These calculations are restricted to binary mixtures at constant temperature. We extend Teja's calculation to multicomponent systems at both constant temperature and constant pressure. Computation methods for bubble points and dew points have been presented by Anderson and Prausnitz (1980). They illustrated the calculation by applying a generalized canonical van der Waals partition function. On the basis of Anderson's work, two schemes for azeotropic calculation have been derived in this study: (1)by setting the criterion to be equal the mole fractions of each component in both phases and (2) by locating the temperature (or pressure) extremum. However, the latter metthod can be used to predict binary azeotropes only. In the present investigation, we developed a new algorithm for azeotropic calculation, which we tested for a number of binary systems and one ternary system in comparison with the previous methods. To be fair in the comparison,the computer programs for each method have been written to be as efficient as possible and have the
* Author to whom correspondence should be addressed.
same amount of input and output. Furthermore, the same tolerance of convergence is used for each method. The equation of state used is derived from the generalized van der Waals' partition function as illustrated by Anderson. Pure-componentparameters as well as binary interaction parameters used in the calculation were fit to P-V-T data, vapor-pressure data, and binary vapor-liquid equilibrium data by the maximum-likelihood principle. The computer programs used for fitting parameters are those prepared by Anderson (1978). The new algorithm is also used to calculate the azeotropic loci of several binary systems to show its capability.
New Algorithm for Calculation of Azeotropes We use the similarity of an azeotropic system to a pure component to simplify azeotropic calculations and to make them more robust. As a liquid mixture of azeotropic composition is heated at constant pressure, the temperature increases until the boiling point is reached. The entire liquid then vaporizes at that constant temperature, and only after all is in the vapor phase does the temperature continue to rise. This first-order phase transition is identical with that for pure-component vaporization but distinctly different from vaporization of a nonazeotropic liquid mixture. As efficient and robust algorithms exist for calculation of pure-component phase equilibria, we use these methods as a starting point for azeotrope calculations. For a pure component, the equality of the Gibbs free energies in both phases required for phase equilibrium leads to the following expression given in basic thermodynamics texts AL - A" = P(V -
VL)
(1)
dV = P(V" - VL)
(2)
which is equivalent to
LTP
For azeotropes (unlike other multicomponent systems), eq 1 must also hold. The chemical potentials of each component in the vapor phase and in the liquid phase must be equal. Since the compositions are identical in the two phases, this leads to an equality of the molar Gibbs free energies for the two phases. Equation 1is a necessary, but not sufficient, condition for azeotropy. Clearly then, if the constant-composition
0196-4305/86/1125-0547$01.50/00 1986 American Chemical Society
548
Ind. Eng. Chem. Process Des. Dev., Vol. 25, No. 2, 1986
Equation of State Used We have used the generalized van der Waals equation of state of Anderson and Prausnitz (1980) for our comparisons of the methods for calculating azeotropes (eq 6). However, any equation of state can be used with our algorithm
START
1 Choose
Choose TO
r
(or
PO)
i
Calculate
I
where [ b/4V. We use the simple one-fluid mixing rules and the Lorentz-Berthelot combining rules NeN.
(7)
AV(T,Vv,z) and AL(T,VL,,)
(or P) Check (A~-A')
e
= p(vV-vL)
No
Secant iteration
~
where aiand bj have the following temperature dependence
Calcurate
( T , V V , z ) and q ( T , V L , z ) - I -
Secant iteration
YES STOP
Figure 1. Constant-pressure (or temperature) azeotropic calculation by the new method.
isotherm for a mixture has no van der Waals loop, there can be no azeotrope. The other requirements for azeotropy are
f? = f:
(3)
for i = 1,2, ..., N , - 1. We need not solve eq 2 for i = N , as this is guaranteed by eq 1. Equation 3 can be rewritten in terms of easily calculated properties as YIP# = &P.,
(4)
For an azeotrope, x, = y, for all i; thus, we use the following set of equations to replace eq 4 in our algorithm In [#/&
=O
(5)
for i = 1,2,3, ...,N, - 1. Fugacity coefficients (4,) are easily calculated from equations of state (see, e.g.: Prausnitz, 1969). Figure 1 shows the flow diagram for azeotropic calculation with this new algorithm. From an equation of state, the molar volume, Helmholtz free energy, and fugacity coefficients for each phase are calculated. The secant method is used to compute both the equilibrium temperature (or pressure) and the azeotropic composition, since it is less sensitive than the Newton-Raphson method to rapid changes in the slope of the function at points far from a converged solution. An IMSL subroutine for the secant method known as ZSCNT was used for multicomponent calculations. A primary advantage of our method is that we avoid spurious roots. We first solve eq 2. If there is no van der Waals loop, we know immediately that there can be no azeotrope. If there is a loop, we have realistic densities for the two phases with which we test whether eq 5 is satisfied. Our algorithm never searches in nonphysical regions because we constrain the search to regions in which eq 2 is satisfied.
We used the maximum-likelihood method for all parameter estimation (Anderson et al., 1978) to account properly for the experimental uncertainty in the experimental data used. For pure-component parameters ( a ,0, y,and 6), we used vapor-pressure and dense-fluid P-V-T data. Primarily, binary interaction parameters (k,) were obtained from the correlation of P-T-x-y data. However, the k,'s of some mixtures are so strongly temperature- and pressure-dependent that the prediction of azeotropes for different conditions is impossible with only a single value of k,. To compare the azeotropic calculation methods, we have used the k, esimated from the specific azeotropic point to be calculated for those mixtures for which we experienced difficulty in calculation with the k , estimated from the P-T-x-y data. The resulting parameters are given in Table I. These parameters and the parameters for pure fluids given by Anderson and Prausnitz (1980) were used in our calculation.
Comparison of Azeotropic Calculation Methods The highlights of the four methods that we compared can be summarized as follows: New. The azeotrope is treated as a pure component. Anderson-I. Normal phase-equilibria calculations are done until the two phases have the same composition. Anderson-11. Normal phase-equilibria calculations are done to find the temperature or pressure extremum. Teja. Azeotropic criteria are solved simultaneously by varying phase densities and composition. For every method, the tolerance of convergence is chosen as 1 x lo4. To have a reasonable initial guess of pressure for the constant-temperature azeotropic calculation, we compute the vapor pressures for all the pure compounds and choose the estimated pressure to be the average of the lowest and highest vapor pressures. In an analogous manner, the initial guess of temperature is chosen for the same azeotropic calculation. For a binary system, the initial guess for the azeotropic composition was z1 = 0.5 and z 2 = 0.5. The ternary azeotropic composition is guessed as tl = 0.3,z2 = 0.3, and z 3 = 0.4. The same initial guesses and tolerance of convergence are used to predict azeotropes so that we can compare the computer time consumed and the robustness for all the methods of azeotropic calculation.
Ind. Eng. Chem. Process Des. Dev., Vol. 25, No. 2, 1986 549
Table I. Binary-Interaction Parameters svstem k;; waterlacetic acid -0.05033 -0.027 24 waterlethano1 -0.094 27 water/acetone water / phenol -0.154 48 -0.008 77 water / propanol 0.075 75 water/ethyl acetate water / butanol -0.003 80 waterlpyridine -0.052 89 0.074 67 hydrogen sulfidelethane -0.11174 chloroform/ethanol -0.269 06 chloroform/acetone -0.040 13 chloroform/ethyl acetate -0.424 18 chloroform/ hexane 0.056 74 nitromethanelethanol -0.027 16 nitromethanel benzene -0.041 54 acetic acidlbutanol 0.098 77 acetic acidlbenzene -0.090 38 ethanollethyl acetate 0.062 35 ethanol/ pentane -0.046 98 ethanol/ benzene -0.014 59 ethanol/cyclohexane 0.018 93 ethanol/ hexane 0.015 01 ethanol/ heptane 0.034 60 ethanol/octane 0.083 74 acetonelpentane -0.058 37 acetone/cyclohexane -0.057 62 acetone/ hexane -0.14094 propionic acidloctane 0.11925 propanol/ benzene 0.134 19 propanol/cyclohexane 0.109 59 propanol/ hexane 0.003 16 propanol/ heptane -0.302 58 ethyl acetate/cyclohexane 0.065 98 butanol/ benzene 0.10460 butanol/cyclohexane 0.108 48 butanol/ hexane -0.075 54 butanol/ heptane -0.037 47 butanolldecane -0.017 11 benzene/cyclohexane 0.015 17 benzene/ hexane Benzene/ heptane 0.007 28
no. of Doints 10 12 13 14 23 5 20 13 13 15 9 18 9 18 12 18 8 14 19 12 13 9 23 17 11 9 14 12 20 10 12 21 13 9 16 21 20 19 21 15 30
T , O C , or P. torr PI760 PI745 TI25 TI75 PI760 PI760 PI760 TI50 TI50 TI45 PI760 PI760 PI760 PI730 TI45 TI117 PI758 TI60 T/-10 P/ 180 TI30 TI25 PI760 TI45 PI760 PI760 TI45 PI750 TI65 T/55 PI760 TI75 PI760 TI45 TI70 PI760 P/ 1445 TI100 P/759 PI760 PI720
ref Rivenc, 1953 Svoboda et al., 1968 Beare et al., 1930 Ito and Yoshida, 1963 Murti and van Winkle, 1958 Garner et al., 1954 Orr and Coates, 1960 Ibl et al., 1956 Robinson and Kalra, 1975 Scatchard and Raymond, 1938 Karr et al., 1951 Nagata, 1962 Kudryavtseva and Susarev, 1963 Nakanishi et al., 1969 Brown and Smith, 1955 Rius et al., 1959 Othmer, 1928 Murti and van Winkle, 1958 Isii, 1935 Nielson and Weber, 1959 Scatchard and Satkiewicz, 1964 Smith and Robinson, 1970 Raal et al., 1972 Lu and Boublikova, 1969 Lo et al., 1962 Marinichev and Susarev, 1965 Schaefer and Rall, 1958 Johnson et al., 1957 Brown and Smith, 1959 Strubl et al., 1970 Prabhu and van Winkle, 1963 Lee and Scheller, 1967 Carr and Kropholler, 1962 Brown and Smith, 1959 Vonka et al., 1971 Govindaswamy et al., 1976 Vijayaraghavan et al., 1967 Lee and Scheller, 1967 Richards and Hargreaves, 1944 Ridgway and Bulter, 1967 Hlousek and Hala, 1970
Table 11. Summary of Binary Azeotropic Calculation constant-T calculation constant-P calculation CPU time CPU time av CPU time method no. of failures mean, s re1 no. of failures mean, s re1 total no. of failures mean s re1 0.099 1.00 0 0.098 1.00 0 1.00 new 0 0.097 0.138 1.39 0 0.144 1.47 Anderson-I 0 0.149 1.54 0 0.459 4.64 40 0.472 Anderson-I1 23 0.493 5.08 17 4.82 1.723 17.40 10 1.111 11.34 Teja 2 0.591 6.09 8
Azeotropes were calculated for 36 binary systems and one ternary system with all four methods. Since the same equation of state is used, the calculated azeotropes are the same for all the algorithms. Because of the choice of k,, the calculated results agree with the experimental data within experimental error. Both constant-temperature and constant-pressure azeotropic calculations were executed for each system. The detailed results of each calculation are given elsewhere (Wang, 1984). The CPU execution time required for each method per azeotropic calculation is summarized in Table 11. The number of convergence failures for each method is also shown, which can be used to infer the robustness of each method. A FORTRAN G compiler was used under an M V S operating system on an Amdahl470 V/7A computer. The source code for all computer programs is available from the authors. Only the new method and the Anderson-I method were free from failures, in our experience. For every azeotropic system, the new method took less CPU time. On the average, the Anderson-I method took 54% more CPU time
for constant-temperaturecalculations and 39% more CPU time for constant-pressure calculations.
Azeotropic Locus Calculation One particular aspect of azeotropes which is of interest to chemical engineers is the difficulty encounted in the distillation of azeotropic mixtures. Since most synthetic organic liquids of commerical interest are purified by distillation, a significant concern is whether a mixture will run into its azetropic locus during this operation. There is only one degree of freedom for an azeotropic state; in other words, all the azeotropes of a system fall on a curve. This uniform-composition curve can be plotted on a P-T diagram regardless of the number of components. Below we present several representative azeotropic loci, including the experimental data and the results calculated by the new algorithm for azeotropic calculation. In Figure 2, the predicted results for the carbon dioxidelethane system show good agreement with experimental data, and the end point of the locus is very close to that predicted by Teja and Rowlinson (1973). The kil
Ind. Eng. Chem. Process Des. Dev., Vol. 25, No. 2, 1986
550
75
3
I
P
I
I
A
predicted areocropic
55
P
-
5 220
240
260
280 1,
-
I
J
300
320
K
Figure 4. Azeotropic locus of the water/ethanol system.
Figure 2. Azeotropic locus of the carbon dioxide/ethane system.
4 60
I P s 1 . 0 0 6 6 3 bars
180
420
1
P
\
K
b
'..
I
e x t e n e search
I
p r e d i c t e d azeotropeT
r
120
-
60
-
380
340
1, K
Figure 3. Azeotropic locus of the water/acetone system.
used in the calculations is estimated from a single azeotropic point at 216.45 K and 7.31 bar. Figure 3 shows the azeotropic locus for the water/acetone system calculated with the new algorithm and the critical locus calculated with Heidemann and Khalil's method (1980). The predicted azeotropic end point is not in accord with the experimental data, but it ends at the critical region like the experimental results. The k, used is obtained by fitting the azeotropic point at 397.25 K, 6.87 bar. The azeotropic locus of water/ethanol in Figure 4 ends at 540.00 K, 102.11 bar. Again, the k, used in the calculation is estimated from the single azeotropic point at 351.45 K, 1.0113 bar. This system shows good agreement between predicted results and experimental data. Spurious Roots Most failures in azeotropic calculation with the Anderson-I1 method (locating the temperature or pressure
300 (
I
0.1
I
I
0.4 0.6 Mole Froction
I
1
0.8
1.0
Figure 5. Temperature-composition curve with spurious roots for the binary system water/butanol.
extremum) arise from the encounter of spurious roots. Figure 5 shows a case in which spurious roots instead of the azeotropic points are located as the temperature or pressure extreme. Similar failures are widely reported for traditional vapor-liquid equilibriumcalculations. The new method avoids such problems because it searches for the azeotrope directly. It does not search in regions where vapor and liquid compositions are not identical. Conclusions A new algorithm has been developed to predict homoazeotropes with equations of state. It is robust, converges rapidly, and generates its own initial estimates. Rapid convergence is obtained by regarding the azeotropic point as a pure compound. Here, by way of illustration, the new method has been compared (for 36 binary systems
Ind. Eng. Chem. Process Des. Dev., Vol. 25, No. 2, 1986 551
and one ternary system) to methods derived from Anderson and Prausnitz's (1980) vapor-liquid equilibrium calculation and from Teja and Rowlinson's (1973) binary azeotropic calculation. The new method is clearly superior. It is also used to predict the azeotropic loci and the azeotropic end points of several binary systems with satisfactory results. This new algorithm can be used for constant-temperature as well as constant-pressure azeotropic predictions. For a constant-pressure azeotropic calculation, the initial estimate of temperature is chosen as a value lying among the lowest and the highest of the pure-component saturation temperatures. In like manner, the initial estimate of pressure is chosen for a constant-temperature azeotropic calculation. The initial guesses of azeotropic composition, however, are picked rather arbitrarily as l/Nc for every component of an N,-component system. An approximate model for fluids, the generalized van der Waals equation of state, has been used to demonstrate the azeotropic calculations. It requires only pure-component parameters and binary interaction parameters to predict the azeotropes of multicomponent systems. A parameter-estimation process based on the maximumlikelihood principle is used. The pure-component parameters are reckoned by fitting the equation of state to vapor-pressure data and dense-fluid P-V-T data. Similarly, the binary-interaction parameters are estimated from vapor-liquid equilibrium P-T-x-y data. This simple equation of state is not good enough to describe every kind of fluid mixture; nevertheless, it is sufficient for our present illustrative purposes. The results of azeotropic calculations show that the computer time required by other methods is 47%-1135% greater than that required by the new method. As discussed by Leesley and Heyen (1977), 75% of the computer-aided process-design cost is due to the evaluation of thermophysical properties. This efficient and robust new algorithm, therefore, should be useful in process design that involves azeotropic predictions.
Acknowledgment This work was supported, in part, by a grant from the West Virginia Energy Research Center. Nomenclature
1
A = molar Helmholtz free energy a = equation of state constant for mixture ai = equation of state constant for pure compound i b = equation of state constant for mixture bi = equation of state constant for pure compound i f i = fugacity of component i G = molar Gibbs free energy Ki= equilibrium ratio ( K factor) for component i k. = binary-interaction parameter dC= number of component in a mixture P = pressure R = gas constant T = absolute temperature TR,= reduced temperature of component i V = molar volume x i = liquid mole fraction of component i y i = vapor mole fraction of component i zi = either vapor or liquid mole fraction of component i
Greek Letters ai,Pi,yi, 6i = characteristic parameters in the equation of state
= chemical potential of component i
t: = reduced density
= molar density & = fugacity coefficient of component i
p
Superscripts
L = liquid V = vapor
Literature Cited Anderson, T. F. Ph.D. Dissertatlon, University of California, Berkeley, 1978. Anderson, T. F.; Prausnitz, J. M. Ind. Eng. Chem. Process Des. Dev. 1980, 19. 1, 9. Anderson. T. F.; Abrams, D. S.;Grens, E. A. AIChE J . 1978, 2 4 , 20. Beare, W. G.; McVlcar. G. A.; Fergusson, J. B. J . Phys. Chem. 1930, 34, 1310. Brown, I.; Smith, F. Aust. J . Chem. 1955, 8 , 501. Brown, I.; Smith, F. Aust. J . Chem. 1959, 12, 407. Carr, A. D.; Kropholler, H. W. J. Chem. Eng. Data 1962, 7 , 26. Garner, F. H.; Ellis, S. R. M.; Pearce, C. J. Chem. Eng. Sci. 1954, 3 , 48. Govlndaswamy, S.;AndLappan, An.; Lakshmanan, Sm. J. Chem. Eng. Data 1976, 21,366. Griswold, J.: Wong, S. Y. Chem. Eng. Prog., Symposium Ser. 1952, No. 3 , 48, 18. Gugnoni, R. J.; Eldridge. J. W.; Okay, V. C.; Lee, T. J. AIChE J . 1974, 20,
357. Heidemann, R. A.; Khalil, A. M. AIChE J . 1980, 26, 769. Hlousek, K.; Hala, E. Collect. Czech. Chem. Commun. 1970, 35, 1030. Ibl, von N.; Dandliker. G.; Trumpler, G. Chem. Eng. Sci. 1956, 5 , 193. Isii, N. J . SOC. Chem. Ind. Jpn. 1935, 38, 705. Ito, T.; Yoshida, F. J . Chem. Eng. Data 1963, 8 , 315. Johnson, A. I.; Ward, D. M.; Furter, W. F. Can. J . Techno/. 1957, 34,514. Karr, A. E.; Scheibel, E. G.; Bowes, W. M.; Othmer, D. F. Ind. Eng. Chem. 1951, 43, 961. Kudryavtseva, L. S.;Susarev, M. P. Z h . Prikl. Khim. 1983, 36, 1231. Kuenen, J. P. "Theorie der Verdampfung und Verflussigung von Gemischen und der Fraktlonierten Destillatlon"; Leipzig, Germany, 1906; p 111 Lee, L. L.; Scheller, W. A. J . Chem. Eng. Data 1967, 72, 497. Lo, T. Ch.; Bieber, H. H.; Karr, A. E. J . Chem. Eng. Data 1962, 7 , 327. Leesley, M. E.; Heyen, G. Comput. Chem. Eng. 1977, I , 109. Lu, B. C.-Y.; Boublikova, L. J . Appl. Chem. 1969, 19, 89. Maiesinska, W. Bull. Acad. Pol. Sci., Ser. Sci. Chim. 1964, 12, 853. Marinichev, A. N.; Susarev, M. P. Z h . Prikl. Khim. 1965, 38, 378. Murti, P. S.;van Winkle, M. Chem. Eng. Data Ser. 1958, 3 , 72. Nagata, I. J . Chem. Eng. Data 1962, 7 , 367. Nakanishi, K. R.; Toba, R.; Shirai, H. J. Chem. Eng. Data 1969, 2, 4. Nielson, R. L.; Weber, J. H. J . Chem. Eng. Data 1959, 4 , 145. Orr, V.; Coates, J. Ind. Eng. Chem. 1960, 52,27. Othmer, D. F. Ind. Eng. Chem. 1928, 2 0 , 743. Othmer, D. F.; Chudgar, M. M.; Levy, S. L. Ind. Eng. Chem. 1952, 44, 1872. Otsukl, H.; Williams, F. C. Chem. Eng. Prog., Symp. Ser. 1953, No. 6 ,49,
55. Powell, M. J. D. Comput. J . 1965, 7 , 303. Prabhu, P. S.;van Winkle, M. J . Chem. Eng. Data 1963, 8 , 210. Prausnitz. J. M. "Molecular Thermodynamics of Fluid-Phase Equilibria"; Prentice-Hall: Englewood Cliffs, NJ, 1969. Raal, J. D.: Code, R. K.; Best, D. A. J . Chem. Eng. Data 1972, 1 7 , 211. Richards, A. R.; Hargreaves, E. Ind. Eng. Chem. 1944, 35, 805. Ridgway, K.; Butler, P. A. J . Chem. Eng. Data 1967, 72, 509. Rius, A.; Otero, J. L.; Macarron, A. Chem. Eng. Sci. 1959, 10, 288. Rivenc, G. Mem. S e w . Chim. €tat. 1953, 38, 311. Robinson, D. B.; Kaka, H. GPA Res. Rep. 1975, RR-15. Scatchard, G.; Raymond, C. L. J . Am. Chem. SOC. 1938, 60, 1278. Scatchard, G.;Satkiewicz, F. G. J . Am. Chem. SOC.1964, 8 6 , 130. Schaefer. K.; Rall, W. 2.Elektrochem. 1958, 6 2 , 1090. Smith, V. C.; Robinson, R. L., Jr. J . Chem. Eng. Data 1970, 15,391. Strubl, K.; Svoboda, V.; Holub, R.; Pick, J. Collect. Czech. Chem. Commun. 1970, 35, 3004. Svoboda, V.; Hynek, V.; Pick, J. Collect. Czech. Chem. Commun. 1968, 33, 2584. Teja, A. S.; Rowllnson, J. S. Chem. Eng. Sci. 1973, 28,529. Vonka, P.: Svoboda, V.; Strubl, K.; Holub, R. Collect. Czech. Chem. Commun. 1971, 36, 18. Vijayaraghavan, S.V.; Deshpande, P. K.; Kuloor, N. R. J . Chem. Eng. Data 1067 - - - - , 12. 13 Wang, S.-H. M.S.Ch.E. Thesis, West Virginia University, Morgantown, WV, 1984.
--.
Received for review April 22, 1985 Revised manuscript received September 11, 1985 Accepted October 15, 1985
