w‘ .- . .4. , 25.3" .2. .7; m fit. ”I _ 1- 2. fisfgmi f __ é. _ j v Q. ABSTRACT DEVELOPMENT OF PSEUDOCRITICAL MIXTURE RULES APPLICABLE TO THERMODYNAMIC CALCULATIONS IN THE CRITICAL REGION by John Douglas Stevens Phase equilibria and other properties of fluid mixtures in the criti- cal region may be computed thermodynamically if equations of state appli- cable to this region are available. Such equations of state usually have generalized forms involving two or more constants and are used in con- Junction with mixture rules which express these constants as functions of composition. The eight constant Benedict-Webb-Rubin equation of state1 has yielded reasonably satisfactory computations of this type. This research work was directed toward establishing satisfactory mixture rules which may be used with a two constant reduced equation of state. The derived reduced equation of state was based on the Benedict- Webb—Rubin equation for propane. Pseudocritical mixture rules express these quantities as a function of composition. The criterion used to evaluate a mixture rule was its ability to predict conditions along the critical envelope curve. This is the boundary of the two phase region, and also represents the pressure at which the equilibrium ratios converge to unity. Eecause of the complexity of the relations involved, computations were prOgrammed on a CDC 3600 computer. Among the programswritten and tested were the following: (a) Computation of compressibility factor and reduced pressure from reduced temperature and density using the correspond. ing states principle with propane, butane, and other hydrocarbons as reference materials. (b) Computation of derived thermodynamx:properties . 5+ .SS:E\.U e be. I ILL; L a! stir ethane- :eztal 338828 rules: here ' 3- John Douglas Stevens including the second and third partial derivatives of free energy with respect to mole fraction at constant temperature and pressure. (c) Com- putation of critical envelope curves by simultaneously converging the second and third free energy derivatives to zero. The eight constant Benedict-Webb-Rubin equation of state gave values of critical pressure, temperature, and compressibility factor for the ethane-n-heptane system which were in good agreement with Kay'32 experi- mental data. With the two constant reduced equation of state, good agreement was obtained for a binary mixture using the following mixture rules: 2 2 1/2 2 Pcm I [1E1x1(Tc/Pc )1/151x1(T6/Pc)1] 2 2 Tcm - [iflxi(Tc/Pol/2)112/1£131(Tc/Pc)1 where Tc1 and Pc1 are the critical temperature and pressure, respectively, of pure component i, and Tcm and Pcm are the pseudocritical temperature and pressure, respectively, of the mixture. x1 is the mole fraction of component i in the mixture. Attempts were made to develop programs for generating pseudocritical pressures and temperatures from experimental data on true critical tQMh perature, pressure, and density. All approaches used the critical prOper- ties of the pure components as known conditions. Using an approach which assumed cubic forms for the pseudocritical curves, it was possible to gen- erate pairs of points on these curves. The generated points were in good agreement with the mixture rule suggested above. 1Benedict, M., G. B. Webb, and L. C. Rubin, J. Chem. Phys., Q, 33“ (l9h0). 2qu, w. 3., Ind. Eng. Chem., 29, h59 (1938). DEVELOPMENT OF PSEUDOCRITICAL MIXTURE RULES APPLICABLE TO THERMODYNAMIC CALCULATIONS IN THE CRITICAL REGION By John Douglas Stevens A THESIS Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Chemical Engineering 1965 ACKNOWLEDGEMENTS The author wishes to express gratitude to the many people who gave help and encouragement: to Dr. Carl M. Cooper for his valuable assistance and guidance, to the members of the guidance committee, to the Division of Engineering Research of the College of Engineering at Michigan State University and to the National Science Foundation for providing financial support, and to his wife for her encouragement and patience. ii TABLE OF CONTENTS CHAPTER PAGE ACKNOWLEDGEMENTS................................................ 11 LIST OF TABLES.................................................. V LIST OF FIGURES................................................. Vi LIST OF APPENDICESOCOOOOOOOOOOOOO.0.0000OOOOOIOOOOOOOOOOOOOOOOO. V111 INTRODUCTION...OOOOOOOIIOOOOOI..00......OOOOOOOOIOOOOOOOOOOOOOOO 1 BACKGROUNDOOOOOOOIOOO00......0....OO.O.OOOOOOOOOOOOOOOOOOOOOO... h General Equilibrium Theory h Corresponding States Principle 5 Low and Moderate Pressure Vapor-Liquid Equilibria 7 High Pressure Vapor-Liquid Equilibria 11 Comparison of Vapor-Liquid Equilibria Prediction Methods 1h Mixture Rules 15 Pseudocritical Mixture Rules , 16 THEORYIOOOOOOOOOOOOOOOOOOOOOOOI00.......0...OIOIIOOOOOOOOOOOOOOO 20 Definitions of Critical Point and Convergence Pressure 22 Thermodynamic Relations at the Critical Point 23 Expression for Free Energy 25 Corresponding States Theorem for Mixtures 26 Equation of State 23 STATEMENT OF THE MATHEMATICAL PROBLEM........o............o..ooo 31 Evaluation of Equations of State 31 Evaluation of Published Pseudocritical Mixture Rules 33 Generation of a Pseudocritical Mixture Rule 35 Numerical Check on Analytical Derivatives 38 EXPERIMENTAL WORKOOOOOO0.00.000OOOOOOOOOOOOOOOOOOOOOOOO00.0.0... 1‘0 Experimental Procedure ”0 Experimental Critical Point Vapor-Liquid Equilibrium Data kl Evaluation of Equations of State A2 Derivation of the Reduced Equation of State hh Derivation of the Analytical Expressions for the Critical Point Free Energy Relationships hS Evaluation of Available Pseudocritical Mixture Rules AS Generation of Pseudocritical Mixture Rules h6 Check for the Analytical Derivatives 50 Chronological Order of the Research Work 50 iii CHAPTER PAGE RESULTS AND DISCUSSION.......................................... 51 Fitting of Curves to Experimental Data 51 Equations of State for Mixtures 53 Generalized Equation of State 60 Generalized Equations for Free Energy and its First Three Derivatives 60 Evaluation of Published Mixture Rules 61 Generation of Pseudocritical Mixture Rules 68 Differentiation Check 69 CONCLUSIONSOO0.000000000000000000000eee.OOOOOOOOOOOOOOOOIOOOOOOO 7h SUGGESTIONS FOR FUTURE STUDY.................................... 76 NOWCLATUREOOOOOOOOOOOIO00600000000000.eOOOeOOeeeOOOOOOeOOOOOOO 78 BIBLIOGRAPHYOOOO.0.00.00.00.00...000.00.OOOOOOOOOOOOOOOOOOOOOOOO 81 iv .l-.': .r-JU TABLE 1. 2. 3. 5. 7. 9. 10. ll. 12. 13. LIST OF TABLES Comparison of vapor-liquid equilibrium prediction methods...... Experimental data for ethane-n-heptane system at the true critical paint (23:h63)00000000000000000000000.000000.00000 Critical temperature, pressure, compressibility factor, and reduced compressibility factor as computed from a least squares fit Of Kay's data........nu...................... Values of critical density, temperature, pressure, and compressiblity factor as predicted by the van der Waals equation of state for mixtures....................... Values of critical density, temperature, pressure, and compressibility factor as predicted by the Benedict- Webb-Rubin equation of state for mixtures.................. Critical temperatures, pressures, and compressibility factors and critical reduced densities obtained from the general- ized equation of state using Kay's mixture rule............ Critical temperatures, pressures, and compressibility factors and critical reduced densities obtained from the general- ized equation of state using van der Waals' combinations... Critical temperatures, pressures, and compressibility factors and critical reduced densities obtained from the general- ized equation of state using Joffe's mixture rule.......... Values of pseudocritical temperature and pressure as computed by Program TOTRANGEeeeeeeeeeeeeeeseeeeeeeeeeeeeeeeeeeeeeeee Results from program CHECK for three values of compositions.... Generated values of pseudocritical temperature (°K) in APPrOQCh Ieeeeeeeeeeeeeeeoeeeoeeeseeeeoeeeeeeeseeeeeeeeeeoe Generated values of pseudocritical pressure (atm.) in Approach 10......I.0IO.ICC...I0.0.000.0.0.O.OOIOCOOOOOOOOOOOOOOOOOCI Generated Tcm and Pcm values from program LINEAR using a set of starting values from program DECLINE PAGE 15 A2 52 55 56 62 63 6h 68 73 138 139 lhh «can. I v .e‘qu ('1') LIST OF FIGURES FIGURE PAGE 1. Illustration of a critical envelope curve...................... 23 2. Free energy at constant pressure and temperature as a function or CGMPOBitioneeeocoeooeeoooeeoooeeeeaoeoeeeeoooeeeoeoeeeeo 2h 3. Path for calculation of the free energy of a binary mixture at temperature T and pressure Poeeoooeeeoooeooeooseeeeeoeeeeee 26 h. Flow diagram for calculation of To , Pen, and Dr using experi- mental data Tk. Pk, and zr dat oceeeeeeoeeeoeooeeoeeeeeeeeo 36 5. Flow diagram for calculation of Tcm and Pcm using experimental Tk and Pk data............................................. 37 6. Flow diagram for calculation of To and Pcm assuming the form or these curves to be CUbiceeeeeeeeeeeeeeeeeeeseeeeeeeeeess 38 7. Flow diagram for calculations in the program for checking analytical derivative8........o..........e................. 39 8. Convergence pressure data - ethane lightest component.......... 53 9. Critical temperature as a function of mole fraction ethane..... 57 10. Critical pressure as a function of mole fraction ethane........ 58 11. Critical compressibility factor as a function of mole fraction ethane.....o............................................... 59 12. Critical temperature as a function of mole fraction ethane..... 65 13. Critical pressure as a function of mole fraction ethane........ 66 1h. Critical compressibility factor as a function of mole fraction ethweOOOOOOOOIOOOOOO0.00000.000.00000...OOOOOOOOOOOOOOOOOO 67 15. Pseudocritical temperatures computed from program TOTRANGE as a function Of mole fraction ethane........................... To 16. Pseudocritical pressures computed from program.TOTRANGE as a function Of mole fraction ethane........................... 71 17. Critical compressibility factors as computed from program TOTRANGE as a function of mole fraction ethane............. 72 18. Formation of a two phase mixture at temperature T and pressure 8% P00000000.0IOOOOOOOOOOOOOO00.0000000000.0000000000000000... vi .gaefl' ..;Ueill FIGURE PAGE 19. Path for calculation of the free energy of a binary mixture at temperature T and pressure Peeeeseeeeeeeeeeeeeeeeeeeeeeeeee 86 20. Pseudocritical temperature as a function of composition........ lhO vii est-1'" “A‘s ed‘v . 'vv 5“ h ‘fvs ea 0-4 LIST OF APPENDICES APPENDIX PAGE I Derivation of an Equilibrium Ratio .Expression for a Two Phase Mixture........................................... 8" II Derivation of Relationships for Calculating Free Energy of a Binary Mixture at Temperature T and Pressure P.......... 86 III Derivation of the Equations for G" and G"'................ 39 IV Derivation of Equations for Dr', Dr", Dr"'...............o 93 V Derivation of Equations for Z, Z', Z", and Z"'............ 95 VI Program LINEAR..........................................o..o 97 VII Program TOTRANGE............................................ 125 VIII Program CHECK.................................e............. 130 IX Generated pseudocritical curve results - Approaches I and II 137 viii :5- O. . me}. ., ‘I A. .. :1! I “A! ' DU- 62‘”. a“ INTRODUCTION The chemical engineer often performs calculations requiring a know- ledge of the values of thermodynamic quantities. This is especially true with mass transfer processes such as distillation and absorption where successful design of equipment is dependent on accurate vapor-liquid equilibrium data. Because vapor-liquid equilibrium properties are one of the most important thermodynamic quantities to a chemical engineer, this thesis is oriented towards obtaining a mixture rule which could be used to compute more accurate vapor-liquid equilibria in the most sensitive region, the critical region. In the critical region, the equilibrium ratio (vapor mole fraction divided by liquid mole fraction) approaches one at the critical envelope curve (convergence pressure curve). This curve, which represents the border of the two phase regiOn, is the locus of true critical temperatures and pressures for a multicomponent system. In order for a method of com- putation to be satisfactory in the critical region, it must predict the equilibrium ratio to be one, at temperatures and pressures corresponding to the critical envelope curve. Or conversely, if the equilibrium ratio is equal to one, the method must give the temperature and pressure to be equal to the critical temperature and pressure. There are two methods currently favored for calculating vapor-liquid equilibria at high pressures. One of these methods, the convergence pres- sure method, involves estimation of the shape of the critical envelope curve. This is generally done for the system under consideration on the basis of data for other systems. Equilibrium ratios are then obtained from empirical equilibrium ratio charts (K-charts) using convergence pres- sure as a parameter. The second method, the equation of state method, is l 2 based on equations of state for the pure components plus a mixture rule. Equations of state are approximations of actual physical behavior and have generalized forms involving pressure, temperature, and volume (or density) plus two or more constants. Mixture rules express these constants as functions of composition. An equation of state used in conjunction with a mixture rule provides a means of computing thermodynamic quantities, including equilibrium relationships, for mixtures. The convergence pressure method has its greatest accuracy in the critical region, but this is because it is based on experimental data in that region. However, the convergence pressure method becomes less accu- rate at lower pressures; it is not amenable to the calculation of other thermodynamic properties, and if extended very far from the region of the experimental data, it can lead to results which are thermodynamically im- possible. The equation of state method has all of the advantages of ther- modynamic consistency. The equilibrium ratios computed at low and moderate pressures are very satisfactory, but considerable errors are often encoun- tered at high pressures near'the critical region. If a more accurate equation of state could be developed for mixtures, through the development of better mixture rules, then the equation of state method would be a much more satisfactory method for calculating equilibria. This method would then permit the calculation of other ther- modynamic data, and it would require less specific experimental data than the convergence pressure technique. The objective of this research was to obtain a better equation of state for mixtures at high pressures. The method used was to derive a reduced equation of state and then to attempt to find or develop pseudo- critical mixture rules which when used with the reduced equation of state, would yield accurate thermodynamic predictions in the critical region. . ~4- ‘ 9 v I1. 6‘ 3'9 ’0 vvehfi e ales AFQIAC bU-‘y‘l. 3 A reduced equation of state has only two constants; in this work these were pseudocritical temperature and pressure. Pseudocritical mixture rules express pseudocritical temperature and pressure as a function of composition. The criterion used for evaluating mixture rules was to use them in conjunction with an equation of state to predict the critical -envelope curve (convergence pressure curve) for a particular multicomp ponent system and compare these predictions with experimental data. f!"- ;se, a 'r’ ht bulk q :ezer: .- J e 9‘ ~ .0'44 t C“ BACKGROUND Accurate methods for predicting vapor-liquid equilibria have always been a prime requisite of the design chemical engineer. However, most methods developed to date have limitations which greatly restrict their use, and the engineer often must use a combination of methods or may even not have any acceptable method at his disposal. General Equilibrium Theory A system is said to be in equilibrium if there is no apparent change in the intensive properties with respect to time. However, the intensive properties themselves, such as concentration, partial molal enthalpy, den- sity, refractive index, etc., may be different in different phases of the system. The only requirement for equilibrium is that all potentials which cause changes should be in a well balanced state. Therefore, the equilibrium state requires that the temperature (thermal potential) and pressure (mechanical potential) be the same everywhere. In addition, a third potential, the chemical potential, should also be in a balanced state. This is equivalent to stating that the partial molal free energy be the same in all phases for all components. The balance of the above three potentials is often referred to as the conditions for equilibrium and they may be written symbolically as TV - TL (1) PV 8 PL (2) Fiv - Ru (3) where TV, T 8 temperature in the vapor and liquid L phases, respectively. P P 8 pressure in the vapor and liquid V’ L phases, respectively. b FiV’ FiL I free energy/mole for component i in the vapor and liquid phases, respectively. Condition (3) may also be written fiv ' f1L (h) where f I fugacity of component i in the vapor and iV’ fiL liquid phases, respectively. The most practical means of expressing vapor-liquid equilibria is in terms of the equilibrium ratio or K-factor which is defined for component i as Ki ' yi/xi (S) where yi I mole fraction of componenti in the vapor phase. x1 I mole fraction of component i in the liquid phase. Another constant called the vaporization equilibrium constant is sometimes used. This is defined as Kv ’ aiV/aiL ' (YiVyi)/(Yiin) (6) where a I the activities of component i in the vapor iV’ aiL and liquid phases, respectively. YiV’ YiL I the activity coefficients of component i in the vapor and liquid phases, respectively. If ideal solutions are formed in the two phases, then both activity coefficients are unity, and equation (6) reduces to equation (5). Other- wise, the relationship between the two equilibrium constants is Ki ‘ (YiL/Viv)xv ' (7) The discussion which follows refers to the equilibrium state in terms of the K-factor. Corresponding States Principle A3 Van der Waals first defined the term "reduced condition" and 6 presented the corresponding states theorem for all pure gases in 1873. Young extended this concept to liquids in 1899. The theorem of corresponding states claims that all pure gases, when compared at the same reduced temperature and reduced pressure, will have the same compressibility factor or, in other words, deviate from perfect gas behavior to the same degree. Reduced temperature and reduced pressure are defined as Tr I T/Tc (8) Pr I P/Pc (9) where T temperature of the fluid P I pressure of the fluid Tc I critical temperature of the fluid Tr I reduced temperature of the fluid Pc I critical pressure of the fluid Pr I reduced pressure of the fluid? The compressibility factor may be written symbolically as z - PV/R'r (10) where Z I compressibility factor V I volume/mole R I gas constant Graphs of compressibility factor as a function of reduced terperature and reduced pressure are found in most thermodynamic textsoflo When the corresponding states concept is extended to rixiwres, the reduced preperties of mixtures based on the true critical p? ,2rties of the mixture do not give the same functional relations for czc~wessibility factors as for the pure components. In order that compressinility factors for mixtures will follow the same functional relationships as do pure 7 components, it is necessary to use hypothetical values for the critical properties of the mixture. These hypothetical critical properties are called pseudocritical properties. The pseudocritical temperature and pseudocritical pressure are denoted by Tom and Pcm, respectively. Any preperty which is put into reduced form by using the pseudocritical prep- erty ;s said to be in the pseudoreduced condition. Low and Moderate Pressure Vapor-Liquid Equilibria If two pure substances, originally at their respective vapor pressures, are isothermally changed to pressure P and then mixed to form a two phase mixture, the thermodynamic expression for the equilibrium ratio is (Appendix I) yA/XA = PAVPvAVPYAL/(PVAPYAV) L . (11) , P where L I exp[(l/RT) I VLdP] ' PAVP PAVP I vapor pressure of pure component A P I total pressure of mixture VAVP I fugacity coefficient of pure component A at its vapor pressure I f /P AVP AVP VAP I fugaCity coefficient of pure component A at total pressure P I fAP/P YAL I activity coefficient of component A in the liquid phase of the mixture YAV I activity coefficient of component A in the vapor phase of the mixture V; I molal liquid volume of pure component A. A similar expression may be written for component B. ‘o e, '0 E II- nu t" (I) 'rl [It ('3 ( 7 I e.‘ ‘ (I) 8 Relationship (11) involves the use of experimental vapor pressures and correction for deviations from ideal behavior. It is this relation- ship which is the basis for many methods of vapor-liquid equilibrium prediction. If there are no deviations from ideal behavior in either the liquid or vapor phase and if the liquid correction factor (L) may be neglected, equation (11) reduces to Raoult's Law: yA/xA = PAVP/P (12) Raoult's Law is limited to the case of a mixture of perfect gases in equilibrium with a liquid phase which is an ideal solution. That is, Raoult's Law is valid when molecules of each component are of approximate- ly the same size and when the pure components mix in both the liquid and gas phases without the complicating effects of molecular association, chemical combination, and the like. Raoult's Law neglects the effect of composition and to some extent, the effect of pressure on the behavior of a component in the vapor and liquid phases. Fugacity coefficients vAVP and ”AP are corrections for imperfect gas behavior due to the effect of volume of and the attraction between the molecules of the vapor. Fugacity coefficient values may be read from the fugacity charts. These are charts of f/P as a function of reduced pres- sure and reduced temperature. They were constructed using compressibility factor charts and should be valid as long as the corresponding states principle holds. 38 used the fugacity charts to construct Souders, Selheimer, and Brown a plot of K-factor as a function of temperature and pressure. They assumed that the liquid and vapor mixtures were ideal solutions, thereby reducing equation (11) to 9 P v yA/XA :s M L (13) AP Since v I f/P, equation (13) becomes yA/xA ' fAVPL/fAP ‘ fAL/fAv (1”) where fAv I fugacity of pure vapor A at the pressure and temperature of the system. fAL I fAVPL I fugacity of pure liquid A at the pressure and temperature of the system. Kecharts were constructed for several of the hydrocarbons. The quantities YAL and YAV are called activity coefficients and are corrections for deviations from Amagat's Law in the liquid and vapor phases, respectively. Because of the dense nature of liquid solutions, YAL is much more important than YAV at low and moderate pressures. At pressures near the critical where both phases are quite similar, YAV becomes a factor. The deviations from Amagat's Law behavior are due to the fact that attractions between similar molecules are different than attractions between dissimilar molecules. Activity coefficients are composition dependent and hence, whenever they are equal to anything except unity, the K-factor is composition de- pendent. The methods that have been used for predicting activity coef- ficients are either empirical or semitheoretical, the theoretical part being the use of thermodynamic equations to direct the development of empirical rules. Many rules have been proposed, but more commonly the activity coefficients are estimated by such equations as the Wohl, Margules, van Laar, or Redlich and Kister equations.20 37 It should be noted that some authors express activity coefficients as a product of two activity coefficients, each of which corrects for ‘ee. V “..EA V - ~ A. A.. .v a, '.»k . I _:::-~.. 0 'i-b“ c I ._.~ . Qg h ..-.. a. *3. h A j. V. A‘Q IV) in 10 different types of deviations. This could be expressed Y . dw' where y I overall activity coefficient 7' I activity coefficient expressing deviations from Amagat's Law as a result of chemical dissimilarity ¢ I activity coefficient expressing deviations from Amagat's Law as a result of molecular size or volatility. If a system is composed of a homologous series, then 7' I l. The liquid volume correction term (L) in equation (11) is a correc- tion for the difference between the fugacity of the liquid at the system pressure and the fugacity at the pure component's vapor pressure. This term is generally not important at low and moderate pressures although it may be calculated any time the liquid molal volume is known or an equation of state is available. By introducing additional parameters into the correSponding states 30 has principle, more accurate fugacity predictions can be made. Pitzer extended the correSponding states principle by using a third parameter. Although others had previously suggested a third parameter, none of the methods prOposed were as successful as Pitzer's. He called the new param- eter the acentric factor, which is supposedly a measure of the deviation of the intermolecular potential function from that of simple spherical molecules. The compressibility factor Z may then be written Z I Z(Pr, Tr,w) where w I acentric factor. Hougen, Watson, and Ragatzao also used a three parameter approach. They calculated equilibrium ratios as a function of the critical ll compressibility factor, reduced temperature, and reduced pressure. A value of 0.27 was used for the critical compressibility factor in con- struction of a table and a chart, but corrections were given for values of critical compressibility factor other than 0.27. The table and chart are not accurate at high pressure and temperature because the assumption of ideal solutions in the liquid and vapor phase is no longer valid. Also, both the Pitzer and Hougen methods are limited by the accuracy of the rules used to calculate the pseudocritical properties. Bloomer and Peck7 also modified the corresponding states principle by introducing a third parameter S. This parameter takes into account the aspherical factor for nonpolar molecules. 8 is determined from P-V-T data and can be expressed in a simple relationship to the critical com- pressibility factor. High Pressure Vapor-Liquid Equilibria Two basically different approaches can be used for prediction of vapor-liquid equilibria near the critical region. These are (l) equilibria from an equation of state and (2) convergence pressure methods. The most commonly used equation of state for prediction of high pres- 3 sure equilibria is the Benedict-Webb-Rubin equation. This equation uses eight constants to describe the vapor and liquid phase behavior of a given compound. The authors published a set of constants for each of twelve light hydrocarbons. They alSo Suggested mixture rules for combining the eight constants to extend the eduation of state to mixtures.h However, 32 Price 23 2.}, indicate substantial deviations between experimental K- factors and those predicted by the Benedict-Webb-Rubin mixture equation. 17 Hester has shown that the Benedict-Webb-Rubin mixture equation also may yield highly improbable values for the pseudocritical compressibility 12 factor. Benedict, Webb, Rubin, and Friend6 recorded K-values of twelve hydrocarbons in a series of 32h charts called the Kellogg charts. In order to reduce the number of composition variables, these authors ex- pressed all vapor and liquid compositions in terms of two variables, the molal average boiling points of the vapor and liquid, respectively. Each chart refers to a particular'compohent at a particular pressure. The Kellogg charts are reliable for many applications, but they re- quire successive corrections for the compositions of both phases. Another disadvantage is the necessity for interpolation between charts with respect to pressure. These 32h charts have been reduced and.made easier to apply by DePriesterll and Edmister and Ruby.13 The DePriester charts include two for each hydrocarbon, one giving the fugacity ratio for the vapor phase, and the second, the fugacity'ratio for the liquid phase. Each chart represents the relation between pressure, temperature, and composition, and thus eliminates the necessity of interpolation between charts with respect to pressure, as with the Kellogg charts. DePriester's consolida- tion of the information on the Kellogg charts into twenty-four charts was with little loss of accuracy. However, the charts only cover the pressure range up to 1000 psia while the Kellogg charts include pressures up to 3600 psia. Edmister and Ruby used only six charts plus one table of vapor pres- sure data to cover the same range of pressure, temperature, and composi- tion as the Kellogg charts. 'The Edmister and Ruby charts were derived from the Kellogg charts by correlating fugacity coefficients in terms of four reduced parameters. They are in considerably greater error than the 13 Kellogg and DePriester charts. Gamson and Watson1h developed equations for activity coefficients based on generalized fugacity coefficients and the pseudocritical con- 37 empirically modified the Gamson and Watson cept. Smith and Watson relationships and constructed a graphical correlation of activity coeffi- cients as functions of the pseudocritical temperatures and pressures of the phases. These generalized correlations are applicable to systems of components of no chemical dissimilarity. On the basis of these correla- tions, Smith and Smith36 published a set of K-charts for hydrocarbons. The generalization of the Smith and Smith charts causes some loss of accuracy for Specific data. However, this generalization does allow the charts to be used for many more hydrocarbons than the Kellogg charts. Recently Mehra and Thodos27 have developed an approach which is fundamentally different from existing methods for the development of equilibrium correlations from experimental data. They predict beactors for binary hydrocarbon systems in the critical region by using the normal boiling point ratio, the reduced vapor pressure, and the pseudoreduced vapor pressure. Good accuracy is obtained with the correlations, but the limitations to a binary hydrocarbon system limit the practicality of the method. All of the methods presented above may be classified as theoretical. Correlations of K-factors based on experimental data were developed about the same time as the Kellogg charts and modifications thereof. These empirical correlations have the advantage of being based on experimental data and are, therefore, very accurate in the region where the experi- mental data was obtained. However, the experimental correlations are not amenable to the calculation of other thermodynamic properties, and the P? 52:: 1h correlations are only valid in the range of the experimental data. The most prominent of the empirical approaches is the convergence pressure method. The standard procedure is to predict convergence pres- sure based on experimental data for various binary systems and then ex- tend the method to multicomponent systems by treating the mixture as a 26 developed an improved fictitious binary system.3h’h5 Lenoir and White method using effective boiling points for predicting the convergence pres- sure. Once the convergence pressure is known, the K-factor is found from empirical correlations between equilibrium ratios and the estimated con- vergence pressure. The NCAA (Natural Gasoline Association of America) charts28 were prepared for prediction of K-factors using the convergence pressure technique. The convergence pressure method is relatively con- venient to use and gives acceptable results in the critical region. Comparison of Vapor-Liquid Equilibrium Prediction Methods The methods discussed above for predicting vapor-liquid equilibria are compared in Table l. The construction of such a table is a matter of judgement and should be taken only as a general indication of the strong and weak points of the methods discussed above. Question marks mean that it is difficult to determine if the method is thermodynamically consistent. One of the aims of this research was to devise a method which would exhibit all of the prOperties listed in the table. (D 15 TABLE 1. Comparison of vapor-liquid equilibrium prediction methods Method General- Composition Accuracy Thermo. ized Dependent Low Moderate High Critical Consist. P P P Region Raoult's Law Yes No Yes No No No Yes Souders at 31;. No No " Yes Yes No No Yes Hougen-Watson Yes No Yes Yes No No Yes B-WIR Eqn. No Yes Yes Yes Yes Yes Yes Kellogg charts No Yes Yes Yes Yes No ? DePriester No Yes ' Yes Yes Yes No 7 Edmister-Ruby Yes’ Yes Yes Yes No No No Smith-Smith Yes Yes Yes Yes No No Yes Convergence No Yes Yes Yes No , Yes No Pressure *Separate charts are needed for methane Mixture Rules Success with the equation of state method for calculation of vapor- liquid equilibria is dependent on the accuracy of the mixture rules used. A mixture rule is some arbitrary procedure for converting an equation of state for pure compounds into an equation of state for a mixture of those compounds. Equations of state for pure compounds are generally converted to equations of state for mixtures by some method of combining the equation of state constants. Amagat's, Dalton's, and Bartlett's ruleshh are simple mixture rules, but they are not applicable in the critical region. h Examples of more sophisticated mixture rules are those of Benedict at 9;. They converted their equation of state to mixtures by specifying rules for l6 combining the eight constants (A0, Bo’ Co’ a, b, c, o, v) for the pure com- ponents into corresponding constants for the mixture. A reduced equation of state only has two constants. When the reduced equation is used for a pure compound these constants are the critical temperature and pressure of that compound. (For mixtures the two equation constants are the pseudocritical temperature and pressure of the mixture. These pseudocritical quantities are generally expressed as a function of the critical properties of the pure components and the composition of the mixture. Such relationships are called pseudocritical mixture rules. Pseudocritical Mixture Rules Kay22 defined the pseudocriticals for a mixture of n components as ; n rem =- 2 xiTci (15) i=1 n Pc I Z x Pc (15) m i=1 i i For binary mixtures, Tangho found that a simple additive relationship was a fair approximation for critical temperature, but the pseudocritical pressure usually deviated widely from linearity. Case and Weber9 also eb- served large deviations from Kay's rule. Van der Waals20 suggested the following pseudocritical mixture rules for use in the reduced form of his equation of state. . n n p.m . [iilxi(Tc/Pcl/2)1/1:1xi(Tc/Pc)112 (17) D i A n Tcm . [1:1xi(Tc/Pc1/2)1]2/1:lxi(Tc/Pc)1 (18) Joffe21 used the van der Waals equation of state to show that the assignment of a pseudocritical temperature and pressure to a mixture is qr.- vie. s ‘o ;, I I‘- 17 entirely consistent with the method of combination of equation of state constants. As a result, Joffe prOposed the following relations for com- puting pseudocritical constants of a gas mixture: n (Tc/Pcl/e)m I i£1x1(Tc/Pcl/2)1 (19) n n 1/3 1/3 3 (Tc/Po) I (1/8) X 2 x x [(Tc/Pc) + (Tc/Pc) ] (20) m 1'1 3'1 .1 j 1 J where the double summation centains one term for each possible permutation in pairs of like and unlike components of the mixture. Joffe claimed that for five different binary gas mixtures, his proposed mixture rules yielded more accurate rules than Kay's rules. However, although Kay's rule and Joffe's rule are effective for non polar gaseous mixtures, they are in error for saturated vapors and liquids. Leland, Chappelear, and Gamsonzh proposed the following mixture rules. -'n n mile 2 Z x x a a . . 131: T..l_1_ll (21) m n n 3 z z x x (b + b ) i=1 3-1 1 J 1 J ‘ n 1 Tom 151 xi(Zc)i Pc I *: (22) m n n , 3 t- 2 x x (b + b ) i-1 5-1 1 3 1 3 L I. I (ZcTcl'm/Pc)';ll'/2 9.1 b1 I (l/2)(ZcTc/Pc)1/3 The parameter a is an empirically determined function of pressure and the pseudocritical pressure as determined by Kay's rule. This mixture rule is applicable to liquids which may be approximated by simple spherical 18 molecules. Stewart, Burkhart, and V0039 simplified equation (20) to 1/2 2 n n (Tc/Po)m I l/3 2 x1 (Tc/Pc)i + 2/3[1£ x1(Tc/Pc) , 131 .1 which is used in conjunction with (19) to compute pseudocritical tempera- ture and pressure. Hougen 53.21920 claim that the Stewart method reduces the average deviations by over fifty percent as compared to the results obtained with Kay's rule. 15 Guggenheim and.McG1ashan 'used the reduced second virial coefficient to propose mixture rules for the pseudocritical temperature and volume. However, because use of the second virial coefficient is sufficient only where deviations from ideality are small, these mixture rules are of little use in vapor-liquid equilihrium calculations. 31 Prausnitz and Gunn adjusted the mixture rules suggested by Guggen- heim and McGlashan to fit experimental data for the second virial coeffi- cient. They suggested the following rules for calculating pseudocritical temperature and pressure. Tom I [B + (82 + ermyN/2]/(23ch) (2h) Pcm - [RTcm/chht yizc (25) 1 i where B I E y (VcTc) y 1,313 13 Y I 2 yiyJ(VcTc2) 13., i: ch 2 yiijcij 1: y1 I composition of component i ch I pseudocritical volume of mixture l9 r,s I parameters which are read from a table published by Prausnitz and Gunn. The quantities Tcij and VciJ are called the characteristic critical tame P9313119 and volume, respectively, and are computed from 1/2 - ATc Tc I (TciTcJ) 13 id Vcid I 0.5(VcchJ) ~ AVc where ATc 1: 13, Men I corrections to the characteristic temperature and volume, respectively. THEORY An equation of state, such as the eight constant Benedict-Webb-Rubin equation,3 can be used to predict thermodynamic properties in the critical region if accurate mixture rules are available for converting the equation of state to mixtures. The factor which prohibits extensive use of such an equation is the fact that values for the eight constants must be known for each of the components in the mixture. Because these values are difficult to determine,they have been found for only a very limited number of com- pounds. Thus, although the equation may be very accurate for systems where all the constants are known, the equation is not completely satis- factory because it can not be generalized to mixtures of other compounds. If an equation of state for pure compounds is put into a reduced form, it is changed to an equation with two constants, for example, critical temperature and pressure. If this reduced equation of state is applied to mixtures, the constants become the pseudocritical temperature and pres- sure since the true critical temperature and pressure of mixtures are known to give inaccurate reduced relationships. Pseudocritical tempera- ture and pressure are.usually expressed as a function of the critical properties of the pure components and the composition of the mixture. These relationships, called pseudocritical mixture rules, are difficult to derive because pseudocritical temperature and pressure have no actual phys- ical significance. However, if a reduced equation of state and pseudo- critical mixture rules could be found which together would give accurate thermodynamic predictions in the critical region, one could apply the re- duced equation of state to many systems with only knowledge about the critical prOperties of the pure components. This research concerned the derivation of pseudocritical mixture rules for light hydrocarbon systems, 20 ‘i... 0“. 21 but most likely such derived mixture rules could be used for many other classes of compounds. A big advantage of making thermodynamic calculations with a reduced equation of state is that the calculations only need to be performed once for a given reduced temperature and pressure.* Any other mixture with the same reduced temperature and pressure will have the same value for the thermodynamic quantity. Thus generalized tables or graphs can be pre- pared with reduced temperature and pressure as parameters. Vapor-liquid equilibrium calculations can be made with a reduced (generalized) equation of state if the equation is applicable in both the liquid and vapor phases. A plot of the quantity, fugacity divided by mole fraction, as a function of reduced temperature and pressure can be made by using the reduced equation of State to express the fugacity quantity as a function of reduced temperature and pressure. To calculate equilibrium ratios with such a chart, the unknown quantities, temperature, pressure, or compositions, are adjusted by trial and error until the sum of composi- tions in the vapor and liquid phases are each unity, and the condition, fiv ' fiL is satisfied for all components in the mixture. The equilibrium ratio is then computed from Ki ‘ yi/xi Although the construction of generalized charts has been possible for some time, such charts are of little ude without accurate pseudocritical mixture rules. a If a variable is put into a reduced form using a pseudocritical quantity, then it is more proper to speak of the variable as being in a pseudore- duced form. However, the more common practice of denoting pseudoreduced variables as simply reduced variables is used throughout this thesis. 9“ Vs 93‘ es la (.3 "h '~ RN e._"1 To“: U 5“ t ‘F‘ e‘htt ‘a a..' v"“: 9.3 i M. rare 22 One approach to obtaining a more accurate method for making thermo- dynamic calculations in the critical region would be to derive a reduced equation of state and then to find or develop accurate pseudocritical mixture rules for use with the reduced equation of state. The reduced equation of state could be derived using the eight constant Benedict- Webb-Rubin equation for a particular compound, for example, prOpane. If the corresponding states principle is valid, the derived reduced equation of state should be the same no matter what reference compound is used. Currently available mixture rules could be evaluated by using them with the reduced equation of state to predict vapor-liquid equilibria (convergence pressure curve) in the critical region and then comparing the results with experimental data. Or conversely, pseudocritical mixture rules could be derived by beginning with experimental critical region vapor-liquid equilibrium data (i.e. an experimental convergence pressure curve) and then computing pseudocritical temperature and pressure. Definitions of Critical Point and Convergence Pressure The critical point of a binary mixture is defined as the temperature and pressure at which the vapor and liquid phases of the mixture become indistinguishable. Or, in other words, the temperature and pressure at which the K-factor equals one. If the critical points are measured for several different mixture compositions, and these points are plotted on a pressure- temperature diagram, the curve which passes through all of the critical points is called the critical envelope curve° Such a curve is illustrat- ed in Figure 1. Points A and B represent the critical points for the two pure components, while the intermediate points represent critical points of mixtures. Curve AB, which connects all of the points, is the critical envelope curve . 23 Pressure i-3 I'-‘ Temperature Figure 1. Illustration of.a critical.envelope curve Convergence pressure is defined as the critical pressure at the temp- erature of the system. Thus, in Figure l, the convergence pressure of any Therefore, curve AB is also system at temperature T is the pressure P 1 1° called the convergence pressure curve. Thermodynamic Relations at the Critical Point Because prediction of the critical enve10pe curve (convergence pres- sure curve) is equivalent to predicting the points where the K-factors equal one, it would be very desirable to have a means of calculating the critical envelope curve. This would be possible if there were thermody- namic conditions which held only at the critical point and in fact, pro- vided a basis of definition of the true critical point. Such relationships do exist,12 and they are: 2 2 _ (a F/ax )T,P - O (l) (33F/3x o (2) 3 )T,P a 2h when T and P are the true critical temperature (Tk) and pressure (Pk), respectively. Relations (1) and (2) can be illustrated with a F-composition curve. /7fib D / / \\ / I, \ \ / I | \ / I ' \ \l ' \ T I J \ 4;; H K B Figure 2. Free energy at constant pressure and temperature as a function of composition In Figure 2, the F-curve at T1(T1>Tk) is concave upwards everywhere. That is, (32F/3x2)T P > o . 1 ’ thus, all honogenous phases on curve EG are stable. On curve AB (T3b Mixture Std. State / po P Pure B Pure B 5 6 (l-yA): r=1 3' P 3 h Figure 3. Path for calculation of the free energy of a binary mixture at temperature T and pressure P- Using the path in Figure 3, it is shown in Appendix II that the expression for G is G = In P + yA ln yA + (l - yA)ln(l - yA) + Z - l - 1n Z D 2 + I [(P/RT - D)/D lab (6) 0 where D a density.‘ The first three terms represent the free energy of a mixture of perfect gases at pressure P. The last four terms are a measure of the non-ideality of the imperfect gas mixture. It should be observed that by following the path of Figure 3, it is not necessary to be concerned with nonideal solu- tions since the gases are.in the perfect gas state when they are mixed. Corresponding States Theorem for Mixtures One of the basic assumptions made in this work was that the corre- sponding states theorem was applicable to mixtures. That is, it was assmn- ed that two parameters (reduced temperature and pressure) were sufficient 27 to represent the state of the mixture. The following discussion concerns the theoretical support available for applying this theory to pure com- pounds and to mixtures. After van der Waalsh3 prOposed the corresponding states principle in the late nineteenth century, many years elapsed before enough knowledge was gained about molecular behavior to support the hypothesis theoretically. During this period, the belief that the state of a system could be repre- sented by Just two parameters had only empirical Justification. With the advent of classical statistical mechanics, a means was obtained to give the proposal theoretical support. Hirschfelder, Curtiss, and Bird18 applied classical mechanics to spherical nonpolar gases to show that reduced pressure is a unique func- 29 showed a similar result tion of reduced temperature and volume. Pitzer for liquids while Guggenheim and McGlashan15 extended the treatment to mixtures of slightly imperfect gases. In these three works the primary assumption was that the energy of interaction between any two molecules was representable by the same general function of two parameters. That is, the expression for the potential energy was of the form v(r) - ef(r/0) where s and a are characteristic of the molecular species and f is the same general function for all molecules. The Lennard-Jones potential is of this two parameter form and represents an approximation to this general function. In each of the works mentioned above the authors were able to show that reduced temperature was a function of c and reduced volume a function of a. By establishing that the reduced pressure of a mixture of slightly imperfect gases was representable by a single reduced temperature and 28 15 provided support to the hypothesis that volume, Guggenheim and McGlashan the corresponding states principle could be extended to mixtures by treat- ing a mixture as if it were a pure compound. Rules which use critical pro- perties of the pure components of a mixture to compute values of reduced temperature and volume (or any other pair of the variables, pressure, temperature, volume) for a mixture are called pseudocritical mixture rules. The development of pseudocritical mixture rules for pressure and temperature was one of the primary goals of this research. The theoretical support for the corresponding states theorem is very limited. However, it is reasonable to assume that certain classes of compounds (e.g. hydrocarbons) will have a common function f. If this is the case, the corresponding states theorem appears to have reasonable theoretical Justification for its use with mixtures within that class of compounds. A third parameter approach as discussed under Background is one method of extending the corresponding states principle to a variety of classes of compounds. This third parameter serves to make adJustments for differences in the function f between classes of compounds. For mixtures, the expressions for reduced temperature, density, and pressure become: Tr I T/Tc m Dr I D/Dc m Pr I P/Pc m where Ten, Dem, Pcm I pseudocritical temperature, density, and pressure, respectively. Equation of State The need for an equation of state is apparent from the last term in 3 equation (6). The Benedict-Webb-Rubin equation of state gives excellent 29 representation of the critical prOperties of pure light hydrocarbons. This equation is applicable to both the liquid and vapor phases with the same set of equation constants. The form of the Benedict-Webb-Rubin equation for pure components is 2 3 P e RTD + (BORT - A0 - Co/T2)D + (bRT - a)D + aoD6 + (cD3/T2) (l + yD2)exp(-1D2) (7) where A0, B0’ C , a, b, c, a, y = Benedict-Webb-Rubin equation constants 0 which are dependent on the substance the equation is to describe. Benedict gt_2l,3 suggested that Co be expressed as a function of temp perature,:tf it was hoped that the equation of state would represent vapor pressures at the subatmospheric level. In this way the equation of state was made extremely accurate at very low pressure as well as in the criti- cal pressure region. The importance of this prOperty to this research was that any mixture rule which was devised could be used at lower pressures if the equation of state was also applicable in that region. Therefore, Co was expressed as a function of temperature. This expression is equa- tion (1) on page 115. Assuming the corresponding states theorem is valid, the expressions for reduced temperature, density and pressure can be substituted into equation (7) to give a reduced equation of state. Similarly, equation (7) can be substituted into equation (6) and the resulting expression put into reduced form by the reduced relationships listed above. The analyt- and (336/3x3) ical expressions for (32G/3x2) in terms of reduced T,P T,P’ parameters, can then be used to evaluate or derive pseudocritical mixture rules in the critical region by using the thermodynamic conditions that 2 (3 G/axz)T P and (33G/3x3)T P equal zero at the critical point. D ’ 30 Experimental critical point data can be used as a check on the results computed from available mixture rules or it can be used as the basis for deriving pseudocritical mixture rules. STATEMENT OF THE MATHEMATICAL PROBLEM This chapter summarizes, in mathematical form, relationships used in this thesis. The discussion is divided into four sections: (1) the search for an equation of state on which to base the derivation of a reduced equa- tion of state, (2) evaluation of published pseudocritical mixture rules, (3) generation of pseudocritical temperature and pressure values on which a new pseudocritical mixture rule could be based, and (h) a numerical method for checking analytically derived derivatives. The method of solu- tion, as attempted in this research, is also discussed. Evaluation of Equations of State In this section two equations of state were evaluated for their ability to predict convergence pressure curves. Each equation was used to compute critical temperature, pressure, and density values for a particular binary system. These computed values were then compared with experimental data. The equation which made the best predictions was used as the basis for the derivation of the reduced equation of state. h3 and Benedict-Webb-Rubin3 equations of state were The van der Waals combined with mixture rules recommended by van der Waals and Benedict 33.5%, to form equations of state for mixtures. That is, van der Waals: Equation of state for pure compounds P I P(T,D,a,b) Equation of state + -————+> for mixtures P = P(T,D,x) Mixture rules a I a(x) b I b(x) 31 32 B-W-R: Equation of state for pure compounds P I P(T,D,a,b,c,A B C a y o 090’! Equation of state + -————> for mixtures P I P(T,D,x) Mixture rules a I a(x) O O ; I y(x) The general expression for dimensionless-free.energy G is (page 26): D G = In P + x In x + (l-x)ln(l-x) + z - 1 - ln 2 + [ [(P/RT - DyDZJdD (1) o By noting that P I P(T,D,x) Z I Z(P,T,D,x) then G = G(T,D,x) (2) for both equations of state. From the expression for G, the derivatives (aQG/ax2)T.P and (336/3x3)T.P may be derived analytically. This gives equations of the form 0" = G°°(T,D,x) (3) G'99 = G°'"(T,D,x) (h) where G"‘ I (326/8x2)T’P 3 3 It has been shown that at the critical point, GH and G"' equal zero. Since the obJective in this part of the thesis was to calculate the criti- cal temperature, pressure, and density, equations (3) and (h) were solved simultaneously for the value of T and D which made G" and G"' equal zero for a given composition (x). Because the computed values were the 33 critical temperature and density, T and D were actually Tk and Dk, respec- tively. The critical pressure was then computed from Pk = Pk(Tk,Dk,x) This procedure was repeated for several values of x and the results were compared with the experimental critical point data. Because equations (3) and (h) were complex expressions, computer pro- grams were written which used the NewtonuRaphsonh2 convergence method to find a value of T and D which made G“' and G"° equal to zero. The partial derivatives required in the Newton-Raphson method were estimated numerically. Evaluation of Published Pseudocritical Mixture Rules This section of the research used a reduced equation of state in conJunction with available pseudocritical mixture rules to compute criti- cal temperature, pressure, and density values for a particular binary sys- tem. The mixture rules were evaluated on the basis of how well their critical point predictions compared with experimental data. The Benedict-Webb-Rubin equation of state for propane (i.e. the con- stants a, b, ..., V were those for propane) was used as the basis for the reduced equation of state. By making the substitutions T I TcTr P I PcPr D I DcDr where Tc, Pc, Dc I critical temperature, pressure, and density of the reference compound (prOpane). Tr, Pr, Dr I reduced temperature, pressure, and density, respectively, the equation of state becomes Pr I Pr(Tr,Dr) (5) 3h By the corresponding states theorem, the function which expresses Pr in terms of Tr and Dr is the same general function for all pure compounds and mixtures of those compounds. Thus, although equation (5) was derived using the Benedict-Webb-Rubin constants and the critical constants for prOpane, this equation is valid for other compounds as well. When equation (5) is applied to mixtures the reduced pressure, temp perature, and density are computed from Pr I P/Pcm Tr I T/Tc m Dr I D/Dcm where Pcm, Tcm, Dcm I pseudocritical pressure, temperature, and density, respectively. Therefore, equation (5) written for mixtures is P/Pcm = Pr(T/Tcm, D/Dcm) (6) If the Benedict-Webb-Rubin equation of state is substituted into equation (1) and the resulting expression put into reduced form, the expression for computing G of a mixture becomes: G I G(T/Tcm, D/Dcm, x) "I from which G”' I G"°(T/Tcm, D/Dcm, x) (T) G"'0 I G'°°(T/Tcm, D/Dcm, x) - (8) Equations (6), (T), and (8) are now equations with only Tom and Pcm as constants since De I Pc /RTc Zc m m m where Zc I critical compressibility factor (equal to the same value for all compounds). Pseudocritical mixture rules generally express Tom and Pcm as func- tions of composition and the critical prOperties of the pure components. 35 Using published expressions for Pcm and Tcm, equations (7) and (8) were solved simultaneously for the values of T (i.e. Tk) and D (i.e. Dk) which made G" and G"' equal to zero. The method of computation and evaluation of results was the same as used.in the evaluation of equations of state. Generation of a Pseudocritical Mixture Rule This part of the research concerns the attempt to generate curves of Pcm as a function of x and Tom as a function of x from experimental Tk, Pk, and Zr data. Three different approaches were used. I. This approach used expressions for Tk(x), Zr(x), and Pk(x) where, for the pure components, (x I O, 1.0), the following was true Tk I Tc Pk I Pc Zk I Zc In addition, the fact that G" equals zero at the critical point was also used. The use of experimental data and the free energy condition to compute Tcm, Pen, and Dr is illustrated by the flow diagram of Figure A. For all of the flow diagrams in this chapter, it is assumed that Tc and Pc are known for the pure components (x I O, 1.0), but this fact is not shown on the diagrams. Also, Tk, Pk, and Zr are known functions of x. The diagrams illustrate only the general overall calculation procedure, excluding many of the minor details. Linear interpolation or extrapola- tion was used to adJust the variables in Figure h. 36 Start T At x I Ax At x I AxI 2Ax Use reduced eqn. of state to Assume values for Tr, Dr compute Pr, Zr. E ' Use Tk, Pk to compute Tcm, Pen. 1 .. At x I Ax No At x I Ax - Is calculated Zr I AdJust Dr experimental Zr? Ye“ At x . 2Ax - ‘ Ad ust Dr At x 2Ax ]N° Use reduced eqn. of state to _At x I 2Ax compute Pr; Zr. Is calculated Zr I Use Tk, Pk to compute Tom, Pcm. experimental Zr? I» (I Yes At x’I 2Ax No Is At x I Ax AdJust Tr G" I 0? Calculate first and second derivaties of Tcm and Pen. Yes Compute G". Repeat 1 through 2 for x I 3Ax, ..., 1.0 Yes t x I 1.0 No At x I Ax St°P ‘ Is To - Tk? AdJnu' at'T'r” Figure A. Flow diagram for calculation of Ten, Pen, and Dr using experi- mental Tk, Pk, and Zr data II. In this approach the experimental Zr data was not used. In- stead, the second free energy condition at the critical point, that (a3G/ax3)T P is equal to zero, was employed. Figure 5 shows a flow 0 diagram for the calculations in this approach. 37 Start l At x I Ax At x I AxI 2AxI §Ax Use the reduced eqn. of Assume values for Tr, Dr state to compute Pr. Use Tk, Pk to compute T07 . PC 0 . m m At x I Ax At x I 2Ax - Calculate first and second Use the reduced eqn. of derivatives of To and Po . state to compute Pr. . Compute G". m m Use Tk, Pk to compute Tc m’ . PcIn . ,, L No At x I 2Ax [13*9 ' 031 AdJust Tr Yes At x I §Ax At x I 2Ax Use the reduced eqn. of Calculate first, second, @ state to cOmpute Pr. and third derivatives of Use Tk, Pk to compute To and Fe . Compute 6" m m . TC m.Pcm e andG"'° T At x . :Ax +} No Are 0" and G"' AdJust Tr, Dr equal zero? Yes of x and compare with vAre Tc I Tk x I hAx,..., 1.0 C) experimental Zr curves. and Pa I Pk? AdJust Dr at x I Ax until a fairly good Compute Zr as a function’s,” At x I 1.0 Repeat®through®for corresponding curve is No AdJust Tr at x I 2Ax, obtained. 5 Dr at x I Ax. StOp Figure 5. Flow diagram for calculation of Tcm and Pcm using experimental Tk and Pk data 38 The maJor problem with this approach was in finding a convergence method which would find values of Tr and Dr that would make G " x+Ax x+Ax x and Gx"' equal to zero. The two methods which were used were the Newtone Raphson method!)2 and a method which assumed G" and G"' to be linear func- tions of Tr and Dr. The linear method is described in detail on page 107. III. In this part of the research, the form of the Tcm and Pcm curves was assumed to be cubic. Using the values of Tc and Po for the pure com. pounds, values of Tom and Pcm were computed at two intermediate values of x. The calculation flow diagram appears in Figure 6. The method used to adJust Tr and Dr was the Newton-Raphson method. Start At x I Ax, 2Ax At x I Ax, 2Ax At x I Ax, 2Ax Assume values of a-Compute Pr. i Compute first, second, Dr, Tr. Use Tk, Pk to and third derivatives compute Ten, Pcm. of Ten, Pom. ' Compute G", G"'. l At x I Ax, 2Ax No At x I Ax 2Ax AdJust Tr, Dr Are G", 0"! equal zero? tes Stop Figure 6. Flow diagram for calculation of To and Po assuming the form m m of these curves to be cubic. Numerical Check on Analytical Derivatives Because of the complexity of the equations for G" and G"', a pro- gram.was written which numerically checked the analytical expressions. For a function P(T,D,x), the derivative with respect to x at constant T fl 'a‘c w I and c lanalv \ U . r 4— 39 and P can be derived as follows: (3F/3‘)T,P . (aF/ax)T’D + (3F/3D)T,x(3D/aX)T,P (39/31)T’P . '(3P/3‘)T,n/(3P/3D)T,x e - (aw/3:0,“;zap/31>)...x (arr/3:0,.P - (3F/ax)T’D - (3F/3P)T,x(3P/3‘)T,D (9) where F can be 0, G', or G". The flow diagram used to compute these derivatives numerically is shown in Figure 7. Start At constant T, D Calculate numerical compute F and P. 4 values for (3F/3x)T D [Let F I G} 3' Then increase x and and (GP/ax), D ’ 9 compute F and P again. - ,_______nu________ Calculate (3F/3x) P Calculate At constant T,x and compare with ’ numerical values+--'compute F and P. analytical G'. for (aF/BP)T . Increment D and 7 ‘ ,x compute F and P _again. Repeat all steps for, G' and GM ; StOP Figure 7. Flow diagram for calculations in the program for checking analytical derivatives '11 ‘ 3:21:11 5 3829? I e :achiz: >- EXPERIMENTAL WORK Experimental Procedure This research involved extensive use of automatic digital computers. During the course of the work two computers were used. For the first year of work, the Michigan State University MISTIC computer was used. Input and output for this computer was by means of paper tape. Printed copies of the tapes were obtained from a teletype machine. After the first year of work on this problem, Michigan State Univer- sity purchased a Control Data Corporation 3600 Computer. This computer used Fortran programming language and was approximately eighty times faster than the MISTIC computer. The procedure for running a program.on this machine was (1) to construct the Fortran program, (2) to punch the program.on IBM cards using a card punch, and (3) to submit the program at the Computer Center for running. For both computers, the basic steps which eventually led to a cor- rectly Operating program were the same. The first step was to decide what the program was to accomplish. Next, the mathematical fermulas and techniques which were needed to obtain this goal were derived. The pro- gram.was written in such a way that the computer would perform exactly the desired mathematical Operations. The computer program was then run, and results were obtained. If these results were correct, or at least appeared correct, this whole process was a relatively simple one. However, if the program did not work properly, the process of removing errors was often very time consuming, especially if the program was long. For long programs, the procedure was to construct the program by parts. This hO all m hl allowed a constant check on each part of the program as it was constructed. With programs such as were written for this proJect, another element was involved in the correction process. Because many of the convergence methods were derived specifically for this proJect, there was always the chance that the basic equations had been derived incorrectly, or that even the method was not applicable to this research. If a method did fail to work, considerable time was spent trying to evaluate why the method failed, and whether there were any changes which could be made to correct the method. Often it was possible to alter the method to obtain a working program, but often no usable results were obtained. Experimental Critical Point Vapor-Liquid Equilibrium Data W. B. Kay's data for the ethane-n-heptane system23 were the specific experimental equilibrium data used. These data appear in Table 2. Values for critical compressibility factor (Zk) were computed from Zk I Pk/RTka. The pseudocritical compressiblity factor (Zcm) was assumed to be a straight line between the critical compressibility factor of the two pure components. Prausnitz and Gunn31 made a similar assumption for the acentric factor a and obtained results that were in good agreement with experimental data. Because Zcm is a linear function of w, the straight line assumption for Zcm has some Justification. Reduced compressibility (Zr) at the true critical point was then calculated from Zr 3 Zk/ZC o m Because the MISTIC library routine for a least squares fit of the data did not yield satisfactory results, a special program was written for fitting the data. In this program the method of least squares was A2 TABLE 2. Experimental data for ethane-n-heptane system at the true critical point (232h63) Weight Percent Mole Percent Tk Pk Dk 3 02 c7 c2 07 0F psi lbs./ft. o 100 0 100 513.3 396 1h.653 9.78 90.22 26.5h 73.h6 h68.2 682 16.22 29.91 70.09 58.71 h1.29 373.9 1106 16.62 50.2h h9.76 77.09 22.91 276.8 1263 17.h3 70.22 29.78 88.71 11.29 189.8 1132 16.97 90.22 9.78 96.83 3.15 120.3 850 15.h8 100 0 100 0 90.1 712 13.736 modified in order that the curve would exactly fit the data for pure ethane and n-heptane. This exact fit for the properties of the pure com- ponents was desirable because, by definition, the reduced temperature and pressure are unity at the critical point. That is, Tr I Tk/Tc I l Pr I Pk/Pc I 1 All degrees of polynomials up to six were fitted to the data and then evaluated on the basis of how well they predicted the data points. In all cases the fourth degree polynomials gave the best fit. The tem- Peramre curve was expressed as a function of mole fraction ethane. How- ever, the pressure and reduced compressibility factor curves were expressed as a function of weight fraction ethane because better fits were obtained using this variable. A computer program was written which computed Tk, Pk, and Zr for any value of x. Evaluation of Equations of State Because of their well known capabilities in both the vapor and liquid phases, the van der Waalsh3 and Benedict-Webb-Rubin3 equations of state were selected to be evaluated for their ability to predict vapor-liquid h3 equilibria in the critical region. The form of the van der Waals equation for a binary system (ethane-n-heptane) is: P = RTD/(l - bD) - aD2 (1) 1/2 al/2]2 1 + (l - x) 2 I van der Waals constant a where a I [xa for the mixture, 0' ll xb1 + (1 - x)b2 I van der Waals constant b for the mixture, 1’ b1 I van der Waals constants for pure ethane, a2, b2 I van der Waals constants for pure n-heptane, x I mole fraction ethane. a The Benedict-Webb-Rubin mixture equation is the same as equation (7) on page 29 except that the constants for the mixture are found by the following mixture rules: B I xB + (1 - x)B o o 1 °2 A a [xAl/2 + (l - x)Al/2]2 O 01 02 and likewise for Co and 7. 1/3 + ( 1 a1/3]3 1 - x) 2 a I [xa and likewise for b, c, and o. where subscript 1 represents ethane subscript 2 represents n-heptane. Each equation of state was substituted into the general expression for dimensionless free energy G (equation (6), page 26). The derivatives (32G/ax2)T P and (330/3: were analytically derived for both equations 9 3 )T,P of state. Computer programs were written which simultaneously solved the G" and G"' equations for the value of T and D which made the expressions Zerix State I d”sity 305.31): e, we.” 1., Pef uh equal to zero. The Newton-Raphson convergence method".2 was used in the calculations. These values of T and D were actually Tk and Dk because 0" and G"' are equal to zero only at the critical point. The critical pressure Pk was computed from the equation of state. These data from the two equations were compared with Kay's experimental data for the ethane-n-heptane sys- tem.2§‘ As a result of this comparison, the Benedict-Webb-Rubin equation of state was selected as the basis for the reduced equation of state. Derivation of the Reduced Equation of State The reduced equation of state was derived using the Benedict-Webb- Rubin3 equation of state for propane. This equation would be equation (7) on page 29 with the constants (a, b, ..., y) those for propane. By making the following substitutions: T I TcTr D I DcDr P I PcPr where Tc, Dc, Pc I critical temperature, density, and pressure of propane, an equation with only two constants can be derived. If Tom and Pcm are selected to be those two constants, then Dcm can be computed from DcIn I Pen/RTcch where Zc has the same value for all compounds, according to the corresponding states principle. Thus, the reduced equation of state expresses reduced pressure as a function of reduced temperature and density and has only pseudocritical temperature and pressure as constants. Computer programs were written which computed Pr and Zr given Tr, Dr, and the Benedict-Webb-Rubin equation constants and the critical constants for the reference compound (e.g. prepane). hS Derivation of the Analytical Expressions for the Critical Point Free Energy Relationships Using the general expression for dimensionless free energy G (equa- tion (6), page 26), the equation for G, when P is expressed by the Benedict-Webb-Rubin equation of state, was derived. The expression for G was then put into reduced form (i.e. made a function of two constants) by making the above listed substitutions for T, D, and P. The analytical equations for (320/3x2)T P and (330/3x were then derived. Computer ’ 3 )T,P programs were written in order that G" and G"' could be calculated on a digital computer. Evaluation of Available Pseudocritical Mixture Rules Three available pseudocritical mixture rules were evaluated by using them in conJunction with the reduced equation of state to compute the temperature, pressure, and compressibility factor of the ethane-n-heptane system at the point where the K-factor was unity (convergence pressure curve). The mixture rules tested were: 1. Kay's rule22 Tom I xTc + (l - x)Tc2 l Pcm I xPcl + (l - x)Pc2 where x I mole fraction ethane subscript 1 represents ethane subscript 2 represents n-heptane 2. Van der Waals' combinations”3 (Tc/Pc)In I x(Tc/Pc)1 + (l - x)(Tc/Pc)2 (Tc/Pcl/2)m . x(Tc/Pc1/2) + (1 - x)(Tc/Pcl/2) 1 2 3. Joffe's rule?1 ”In H. "c A6 1/2) 1/2) 1/2) (Tc/Pc m I x(Tc/Pc + (l - x)(Tc/Pc 1 2 (Tc/Pc)m = x2(Tc/Pc)l + (1 - x)2(Tc/Pc)2 + l/lt[(Tc/Pc)i/3 + (Tc/Pc):/3]3 x (1 - x) The procedure for each x was to use the expressions for Tcm and Pcm in the generalized equations for G" and G"' and then solve G" and G"' simultaneously for the value of T (i.e. Tk) and D (i.e. Dk) which made 0" and G"' equal zero. Pk was computed from the reduced equation of state. The mixture rules were evaluated on the basis of how well the values of Tk, Pk, and Dk, for several values of x, compared with Kay's experimental data.23 Computer programs were written for the evaluation of the three mix- ture rules. The Newton-Raphson convergence methodh2 was used in the three programs. Generation of Pseudocritical Mixture Rules Three different approaches were used in an attempt to generate curves of Tcm as a function of x and Pcm as a function of x. I. In this approach Kay's experimental critical point data23 were used along with the condition that 0" equals zero at the critical point, to compute curves of Tom and Pcm as a function of composition. A computer program was written for use on the MISTIC computer, and it followed the calculation flow diagram.in Figure h, page 36. In addition, several varia- tions of this approach were tried. The size of the increment Ax was varied, although 0.05 was the most commonly used value. Also, the computer program was changed to begin at the other end of the composition range (i.e. pure 1+7 n-heptane). Finally, butane was used as a reference compound instead of propane. This meant that the Benedict-Webb-Rubin constants were those for butane. To aid in obtaining convergence, expressions for the maximum possible values of Dr, given Trmax and Zr, were used to put some bound on new guesses of Tr and Dr. A check program was written by another person for checking results from this approach. II. In this part of the research the experimental reduced compress- ibility factor data for the true critical point were not used. Instead, the condition used was that the third partial derivative of free energy with respect to composition, at constant temperature and pressure, was equal to zero at the critical point. This part was begun on the MISTIC computer, but the arrival of the Control Data 3600 Computer made it necessary to rewrite all programs in Fortran programming language. Fol- lowing the flow diagram in Figure 5, page 37, computer programs were written for computing Tcm and Pcm at successive values of composition. The Fortran program was called Program LINEAR and is listed and described in Appendix VI. Because the equations for G" and G"' are functions of both Tr and Dr, the generation of values for Tr and Dr required a simultaneous con- vergence process for the two variables. The maJor problem connected with this part of the research was the development of such a convergence jprOcess. The first convergence scheme used was a process which assumed that G" and G"' were linear functions of Tr and Dr. The details of this convergence scheme are discussed in the description of subroutine CONVERG h8 in Appendix VI. A variation of the linear convergence method was also used. Since three sets of guesses for (Tr,Dr) were needed to predict a new set of guesses, the convergence scheme was changed so that the three sets of guesses included the last set and the best two previous sets. Before, the immediate past three sets had been used. Since the programs written for this approach had failed because of lack of convergence, two methods were used to obtain better initial s guesses for Tr and Dr at x2, x and x, , for use in Program LINEAR. The 3. as first method changed six variables (Trz, Tr3, Trh, Dre, Dr3, Drh) simultaneously until 0" was equal to zero at x and x and G"' was 2 3’ equal to zero at x These six variables were changed to meet these 3. three conditions by a method of steepest descent. The method of steepest descent involved calculating the partial derivative of a function e (defined as the sum of the squares of G2", 03", and 03"') with respect to each variable while the other variables were held constant. These slopes were used to compute the estimated changes in each of the six variables necessary to make 0 equal to zero. The equations used were a . -¢/[(a¢/axx1)2 + (8¢/axx2)2 + ... + (a¢/axx6)21 (1) AXXl I (at/axx1)p , etc. where XXl represents a dummy variable used for interpolation and “Subscripts on x represent a numbering system on successive values of x beginning with x1 I 0, x2 I Ax, etc. as Subscripts on variables other than x represent the values of these variables at the value of x which has the same subscript. 119 extraplation instead of Tr or Dr. There is one XX variable for each of Tr, Tr3 ,Trh, Dr2 Dru. p is a proportionality factor. ,Dr3, To speed convergence, a minimization procedure was used. This pro- cedure used the slope of the p versus ¢ curve at p I O and the value of ¢ at some other p value to fit a quadratic equation to these data points. The calculated minimum of this equation was used to compute the next values of 0. Except for a variation in the minimization procedure, the above convergence method is approximately the same as one described by Booth8. Because the speed of convergence was extremely slow with this method, several programs were prepared in an attempt to find the reason for this slowness. These included (1) programs for learning more about the ¢ versus 9 curve and (2) a program for determining the best increment of XX to use in the numerical differentiation calculations. The second method used an approach similar to the Newton-Raphson h2 method to simultaneously vary the variables (Tr, Tr3 ,Trh, Dr2 ,Dr3, Dru) until 0" and G"' were zero at x2 and x3. Since the fact that Tr I l and Pr I 1 at x I O I x1 was also used, four points were available to compute up to third derivatives of Tr and Pr. Because the six variables (Tre, Tr Trh, Dr2, Dr3, Drh) were varied to satisfy only 3’ four requirements (0" I O and G"' I 0 at x and x3), an infinite number 2 of solutions was possible. This permitted other requirements, such as non-negativity and well behaved curves, to be satisfied. The computed sets of starting data were then used as initial values in Program LINEAR. III. The third approach assumed a cubic form for the Tom and Pcm curves as a function of x. The computer program written for this section follows the calculation flow diagram in Figure 6, page 38. This program 50 was called Program TOTRANGE and is listed in Appendix VII. The Newton- Raphson convergence methodhz was used for solving a system of four equa- tions and four unknowns in the program. Check for the Analytical Derivatives A check program.was written to check the equations for G" and G"'. This program, called Program CHECK, follows the flow diagram in Figure 7 on page 39. It is listed and described in Appendix VIII. Chronological Order of the Research Work To provide greater continuity to the description of the research work, the discussion does not necessarily follow the same chronological pattern in which the work was done. The work was actually performed in the fol- lowing order: 1. Generation of pseudocritical mixture rules. 2. Evaluation of published pseudocritical mixture rules. 3. Evaluation of equations of state. Also, the placing of the computer programs in the appendices should not detract from the fact that they constituted the bulk of the research work. They have been included in the appendices to give continuity to the re- mainder of the report. RESULTS AND DISCUSSION Fitting of Curves to Experimental Data Early in the research work it was necessary to fit polynomials to the 23 to adapt this data for use on a experimental critical point data of Kay digital computer. Critical temperatures were expressed as a function of mole fraction, but critical pressures and reduced compressibility factors were expressed as functions of weight fraction because this variable yielded better fits. The polynomials which resulted are: Tk = (973.0 — 186.h392126r + 200.h13825x2 - 558.370h113x3 + 121.195736xh)/1.8 (1) Pk a (396. + 3051.927387w - 1199.7h6006w2 - 1.263.827785W3 + 2727.6h6362wh)/1h.696 (2) Zr 2 1.0 + h.03061128hw - 6.695218637w2 + 1.855922080w3 I + 0.808685291swh (3) where x I mole fraction ethane W weight fraction ethane I 30.06x/[30.06x + 100.17(l-x)] Tk critical temperature, °K Pk I critical pressure, atm. Zr critical reduced compressibility factor. Table 3 lists computed values of Tk, Pk, Zk, and Zr at .05 increments of mole fraction ethane. One-half weight fraction ethane is approximately 0.769 mole fraction ethane. TABLE 3. 52 Critical temperature, pressure, compressibility factor, and reduced compressibility factor as computed from a least squares fit of Kay's data .00 .05 .10 .15 .20 .70 .75 .80 .85 .90 .95 1.00 Tk (°K) 5h0.53888 535.60000 530.9910? 526.h9h58 521.90313 517.019h0 511.65618 505.63639 h98.79300 h90.96911 h82.01792 h71.80272 h60.19689 hh7.08393 h32.357h2 h15.92108 397.6886? 377.58h10 355.5h13h 331.50h50 305.h2776 Pk (Atm.) 26.9h6107 30.15h282 33.552553 37.1h8395 no.9h6990 uu.9u9aeu h9.152915 53.5h3590 58.096802 62.76911h 67.h903h8 72.151619 ’76.588316 80.556178 83.698999 85.507h51 85.275775 82.08393h 7h.900183 63.1088h5 h8.hh8557 Zk .25919000 .27526831 .2916h371 .30825131 .3250022h .3h177661 .358h1h01 .37h7016h .39035853 .h0501526 .h1818773 .h292hh55 .h3736661 .hh150205 .hh032381 .h3221376 .h153365h .38796h87 .3h9h7786 .30315790 .26395000 Zr 1.0000000 1.0610586 1.1231h9h 1.1860199 1.2h93263 1.3126070 1.3752h66 1.h36h310 1.h950881 1.5h98112 1.5987603 1.6395395 1.6690h52 1.6832977 1.6772835 1.6hh899h 1.579238h 1.h738291 1.326h226 1.1h95795 1.0000000 53 Equations of State for Mixtures 3 equations of state for The van der Waalsh3 and Benedict-Webb-Rubin mixtures were evaluated on the basis of their ability to predict conver- gence pressure curves and critical compressibility factors. Prediction of these quantities was used as a basis for evaluation because good predic- tion indicates computational accuracy in the critical region. The prediction of convergence pressure curves is of particular inter- est because the curves are used in the convergence pressure method of pre- dicting vapor-liquid equilibria. Because of a shortage of critical point data, the shape of convergence pressure curves is generally estimated from data for other systems. This is illustrated in Figure 8 which is a typical convergence pressure diagram published by Hadden.16 In this diagram only three of the curves are substantiated by experimental data. mo , , , 7“ '3’ ' T'!I"{'172.“7 .‘ , .. , , 7 .1 . . _ . ‘7 , ‘1 ,r‘ H , .- 3” )3" it; hz' 7‘72! A 7‘) L {(1 v.:‘ 'i ' 3'1”” 17v) ( . (1:1- ' wt 7. 012:: 1.1 “"78" 77-, 7:: 71503)”)? , '1' an v 1. “7 7 H 1. . 7 7 . ' '7 T3 (9’ 7+ ’ ' ‘ ‘1 ‘7 7 M77 7 7.. . u 7. . : ' 1 )9 v :I 71W) ’1 '1 r 7 7 s 7 a, _ < 1 Mess ' ' __ ‘ h i see - g in . - m .. .’ _ c ’ us ' ' . . Mo ‘ ‘ - Q . . . ‘5‘ . . “' r ' _ ”f" ~'~. ‘ .1 ,L, A c7 m ' , A 1§,YZ? 7° 3 4. Ollperlssntsl 7- - - - n - '7“ w" 7.. J; rm 1 6:73:27 , ‘1 t H. ; ': L r r; 1+1 ‘ A 1 ~ ' 7 -' ' ‘.' g x ' ‘ ; .1: 'EYL :119- (’7' “::7.- (1L4) VI»??- 77+) ‘71, 'H " if'7'Jf77rl’v Ir 71711-1 -' ‘4);3g-11'n]. if}; 7") 7 I” 100 no ‘00 SW momma °r Figure 8. Convergence pressure data - ethane lightest component 5h The convergence pressure curves shown in Figure 8 are plots of pres- sure as a function of temperature. Throughout the remainder of this dis- cussion, convergence pressure data will be plotted on two separate graphs. These graphs will be temperature and pressure as a function of composition. Because composition is a common parameter, the two separate plots are equivalent to a single plot such as in Figure 8. Using the van der Waals equation plus the rules recommended by van der Waals for computing equation constants for a mixture, the critical state properties shown in Table h were obtained. Table 5 lists similar results obtained by using the Benedict-Webb-Rubin equation of state with the combination rules recommended by Benedict £3.2l9 for mixtures. Figures 9, 10, and 11 are plots of the critical temperature, pressure, and compressibility factors, respectively, as a function of composition. The critical temperatures predicted by both equations of state agree very well with Kay's data (Table 2, page A2). The van der Waals equation was in considerably greater error in predicting pressures than the Benedict-Webb- Rubin equation. Table h shows the maximum percent pressure deviation for the van der Waals equation to be over 3h percent. The Benedict-Webb-Rubin equation was also in good agreement with known compressibility factors while the van der Waals agreement was relatively poor. Overall, the Benedict-Webb-Rubin equation of state for mixtures was in good agreement with Kay's experimental data. However, it should be realized that Kay's data were probably used in the empirical derivation of the Benedict-Webb-Rubin equation of state, and this could account for the excellent agreement. Nevertheless, the excellent results obtained with the Benedict-Webb-Rubin equation indicate that this equation of state is capable of being applied to pure compounds and mixtures with TABLE h. 55 Values of critical density, temperature, pressure, and compressibility factor as predicted by the van der Waals equation of state for mixtures 0.05 0.10 0.15 0.20 0.25 0.30 0.35 O.hO 0.h5 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 Dk (moles/liter) 1.682h28h 1.7h8709h 1.8202682 1.8977382 1.9818530 2.073h668 2.1735777 2.2833569 2.h0h1839 2.5376893 2.6858055 2.850825h 3.035h679 3.2h29h07 3.h769817 3.7h182h5 h.0h19h09 h.381llh8 b.7592177 Tk (°K) 535.07767 529.26780 523.07300 516.h5669 509.37809 501.79175 h93.6h71h h8h.88813 h75.h5271 A65.27276 h5h.27h19 hh2.37753 A29.h9938 h15.55h8h h00.h6158 38h.1h5h7 366.5%679 3h?.62230 327.32166 Dev. From Data (5) -0.10 -0.32 -0.65 -1.0h -1.h8 —l.93 -2.37 -2.79 -3.16 -3.h7 -3.72 -3.87 -3.93 -3.89 -3.72 -3.h1 -2.92 -2.23 -1526 Pk (fitne) 28.771692 30.6h2235 32.575061 3h.567900 36.616655 38.71h76h h0.852363 h3.015181 h5.18308h h7.328168 h9.h1228h 51.3838h0 53.173772 5h.690577 55.81hhh3 56.390708 56.223212 55.067867 52.623hh5 Dev 0 From Data (1) _ h.59 - 8.67 -12.31 ~15.58 -18.5h -21.2h -23.70 -25.96 -28.02 -29.87 -31.52 -32.91 -33.99 -3h.66 -3h.73 -33.87 -31.51 -26.h8 -16.6l Zk 0.389h2833 0.h03h0651 o.h1687237 Ooh29752h5 o.hh195956 0.h5338962 0.h639177l 0.h733931h 0.h816332o 0.h88h1521 0.h93h6673 0.h96h5332 0.h9696389 0.h9hh9366 o.h88h2613 0.h7801771 0.h6239363 0.hh05767h 0.h1160868 55 TABLE h. Values of critical density, temperature, pressure, and compressibility factor as predicted by the van der Waals equation of state for mixtures x Dk Tk Dev. From Pk Dev. From Zk (moles/liter) (°K) Data (1) (atm.) Data (1) 0.05 1.682h28h 535.07767 -0.10 28.771692 - h.59 0.389h2833 0.10 1.7h8709h 529.26780 -0.32 30.6h2235 - 8.67 o.ho3h0651 0.15 1.8202682 523.07300 -0.65 32.575061 -12.31 0.h1687237 0.20 1.8977382 516.h5669 -1.0h 3h.567900 -15.58 0.h29752h5 0.25 1.9818530 509.37809 -1.h8 36.616655 -18.5h 0.hh195956 0.30 2.073h668 501.79175 -1.93 38.71h76h -21.2h 0.h5338962 0.35 2.1735777 h93.6h71h -2.37 h0.852363 -23.70 0.h6391771 0.h0 2.2833569 h8h.88813 -2.79 h3.015181 -25.96 0.9733931h 0.h5 2.h0h1839 h75.h5271 -3.16 h5.18308h -28.02 0.h8163320 0.50 2.5376893 h65.27276 -3.h7 h7.328168 -29.87 0.h88h1521 0.55 2.6858055 h5h.27h19 -3.72 h9.h1228h -31.52 0.h93h6673 0.60 2.850825h hh2.37753 -3.87 51.3838h0 -32.91 0.h96h5332 0.65 3.035h679 h29.h9938 -3.93 53.173772 -33.99 0.h9696389 0.70 3.2h29h07 h15.55h8h -3.89 5h.690577 -3h.66 0.h9hh9366 0.75 3.h769817 h00.h6158 -3.72 55.81uhh3 -3h.73 o.h88h2613 0.80 3.7h182h5 389.1h5h7 -3.h1 56.390708 —33.87 0.h7801771 0.85 h.0h19h09 366.5h679 -2.92 56.223212 -31.51 0.h6239363 0.90 h.38111h8 3h7.62230 -2.23 55.067867 -26.h8 0.hh05767h 0.95 h.7592177 327.32166 -1.26 52.623hh5 -16.61 0.h1160868 56 TABLE 5. Values of critical density, temperature, pressure, and compressibility factor as predicted by the Benedict-Webb- Rubin equation of state for mixtures 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.h0 0.h5 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 Dk (moles/liter) 2.198882h 2.3lh9502 2.h500561 2.6163878 2.8679907 3.3h78922 3.6516231 3.9338281 h.22h0173 n.5331285 h.868h388 5.2362520 5.6h25121 6.092%091 6.588h252 7.12h3169 7.6667710 8.0986986 8.0519223 A Tk (°K) 537.15117 532.15117 528.36886 523.1931h 516.98032 508.h5078 500.59817 992.32922 h83.37107 h73.5680h h62.77669 h50.8h712 h37.61983 122.93209 h06.63972 388.67510 369.19h15 3h8.88235 328.97295 Dev. From Data (2) +0.29 +0.37 +0.36 +0.25 -0.0l -0.63 -1.00 -l.30 -l.55 -1.75 -1.91 -2.03 -2.12 -2.18 -2.23 -2.27 -2.22 -1.87 -O.76 Pk (atm.) 30.070395 33.15158h 36.h81h3h 10.1h0257 hh.hh5773 50.h2855h 55.910509 61.h27057 67.009787 72.55h512 77.91885h 82.860h62 87.03h909 89.953776 90.9h53au 89.1hh539 83.63009h 7h.057359 61.997280 Dev. From Data (1) -0.28 -l.20 -l.80 -l.97 -1.12 +2.60 +h.h2 .+5.73 +6.75 +7.50 +7.99 +8.19 +8.0h +7.h7 +6.36 +h.5h +1.88 -l.13 -1076 Zk 0.3102102? 0.3273938? 0.3h338080 0.35729879 0.365253h1 0.360970h9 0.3726780h 0.386h5976 0.3998665h 0.h1181255 0.h21h0233 0.h2767h23 0.h29h766h 0.h2537825 0.h1362313 0.39226573 0.36000691 0.3193666h 0.28518659 57 o.H m.o w.o ~.o scenes sowveeuu eao: 6.0 m.o 4.0 _ — — «.0 H.o m p- O 7 <3 J— 0mm H _ fi— epew a.»ex .eem m-:-m cusses scupeenm eaos q . A .mom .qudz new as» «o scavessh e as suspeaeaaeu Heowawuo .a shaman 1700M 1 70mm 1 foo: e..se Iowa 1 4100a .-emm enemas mowuoenu eaoz 58 0.0 h.0 0.0 m.0 3.0 m.0 N.0 H.0 0 _ b — P P P b — b O a . _ I _ . a a « LTON .cvm .naeez you es> . . 110: . . .. .3. ..E H .IYIOW seep e.bdM . . 1700 . . . .eem muzum 1700a enemas no«voean eaoa no mauvemmh s as ease-sun Heoukuo .OH skewna enemas sewaeeau sac: 0.0 ~.0 0.0 m.0 :.0 m.o «.0 H.0 0 _ r _ 0.H 4 O‘ ,_ e O q — A I Lr~.o .eam muzum JT:.O «use ..aae .svm .oaee: use se> 1rm 0 enemas nowvoeum sacs no ceavoesu e ea haunanweeeuaaoo Heowuuuo .HH shaman 60 approximately the same accuracy. Generalized Equation of State A generalized equation of state expressing reduced pressure in terms of reduced temperature and density was derived. The form of the equation was based on the Benedict-Webb-Rubin equation of state. The resulting reduced equation is Pr I fTrDr + (gTr - h - i/Tr2)Dr2 + (JTr- k)Dr3 + lDr6 + (mDr3/Tr2)(l. + nDr2)exp(-nDr2) (h) where constants f through n are the reduced equation of state constants ' and are defined in Appendix III. According to the corresponding states principles, these constants should be the same for all compounds. How- ever, in reality, these constants do vary slightly depending on the refer- ence compound used to compute them. In this research propane was used as the reference compound unless otherwise noted. Generalized Equations for Free Energy and its First Three Derivatives Because the second and third derivatives of free energy (G) with respect to composition, are equal to zero at the critical point, it was desirable to have generalized expressions for these quantities. The generalized expression for G was found to be G I 1n PcmPr + x 1n x + (1. - x)ln(l. - x) + Z - l. - ln Z +-[(g - h/Tr - i/Tr3)Dr + (J - k/Tr)Dr2/2 + (l/Tr)Dr5/5 2 + (mDr /2Tr3)[(2 - 2exp(-nDr2))/nDr2 - exp(-nDr2)]]/f (5) Details of the above generalization appears in Appendix III. The equations for G" and G"' and the procedure for obtaining them are shown in Appendices III, IV, and V. 61 Evaluation of Published Mixture Rules Using the derived generalized equation of state and generalized ex- 22 van der Waals',“3 and Joffe's21 mix- pressions for GM and G'°', Kay’s, ture rules were evaluated for their ability to predict critical state properties. Tables 6, 7, and 8 list the predicted values of critical temperature, pressure, and compressibility factor obtained with each of the three mixture rules. These results are plotted as functions of comp position in Figures 12, 13, and lb. The critical temperature agreement was very good in all cases. The maximum temperature percent deviation was h.67 percent; this occurred with Kay's rule. The predicted critical pressures were best with van der Waals0 mixture rule and poorest with Kay's rule. Joffe's rule, although slightly poorer than van der Waals' combinations, was nevertheless in good agreement with Kay's data. Figure 1h shows only fair agreement with the experimental compressibility factors. Van der Waals9 and Joffe"s rules are approximately equivalent in their ability to match the experimental compressibility factors. Overall, the van der Waals9 mixture rule appears to be the best of the three mixture rules evaluated. In fact, the critical temperature and pressure predictions made by using the van der Waals mixture rule in the generalized equation of state were about as good as the predictions made from the Benedict-Webb-Rubin equation of state for mixtures. The average critical temperature error was 1.16 percent using the van der Waals rule with the generalized equation and 1.32 percent using the Benedict-Webb- Rubin equation for mixtures. The average critical pressure error was 5.hl percent with the reduced equation of state plus van der Waals' rule and h.25 percent with the eight constant Benedict-Webb-Rubin equation. The relatively good temperature and pressure agreement obtained with TABLE 6. 62 Critical temperatures, pressures, and compressibility factors and critical reduced densities obtained from the generalized equation of state using Kay’s mixture rule 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.h0 0.h5 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0°95 Dr 0.9231651 0.8836010 0.8713950 0.8769hh2 0.89h6219 0.9206883 0.952hh9h 0.9878553 1.02529h5 1.063h892 1.101h103 1.1381hh7 1.1726882 1.2036582 1.22892h1 1.2h51693 1.2h73h05 1.2277h00 1.17153h5 Tk (°K) 5h0.881h6 5h1.739h2 5h2.17612 5h1.5h189 539.36397 535.32072 529.2332h 521.03687 510.75hh9 h98.h6689 h8h.29001 h68.36h58 h50.8620h h32.00hh6 A12.08810 391.h8666 370.59239 3A9.6h951 328.A573h Dev. From Data (%) +0.99 +2.02 +2.98 +3.76 +h.32 +h.63 +h.67 +h.h6 +h.03 +3.hl +2.65 +1.77 +0.85 -0.08 -0.92 -1.56 n1.85 -1.66 -0.92 Pk (atm.) 32.081979 37.h29280 h3.160219 h9.3h8115 55.981h93 62.97501A 70.179159 77.385910 8h.331900 90.701008 96.126566 100.19h502 102.h5hhh3 102.h51997 99.7975h8 9h.2750h3 85.955295 75.221599 62.618319 Dev. From Data (%) +6.39 +1l.55 +16.18 +20.52 +2h.5h +28.12 +3l.O7 +33.20 +3h.35 +3h.39 +33.23 +30.82 +27.18 +22.h1 +16.71 +10.55 + h.72 + 0.h3 - 0.78 Zk 0.32837705 0.37630737 0.h1h3h263 o.hhhh9636 0.h683h217 0.h8701h30 0.50125h31 0.5111635h 0.5177h698 0.5199h559 0.51766700 0.51032823 0.99722998 0.h776917h 0.h512833h 0.h1815902 0.379h0017 0.33717536 0.29508380 63 TABLE 7. Critical temperatures, pressures, and compressibility factors and critical reduced densities obtained from the generalized equation of state using van der Waals' combinations x Dr Tk Dev. From Pk Dev. From Zk (°K) Data (1) (atm.) Data (5) 0.05 0.9768298 537.h5599 +0.35 30.6hhhh0 +1.63 0.3062h680 0.10 0.963326h 531.3691h +0.6h 3h.h88729 +2.79 0.33903h93 0.15 0.958102h 531.05083 +0.87 38.532570 +3.73 0.36912509 0.20 0.9595236 527.30369 +1.03 h2.815067 +h.h6 0.39668929 0.25 0.966ho8h 522.9536h +1.15 17.359121 +5.36 0.h218h106 0.30 0.9779152 517.78913 +1.20 52.168630 +6.1h 0.hhh59261 0.35 0.9933872 511.65693 +1.19 57.221872 +6.87 0.h6h8hhhs 0.h0 0.0122h75 50h.36183 +1.12 62.h61h97 ,+7.51 0.h8236720 0.95 1.0339289 “95.71973 +0.97 67.781h20 +7.99 0.h9677399 0.50 1.0578205 h85.55h37 +0.73 73.011208 +8.18 0.50798858 0.55 1.0832113 h73.7093h +o.h0 77.899690 +7.97 0.51371798 0.60 1.1092170 h60.06761 -0.03 82.102129 +7.20 0.51hhh616 0.65 1.13h6673 hhh.58183 -0.56 85.179809 +5.7h 0.508h7929 0.70 1.1579352 h27.31567 -1.17 86.62713h +3.50 0.19h59062 0.75 1.1766970 h08.h8h2 -1.79 85.9hh862 +0.51 0.h7181998 0.80 1.1876393 388.h9618 -2.31 82.766929 -2.9h 0.h3993903 0.85 1.1861305 367.86hhh -2.57 77.0082h8 -6.18 0.39996076 0.90 1.1656h61 3h7.09781 -2.37 68.9h9631 -7.9h 0.35hh5651 0.95 1.1152699 326.h5019 -1.52 59.191060 -6.21 0.30786518 61 TABLE 8. Critical temperatures, pressures, and compressibility factors and critical reduced densities obtained from the generalized equation of state using Joffe°s mixture rule x Dr Tk Dev. From Pk Dev. From Zk (°K) Data (5) (atm.) Data (1) 0.05 0.9727908 537.69519 +0.39 29.889297 -0.88 0.29790306 0.10 0.9511113 531.81787 +0.73 32.996163 -1.66 0.32103773 0.15 0.9362128 531.88385 +1.02 36.302932 -2.28 0.31880851 0.20 0.9265177 528.69177 +1.30 39.813636 -2.69 0.37217618 0.25 0.9211978 525.16053 +1.57 13.650313 -2.89 0.39117022 0.30 0.9206871 521.17216 +1.86 17.751517 -2.85 0.11177931 0.35 0.9238190 516.59191 +2.17 52.171167 -2.56 0.1339313? 0.10 0.9308000 511.27712 +2.50 56.926158 .-2.01 0.15118313 0.15 0.9116617 505.01261 +2.87 62.016929 -1.20 0.16720818 0.50 0.9565150 197.68220 +3.25 67.118917 -0.11V 0.18076911 0.55 0.9756109 188.95109 +3.63 73.062511 +1.26 0.19167708 0.60 0.9991835 178.56353 +3.99 78.803787 +2.89 0.19922581 0.65 1.0273678 166.19395 +1.27 81.379008 +1.75 0.50210232 0.70 1.0602171 151.18901 +1.12 89.338656 +6.71 0.19976137 0.75 1.0973090 131.10251 +1.37 92.962651 +8.72 0.18930956 0.80 1.1371929 113.77168 +1.01 91.17901 +10.11 0.16811056 0.85 1.1761191 390.15139 +3.11 91.567138 +11.25 0.13101523 0.90 1.2018912 361.18166 +2.52 83.588113 +11.60 0.38386639 0.95 1.198557% 336.58129 +1.53 69.258505 +9.7h 0.31938622 65 0.H escape souvoeuh eHoz 0.0 0.0 v.0 0.0 m.0 3.0 m.0 «.0 H.0 0 _ _ a P . a I r. _ _ I 1 d 3 1 4 J) d -700m [roam an a so he nor 3 .ca. on some 0H . a 3 0 new m 0 was 0 1700: e..ea 1T0ma ease n.000o0 new: .:00 ceuwaeseseo 60.0.0 wing 700m ease e.hex new: .sdm denudeaeeeo 1mm escape soHIosuu «Hos Mo scavenge e as suaveuenae» Hsowawuo .NH saunas ..( N‘llluls‘ nsl‘v‘vtlsk 70H enemas soauoeam «Ho: ease «.00000 new: .eum veuaaehoeec spec u.heg hem 0H3 0., gurus» .svm confideaeaeo ease .maeez new as» new: .sdm veaudsueseo enemas soHuueum eaoa mo ecuuommh e as ens-seam Heowvuho .ma shaman oJ - To m6 To 06 To ...o To To to o _ _ a _ p _ I I. . a J. a . 47 I 4. 7e a 0 .8. Tom .3. ..E 700 woos '((“‘IC 67 enemas segues.“ sac: 0.H 0.0 0.0 e.0 0.0 “.0 3.0 0.0 «.0 H.0 0 . _ a _ _ . p 7 p a H.0 a q d d I q 4 q # iTN.0. A. case ...0000 near .eem uouaasuoeoo . 0.0 spec s.hiu 1T..0 . . 170.0 sass Lassa use as» . eves .eua eased-souso .Hsh ..aee ease .eum 0ouaastoeo0 110.0 enemas segues.» sacs no 8305C s .e.e ham—«3:93.00 133.6 ...H 93»: Wine tian: Value Vi: 68 the van der Waals9 combinations indicates that this mixture rule may be very satisfactory for vapor-liquid equilibrium calculations. However, the poorer compressibility factor agreement indicates that the predicted densities are in error from the eXperimental densities. This lack of agreement means that the mixture rule probably would not be entirely satisfactory for predicting other thermodynamic quantities from a gen- eralized equation of state. Generation of Pseudocritical Mixture Rules The results from the first two approaches which were used to compute Tcm and Pcm values were not completely satisfactory. They are listed along with a discussion of the data in Appendix IX.. By assuming the curves of Tcm and Pcm to be cubic (Program TOTRANGE), values of Tom and Pcm were computed at two intermediate values of composi- tion. Table 9 lists the results from program TOTRANGE for two pairs of x values. The pairs used were x 8 0.6215, 0.8723 (approximately symmetrical with respect to weight fraction) and x = 0.333, 0.667 (symmetrical with respect to mole fraction). TABLE 9. Values of pseudocritical temperature and pressure as computed by program TOTRANGE x Tc Pc Graph Symbol 0.333 h5h.11823 °K 28.105853 atm. A 0.667 378.39011 35.058287 15 0.6215 39h.77h62 35.821393 0 0.8723 336.30273 h3.010320 D Figures 15 and 16 are plots of pseudocritical temperatures and 0“ ‘ 69 pressures as a function of composition. Figure 15 shows the Tcm curves from program TOTRANGE and van der Waals9 mixture rule to be nearly linear. The generated pseudocritical pressure curve (Figure 16) is reasonably close to the curve predicted by the van der Waals' mixture rule, but both of these curves deviate considerably from linearity (Kay's rule). Using the computed reduced compressibility factors from program TOTRANGE, the values of critical compressibility factor were calculated and plotted in Figure 17. The agreement with Kay°s data is only fair but it is better than the agreement obtained using van der Waals' mixture rule in the generalized equation of state. Errors were undoubtedly introduced because the computations in program TOTRANGE involved fitting a third degree polynomial to only four points and then taking up to third derivatives of this polynomial. This numerical procedure would explain the variation of the values at x = 0.6215 and x = 0.667 in Figures 15, 16, and 17. If more accurate pseudocritical curves could be generated by using smaller increments of composition then, theoretically, the computed curve should come closer to the data curve in Figure 17. Differentiation Check The results of a program written to numerically check the expressions for GM and G“9 are shown in Table 10. Agreement to at least three sig- nificant figures was obtained between the numerical and mathematical cal- culations of the derivatives. .eeenae scanning ode: m.o m.o v.0 v.0 m.o 3.0 m.o «.0 H.o o P a _ _ n .4 47— TO -4|- — L _ n7 877 . . 10mm J J1oom . , sash assuage .uHsdz you ns> soak Boa 70mm :— moz I. > iron a O 8 o .c m. .\ D I can» e.»iu luau um 1102 I . - . e. o 5352. so...“ and 5 a 110m ends». scuooduu one: he scavonsu s as aoa<2/Yn2 - [2 exp<~yn2)1/Yn2 - ...(...2))] (5) Putting (5) into reduced form G I 1n PcmPr + x In x + (1 - x)ln(1 - x) + Z - 1 - 1n Z + 3 (Bo - Ao/RTcTr - Co/RTr Tc3)DcDr + (b - a/RTcTr)D02Dr2/2 + aoDcSDrS/SRTcTr + (c0r20c2/2RTc3Tr3)[(2-2exp(-y0e20r2))/ yDc2Dr2 - exp(-yDc2Dr2)] (5) Note: BoDc I g/f AoDc/RTC I h/f C Dc/R'I'c3 O I i/f 2 ch /2 = J/2f aDc2/2RTc -Jh/2r aoDcs/SRTc I 1/5f 3 ch2/2RTc I m/2f c/ZYRTC3 I m/2nf yDc2 I n Then, G I 1n PcmPr + x 1n x + (1. - x)ln(1. - x) + Z - l. - 1n 2 +[(g - h/Tr - i/Tr3)Dr + (J - k/Tr)Dr2/2 + (l/Tr)Dr5/5 +(mDr2/2Tr3)[(2-2exp(-nDr2))/nDr2 - exp(-nDr2)]]/f (7) Let 91 TO ' le/Tr (8) To' I -Tr'/Tr2 (9) To" I 2(Tr')2/Tr3 - Tr"/Tr2 (10) To"' a -To[3(To'Tr" + To"Tr’) + ToTr"'] (11) E T . 1./Tr3 = T 3 (12) 1 n O .J 2 . Tn' . 3To To' (13) : 2 2 . Tn” . 6T°(T°') + 3T0 To“ (1") :. T "' - 3[2T ' + 6T T 'T v' + T 2T m] (15) g n O O O O O O - G I x In x + (l.-x)ln(l.-x) + Z - 1. - 1n Z + 1n PcmPr 2 S +[(g - hTo - iTn)Dr + (J - kT°)Dr /2 + 1T°Dr /5 + an/n + (m/2n)Tnexp(-nDr2)(-2 - nDr2)]/f (16) 0' . (%%.T,P = 1n[x/(1, - x)] + 2' - z'/z+[(g - hTo - iTn)Dr' 2 h -(hT°' + iTn')Dr + (J - kTo)DrDr' - kTO'Dr /2 + 1T°Dr Dr' + lTo'DrS/S + an'/n + (m/2n)°exp(-nDr2) (~2Tn' - nDreTn' + 2nDrDr'Tn + 2n20r30r'Tn)][f (17) 2 G" ’ (%)T p = 1./(x - x2) + z“ - z"/z + mm2 + [sDr" 3x ’ -h(T Dr°° + T "Dr + 2T 'Dr') - i(T Dr" + T "Dr + 2T 'Dr') 0 O O n n n + JEDrDr': + (Dr')2] - k[ToDrDr" + T°"Dr2/2 + 2To'DrDr' + T°(Dr')2] + 1[2T°'Dthr' + To"Dr5/S + hToDr3(Dr')2 h + ToDr Dr") + an"/n + (m/2n)exp(-nDr2)(-2Tn" - nDr2 II!!! n + hnDrDr'Tn' + 2nDrDr"Tn + 2n(Dr')2Tn + hnaDrBDr'Tn' 3 + 2n2Dr Dr"Tn + 2n2Dr2(Dr')2Tn - hn3Drh'(Dr')2Tn]/f (18) 92 GH: ,(___, 3x 3 T, P = (2x - 1.)/(x - x2)2 + Z9“ - Z°"/Z + 3Z°Z"/Z2 - 2(Z°/Z)3 + [gDr"' - h(ToDr'°' + To"'Dr + 3Dr'T°" + O! O _ '90 '9? + ! t0... !! ! 3Dr To ) i(TnDr + Tn Dr 3Tn Dr 3Tn Dr ) + J(DrDrHo + 3Dr°Dr°°) - k(TODrDr"' + To"'Dr2/2 + 3T0Dr°Dr" + 3To'DrDr'° + 3To°°DrDr' + 3To'(Dr')2) 1(3To"Dthr' + 12To'Dr3(Dr')2 + 3To'Dthr" 3 + + 12 ToDr2(Dr')3 + 12 ToDr Dr°DrH + ToDthr"' To°'°DrS/5 + an°"/n + (m/2n)exp(-nDr2)(-2Tn"' + nDr2Tn "9 + 6nDrDr°Tn V” + 6nDrDr“Tn ' + 6n(Dr')2Tn' 3 + 2nDrDr"'Tn + 6nDr°Dr"Tn + 6n2 Dr Dr'Tn" + 6n2' 3 Dr3Dr'°Tn° + 6n2 Dr 2(Dr9) 2Tn° + 2n 2Dr Dr"'Tn 3 h 2Dr°Dr"Tn - 12 n3Dr"(DrV)2Tn ° - 12 n Dr . + 6n2Dr Dr'Dr°°Tn n20n3Dr3(Dr"3 ') Tn + 8n hDr5(Dr° )3'1‘n )]/f (19) See Appendix IV for derivation of equations for Dr', Dr", and Dr"'. See Appendix V for derivation of equations for Z, Z', Z'“, and Z"'. APPENDIX IV Derivation of Equations for Dr', Dr°°, and Dr'°' Taking the derivative with respect to x of equation III-3: Pr' I [fDr + (g + 2i/Tr3(Dr2 + JDr3 - (2mDr3/Tr3)(l. + nDr2) exp(-nDr2)]oTr° + [fTr + (gTr - h - i/Tr2)2Dr + (JTr - k)3Dr2 + 6(1)Dr5 + (3mDr2/Tr2)(l. + nDr2)Oxp(-nDr2) + (mDr3/Tr2)(2nDr)exp(-nDr2) + (mDr3/Tr2)(1. + nDr2) (-2nDr)exp(-nDr2)]Dr' Let U I Pr' -[fDr + gDr2 + 2iDr2/Tr3 + JDr3 - 2mDr3 (l. + nDr2) exp(-nDr2)/Tr3]Tr° (1) V I fTr + 2gTrDr - 2hDr - 2iDr/Tr2 + 3JTrDr2 - 3kDr2 + 6(1) DrS + (m/Tr2)exp(-nDr2)(3Dr2 + 3nDrh - 2n2Dr6) (2) therefore, Dr' 2 U/V . (3) also, Dr°° = U°/V - UVV/Vz a (U0 - Dr°V°)/V (A) Dr“0 = (UM - Dr'VH - Dr°“V')/V - (Uo - Dr'V')V°/V2 I (U°° - Dr°V°° - Dr°°V°)/V - Dr'°V'/V - (U°° - V°"Dr° - 2V°Dr°°)/V (5) Let It a Tr'/Tr3 a TrVTn (6) v 3 0! v o Tt Tr Tn + Tr TD (7) T '9 I Tr°°'T + Tr"T ° + Tr'VT ' + Tr'T '° t n n n n = Tr'°°T + Tr'T 0° + 2Tr"T ° (8) n n n Therefore 93 01‘ U I Pr0 - fDrTro - gDr2Tr' - 2iDr2TnTr" - JDr 9h 3Tr' + 2mDr3(1. + nDr2) exp(-nDr2)TnTr' U' I Pr" - fDrTr'° - fDr'Tr' - gDr2Tr" - 2gDrDr'Tr' - 2iDr2TnTr" U9 V! - 2iDr2Tn'Tr° - hiDrDr'Tt - JDr3Tr" - 3JDr2Dr'Tr' + 2mDr30 (1. + nDr2)exp(-nDr2)Tt° + 2mDr3(1. + nDr2)(-2nDr)Dr'exp(-nDr2)Tt + 2mDr3(2nDrDr")exp(-nDr2)Tt + 6mDr2Dr'(1. + nDr2)exp(-nnr2)Tt Pr'° - fDrTr°° - fDr°Tr' - gDrZTr°° - 2gDrDr'Tr' - 2iDr2Tt' - hiDrDr°Tt - JDr3TrH - 3JDr2DrVTr° + 2mDr3(1. + nDr2) 2 2 h exp(-nDr2)oT ' + 2moexp(-nDr2)Tt(3 + 3nDr - 2n Dr )Dr2Dr' t (9) fTr' + 2gDr'Tr + 2gDrTr' - 2hDr' - 2iDr'/Tr2 + hiDrT t + 6JTrDr'Dr+3JTr°Dr2 - 6kDrDr° + 30(l)Dthr' + moexp(-nDr2) (6DrDr9/Tr2 + 6nDr3Dr9/Tr2 - l8n2DrSDr'/Tr2 + hn3Dr7Dr'/Tr2 2 h 2 6 . - 6Dr Tt - 6nDr Tt + hn Dr Tt) (10) UN I Pr°°' - fDrTr°'° - 2fDr°Tr°° - fDr"Tr° - gDrzTr"' thrDr'Trflo - 2gDrDr°°Tr' - 23(Dr°)2Tr9 - 2iDr2Tt" 81DrDr°Tt° - hiDrDrWTt - hi(Dr°)2Tt - JDrBTr"' 6JDr2Dr°Tr°° - 3JDr2Dr'VTr° - 6JDr(Dr')2Tr' . + [2m°eXP(-nDr2)(Dr3Tt°' + nDrsTt'° + 6Dr2Dr'Tt' + 6DrhnDr°Tt' - hnzDr6Dr°Tt' + 3Dr2Dr"Tt + 3nDthr"Tt 2 6 2n Dr DrV'Tt + 6nDr3(Dr“)2Tt + 6Dr(Dr')2Tt 18n2Dr5(Dr°)2Tt + hn3Dr7(Dr°)2Tt] (11) 95 V99 I fTrH + 2gDr°°Tr + thr°Tr° + 2gDrTr°° - 2hDr" - 2iDr"/Tr2 + 81Dr'Tt + hiDth' + 12JTr°DrDr' + 6JTr(Dr')2 + 6JDrDr"Tr + 3JTr°'Dr2 - 6k(Dr°)2 - 6kDrDr" + 30(1)Dthr'{. 2 + 120(1)Dr3o(Dr°)2 + 2moexp(-nDr2)(3DrDr"/Tr2 + 3nDr3Dr"/Tr 9n2Dr5Dr°°/Tr2 + 2n3Dr7Dr”/Tr2 + 3(Dr')2/Tr2 + 3nDr2(Dr°)2/Tr2 - 51n2Drh(Dr')2/Tr2 + 32n3Dr6(Dr')2/Tr2 hnhDr8(Dr°)2/Tr2 - 3Dr2T ' - 3nDrhT ° + 2n2Dr6T ' t t t . 3 7 12DrDr Tt - l2nDr t Dr°T + 36n2DrSDr°Tt - 8n3Dr Dr'Tt) (12) Let Then APPENDIX V Derivation of Equations for Z, Z', 29', and Z"' Z I P/RTD I PrPc/RTrTcDrDc I (Pr/TrDr)(Pc/RTch) I Pr/fTrDr Y I Pr/fTr Y' :- (Pr' - fYTr")/(fTr) Y" I (Pr" - fYTr°' - 2fY°Tr°)/(fTr) Y°°' = (Prm - fYTr“9 - 3fY'Tr°° - 3fY°°Tr')/(fTr) Z I Y/Dr Z' I (Y' - ZDr')/Dr Z" I (Y" - ZDrH - 2Z'Dr')/Dr zone = (Yooo _ ZDr°'° _ 3Z°Dr'° _ BZ'°Dr°)/Dr 96 (1) (2) (3) (h) (S) (6) (7) (8) APPENDIX VI Program LINEAR The objective of program LINEAR is'to vary simultaneously the value of reduced temperature (Tr) and reduced density (Dr) at successive values of composition (x) until the second and third derivatives of the free energy, 0" and G"'“, become equal to zero. This program follows the flow diagram in Figure 5, page 37. Starting at one end of the composition range (x1 I 0), the boundary conditions Tr I 1 and Pr I l are used. The value of composition is incre- mented twice and guesses are made for Tr and Dr at these compositions (x2, x3). Then by simple linear interpolation the value of Dr3 (Dre, ".is equal to Tr2,or Tr could also have been chosen) is changed until G 3 2 zero. Guesses for Tr and Dr are then made for the next composition (xh). These guesses (Dru, Trh) are changed simultaneously by a linear convergence method until 03" and G3"° become equal to zero. This process is repeated until the entire composition range has been covered. The only restriction on the values of Tr and Dr which are selected is that they be positive since negative values have no physical signifi- cance. The main divisions of program LINEAR are: 1. Program LINEAR-Main Body 2. Subroutine CONVERG 3. Subroutine POLY h. Subroutine WBKAY 5. Subroutine PRCALC 6. subroutine G2CALC 7. Subroutine G3CALC 80 Data 97 98 PROGRAM LINEAR DIMENSION X(25)e W(25)9 PC‘ZS). TC(25)9 PK(25)9 TK!25)9 PR‘ZS). 1 TR(25)0 DR(25)0 ZR(25)9 2(5). FUNC(5)9 DELTAD(4)0 DELTAT‘4).G(6) 2 CPZRWRPR(25) COMMON I oDCBt TCB oAC.BC .CC .DC .EC . C .Cl 9C2 .ij .C 1+ 9 Cli- . Cb .(LY . C‘fi -, C'F' .0253 . 1 DR 1 DER o DRZDER 0 TR 9 TR I DER t TRBDER 9 TRIJDER . PR 9 PR 1 DER tPHZDhR t Fifi-{37.25.72 , 2 pCOpKOTCOTKOQQOTOQTOIDEROTOZDERQTNQTNIDERQTNEDER'TTQTTIDERQVQ 3 VlDEReYeYlDERoYZDEReZKtZKlOERoZKZOERoZRQZeXoWeGZDEReGBDERQFUNC 4 ePZRWRPR 1 FORMAT (513.11: '- 2 FORMAT (ZSXQZSHDOUBLE LINEAR CONVERGENCEt4(/)) 3 FORMAT (I2) 16 FORMAT (6X05E21e10) 4O FORMAT(20H NO OF ITERATIONS I 0I29/) 41 FORMAT(IHI) 42 FORMAT (2/915X08HFUNCTION013X09HFIRST DEROIIXOIOHSECOND DEROIZXO 1 9HTHIRD DERo/93H TC93X04EZIeIOo/93H PC03X04E21eIOe/03H TRQSXe 2 4EZIeIOs/03H PR03X04EZIelO) 43 FORMAT (leIlXeIHXelBXoZHTRQlBXoZHPRsIBXoZHDRo18X92H2R015X9 l 7HPZRWRPRe/06E20eIOQ/eéEZOeIO) 44 FORMAT (9H GZDER I eElBe1003X98H63DER I vElBelO) 45 FORMAT (25H DOES NOT CONVERGE AFTER QIZOIIH ITERATIONS) 46 FORMAT(/ZOH CALCULATION AT X I 0F5039//17X92HTR919X02HDR917X0 I SHGZDER916X05HG3DER) 50 FORMAT (4H PK(01294H) I 9E18e1004X03HTK(01204H) I 0E18e10004X0 1 4H PC(91294H) I cEIBelOo4XoBHTC‘01294H) = oElBelO) 51 FORMAT (11H INPUT DATAe/ogH TR(2) x oEIBelOsSXeBHTR(3) 8 9E18.10. l 5X08HTR(4) I eElBeIOs/09H DR(2) I 0E1801005X08H0R(3) I tEIBelO! Z 5X08HOR(4) 3 OEIBeIOel) PRINT 2 READ lo DCBOTCBQPCBIRQBOQAOoCOoBerCoALPHAeGAMMAoACoBCtCCQDCoEC PRELIMINARY CALCULATIONS C! I (DCBITCB*R)/PCB C2 I C1*80*DCB C3 I (AO*DCB**2)/PCB C5 I C1*DCB§*2*B C6 I (DCB**3/PCB)*A C? I C6*DCB*53*ALPHA C8 I (Cé/(AiTCB**2))*C C9 I DCB**2*GAMMA C13F9 C2369 C3IH0 C4=Is C5=Jo C6=Ke C7=Lo C8=Mo C9=N TOL I leE‘IZ X(l) = O X(2) I 0e025 X(3) = Oe05 TR(1) I 100 DR‘I) ' IeO PR(1) 3 leO ZRII) I leO DELTAX 3 OeO5 NUM I (Is - X(3))/DELTAX + 3e DO 6 II4ONUM 6 X(I) 3 X(I-1) + DELTAX DO 5 IIIONUM 99 5 CALL WBKAY I = I PC(I) I PK(I)/PR(I) TCII) I TKII)/TR(I) READ 39 NUMDATA DO 25 NN=I¢NUMDATA 4 READ Io (TR(I)oI I 204). (DR(I)9 I =294) PRINT 510 (TR(I)01:294)9(DR(I)9I=204) PRINT 500 Is PKCI)9 Io TK(I)9 Io PC(I)o Io TCII) 104 I = 2 CALL PRCALC TC(I) I TK(I)/TR(I) PC(I) I PK(I)/PR(I) C THIS PROGRAM USES DR(3) T0 CONVERSE G£DER(2) DELTADR = 0.02 GZDERI = 0 DO 8 LIlo4O I 8 I + 1 URI!) 8 DR(I) + DELTADR CALL PRCALC TC(I) I TK(I)/TR(I) PCII) I PK(I)/PR(I) I 8 I-I DELTAX I X(I+I) - XII) TCIDER I (TC(I+I) - TC(I-I))/(2.*DELTAX) TCZDER = (TC(I+I) - 2.*TC(I) + TC(I-1))/DELTAX**Z TRIDER I ~TRII)*TCIDER/TC(I) TRzDER I -2.*TRIDER*TCIDER/TC(I) ~ TRII)*TCZDER/TC(I) PCIDER I (PCII+I)-PC(I~I))/(2.*DELTAX) PCZDER I (PC(I+1) - 2o*PC(I) + PC(I-I))/DELTAX**2 PRIDER I -PR(I)*PCIDER/PC(I) PRzDER I -2.*PRIDER*PCIDER/PC(I) - PRII)*PCZDER/PC(I) CALL GZCALC IF (ABSF(DELTADR) - TOLI90907 7 DELTADR I -GZDER*DELTADR/(GZDER - GZDERI) 8 GZDERI 3 GZDER STOP 0007 9 PRINT 420 TC‘I)! TCIDERQ TCZDER. XII).PC(I)1 PCIDEROPCZDERQ XII)! I TR‘I). TRIOER. TRZDER. X(I)o PR(I)9 PRIDERO PREDERQ XII) PRINT 430XIIIOTR(I)9PR(I)ODR(I)QZR(I)0PZRWRPR(I)9XII+I)OTR(I+l)o 1 PR(I+IIODR(I+I)OZR(I+I)oPZRWRPRII+l) I = I + I 29 NUMBER 3 0 PRINT 469 X(I) (3 TO STATEMENT 31-12 SETS UP 3 SETS OF GUESSES (TRODR) AND C CALCULATES GZDER AND G3DER DELTATII) ‘ 0.1 DELTATI2) = O DELTADII) I O DELTAD(2) = 0.1 M = 1 GO TO 10 28 M = M + 1 TR(I+I) = TR(I+1) + DELTATIM-l) lOO DR(I+I) I DR(I+I) + DELTAD(M-I) 10 I I I + I CALL PRCALC PC(I) I PK(I)/PR(II TC(I) I TK(I)/TR(I) DO 30 KIIO4 Z(K) I X(I“4+K) 30 FUNC(K) I TC(I-4+K) CALL POLY (TCIDERQTCZDERoTC3DER) DO 31 K 2104 31 FUNC(K) I PC(I‘4+K) CALL POLY (PCIDERoPCZDERvPCBDER) I I I - 1 TRIDER : -TR(I)*TCIDER/TC(I) TRZDER = ~2o*TRlDER*TClDER/TC(I) - TR(I)*TCZDER/TC(I) TR3DER : 6o*TK(I)*TCIDER*TCZDER/TC(I)**3-6.*TK(I)*TCIDER**3/TC(I) I **4 - TK(I)*TCBDER/TC(I)**2 PRIDER -PR(I)*PCIDER/PC(I) PR2DER = -2.*PRIDER*PCIDER/PC(I) - PRII)*PC2DER/PC(I) PR3DER = 6.*PK(I)*PCIDER*PCZDER/PC(I)**3-6.*PK(I)*PCIDER**3/PC(I) I **4 - PKII)*PC3DER/PC(I)**2 CALL GBCALC G(M) I GEDER G(M+3) = GSDER IF'IM-3)28c32035 32 IF (G(3)**2 + G(6)**2 - TOL)15qISoI4 THIS CONVERGENCE PROCEDURE USES A LINEAR INTERPOLATION 14 CALL CONVERG (DELTAT(I)0 DELTAT(2)0 DELTAT(3)9 DELTAD(1)9 DELTAD 1 (2)0 DELTAD(3)9 G(l)o 6(2). 6(3). 6(4). 6(5). G(6)) I7 IF(DR(I+I) + DELTAD(3))180180750 750 IFIDRII+1)+DELTAD(3)-2.)I9018v18 18 DELTAD(3) = DELTAD(3)/Zo GO TO 17 I9 IF(TR(I+1) + DELTAT(3))ZOo20v751 751 IF(TR(I+1)+DELTAT(3)‘4.)21.20.20 20 DELTAT(3) = DELTATI3)/Zo GO TO 19 21 DO 23 LIIQZ DELTAT(L) = DELTAT(L+1) DELTAD(L) = DELTAD(L+1) G(L) I G(L+I) 23 G(L+3) = G(L+4) PRINT l6. TR(I+I)0 DR(I+1)9 GZDER. GBDER TR(I+I) I TR(I+1) + DELTAT(3) DR(I+I) I DR(I+1) + DELTAD(3) NUMBER I NUMBER + 1 IF (NUMBER-30)10910024 15 PRINT 42.TC(I)¢ TCIDER. TCZDER. TC30ER9PC(I)9 pCIDEROPCZDERO I PCBDERO TR(I)9 TRIDER. TRZPERo TRBDERo PR(I)0 FRIDERO PRZDERc 2 PRBDER NNN = 1+1 PRINT 439 (X(J)0TR(J);PR(J)oDR(J)oZRIJ)oPZRURPR(J)9 J-IONNN) PRINT 44o GZDERQ GSDER PRINT 409 NUMBER lOl 13I+I DRII+1) = 1.1*DR(I) TR(I+1) = 1.1*TR(I) [P(I-NUM)29Q35035 24 PRINT ‘5. NUMBER 25 PRINT 41 35 STOP 5555 END 102 Program LINEAR-Main Bogl Purpose: Input: Details: To generate values of Tr and Dr at x (1 equal to or greater 1+1 than 2) such that 61'“ and G1°°' equals zero. P-V-T data on a reference compound (in this case, prOpane). This data includes the critical properties, the gas constant R, eight Benedict-Webb-Rubin constants, and five constants for computing Co as a function of temperature° Also input are the first 2 and x30 Using the values of Tr and Dr at x1, x guesses of Tr and Dr at x 2, and x3 (Tr I l and Dr I l at x1 I 0), this program changes Dr3 until the second derivative of free energy, G 9°, is equal to zero at x For 2 2° x1+1 (1 equal to or greater than 3) the guesses of Tr1+1 and Dr1+1 are adjusted until Gi process is repeated until the entire composition range has been 9' and Gi'°' equal zero. This latter coveredo A linear convergence scheme is employed to adjust Dr3 until G2" equals zeroo Two guesses are made for Dr then a new guess is calculated from: 3, and (new Dr - present Dr3)/(present Dr - previous Dr3) I(desired 3 62°" - present G 3 '°)/(present G 0° - previous G2") 2 Since desired 62°“ I O, 2 new Dr I pres Dr + [-pres G2°°/(pres G “9 - prev 62")] 3 3 2 (pres Dr3 - prev Dr3)o (1) Calculation of 62"“ requires first and second derivatives of Tr and Pr at :20 These derivatives are computed by fitting a quadratic equation through the values of To and Fe at :1, :2, and x and then computing Tc ' and To "' from: 3’ 2 2 103 9 g - Tc2 (Tc3 Tcl)/(2Ax) (2) 2 00 a - Tc2 (Tc3 2Tc2 + Tcl)/(Ax) (3) Similar equations may be written for Pea9 and Pc2"o First and second derivatives of Tr and Pr are computed from: Tr" = [a(T/Tc)/3x] = -T°Tc°/Tc2 a -Tr-Tc'/Tc (u) T,P Tr°' I [32(T/Tc)/3x2] I —2Tr°°Tc°/Tc - TroTc'°/Tc (5) T,P Values of Tr and Dr at x1+1 (1 equal to or greater than 3) are computed using a convergence process which assumes G" and G"'9 to be linear functions of Tr and Dre The values of Tr and Dr at x1, x , and x must be known because third derivatives i-l i-2 are needed to calculate Gi°"°o Three sets of guesses for (T D ) are needed before the linear interpolation r1+1 method can be appliedo Subroutine CONVERG is called to compute r1+1’ the next increment of Tr“l and Dr1+1 from the previous three sets of guesses for Tr1+1 and Dr1+1, and the values of 61" and 61"“ calculated from those sets of guesses. Whenever the double linear convergence process is used, all derivatives of To and Fe are computed using subroutine POLYs This subroutine fits a third degree polynomial to the values of Tc and Po at x , and x Subscript i is necessari- 1+1' ‘1' x1.1 1-2° 1y at least 3. The value of Tri000 is calculated from: Tr'V' . [33(T/Tc)/3x3]T9P= 6T°Tc°°Tc°°/Tc3 - 6T(Tc')3/Tch - T°Tc°"°/Tc2 (6) and similarly for Pr’Vgo The only restrictions imposed on the calculated values of Tr and Dr are that they must remain positive since negative values do not have physical significanceo This program also Method of Checking: Fortran Nomen- clature: 10h limits the number of iterations to thirty before the convergence process is interruped because.of lack of convergenceo Since the primary function of the main body of this program is to provide convergence, the program must work prOperly if such convergence is obtainedo AC, BC, CC, DC, EC I five constants to compute Co as a function of temperature A0, B0, CO, A, B, C, ALPHA, GAMMA I the eight Benedict-Webb- Rubin constants Bo, A 59 b0 cs “9T C o’ o’ 0000 \anH llllllllllllllll 0000 \00340\ f 8 h J k l m n T DCB, CB, PCB, I critical density, temperature and pressure of propane DELTAD I variable used to represent differences in Dr DELTADR I difference between the most recent value of Dr DELTAT I variable used to represent differences in Tr DELTAX I increment of x after x DR I reduced density (Dr) 3 NUM I highest value of the subscript on x NUMBER I a variable used to count iterations NUMDATA I a variable used to represent the total number of sets 3 of data PCI pseudocritical pressure (Pc) PClDER I Pc“ PCQDER I PcH PC3DER I Pc°°° PK I critical pressuve of mixture (Pk) PR I reduced pressure (Pr) PRIDER I Pr' PRZDER I Pr°° PR3DER I Pr°°° R I gas constant (literoatmo)/(8-moleo°K) TC I pseudocritical temperature (Tc) TClDER I ch TCZDER I Tc°° TC3DER I TC°" TX I critical temperature of mixture (Tk) TR I reduced temperature (Tr) 105 TRlDER I Tr0 TRZDER I Tr°° TR3DER I Tr'M X I mole fraction ethane ZR I reduced compressibility factor (Zr) l SUBROUTINE CONVERG (DELTAXle DELTAXZQ DELTAXQ DELTAYXQ DELTAYO DENOMI COEFFI COEFFZ COEFF3 COEFF4 DENOMZ DELTAX DELTAY END Illflflflfl 106 F19 F20 F39 Cl. 629 G3) DELTAX1*DELTAY2 - DELTAX2*DELTAY1 ((FZ”F1)*DELTAY2 I (F3-F2)*DELTAY1)/DENOM1 ((F3‘F2)*DELTAX1 ‘ (F2-Fl)*DELTAX2)/DENOM1 ((GZ‘GI)*DELTAY2 ‘ (GB‘GZ)*DELTAY1)/DENOM1 ((63-GZ)*DELTAX1 ‘ (GZ‘GI)*DELTAX2)/DENOM1 COEFF1*COEFF4 ’ COEFF2*COEFF3 (53*COEFF2 ‘ F3*COEFF4)/DENOM2 (F3*COEFF3 “ G3*COEFFl)/DENOM2 DELTAYZQ 107 Subroutine CONVERG Purpose: Input: Details: To compute the estimated increment for each of two variables which will make two functions of those two variables equal zero. The differences between the previous three guesses of both variables and the previous three values of both functions computed from those variables. The equations for a double linear convergence process are developed for two general functions, F(X,Y) and G(X,Y). This procedure is called a double linear convergence process because two functions, F(X,Y) and G(X,Y) are assumed to be linear functions of the two variables, X and Y. Required initially are three different sets of guesses for X and Y [(X1,Y1), (X2,Y2), (X3.Y3)] and the values of F and G for each set of X and Y. Values of X1‘ and Yh are de- sired such that Fh I O, Gh I 0. Letting AX1 I X2 - X1, 2 3 2 and six unknowns (a, b, c, d, Xh, Y“) can be written: aAX + bAY I AF 1 1 l 8AX2 + bAY2 I AF2 an3 + bAY3 AFB - -F3 \ cAXl + dAYl I AGl cAX2 + dAY2 I AG2 cAX3 + dAY3 I AG3 I -G3 AX I X - X , etc., six equations (1) (2) (3) (h) (5) (6) Equations (1) and (2) are solved for a and b, equations (1:) and (5) for c and d, and equations (3) and (6) for AX3 and AY3. The expressions for AX3 and AY3 are: Method of Checking: Fortran Nomen- clature: 108 AX3 a x,4 - x . (G3b - F3d)/(ad - bc) (7) 3 AY3 I Yh - Y3 The last two mathematical statements of the subroutine are I (F3c - 63a)/(ad - bc) (8) equivalent to equations.(7) and (8). This subroutine was checked by using it on a systcm of two linear functions of two unknowns. Three sets of guesses were made for the two unknowns, and the subroutine was employed to compute the increments necessary to make the functions equal to zero. Only one step was needed which indicated that the subroutine was extrapolating properly. COEFFl I coefficient a in equations (1) to (3) COEFF2 I coefficient b in equations (1) to (3) COEFF3 I coefficient c in equations (h) to (6) COEFFh I coefficient d in equations (h) to (6) DELTAXl I X2 - X1 DELTAX2 I X3 - X2 DELTAX I Xh - X3 I difference between the next value of X and the most recent value of X DELTAYl I Y2 - Y1 DELTAYZ I Y3 - Y2 DELTAY I Yh - Y3 I difference between the next value of Y and the most recent value of Y DENOMl - (Y3 - Y2)(X2 - x1) - (Y2 - Y1)(X3 - x2) DENOM2 I ad - be bUBROUTINE POLY 109 DIMENSION X(25)o W(25)9 PC‘ZS). 2(5)e l TR(25)¢ DEN(1) DEN(2) DEN(3) DEN(4) S I 2(3) FSTDER = (fibUNI‘ SECDER um.- +++u THRDER I END DR(25)0 ZR(25)9 (FSTDERoSECDERqTHRDER) TC(25)§ PK(25)9 TK(25)e PR(25)0 FUNC(5)O DEN(5)O PZRWRPR(25) COMMON IeDCBeTCBeACeBCeCCoDCoECOCOCI9C20C3QC49C59C61C79C89C990Rc 1 DRIDERODRZDERsTReTRlDERvTRaDEReTRBDERQPRePRlDERePRZDERopflaDERs 2 PCePKeTCeTKoQQeTOeTOlDERsTOZDEReTNeTNlDERQTNZDEReTTsTTIDER9V0 3 VIDERQY'YIDEROYZDERQZKOZKIDERQZKzoEROZROZ!XOWOGZDEROGBDEROFUNC 4 OPZRWRPR (Z(l)“Z(2))*(Z(1)‘Z(3))*(Z(l)‘Z(4)) (2(2)‘Z(1))*(Z(2)'Z(3))*(Z(2)‘Z(4)) (Z(3)‘Z(1))*(Z(3)‘Z(2))*(Z(3)‘Z(4)) (Z(4)‘Z(1))*(Z(4)‘Z(2))*(Z(4)’Z(3)) ((30*S**2-20*S*(Z(2)+Z(3) 2(4))+Z(3)*(Z(2)+Z(4))+Z(2)* (6.0*S -2o*(Z(4) (600*5 -2.*(Z(l) (6.0*$ -2.*(Z(1) (6eO*S -2.*(Z(1) + + + + 2(2) 2(3) 2(2) 2(2) + + + + 2(4))*FUNC(1))/DEN(1’+((3e*5**2“20*$*(2(1)+Z‘3)+Z(4))+Z(l)* (Z(3)+Z(4))+Z(3)*Z(4))*FUNC(2))/DEN(2)+((3.*S**2-20*$*(Z(1)+ Z(2)+Z(4))+Z(1)*(Z(2)+Z(4))+Z(2)*Z(4))*FUNC(3))/DEN(3)+((3.*S**2 -2.*S*(Z(1)+Z(2)+Z(3))+Z(1)*(Z(2)+Z(3))+Z(2)*Z(3))*FUNC(4)) /DEN(4) Z(3)!)*FUNC(1)/DEN(1) Z(4)))*FUNC(2)/DEN(2) Z(4)))*FUNC(3)/DEN(3) Z(3)))*FUNC(4)/DEN(4) 6o*(FUNC(1)/DEN(1)) + 6o*(FUNC(2)/DEN(2)) 1 + 60*(FUNC(3)/DEN(3)) + 6o*(FUNC(4)/DEN(4)) llO Subroutine POLY Purpose: To calculate numerically first, second and third derivatives of a function at a given value of the independent variable. Input: Four values of the function, f(Zl), f(22), f(23), and f(Zh), and four values of the independent variable, 21, 22. 23, and 2h are placed in memory prior to the calling of the subroutine. Details: A third degree polynomial (four points) may be writtenhl: r(z) — a23 + bz2 + cZ + d (1) or (z-zz)(z-z )(z-zh) (z-z1)(z-z )(zazh) “2’ ' T7477 “21’ * r—T'Lr—J “22’ * (Z-Zl)(Z-22)(Z-Zh) r(z ) + (Z-Zl)(Z-22)(Z-Z ) {(2 ) W> a W. m ..-z. .. 2 (2) 32 - 2(22+z +zh)z + z z + zazh + 2 Zn 21-22 Zl-Z3 ZlIZh 1 + see (3) f"(Z) I W f(Z1)+ u. (h) 3 :vv 3 6 I f (Z) 133:3-7T§::§;Tr2;:2:7‘f(Zl) + 00° ' (5) Equations (3). (h), and (5) can be used to calculate numerically first, second, and third derivatives, respectively, at any value Z. 'The value of Z is called S in the Fortran program. In this particular program, the derivatives are computed at 23. Method of Checking: A third degree polynomial was constructed and values of {(z) cal— culated at four arbitrary values of Z. These results were then Fortran Nomen- clature: 111 input into the subroutine to see if the generated values agreed with the original polynomial and its derivatives, for specific values of Z. DEN I the.denominators_of the terms in equations.(2), (3), (h), and (5). DEN(1) represents the denominator of the first term,.etc. FSTDER I calculated value of the first derivative FUNC I value of the function . S I the value of 2 at which the derivatives are being calculated SECDER I calculated value of the second derivative THRDER I calculated value of the third derivative Z I specific values of the independent variable (2) for which the values of the function are known 112 SUBROUTINE WBKAY DIMENSION X(25)v W(25)9 PC‘25)9 TC‘25)9 PK(25)0 TK(25)0 PR‘ZS). 1 TR(25)9 DR‘ZS)! 29(25)e 2(5)s FUNC(5)QPZRWRPR¢25) COMMON IODCBOTCBsACeBCsCCeDCsEGeCeCl9C29C39C49C59C69C79C80C990Rs l DRIDERODRZDEReTReTRlDEReTRZDERoTRBDERepRePRlDERQPRZDERQPR3DERQ 2 PCQPKeTCeTKQQQeTOQTOlDERsTOZDEReTNeTNIDEReTNaDERsTTeTTlDEReVe 3 VIDERQYeYlDERsYEDERoZKoZKlDEReZKZDEReZRsZeXoWoGZDEReGBDERoFUNC 4 ePZRWRPR W‘I) = 30e06*X(I)/(30e06*X(I)+(1.‘X(I))*100e17) TKCI) I (51303-18604392126*X(I)+20 e413825*X(I)**2-558e3704113* 1 X(I)**3+12101957936*X(I)**4+459e67)/1e8 PK(I) I (396e+3051e927387*W(I)-1199e746006*W(I)**2-4263e827785* 1 W(I)**3+2727.646362*W(I)**4)/14e696 END 113 Subroutine WBKAY Purpose: To calculate Tk and Pk at a given composition. Input: The value of composition (x1) and the value of subscript i must be placed in memory before this subroutine is called. Details: The experimental Tk and Pk values for the ethane-n-heptane system, as published by W. B. Kay,23 were fitted to a fourth degree polynomial by a modification of the method of least squares. This modification was necessary in order to obtain an exact fit for pure ethane and pure n-heptane. The critical temperature envelOpe curve is expressed as a function of mole fraction ethane. However, the critical pres- sure envelope curve is expressed as a function of weight frac- tion ethane because a better fit of the pressure data was ob- tained using this variable. Method of , Checking: The compositions as recorded in the published article were input, and the results were compared with the published values at these compositions. Fortran Nomen- clature: I I subscript i on composition (x) PK I true critical pressure of mixture (Pk), atm. TX I true critical temperature of mixture (Tk), °K W I weight fraction ethane X I mole fraction ethane (x) 11h SUBROUTINE PRCALC DIMENSION X(25)9 W‘25)! PC‘25)0 TC‘ZS)! PK‘ZS). TK‘ZS). PR‘25)9 1 TR(25)9 DR(25)e ZR‘25)0 2(5)! FUNC‘S)OPZRWRPR(25)QSAVE(4) COMMON IQDCBOTCBeACeBCeCCeOCoECeCQCX9C2eC3sC4eC59C69C7eC89C9eDRe DRIDERQDRZDEReTReTRlDEReTRBDERsTR3DERePRePRIDERePRZDEROPR3DER9 PCQPKQTCeTKeQQeTOeTOlDERoTOZDEReTNeTNIDEReTNzDERsTTeTTIDERQVQ VIDEReYeYlDEReYZDEReZKeZKlDERsZKZDERsZReZonWeGZDEReGSDEReFUNC QPZRWRPR Q I (TR(I)*TCB)/IOOe «bUNH ox = (10.*Q**6)**4 co = (EC/(Q**2*Ql) + DC/G/Ol + BC/Ql + AC)/(CC/Ol + 1.) C4 = (cos/C)*c0)/oca DO 5 KIloZ DR(I) = loE-05*(-le)**K + DR(I) 00 = EXPF(-C9*DR(I)**2) PR(I) I CI*DR(I)*TR(I)+C2*TR(I)*DR(I)**2-C3*DR(I)**2-C4*DR(I)**2 1 /TR(I)**2+C5*TR(I)*DR(I)**3-C6*DR(I)**3+C7*DR(I)**6+(le+C9*DR(I) 2 **2)*QQ*C8*DR(I)**3/TR(I)**2 29(1) I PR(I)/TR(I)/DR(I) SAVE(K) = PR(I) 5 SAVE(K+2) I ZR(I) PZRWRPR(I) = (SAVE(4)-SAVE(3))/(SAVE(2)“SAVE(1)) END 115 Subroutine PRCALC Purpose: Input: Details: Method of To calculate Pr and Zr at a given composition. Values for constants f, g, h, J, k, l, m, n, c, aCO, bCO, cCo, dCo, eCO, critical temperature (TCB), critical density (DOB), Tri, Dri, and the value of the subscript i must be placed in memory before this subroutine is called. The expression for the Benedict-Webb-Rubin constant Co as a function of temperature is: co - [ICC/Qz-Ql) + dC°/(Q-Ql) + bCo/Ql + acol/(cco/Ql + 1) (1) where Q I Tr-TCB/lOO Q1 I (10°Q6)h From the value of Co the constant 1 is computed from: 1 - m-Colc/DCB (2) For use in the subroutine and subroutines GZCALC and G3CALC, a quantity called QQ is defined: QQ I exP(-n-Dr12) (3) Pri is calculated according to equation (3), Appendix III. Zr1 is computed from: Zr1 I Pr1/(Tr1 Dri) f (h) As added information, the slope of the reduced compressibi- lity factor curve as a function of reduced pressure at constant temperature is computed numerically. This quantity is called PZRWRPR. Checking: This expression for reduced pressure was used to calculate fugacities. These calculated fugacities checked very well 116 with known fugacities. Fortran Nomen- clature: AC I aCO BC I b0 0 O r:- I Benedict-Webb-Rubin reduced equation of state constant 1 DC I dCo EC I eCo I I subscript i PZRWRPR - (azr/aPr)Tr 117 SUBROUTINE GaCALC DIMENSION X(25)9 W‘25)9 PC‘ZSII TC‘25). PK(25)0 TK‘25)9 PR‘ZS)O l TR(25)0 09(25)e 29(25)e 2(5)! FUNC(5)OPZRWRPR(25) COMMON IQDCBOTCBQACQBCQCCsDCeECszC)9C29C3sC40C5eC6eC7eCBeC990Re 1 DRIOERQDRZDER.TReTRlDEReTREDEReTRBDERoPRePRIDERoPRZOERepR3DERe 2 PC.PK.TC.TK.OO.TO.T01DER.TOEDER.TN.TNIDER.TNaoERoTToTTIDER.V. 3 VIDERQYQYIDERoYZDEReZKeZKlDEReZKZOERe29eZeXeWeroeReG3DERQFUNC 4 ePZRWRPR CALL PRCALC TO I le/TR‘I) TOIDER : —TO**2*TRIDER TOZDER I (2e*TRIDER**2/TR(I) - TRZDER)/TR(I)**2 TN I TO**3 TNIDER I 3e*TO**2*TOlDER TNZDER I 6e*TOIDER**2*TO + 3e*TO**2*TOZDER TT I TRIDER*TN TTIDER I TREDER*TN + TRIDER*TNIDER U = PRIDER-(C1*DR(I)+C2*DR(I)**2+20*C4*DR(I)**2/TR(I)**3)*TRIDER 1 + (-C5*DR(I)**3+2e*C8*DR(I)**3*(le+C9*OR(I)**2)*OQ/TR(I)**3)* 2 TRIDER V = C1*TR(I)+2e*C2*DR(I)*TR(i)-2e*C3*DR(l)‘2e*C4*DR(I)/TR(I)**2 l +3.*C5*TR(I)*DR(I)**2-3e*C6*OR(I)**2+6e*C7*DR(I)**5+(C8/TR(I)**2) 2 *00*(3e*DR(I)**2+3e*C9*DR(I)**4-2e*C9**2*DR(I)**6) DRIDER I U/V UIDER I PRZOER‘C1*OR(I)*TR2059-C1*DR1DER*TRIDER-C2*DR(I)**2*TRZDER 1 -2.*C2*DR(I)*DRIDER*TRIDER-2.*C4*DR(I)**2*TTIDER-4.*C4*DR(I)* 2 ORIDER*TT-C5*DR(I)**3*TRZDER-3e*C5*DR(I)**2*DRIDER*TRIDER 3 +20*C8*DR(I)**3*(10+C9*DR(I)**2)*QQ*TT1DER+20*C8*TT*QQ* 4 (3e+3e*C9*OR(I)**2‘2e*C9**2*DR(I)**4)*OR(I)**2*DRIDER VIDER I C1*TRIOER+2e*C2*DRlDER*TR(I)+2e*C2*DR(I)*TRIDER-2s*C3* 1 DRIDER-2e*C4*DRlDER/TR(I)**2+4.*C4*DR(I)*TT+6e*C5*TR(I)*DR(I)* 2 DRIDER+3e*C5*TRIDER*DR(I)**2'6e*C6*DR(I)*DRIOER+30e*C7*DR(1)**4* 3DRIDER+C8*OO*(6e*DR(I)*DRIOER/TR(I)**2+6e*C9*DR(I)**3*DRIDER/TR(I) 4 **2’IBeO*C9**2*DR(I)**5*DRIDER/TR(I)**2+4e*C9**3*DR(I)**7*DRIOER 5 /TR(I)**2-6e*DR(I)**2*TT-6e*C9*DR(I)**4*TT+4e*C9**2*DR(I)**6*TT) DREDER I (UIDER - VIDER*DRIDER)/V Y = PR(I)/(CI*TR(I)) ZK I Y/DR(I) YIOER I (PRIDER - Y*C1*TRIDER)/(CI*TR(X)) Y2DER m (PREDER - Y*C1*TRZDER - 2e*YlDER*Cl*TRlDER)/(Cl*TR(1)) ZKIOER I (YIOER - ZK*DRIDER)/DR(I) ZKZDER I (YZOER-ZK*DRZDER-2e*ZK1DER*DR1DER)/OR(I) GEDER I leO/(X(I)‘X(I)**2)+ZKZDER-ZKZDER/ZK+(ZKIDER/ZK)**2 + (C2*DRZDERIC3*(TO*DPEDER+TOZDER*DR(I)+2e*TOIOER*DRIDER)‘C4* (TN*DRZDER+TNZDER*DR(I)+2e*TNlDER*DRlDER)+(2e*C5/2e)*(DR(1)* DRZDER+DRIDER**2)‘(C6/Ze)*(2e*TO*DR(I)*DRZDER+TOZDER*DR(I)**2+ 4e*TOIDER*DR(I)*DRIOER+2e*TO*DRlDER**2)+C7*(2e*TOIDER*DR(l)**4* DRIDER+TOZDER*DR(I)**5/5.+4.*TO*DR(I)**3*DRIDER**2+TO*DR(1)**4* DRZDER)+C8/(2.*C9)*(2.*TN20ER+QG*(-2.*TN20ER-C9*DR(I)**2*TNZDER+ 4e*C9*DR(I)*DRIDER*TNIDER+2e*C9*DR(I)*DREDER*TN+2e*C9*ORIDER**2* TN+4.*C9**2*DR(I)**3*DR10€R*TN10ER+2.*C9**2*DR(I)**3*DnzoER*TN+2. *C9**2*OR(I)**2*DRIDER**2*TN“4e*C9**3*OR(I)**4*DRIDER**2*TN)))/Cl END \OCDIJO‘UIbUNH 118 Subroutine 62CALC Purpose: Input: Details: Method of Checking: Fortran Nomen- clature: To calculate G". Derivatives Tr', Tr", Pr', and Pr" at x must be in memory 1 when GZCALC is called. In addition, the input requirements for subroutine PRCALC must be satisfied because PRCALC is called in this subroutine. Subroutine PRCALC is called to compute the constant 1 and the quantity QQ. The derivations of all equations used in this subroutine are in Appendices III, IV, and V. The equations in the order used are: III 9, 10, 11 III - 12, 13, 1h Iv-6,7 IV - 1, 2, 3 IV - 9, 10 :40.»er NU) III - Program CHECK was used to check the equations used in this subroutine. DR I reduced density (Dr) DRlDER I Dr' DRZDER I Dr" GZDER I G" PR I reduced pressure (Pr) PRIDER I Pr' PR2DER I Pr" TN I Tn TNlDER I Tn. TN2DER I Tn" 119 TO I To T01DER I To. TOZDER I To" TR I reduced temperature (Tr) TRlDER I Tr' TRZDER I Tr" TT I Tt’ TTlDER I Tt' U I U UlDER I U' V I V VlDER I V' X I.mole fraction ethane (x) Y I Y YlDER I Y': Y2DER I Y" ZK I compressibility factor (Z) ZKlDER I Z' ZKZDER I Z" 120 SUBROUTINE G3CALC DIMENSION X(25)e W(25)9 PC‘25). TC(25)9 PK(25)O TK(25)0 PR(25)O 1 TR(25)9 09(25)e ZR(25)1 2(5). FUNC(5)OPZRWRPR(25) COMMON IoDCBoTCBeACoBCoCCoOCoECQCoCI9C29C39C49C59C60C71C80C9QDRO 1 DRIDERQDR2DERsTReTRIDERsTR£DERQTR3DERQPRQPRIDERQPRZDEROPR3DERO 2 PCQPKOTCeTKeQQeTOsTOlDERQTOZDEReTNeTNIDEReTN2DEReTTeTTlDERsVe 3 VlDEReYeYlDEROY2DEReZKQZKIDEROZK2DEReZReZeXeWeG2DEReG3DEROFUNC 4 sPZRWRPR CALL G2CALC TO3DER =-TO*(3e*(TOIDER*TR2DER+TO2DER*TR1DER)+TO*TRBDER) TNBDER I 3e*(2e*TOlDER**3+6e*TO*TOIDER*TO2DER+TO3DER*TO**2) TT2DER I TRBDER*TN+TRIDER*TN20ER+2e*TRZDER*TN1DER UZDER I PRBDER-Cl*(DR(I)*TRBDER+2e*DRIDER*TR2DER+DRZDER*TRIDER) 1 ‘C2*(DR(I)**2*TR3DER+4e*DR(I)*DRIDER*TRZDER+2e*DR(I)*DRZDER* 2 TRIDER+20*DR1DER**2*TRIDER)’C4*(20*DR(I)**2*TTZDER+80*DR(I)* 3 DRIDER*TTIDER+4e*DR(I)*DR2DER*TT+4e*DRIDER**2*TT)”C5*(DR(I)**3 4 *TR3DER+6e*DR(I)**2*DRIDER*TRZDER 3e*DR(I)**2*DRZDER*TRIDER+ 5 6e*DR(I)*DRIDER**2*TRIDER)+2e*C8*QQ*(DR(I)**3*TT2DER+C9*DR(I)**5 6 *TTZDER+6e*DR(I)**2*DRIDER*TTIDER 6e*DR(I)**4*C9*DRIDER*TTIDER 7-4e*C9**2*DR(I)**6*DRIDER*TTIDER+3e*DR(I)**2*DR2DER*TT+3.*C9*DR(I) 8 **4*DRZDER*TT"2e*C9**2*DR(I)**6*DR2DER*TT+6e*C9*DR(I)**3*DRIDER** 9 2*TT+6.*DR(I)*DRlDER**2*TT-18.*C9**2*DR(l)**5*DRIDER**2*TT) U2DER I UZDER+2e*C8*QQ*4e*C9**3*DR(I)**7*DRIDER**2*TT V2059 I C1*TR2DER+C2*2e*DR2DER*TR(I)+4e*C2*DRlDER*TRIDER + 2e*C2*DR(I)*TREDER-2e*C3*DRZDER+C4*(-2e*DRZDER/TR(I)**2+8s*DRIDER *TT+4e*DR(I)*TTIDER)+C5*(12e*TRIDER*DR(I)*DRIDER+60*TR(I)*DRIDER **2+6e*DR(I)*DRZDER*TR(I)+3e*TR2DER*DR(I)**2)‘6e*C6*(DRIDER**2+ DR(I)*DRZDER)+30e*C7*(DR(I)**4*DRZDER+4e*DR(I)**3*DRIDER**2) +2s*C8*QQ*((3e*DR(I)*DRZOER+3.*C9*DR(I)**3*DREDER~9e*C9**2*DR(I) **5*DR2DER+2.*C9**3*DR(I)**7*DRZDER)/TR(I)**2+(3e*DR)DER**2+ 3e*C9*DR(I)**2*DRIDEQ**2'51e*C9**2*DR(I)**4*DRIDER**2+32e*C9**3* DR(I)**6*DRIDER**2‘4e*C9**4*DR(I)**8*DRIDER**2)/TR(I)**2-30*DR(I) **2*TTIDER-3e*C9*DR(I)**4*TTIDER+2e*C9**2*DR(I)**6*TTIDER) V2DER : V2DER + 2e*C8*QQ*('12e*DR(I)*DRIDER*TT-12e*C9*DR(I)**3* 1 DRIDER*TT+36e*C9**2*DR(I)**5*DRlDER*TT“8e*C9**3*DR(I)**7*DRIDER 2 *TT) DRBDER I (UZDER-(2e*VlDER*DR2DER+V2DER*DR1DER))/V Y3DER I (pRBDER‘Y*CI*TR3DER-3e*YIDER*Cl*TRZDER-3e*YZDER*C1*TRIDER) I /(C1*TR(I)) ZK3DER I (Y3DER“ZK*DRBDER-3e*ZK1DER*DRZDER“3e*ZKZDER*DR1DER)/DR(I) II G3DER = (2e*X(I)‘Ie)/((X(I)“X(I)**2)**2)+ZK3DER-ZKBDER/ZK+3e* ZK2DER*ZKIDER/ZK**2‘20*(ZKIDER/ZK)**3+(C2*DRBDER-C3*(TO*DRBDER+ TOBDER*DR(I)+3e*DQIDER*TO2DER+3e*DRZDER*TOlDER)“C4*(TN*DR3059+ TN3DER*DR(I)+3e*DRlDER*TNZDER+3e*DRZDER*TNIDER)+C5*(DR(I)*DRSDER +3e*DRIDER*DRZDER)‘(C6/2e)*(2e*TO*DR())*DR3DER+TO3DER*DR(I)**2) -(C6/2e)*(6e*TO*DRIDER*DR2DER+6e*TOIDER*DR(I)*DR2DER+6e*TO2DER* DR(I)*DRIDER+60*TOIDER*DRIDER**2) C7*(3e*TOZDER*DR(I)**4*DRIDER+ 12e*TOlDER*DR(I)**3*DR10ER**2+3.*TOIDER*DR(I)**4*DR2DER+12e*TO *DR(I)**2*DRIDER**3+I2O*TO*DR(I)**3*DRIDER*DR2DER+TO*DR(I)**4* DR3DER+TO3DER*DR(I)**5/50)+C8/C9*TN3DER)/Cl G3DER I G3DER+((C8/(2e*C9))*QQ*('20*TNBDER-C9*DR(I)**2*TN3DER+6e* 1 C9*DR(I)*DRIDER*TN2DER+6.*C9*DR(I)*DR2DER*TN1DER+6e*C9*DRIDER**2* 2 TNIDER+2e*C9*DR(I)*DRSDER*TN+6e*C9*DRIDER*DRZDER*TN+6.*C9**2* 3 DR(I)**3*DR1DER*TNZDER+6e*C9**2*DR(I)**3*DR2DER*TNIDER+6e*C9**2* ‘OCDIJO‘UIbUN-I ‘OCDIJO‘UIbUNt-I 40‘UIA 121 DR(I)**2*DR1DER**2*TN10ER+2.*cg**2*DR(1)**3*DRgoER*TN+5,*c9**2* DR(I)**2*DR10ER*DR20ER*TNI12,*c9**3*on(1)**4*DR1DER**2*TNIDER-12, *C9**3*DR(I)**4*DRIDEQ*DR20ER*TN-20.*C9**3*DR(I)**3*DRIDER**3*TN+ 8.*C9**4*DR(I)**5*DRIDER**3*TN))IC1 END 122 Subroutine G§CALC Purpose: Input: Details: Method of Checking: Fortran Nomen- clature: To calculate G"'. Tr"' and Pr"' at x1. Also, the input requirements of sub- routine G2CALC must be satisfied because GZCALC is called in this subroutine. Subroutine GZCALC is called to make preliminary calculations. The derivations of the equations used in this subroutine are found in Appendices III, IV, and V. The equations in the order used are: III III IV - 8 IV IV - 5 V - h V - 8 III 19 I I P'F' \n04 I p.) H. O ...n (0 Program CHECK was used to check the equations used in this subroutine. DRBDER I Dr"' GBDER I G"' PRBDER I Pr"' TN3DER I Tn"' TO3DER I To... TR3DER I Tr"' TT2DER I Tt" U2DER I U" v V2DER I V" ZKBDER . 2'” 123 DATA 50153462298 37000505321 42e38841074 008207 e097313 6e87225 508256e .0225 e9477 129000e e000607l75 0022 508256.1345 le388457033E-06 2e8445691415‘12 .0 e0 1 1.005 10010 1e017 1e002 1.004 1.007 12h Data The data represents the following: 5.153h62298 I the critical density of propane as computed from the Benedict-Webb-Rubin equation of state, g-moles/liter , 370.0505321 I the critical temperature of propane as computed from the Benedict-Webb-Rubin equation of state, °K h2.388h107h I the critical pressure of propane as computed from the Benedict-Webb-Rubin equation of state, atmt 0.08207 I gas constant B, (liter-atm.)/(8Imole-K°) 0.097313 I Benedict-Webb-Rubin constant B 6.87225 - Benedict-Webb-Rubin constant A ° 508256. I Benedict-Webb-Rubin constant 0° 0.0225 - Benedict-Webb-Rubin constant b 0.9h77 I Benedict-Webb-Rubin constant a 129000. I Benedict Webb-Rubin constant c 0.000607175 I Benedict-Webb—Rubin constant a 0.022 I Benedict-webb-Rubin constant 7 508256.13h5 I aCo 1.388h57o33x10‘6 - bc° 2.8hh5691h1x10'12 - cco 0 I dC o 0 I eCo 1 I number of sets of guesses of Tr2, Tr3, Trh, Dre, Dr3, Drh 1.005 I first guess of Tr2 1.010 I first guess of Tr3 1.017 I first guess of Trh 1.002 I first guess of Dr2 1.00h 3 1.007 I first guess of Drh first guess of Dr val um 1‘81 1101 ab 1i APPENDIX VII Program TOTRANGE The objective of program TOTRANGE is to vary simultaneously the values of reduced temperature (Tr) and reduced density (Dr) at x and x3 2 until G ", G ", 02"', and 03"' become equal to zero. The composition 2 3 range is split into only three intervals. Therefore, possible values of :2, x3, and xh could be 0.33, 0.67, and 1.0, respectively. This program follows the flow diagram in Figure 6, page 38. Program TOTRANGE makes use of the known values of Tr and Pr for both pure components. Because four variables (Trz, Tr3, Dr2, Dr3) are changed simultaneously to satisfy four known conditions, it should be possible to converge to a solution. Negative values of Tr and Dr are not permitted. The convergence method used is similar to Newton's convergence method for systems of equations.he The main body and subroutine PHICALC are listed for reference. The nomenclature is identical to that listed on page 10h except for the vari- able PHI which is defined as the sum of the squares of the four above listed 0 quantities. The main divisions of program TOTRANGE are: 1. Program TOTRANGE-Main Body 2. Subroutine PHICALC 3. Subroutine WBKAY h. Subroutine POLY 5. Subroutine PRCALC 6. Subroutine G2CALC 7. Subroutine G3CALC 8 a Data. 125 PR! DI CCCOCCCC$QDRDXXXXTDDZTDDZDCD s 1 F.» a 2 NU fl: C C 126 PROGRAM TOTRANGE DIMENSION XIZS)! WIZS): PC‘25Ie TCI25)9 PK(25)9 TK(25’9 FRIES). I TRI25)O DR‘25)9 ZR(25)9 2(5)9 FUNC(5)O 6(4). XX(7)9 TINCI(4)9 2 DINCI(4)0TINC(4)ODINC(4)0PHIN(7)ODELXX(7)e PXIIO)OD(794)O 3 PPDXC7IO DELTAXX(6)9 DA(5e4)e PDXI7): 88(5e4)s RZRWRPRIZS) COMMON IODCBQTCBOCOACOBCOCCODCOECOCIOC20C30C4'CSOCO!C7QC8QC90XOWO I DReDRIDEReDRZDEReTReTRIDEReTRZDERQTR3DER9PR0PRIDER0PRZDER0PR3DER 2 eTCePCeTKePKsQQeTOeTOIDEReTOZOERoTNeTNIDEReTNZDEReTTeTTIDER 3 CVOVIDERQYOYIDEROYZDERQZKOZKIDER9ZKZDERQZRQZOSQFUNCQGZDEROGBDERO 4 GePHIePZRWRPR I FORMAT (518.11) 2 FORMAT (/913H CALCULATION eIEs/I 3 FORMAT {/03H Xs4EZIeIOe/93H TR94E21eIOe/03H DRe4E21eIOs/03H PR: 1 4E2IeIOe/93H ZRs4E21eIOe/98H PZRWRPRCIGXQZEZIelO) 4 FORMAT (1H1) 6 FORMAT (I2) 7 FORMAT (8H RESULT598X06HPHI = eElBeIO) 8 FORMAT (10H GZDERZ = eElBelOe3X09HG3DER2 = eE18e1003X99HGZDER3 3 e I EIBeIOe3X99H63DER3 = sEIBelOe/I 9 FORMAT (9H TR(2) = 9E1801094X98HTR(3) = oElBelOe4XeBHTRC4) 3 s 1 ElBeIOe/egH DRIE) = 9E18e1004X98HDRISI = eEIBelOs4XeBHDR‘4) = o 2 ElBeIO) ll FORMAT (25Xo34HFOUR VARIABLES AND FOUR CONDITIONSe4(/)) PRINT 11 READ 19 DCBQTCBQPCBQRQBOeAOeCOeBOAeCeALpHAsGAMMAeACeBCQCCODCOEC PRELIMINARY CALCULATIONS C1 = (DCB*TCB*R)/PCB C2 = C1*BO*DCB C3 = (AO*DCB**2)/PCB C5 = CI*DCB**2*B C6 = (DCB**3/PCB)*A C7 = C6*DCB**3*ALPHA C8 = (C6/(A*TCB**2)I*C C9 = DCB**2*GAMMA CI=F9 C2=Go C3=Ho C4=Io C5=Jo C6=Ko C7=Lo C8=Mo C9=N READ 69 NUMDATA DO 39 NN=19NUMDATA 200 READ Iv (TR(I)oI=293)e(DR(I)9I=293) DELTA = loE-OS X(II = OeO X(2) = Oe6215 X(3) 3 0e8723 XI4) = 1.0 TRII) 3 leO DRIII 8 100 PRII) = 1.0 29(1) 3 IeO TR(4) 3 100 DRI4) 3 leO RRI4) 3 IeO ZR(4) 3 100 00 201 I = 194 201 CALL WBKAY PCII) I PK(1)/PR(1) 10 15 20 66 69 7O 71 77 72 73 74 75 41 43 42 127 TC(1) TK(I)/TR(I) PC(4) PKIQI/PRI4) TCI4) = TKI4)/TR(4I CALL PHICALC PHINIII 3 PHI PRINT 7. PHI PRINT 9. (TR(IIOI=ZQ4IOIDR(I)9I=204) PRINT 89 (G(I)OI=194’ DO 37 M = 1050 PRINT 20 M DO 18 L=194 DI59L) = ’GIL) DO 10 13192 XXII) = TRII+11 XXII+2I = DRII+11 DO 20 J3104 XX(J) = XXIJ) + DELTA DO 15 I=102 TR(I+1) = XX(I) DRII+11 = XXII+21 CALL PHICALC XXIJ) = XXIJ) - DELTA DO 20 1:104 D(J9I) = (6(1) + DI59111/DELTA DO 66 1:195 DO 66 J=le4 DA(I9J) = DIIeJ) DO 74 13194 DO 71 J=194 DIVISOR = DA‘IQJ) IF (DA(IeJ))69971969 DO 70 K = 195 DA(KoJ) = DA(KgJ)/DIVISOR BB(K¢I) = DA(KoJ) CONTINUE DO 74 J=194 IF (DA(I9J))72974072 DO 73 K=I05 DAIKQJ) = DA(K9J) - BB(KQI) CONTINUE PDX(5) = .10 DO 75 MA = 194 I = 5 - MA PDXII) 3 De L = I + 1 DO 75 K=L95 PDXII) = PDXII) - BB(KoI)*pDX(K) DO 41 N=194 DELTAXXIN) = PDX(N) XXIN) = XXIN) + DELTAXXIN) DO 88 KK=1015 DO 42 1:192 TRII+1I = XX(I) DRII+11 = XX(I+2) 87 88 79 34 37 38 1 39 CALL PHICALC PRINT 7. PHI PRINT 99 (TR(1,91329419‘DR‘IIQI3294) IF (PHIN(11-PHI)87987079 DO 88 N = 104 DELTAXX(N) = DELTAXX(N)/2o XX(N) I XX(N) - DELTAXXIN) GO TO 38 PHIN(1) = PHI PRINT 8. (G(I)ol=194) IF (PHI - loE-IO)38938037 CONTINUE PRINT 39 (XII)9I=IO4)0 (TR(I)eI=Ie4)9 (ZRII)OI=1'4)O (PZRWRPRII)9I=203) PRINT 4 END (DR(I)eI=le4)e(PR(I)eI=ls4)9 129 SUBROUTINE PHICALC . DIMENSION XIZS): WI2510 PC1251! TCI25)9 PK(25)0 TK(25)9 PRIZS): 1 TR(2519 DR(25)e ZRIES). 2(5). FUNG(5)9 6(4). PZRWRPRIZS) COMMON IQDCBQTCBoCoACoBCoCCoDCsEC¢C19C29C39€4QC50C60C79C89C90X9We 1 DRoDRIDERoDRZDERoTReTRlDERvTRZDER¢TR3DERQPR9PRIDERQPRBDEROPRBDER 2 eTCePCeTKcPKoOO‘TOoTOIDERsTOEDERoTNeTN1DEReTNBDEReTTQTTIDER 3 9V.VIDERQYQYIDEReYZDERoZKoZKlDERoZKZDERoZRsZeSoFUNCQGZDERoG3DER9 4 GoPHIoPZRWRPR DO 10 I=293 CALL PRCALC TCII) = TKII)/TR(1) PC(I) = PK(I)/PR(I) DO 14 1:203 DO 11 K=ls4 Z(K) = XIK) FUNCIK) I TCIK) S = 2(1) CALL POLY (TCIDERQTCZDEROTCBDER) DO 12 K3194 FUNC(K) = PC(K) CALL POLY (pCIDERQPCZDERQPCBDER) TRIDER I “TR(I)*TCIDER/TC(I) TRZDER *2.*TRIDER*TCIDER/TC(I) - TR(I)*TCZDER/TC(I) TR3DER 60*TK(I)*TCIDER*TCZDER/TC(I)**3‘60*TK(11*TCIDERii3/TCIII 1 **4 - TKII)*TC3DER/TC(I)**2 PRIDER = -PR(I)*PCIDER/PC(I) PRZDER = -2.*PRIDER*PCIDER/PC(I) - PRII)*PCZDER/PC(I) PRBDER = 60*PK(1)*PC1DER*PC2DER/PC(I)**3-6.*PK(I)*PCIDER**3/PC(I1 1 **4 - PKII)*PC30ER/PC(1)**2 CALL G2CALC CALL GBCALC G(2*I-3) = GZDER G(2*I-2) = G3DER PHI = G(I)**2 + G(Z)**2 + G(3)**2 G(4)**2 END 129 SUBROUTINE PHICALC ' DIMENSION X(25)0 W‘ZS). PC1251! TCIZSIQ PK(25)0 TKI25)Q PRI25)9 I TRIESIO DR(2519 ZRIZSI: 2(5). FUNCISIO G(4)e P2RWRPRI25) COMMON IoDCBeTCBeCoACeBCoCCeDCoECOCI9C29C39C40C59C6oC7eCBoC9eXeWe 1 DRQDRIDERQDRZDERQTRsTRIDERQTRZDEROTR3DER9PR9PR1DEROPRBDEROPR3DER 2 oTCePCeTKoPKeOQeTOeTOIDEReTOZDEReTNQTNIDEROTNZDEROTTOTTIDER 3 9V0VIDERQYOYIDERcY2DEReZKOZKlDEReZKZDEReZRQZOSeFUNCeGZDER063DER9 4 GQPHIQPZRWRPR DO 10 I=293 CALL PRCALC TCII) = TK(1)/TR(I) PCII) : PKIII/pR(I) DO 14 1:293 DO 11 K=le4 Z(K) = XIK) FUNC(K) ‘ TC(K) S = Z(I) CALL POLY (TCIDEReTCZDERoTC3DER) DO 12 K3194 FUNC(K) = PCIK) CALL POLY (PCIDERQPCEDER¢PC3DERI TRlDER 8 ”TRII)*TC1DER/TC(II TREDER *2e*TR1DER*TC1DER/TC(I1 - TR(I)*TCZDER/TC(I) TR3DER 6e*TK(I)*TC1DER*TCZDER/TC(I)**3‘6e*TK(I)*TCIDER**3/TC(I) 1 **4 - TKII)*TC3DER/TC(I)**2 PRIDER = -PR(I)*PC1DER/PC(I) PRZDER = -2.*PRIDER*PCIDER/PC(I) - PR(I)*PC20ER/PC(I) PR3DER = 6.*PK(I)*PCIDER*PCZDER/PC(1)**3-6.*PK(I)*PCIDER**3/PC(I1 1 **4 - PKII)*PC3DER/PC(I)**2 CALL GZCALC CALL G3CALC G(2*I'3) = GZDER G(2*I-2) = G3DER PHI = G(I)**2 + 6(2)**2 + G(3)**2 G(4)**2 END eeeeeee nnnnnnn APPENDIX VIII Program CHECK The purpose of program CHECK was to check the equations for G" and G"' and their supporting equations in subroutine GZCALC and GBCALC. Beginning with the expression for G, the derivatives were calculated by numerical differentiation and from the mathematical expressions for the derivatives. These two answers were then compared, and if they were nearly identical, the mathematical equations for the 0 quantity were assumed to be correct. This program follows the flow diagram in Figure 7. page 39. The main division of program CHECK are: 1. Program CHECK-Main Body 2. Subroutine PRCALC 3. Subroutine G2CALC h. Subroutine GSCALC 5 a Data 130 C C ”DUN“ 1 2 3 4 UIDUNF‘ 131 PROGRAM CHECK DIMENSION XIES)! WI25)0 PCI25)0 TCIBS)! PK‘ESII TK‘ZS)! PRIZS). TRIES). DR(25)O ZR(25)O 2(5)! FUNCISIO 6(4). XXI7). TINCII4IO DINCI(4)9TINCI4)9DINCI4)9PHIN‘7IODELXX‘7)e PX‘IO)OD(703)O PPDXI710 DELTAXXI6IO DAI49310 PDX‘7)! 88(493) QDCCIZO)‘ DKI20) COMMON IQDCBOTCBQPCBORQBOQAOeCOeBsAeCQALPHAeGAMMAeACQBCOCCODCQECO C19C20C39C41C5eC69C70C81CQQXeW9DR9DRIDEReDREDEReTReTRIDER.TRZDERe TR3DERePRePRIDERePRZDERePR3DEReTCePCeTKoPKeOQeTOeTOIDEROTOZDERQ TN.TNIDER.TN20ER.TT.TTIDER.V.VIDER.Y.YIDER.Y2DER.ZK.ZKIDER. ZKaDEROZRQZ‘SQFUNCOFSTDEROSECDEROTHRDEROcaoERQGBDEROGOPHI OGGQGIDER FORMAT (EIBeII) FORMAT (12) FORMATI25X913HCHECK PROGRAM94(/)) FORMAT (3EZOOIO) FORMAT I/96H CHECKe3EZOe109/96X03E20e1003(/)I PRINT 3 READ 19 DCBOTCBOPCBOR‘BOOAOQCO0BQAOCIALPHA96AMMA9ACOBCOCCODCOEC PRELIMINARY CALCULATIONS C1 = (DC8*TCB*R)/PCB C2 = C1*BO*DCB C3 = (AO*DCB**2)/PCB C5 = C1*DCB**2*B C6 = (DCB**3/PCB)*A C7 = C6*DCB**3*ALPHA C8 = (C6/(A*TCB**2))*C C9 = DCB**2*GAMMA C1=Fc C2=Go C3=Hc C4=Io C5=Je C6=Ke C7=Lo C8=Me C9=N 59 58 1 I DELTA = 0.00001 READ 29 NUMDATA DO 64 NN=10NUMDATA READ ls XII) ASSIGN 60 TO JOHN I = 1 TK(I) = (513-3-186o4392126*X(I)+20 0413825*X(I)**2-558.37O4113* X(I)**3+12101957936*X(I)**4+459o67)/1o8 DKII) = 2.3 + 20*XII) + 30*XII)**2 TC(I) = 530. + 200.*X(I) - 400.*X(I)**3 TCIDER 2 +200. ~1200.*X(I)**2 TCEDER = ~2400e*X(I) TC3DER = ~2400. PC(I) = 25. + 1400*X(I) - 120.*X(I)**3 PClDER I 140- - 360.*X(I)**2 PCZDER = -720.*X(I) PCBDER = ~720o DCC(I) = C1*PC(I)/TC(I)/R TR(I) = TKII)/TC(I) DR(I) I DK(I)/DCC(I) CALL PRCALC PK(I) = PCII)*PR(I) PRINT 109 TC(I)9TK(I)oTR(I)oDCC(I)qDK(I)oDR(I)oPC(I)oPK(I)ePR(I)e X(I) TRIDER = -TR(I)*TC1DER/TC(I) 132 TRZDER = -2.*TRIDER*TCIDER/TC(I) - TR(I)*TC2DER/TC(I) TRSDER = 60*TK(I)*TCIDER*TCZDER/TC(I)**3*69*TK(I)*TCIDER**3/TC(I) 1 **4 - TKII)*TCBDER/TC(I)**2 PRIDER = *PR(I)*PCIDER/PC(I) PRZDER ~29*PRIDER*PCIDER/PCII) ~ PRII)*PCZDER/PCII) PRBDER = 69*PKII)*PCIDER*PC2DER/PC(I)**3“69*PK(11*PCIDER**3/PC(I) 1 **4 - PK(I)*PC3DER/PC(I)**2 CALL GZCALC CALL GBCALC GO TO JOHN XII) = XII) + DELTA PRINT 109 669 GlDER9 620ER9 GBDER PKl = PKII) 661 = 66 611 GlDER 621 = 62DER G31 = GBDER ASSIGN 61 TO JOHN GO TO 59 XII) = X(I) ' DELTA I! PGWRXDT = (66 - GGl)/DELTA PIWRXDT = (GIDER - 611)/DELTA PFWRXDT = (GZDER - 621)/DELTA PPWRXDT I (PK(I) - PK1)/DELTA PRINT 109 P6WRXDT9 P1WRXDT9 PFWRXDT9 PPWRXDT ASSIGN 62 TO JOHN GO TO 59 DK(I) = DK(I) + DELTA ASSIGN 63 TO JOHN GO TO 58 PFWRPXT = (GZDER - 621)/(PK(I) “ PKI) PGWRPXT = (66 - GGl)/(PK(I) - PKI) PIWRPXT = (GIDER - 611)/(PK(I) . PK1) PRINT 109 P6WRPXT9 PIWRPXT9 PFWRPXT GIDERN = PGWRXDT - PGWRPXT*PPWRXDT GZDERN = PIWRXDT - P1WRPXT*PPWRXDT GBDERN = PFWRXDT - PFWRPXT*PPWRXDT PRINT 119 6119 6219 6319 GIDERN9 GZDERN9 G3DERN END 133 Program CHECKAMain Body Purpose: Input: Details: To check the equations for G" and G"' and all of the support- ing equations in subroutine GZCALC and C3CALC. P-V-T data on a reference compound (in this case propane). This data includes the critical properties, the gas constant R, eight Benedict-Webb-Rubin constants, and five constants for computing Co as a function of temperature. Also, the value of composition (x1) at which the derivatives are to be calculated. This program computes derivatives 0', G", and G"' from the mathematical expressions for these derivatives, and also by numerical differentiation. If each pair of results are nearly identical, the equations are assumed to be correct. The expression for the derivative of'a function F(T,D,P,x) with respect to x at constant T and P was shown on page 39 to be: 2+ ' (BF/ax)T.P - (BF/ax)T’D - (BF/3P)T.x(aP/ax)T’D (1) where F can equal G, G', or 0". [Up to third derivatives are needed for To and Pc. For this program it is sufficient to use simple, arbitrary functions for these variables. The functions were selected such that they approximated the To and Po values for pure ethane and pure n-heptane. The expressions used were: Tc I 530 + 200x - I400x3 (2) Pc - 25 + 1&0: - 120x3 (3) Equations (2) and (3) are readily differentiable with respect to x. This differentiation was done by hand and the resulting formulas put into the program. Method of Checking: Fortran Nomen- clature: 13h For the sake of accuracy, an expression was devised for Dk which approximated the experimental data of Kay23. For purposes of this check program, a constant value of Dk could also have been used. The expression for Dk is: Dk - 2.3 + 2x + 3x2 (A) The expression used for Tk is the same one which is used in sub- routine WBKAY. Extensive use is made of subscripted variables in this pro- gram, although this is unnecessary. This is done for two rea- sons. IFirst, program CHECK was prepared from other programs which were already written, and these had subscripted variables in them. Second, since the equations in this project actually use subscripted variables it is better to check those specific equations. The specific subscript used throughout this program 18 la A second person checked over this program. DCC I pseudocritical density, moles/liter CC I G GlDER I G' G2DER I G". G3DER I G"' GlDERN I numerical value of G' I (aG/ax)T P 9 G2DERN I numerical value of 0" I (aG'/ax)T P G3DERN I numerical value of G"' I (30"/3x)T P ’ 011 I name of variable which is used to save the first value of G' » 021 I name of variable which is used to save the first value of 0" G31 I name of variable which is used to save the first value or G“! PKl I name of variable which is used to save the first value of Pk PGWRXDT I (EC/ax)T D O 3 9 lenan (as /ax)T’D 135 PFWRXDT I (30"/8x)T’D PPWRXDT . (aP/ax),1..D PGWRPXT - (ac/310T.x PlWRPKT - Wen/310.1,,x PFWRPXT - (MN/313),“x 136 Subroutine PRCALC - Sam; as for program LINEAR Subroutine GZCALC - Same as for program LINEAR except that G and G' are also computed. G is calculated from equation III-18 and G' is calculated from equation III-19. Subroutine G3CALC - Same as for program LINEAR. Data - Same as for program LINEAR except that no values for Dr and Tr are put in the data. Instead, after the data card for eCO, values are read for NUMDATA and the compositions at which a check is to be made. The number of values of x must be equal to the value of NUMDATA. APPENDIX IX Generated Pseudocritical Curve Results - Approaches I and II Approach I Using as known conditions, experimental critical temperatures, pres- sures, and compressibility factors and the fact that 0" equals zero at the critical point, an attempt was made to generate pseudocritical tem- perature and pressure curves. Representative data points are listed in Tables 11 and 12. Table 11 contains the generated values of pseudocriti- cal temperature as a function of composition, while Table 12 contains the generated values of pseudocritical pressure as a function of composition. Figure 20 is a plot of the pseudocritical temperature data. The temperature data in Figure 20 show a large amount of instability. The stability of each curve was directly dependent on the original guess made for Tc. * (shown as Tr in Figure 20). In general, Figure 20 OS 05 contains two shapes of curves. The S-shaped curves (Tr I 1.02. 1.0202) 05 exhibit the greatest irregularity. Reference to a typical diagram20 of compressibility factor as a function of reduced pressure and temperature indicates that there are two values of reduced pressure which give the same compressibility factor at the same reduced temperature. Because reduced pressure is a function of reduced density and temperature, there must also be two values of reduced density for a given compressibility factor. If the computer used the high and low values of reduced density in an irregular manner, as the composition range was crossed, erratic results such as the S-shaped curves in Figure 20 could occur. *Tc 05 represents Tc at x I 0.05. 137 138 TABLE 11. Generated values of pseudocritical temperature.(9K) in Approach I Tr.05 Tc.05 Tc.10 Tc.15 Tc.2o Tc.25 Tc.3O 1.035 ' 1.032 518.99 508.52* 1.031875 519.06 508.65* 1.03175 519.12 508.78 199.82 191.53 181.61 185.55“ 1.0315 519.21 509.06 500.10 193.10 190.92* 1.031 519.50 509.69 $02.05*” 1.0305 519.75 510.12 501.16 512.17** 1.03 520.00 511.21 507.76 510.67** 1.0298125 520.10 $11.56 $09.25* ' 1.02925 520.38 512.62 511.63* 1.0285 520.78 511.19 521.11* 1.025 $22.51 $23.56* 1.0202 525.00 510.11 511.12* 1.02 v525.10 510.86 516.83* 0.99 511.01 532.57“ Note: T0 for pure n-heptane I 510.51 °K “Indicates that 0" did not converge to zero at that value of x. "Indicates that G" did not converge to a single value, but that it was near zero. 139 TABLE 12. Generated values of pseudocritical pressure (atm.) in Approach I Tr.05 Pc.05 Pc.10 Pc.15 Pc»20 Pc.25 Pc.30 1.035 ' 1.032 23.28 23.35“ 1.031875 23.32 23.15* 1.03175 23.35 23.56 21.15 26.03 28.95 36.589 1.0315 23.12 23.77 21.96 27.26 32.72“ 1.031 23.56 21.10“ 1.0305 23.69 21.65 27.15 36.88** 1.03 23.82 25.12 29.15 19.37.. 1.0298125 23.86 25.30 29.81! 7 1.02925 21.01 25.85 32.21* 1.0285 21.19 26.62 32.21!I 1.025 25.02 30.62% 1.0202 26.07 37.06 32.25“ 1.02 26.12 37.31 33.31* 0.99 32.12 31.159 Note: Fe for pure n-heptane I 26.916 atm. ”Indicates that G" did not converge to zero at that value of x. ”*Indicates that 0" did not converge to a single value, but that it was near zero. lhO Figure 20. Pseudocritical temperature as a function of composition (Curve parameters are Tr at x 550* ‘ 0005) 510 ' Tc, OK 530‘ 1.03 $20- .‘ 9' 1.032! X ‘ - 510‘” , \: 02982 , 1.0305 ' .031 500 ~1- '3 ‘ 1.0315 190-4— ' 1.02175. 180 . % % f 1r + % 0 0.05 0.10 0.15 0.20 0.25 0.30 Mole fraction ethane 110 Figure 20. Pseudocritical temperature as a function of composition (Curve parameters are Tr at x = 0.05) 560+ 550 ~11— 510 ' Tc, OK 530— 1.03 520-* 1.029 0 510‘ 1.0305 .031 500-7 ‘ 1.0315 hw—U— . 1.02175. 18o . % ~ % f f % + 0 0.05 0.10 0.15 0.20 0.25 0.30 Mole fraction ethane lhl The second type of curve in Figure 20 is parabolic. The reason for this shaped curve is not as obvious as that for the S-shaped curves. Although a thorough mathematical analysis is not possible, it could be speculated that errors in the numerical differentiation eventually became magnified enough to cause unreasonable results. Each of the curves in Figure 20 was terminated before the entire composition range had been traversed because no further convergence could be obtained. All values of Tr which showed any promise of yielding satisfactory results were 05 thoroughly investigated. This was especially true of the transition area between the parabolic and S-shaped curves shown in Figure 20. The programs which began at the pure ethane end of the composition range instead of at the pure n-heptane end, showed no improvement in behavior. The same instability was still indicated. Likewise, the chang- ing of the increment size of x was not an improvement. When butane was used as a reference compound instead of propane, the curves were similar in shape to those of Figure 20, but the butane curves were offset slightly from the propane curves. This approach was eventually abandoned because the condition that G"' also equals zero at the critical point was not being satisfied by the results. Approach II The approach used here was to generate pseudocritical temperature and pressure curves using as known conditions, experimental critical tem- peratures and pressures and the fact that G" and G"' equal zero at the critical point. A mathematical convergence method was required for the calculations. The first convergence method used, the double linear process, was not 112 completely satisfactory. This method sometimes diverged if the quantities (AX1.AY2) and (AxaoAYl) [See page 107 for definitions of these quantities.] became nearly equal. These quantities being nearly equal meant that the process was near an unstable region and thus, behaved inconsistently. The usefulness of this convergence method was limited by its tendency ’ to diverge in certain cases. If only one or some small number of conver- gences are required, then this procedure may be of some use. However, the percentage of failures is too high to use this method without some dis- cretion. Of course, the more linear the equations, the better the method E will work. The success achieved in this project from using the double linear con- vergence process appeared highly dependent on the starting values for the variables Tr and Dr at me and x3. This pointed upagain the possibility of an error becoming magnified in the system. At its best, the program using linear convergence did not proceed higher than 1 I 0.35 before lack of convergence caused termination of the program. Program LINEAR is listed in Appendix VI in lieu of specific results which are meaningful. Another convergence method which was used to solve systems of two equations and two unknowns was Newton's method for systems of equations (also known as the Newton-Raphson method).h2 This procedure extrapolated using the partial derivatives at a single point. In cases where both the linear method and Newton's method worked, Newton's method always required less iterations. However, the linear method often gave convergence when Newton's method did not. For this reason the linear method was considered the better method for use in this research. For solving systems with more variables than conditions, the method 113 of steepest descent as defined in this thesis, was not a good convergence method. There appeared to be two reasons for the difficulties with this method. First, the assumption that d¢/ dz I A0/Ax was not always a good one to make. Even for the simplest of the linear cases, d¢/ dx is twice as large as A0/Ax. ’This means that the value calculated for p was not always as good as it could have been. Second, the amount of interplay between the quantities (0")2, 03")2, and (G3'”)2 in the expression for 0, required the method to make many more iterations than would have been necessary if the method did not use a function such as 0. The convergence method which used an approach similar to the Newton- Raphson methodl‘2 to simultaneously vary the variables (Tr2, Tr3, Tr , Dr2, Dr3, Dru) until 6" and G"' were zero at x2. 33. and 1h gave 8005 303‘ vergence. The program using this convergence method gave results for almost all sets of input tried. The objective was to obtain what appeared to be a good set of starting values (i.e. Tr and Dr at x2, x3.and xh) and then use these values as starting values in program LINEAR. Table 13 lists a typical set of results from this procedure. The generated values are reasonable until x I 0.35. At this point large changes take place in Tc and Pc. Using double linear convergence, no set and Pc of values was found for Tc which made 6" and 0"' equal 0.10 0.10 to zero at x I 0.35. These large changes in To and Pc at x I 0.35 could once again have been caused by the magnification of an error. The values of (BZr/aPr)Tr are quite regular which indicates that the computer did not use the low and high values of reduced pressure in an arbitrary manner 0 111 Table 13. Generated Tcm and.Pcmvalues from program LINEAR.using a set of starting values from program DECLINE ' x Fe Tc (azr/aPr)Tr 0 26.95 atm. 510.51 °K -- 0.025 26.77 533.06 -10.73 0.050 26.61 525.80 -5.12 0.100 26.10 511.55 -2.27 0.150 25.91 196.31 -1.27 0.200 21.60 177.33 -0.70 0.250 21.79 151.15 -0.29 0.300 16.97 111.61 +0.03 0.350 11.63 338.62 - +0.26 “Nan, I ids .52. L13. 3175 3688 “WWW R]! m" L" Y” I" u "I H fl 1293 Ilfl'llmumumum