Evaluation of Equation of State Constants with Digital CornDuters 1
Y
H. W. BROUGH Shell Chemical Corp., Torrance, Calif.
W. C. SCHLINGER
B. H. SAGE
AND
California Instifute of Technology, Pasadena
4, Calif.
Equations of state appear to be useful in predicting the volumetric and phase behavior of multicomponent mixtures of the lighter hydrocarbons. The evaluation of the constants for each of the several components is a tedious and somewhat time-consuming operation. Automatic digital computing equipment offers a possible means of performing formalized calculations of this type and decreasing the time involved. Commercially available digital computing equipment was applied t o the evaluation of the constants for the Benedict equation of state from volumetric data. Conventional least squares techniques have been utilized in determining values of the seven constants in this equation for methane and for propane. These quantities were calculated for each of three different values of one constant
which could not be conveniently established by least squares methods. The details of the procedures used are available elsewhere and the results are tabulated. The calculations described can be controlled automatically to such a n extent t h a t a material saving in calculating time is realized and the procedures are formalized. When applied to this problem, least squares techniques avoid uncertainties in the determination of the weight t o be ascribed t o any part of the experimental data. The methods may perhaps be extended to include heterogeneous equilibrium without undue increase in complexity. The application of automatic digital computing equipment to the evaluation of the constants of equations of state from experimental data should increase the utility of such methods of presenting volumetric measurements.
H E prediction of the volumetric and phase behavior of T h y d rocarbons is a matter of considerable importance to the petroleum industry. The status of thermodynamic data for the industry has been reviewed (18, 19). For many years it has been known that the simplifying assumptions of perfect gases and negligible volume for the liquid phase which form the basis for Raoult’s law are inadequate for describing the behavior of hydrocarbon mixtures a t elevated pressures. Ideal solutions ( l a ) based upon additivity of volumes of the components also failed to yield a eatisfactory generalization of the volumetric and phase behavior of hydrocarbon mixtures. The more complex equations of state which recently have become available primarily aa a result of the efforts of Beattie and Bridgeman ( 9 ) have eatablished another approach to the prediction problem. The “general limit method” discussed by Beattie (1)&or& a means of estimating both the volumetric and phase behavior of mixtures with an accuracy comparable to that with which the equation of state represents the volumetric behavior of the system. The recent work of Benedict (6, 6) resulted in m equation which accurately describes (10)the volumetric behavior of hydrocarbon mixtures in the gas phase and for the bubble point liquid. AvailFble digital computing equipment can be employed to ascertain the volume, enthalpy, and entropy of phases m well as their composition in heterogeneous equilibria with a minimum of manual effort (8). However, the evaluation of the constants for the Benedict equation from experimental data is a laborious and time-consuming process. Benedict (6) employed manual methods which coupled with experience in the field enabled him to attain a consistent set of constants for the lighter paraffin hydrocarbons. More recently Marchman and Prengle (23, 16) have computed the constants for propene and for mixtures of propene and ethane. The determination of the constants for ell the hydrocarbons of practical interest to the petroleum industry is a task of some magnitude. The present discussion considers the application of least squares techniques ( 1 4 ) by means of an automatic digital computer to the evaluation of constants of the Benedict equiitiori of state for the lighter
hydrocarbons yielding a minimum standard deviation from the experimental data. A modified manual form of such methodR was applied earlier ( 4 ) by Benedict to establish constants for an equation of state for nitrogen. These methods permit differing weights to be given to the several sets of experimental data employed. However, the least squares techniques do not necessarily yield the most probable values of the constants since all of the error is assumed to exist in the compressibility factor or pressure and none is ascribed to the specific volume or temperature. MATHEMATICAL CONSIDERATIONS
For present purposes, the Benedict eauation (6) may be rewritten in the following form:
Z. = 1
.
2442
+ (B O -
Ao
-
d
+ ( b - &) d2 +
The problem involved the determination of the constants so that the sum of the squares of the errors was a minimum. This error was taken as the difference between the compressibility factors determined from Equation 1 and the corresponding experimental values. This choice of variable was made because of the relatively small range in 2 as compared to other variables such as pressure or specific weight. The error may be defined for each experimental state in the following way: Zn =
Zo,* -
z,
(2)
The subscript, n, refers to the nth measurement for a particular substance under consideration. When the following summation assumes a minimum value, the corresponding best values of y and Ki are established. n=N
z; = n =1
x2
(Y,
K,)
(3)
INDUSTRIAL AND ENGINEERING CHEMISTRY
November 1951
The quantity Kc in Equation 3 is defined by rewriting Equation 1 in the form Zc = 1
+
i=6 cpi
(d, T ) K ,
+
cp:
K7 (1
+ y dB)
e-72
(4)
i=l
2443
the two derivatives aa identical a t the minimum aa long aa Equation 7 applies. The probable effect of y upon the values of Xz(y) is shown in Figure 1. This information waa estimated from consideration of the possible solutions to Equation 9. Curve A of this diagram corresponds to negative value8 of
The conditions which are desired for the quantities, y and K,, are defined by the following implicit expressions:
Equations 5 and 7 indicate that the sum of the squares with reEpeCt to the variation of y and K i is an extreme value, whereas Equations 6 and 8 ensure that the extreme value is a minimum. The following approach to the problem is necessary because of the yd2) e - y d ’ in Equation 4 . nonlinear nature of the term (1 If such a nonlinear term did not exist, the problem would reduce to the least squares consideration of linear equations which has been treated elsewhere (14). Differentiation of Equation 3 and combination with Equations 2 and 4 yield the following expression:
-+
I
I C 1
I
NEGATIVE
0
POSITIVE
Y Equation 9 which is implicit in both 4 and K i is not suitable for direct calculation. The exponential form of this expression arises from Equation 1 and adds materially to the effort associated with the calculations. A detailed derivation of these expressions is available (7). If it is desired to minimize X2(y, K i ) with respect to y, Equations 9 and 5 may be combined.
Figure 1. Possible Effect of y on the Sum of the Squares of the Deviations
d P X 2 ( yin ) the vicinity of y = 0 and should be the normal behavior
dr2 encountered. Curve B represents the situation for a small psitive value of da X z ( y )in the vicinity of y = 0. In both casea A drZ and B, there is only one minimum for the positive values of y. The third case corresponds to situations where there exist positive
By a choice of trial values of y it was found that the corresponding values of Ki( 7) which would yield a minimum of Xz(y, K i ) for all values of Ki could be astabliahed. The substitution of theee calculated values of K , ( y ) in Equation 10 enables a new trial value of y to be obtained. However, owing to the computational difficulties associated with the solution of Equation 10 for 7 , it waa not employed in the present work. Equation 9 is particularly valuable for investigating the interrelation of X2(y,K i ) and y. However, the computational effort associated with such an investigation is significant. It can be shown along the path described by Equation 7 that
In this equation Ki(y) refers to the values of Ki which satisfy Equations 7 and 8 for a given value of y . It is possible to treat
z X z ( y )near y = zero. The latter situation is the devalues of ddr2 generate case with no niinunum except when y = zero. Conventional least squares methods were employed to evaluate Ki since Equation 7, when combined with Equations 2 and 4, yields the following set of seven linear simultaneous equations corresponding to each of the seven values of cp7:
i=1
n-1
This set of equations may be solved for values of each of the Zci’s corresponding to the chosen value of y. From a knowledge of the values of K i it is possible to determine the magnitude of X z ( y ) . By substitution Xz(y) may be related (7) t o the compressibility factor and the yalues of K , in the following way:
2444
INDUSTRIAL AND ENGINEERING CHEMISTRY
Experience has indicated that Equation 13 is a much simpler means of evaluating Xz(y) than Equation 3. However, in the present work both Equations 3 and 13 were employed for the majority of the states investigated since the actual values of the error, zn, were desired. In order to establish that Equation 8 is satisfied by the values of Ki obtained from Equation 12 the following expression was derived :
Since the value (pi is a real number in all cases the eauation indib W Y , Ki) cated that must be greater than zero. The effect of 3K.T small changes in Ki may be determined from the integrated form of Equation 14. XYY, Ki) = X Y Y )
Vol. 43, No. 11
digits carried in the present calculations inadequate for representing the experimental data with desirable accuracy a t all times. CALCULATED RESULTS
Data are available concerning the volumetric behavior of propane. For present purposes the measurements of Beattie, Kay, and Kaminsky ( 3 ) were employed because this information was utilized by Benedict (6) in an earlier evaluation of constants for the same equation. The constants were obtained for three different values of y, which were 5% above and below as well as equal to that used by Benedict. The results of these calculations and several measures of the over-all agreement of the predicted and experimental results are presented in Table I. The method used to arrive a t each of the several measures of agreement has been indicated in the following expressions:
+ n=l
(18)
CALCULATIONS n=l
The equipment employed for these calculations has been dcscribed by Eckert (9) and King (11). The basic equipment was an electronic calculating punch Type 604 manufactured by the International Business Machines Corp. Auxiliary equipment included a sorter, collator, tabulator, reproducing punch, and a manual key punch. A discussion of the characteristics of the electronic calculating punch was recently made available (8). When using such equipment it is necessary to choose carefully the calculating methods to avoid loss of accuracy from an insufficient number of significant figures in the result. Thip difficulty was particularly important in the present work because of the large number of significant figures required for each of the constants to obtain volumetric and phase equilibrium data of reasonable engineering accuracy. Throughout most of the calculations the data and final resultc were reported to seven significant figures. However, in portions of the intervening details of the calculations '11 figures were carried. The last four were not considered to be significant because of rounding errors and small differences resulting from the relatively large numbers which were often involved. For each state the coefficients of Equation 4 were computed along with the corresponding value of 2, - 1. All the cross - l ) ]and (2, - 1)2were computed. products [rp,,,(nl,,][p,,,(Z, As indicated in Equation 12 their sums were determined and the sets of simultaneous equations, whose coefficients and constant terms were the sums of the above cross-products, were solved. The normal matrix wae evaluated by using the conventional reduction to ' h i t diagonal form." In these latter operations it was necessary to carry as many significant figures as possible to avoid loss of accuracy by rounding errors. The values were checked by the independent calculation of X 2 ( y ) from Equation 13 and Equation 3. Consistency of the values of X z ( y ) up to 10 significant figures was considered a satisfactory indication of agreement and an absence of gross errors in the evaluation of the K's. For subsequent use, the values of K , were maintained to only eight significant figures. In all the calculations, rounding operations were performed upon the sums rather than upon the individual quantities. The average values of 2% with and without regard to the sign of the deviations were computed. I t is again emphasized that in such least squares procedures in order to avoid loss of accuracy in the end result it is necessary t o carry a larger number of significant digits than first appears worth while. The marked variation of such terms as d6/T makes even the eight
n=N
n=l
Table 1. Comparison of Several Values of the Constants of the Benedict Equation for Propane Function Y
B O
Ao X 10-8 co x 10-9 b a X 10-8 a c x 10-8
XP(?)*
42 IZ/Z
Trial 1 5.946290 1.5913065 28.223008 5 0238934 4,6556428 35 088861 3.0273172 17.865402 0.00127724 0,00839 -0,00018 0.00273 0.00483 0.004730
Trial 2 5.653610 1.4480064 26,060274 5.5362459 5,1574946 42,927363 2.6366092 19.713905 0.00096298 0.00295 -0.00011 0,00229 0.00401 0.0O39lc
Trial 3 5.370930 1.2756728 23.513037 6,1228970 5.7401757 51.846881 2.3217096 21,759351 0.00071161 0 00253 -0.00004 0.00191 0.00332 0.00324c
Benedict" 5.653610 1.559992 25.95383 6.219139 5,782101 57,37567 2 . ii01316 25,30414 0.00936535 0,00919 -0.00043 0 00370 0.00617 0.00408C
See (7).
b See Equation 3 or 14. 0
Data a t d of 6.24 pound-moles per cubic foot omitted.
The average deviation of the relative value of z as defined by Equation 19 also has been computed for the situation where all data obtained by Beattie et al. ( 3 ) a t a molal weight of 6.24 pound-moles per cubic foot were eliminated. The elimination of these data made little difference in the fit of the predicted values from the least squares calculations. However, their deletion materially improved the fit of Benedict's constants (6). In Figure 2 is presented the deviation of the sum of the squares of z as a function of the fraction of Benedict's value for y. From a consideration of the nature of this function which is set forth in Figure 1,it is probable that the minimum deviation may occur a t approximately 80% of Benedict's value for y . The minimum value of y employed in these calculations gave a value of X2(y) nearly twice that estimated by extrapolation of Figure 2. However, the values obtained were only a quarter of the root mean square error found by Benedict. For this reason it appears that material improvement in the fit of the Benedict equation may be realized. A detailed comparison of the agreement among the compressibility factors computed from the four sets of constants in Table
I N D U S T R I A L A N D E N G I N E E R I N G CHEMISTRY
November 1951
I indicates only slight differences in the agreement with experimental data except at high densities and low pressures. However, at all temperatures the least squares values yielded smaller average errors in the calculated compressibility factors than were obtained from the constants determined by Benedict. These calculations were baaed only upon the measurements of Beattie (3)for propane and somewhat different values may be obtained if
Table 11. Function
df: 3
16
lKi
IdZI
12
0 *
CI
b
wX 4
I 0.4
L
0.6
I
I
I
08
I .o
I .2
RATIO
p~
Figure 2. Effect of y on the Sum of the Squares of the Deviations of Compressibility Factors for Propane
a
Comparison of Sets of Constants for the Benedict Equation for Methane
Trial 1 Trial 2 Trial 3 1 ,4646060 1.6187740 1.6416900 0.70336500 0.70223723 0.69791469 7.386808 7.333657 7.858538 0.1307256 0.1464787 0.1709262 0.57060456 0.57487211 0.58476176 -0.5152522 -0:4184701 -0.2407344 1.7897648 2.2371713 4 1126911 -0.2134281 -0.1728780 -0.1181582 0.00077752 O.OOO78564 0.000791816 0.00343 0.00345 0.00346 -0.OOO25 -0.00026 -0.00026 0.00278 0,00278 0.00276 0.00243 0,00242 0 00240
Benediot 1.54169oO 0.6828612 7.004694 0.2761352 0.8684956 2.990183 0.5122074 0.4991189 0.00485803 0.00858 0.00176 0.00640 0.00541
surprising sinoe a somewhat smaller number of points than waa used by Benedict (6) was employed by the authors. It is possible that a value of 7 markedly different from that suggested by Benedict would yield still smaller standard deviations for t h e compressibility of methane than those recorded in Table 11. However, the agreement is better than the absolute uncertainty of the experimental data and there seems to be little gain from further refinement of the constants except by the inclusion of additional experimental data. If a correlation including all t h e available experimental data for methane were made, a more extended investigation of the effect of 7 upon the standard deviations of the compressibility factor should be included.
d
x
2445
data for the higher pressures (17)are included. In order to make a direct compariaon with the constants suggested by Benedict (6) other measurements upon this substance were not included in these calculations. It is emphasized that the constants presenteQin Table I may not yield as satisfactory an estimation of vapor pressure and other behavior in the heterogeneous region aa would be obtained from the constants which were proposed by Benedict and which have been based partly on the behavior in the heterogeneous region. In many industrial applications, primary interest is focused upon the volumetric behavior of compounds and little regard is taken of the vapor-liquid equilibria. Methane servea as an example of such a compound. The constants for this hydrocarbon have been computed from 66 measurements of pressure, specific volume, and temperature which were made from 100" to 460' F. and a t pressures from 1500 t o 10,OOO pounds per square inch (16). The unpublished detailed experimental data were employed. These data permitted a study of the influence of experimental information a t elevated pressures upon the values of the constants. The values of y were taken 5% above and below aa well as equal to that employed by Benedict (6). The resulting values of the constants obtained for these three values of y as well as those reported hy Benedict are recorded in Table 11. In this instance marked differences in the constants obtained by the least squarea method and those proposed by Benedict (6) were realized. These differences probably result from the fact that data at higb reduced pressures and temperatures were used in the present calculations. The three values of y employed in the present evaluation were chosen as a result of Benedict's experience with data at somewhat lower pressurea. The standard deviation of compressibility factors computed from the present constants is approximately 40% of that obtained by m e of the constants determined by Benedict (6). Thii daerence is not
DISCUSSION
A method has been proposed for the evaluation of the constants of an equation of state of the form proposed by Benedict (6). The method permits the direct evaluation of all the constants but one of the equation for a condition that corresponds to a minimum in the deviation of the compressibility factor from experimental values. Such calculations for a series of chosen values of one nonlinear constant yield results with much less effort than if the problem were treated in ita entirety by a least squares method. It is possible that other forms of an equation of state would lend dhemselves simply to an extension of the least squares method to all the constants. The method appears t o be promising for the direct evaluation of constants for equations of state from experimental data without the need for intervening graphical or numerical smoothing operations. This discussion is only a preliminary consideration of the problem but serves t o point out the ut8ilityof commercially available automatic digital computing equipment. The constants recorded in Tables I and I1 are not proposed as an improve ment over those originally indicated by Benedict (6)sinae only a limited background of expethental data was employed. However, their application describes the behavior of propane and methane with a somewhat smaller deviation from thf: limited experimental data than was obtained previously. ACKNOWLEDGMENT
The aid of G.D McCann, Jr., and S. P. Frankel in making available the facilities of the computing laboratory at the California Institute of Technology is acknowledged. Olga Strand: vold rendered assistance in the preparation of the manuscript. Manson Benedict kindly consented to review this material while in manuscript form. NOMENCLATURE
Ao, a, BO,b, CO,c = constants of the Benedict equation of state (Eauation 1)
= molal*weight,lb.-mole per cu. foot Ki = derived constants of the Benedict equation of state d
Kt
-
-Ao/R
Kr Ks
= -a/R uu
I N D U S T R I A L A N D E N G I N E E R I N G CHEMISTRY
2446
KJ -Go/R K7 = c / R K, = b K 6 ( r )= values of K , which satisfy Equations 7 and 8
N = number of experimental states employed in the evaluation of a set of constants R = universal gau conutant, (Ib./sq. inch)(cu. foot) per (1b.mole)( R.) 2’ = absolute temperature, a R. Z = compre~sibilit factor 4 = error in calcuited compressibility factor (Equation 2) z 3: average error in compressibility factor (Equation 17) 4 3 = root mean square error of compressibility factor (Equation 16) = constant in the Benedict equation of state (Equation 1) = constant in the Benedict equation of state (Equation 1) y = functions of d and T defined by Esuation 2 cpi d ‘p6 d87T cp2 d/T ‘pe = d6/T ( ~ 3= d / T 3 = d 2 / T * (1 rd2) e - @ = de = d:/Ta = cium of the squares of the errors in 2 for a given value of 7 and values of Ki which satisfy Equations 7 and 8 K i ) = sum of the squares of the errors in Z for a set of values of Ki and y (Y
:!
+
Subscripts R = value used by Benedict c = value calculated from an equation of state i = any one of the seven values of Ki = any one of the seven values of Ki = the nth experimental point LITERATURE CITED (1) Beattie, J. A . , Chem. R m . , 44, 141-92 (1949). (2)
I3eattie. J. A.. and Bridgeman, 0.C., Proc. Am. Actrd. Arts Sci., 63,229-3x3 ( 1 ~ 8 ) .
Val. 43, No. 11
(3) Besttie, J. A.. Kay, W. C., and Kaminsky, J., J . A m . Cht.7n. SOC.,59, 1589-90 (1937). (4) Beiiedict, M.,Ibid., 59,2224-33 (1937). (5) Benedict, M.,Webb, G. R . , and Rubin, L. C., J . C~KWL Ph~p.. 8,334-45 (1940). (0) I b d . , 10,747-58 (1942). (7) Bwugh, H.W‘., Schlinger, W. G., and Sage, B. H., American
Documentation Institute, Washington, D. C., Document
3235 (1951). (8) Connolly, T.J., Frankel, S. P., and Sage, B. H., Elec. Eng., 70, 47 (1951); American Documentation Institute, Document 3036 (1951). (9) Eckert, W.J., “Punched Card Methods in Sojentific Computation,” New York, Columbia University Preea, 1940. (10) Hough, E. W.,and Sage, B. H., Chem. Baa., 44, 193-204 (1949). (11) King, G . W.,J. C h m . Education, 24,61-4 (1947). (12) Lewis, G.N., J. Ant. Chen. Soc., 30,668-83 (1908). Jr., and Motard, R. L., IND. (13) Marchman, H.,Prengle, H. W-., ENG,CEEM.,41,2658-60 (1949). (14) hIilne, W.H., “Numerical CalcuIus,” Princeton, N. J., Princeton University Press, 1949. (15) Olds, R. H.,Reamer, H. H., Sage, B. H., and h e y , W. N., IND. ENO.CEEM., 35, 9224 (1943). (16) Prengle, H. W., Jr., and Marchman, H., Ibid., 42,2371-4 (1950). (17) Reamer, H.H., Sage, B. H., and Lacey, W.N., Zbid., 41,482-4 (1949). (18) Sage. B. H., Ibid., 42,631-9 (1950). (19) Sane, B. H., and Lacey, W. N., Advances in C i m . Series, No. 5, 372 (1951). RECEIVED March 21, 1951. For detailed derivation of t h e equations in thia article, order Uociiinent 3235 from the Amerioan Doounientation Institute, 1719 N St., N.W., Washington 6, D. C., remitting $1.00 for niicrofilm which yield8 irnage- 1 inch high on standard 35-min. motion picture film or $3.00 for photocopies 6 by 8 inches which are readable without optiaal aid.
Calculation of Equilibrium Gas Lompositions Application o f Digital Computers to the Iteration
P rocedare
FREDERICK J. MARTIN AND MORRIS YACHTER The M. W.Kellogg Co., Jersey City, N . J .
I n order to evaluate performance of combustion processes it is generally necessary to calculate equilibrium gas product compositions. Mathematically this requires the solution of a system of simultaneous, nonlinear, algebraic equations and ordinarily this represents a tedious and timeconsuming effort. Consequently, it is necessary to attempt to simplify the solution of such system equations. Because of the nature of the system equations, a simple iteration method for the solution of such systems can be developed. This has been done and the method has proved far less time-consuming than methods commonly used. Moreover, this iteration method has lent itself nicely to solution by the aid of digital computers. This type of equilibrium calculation is practically a universal necessity whenever performance has to be evaluated in advance for rocket or other combustion process designs. A simple, accurate, and fast method of calculation should prove to be of value to workers in these fields.
I
N ROCKET engineering it is often desirable to compute the composition of the product gases in the combustion chamber. Such computations are also required where gas temperatures and rocket performance are to be calculated. In the usual calculation the chamber pressure is BpecifiedLe., it is an independent variable. The product gmes are assumed to be in chemical equilibrium and to be perfect in that Dalton’s law applies, the fugacities are equal t o the partial pressures, and, for any molecular specieb, the enthalpy is a function of temperature. The chamber temperature is then determined in accordance with the first law of thermodynamics. While the necessary equation8 can be set up readily, the solution is usually accomplished only by tedious iterative procedures. The first step is the calculation of the product gas composition for a stated temperature and pressure. This step is by far the most time-consuming in the temperature and performance calculations, and it will be discussed in this paper. The product gas composition for any stated temperature and pressure is obtained by the solution of a series of simultaneoua equations. There is one equation for each chemical element expressing the conservation of mass for that element. Theae
