A generalized computer program for analysis of mixture refrigeration cycles

H . A . Connon

Key words: refrigerant, non-azeotropic mixture, computer program

Un programme d'ordinateur g6n6ralis6 pour I'analyse des cycles frigorifiques des m61anges

L "int#r#t croissant pc ur les m#langes de frigorig~nes non az#otropiques destines aux pompes ~ chaleur et autres appareils a mis en #vidence la n#cessit# de disposer d'une m#thode d'analyse des propri#t#s de ces m#langes pour les

cycles frigorifiques. Les tables de propri#t#s thermodyna- miques ne suffisent pas, parce qu'elles ne donnent pas de renseignements pour la r#gion biphas#e cO la temperature et les compositions de la vapeur et du liquide ne sont pas constantes ~ pression constante, Cet article d#crit le pro- gramme d'ordinateur CYCLE qui peut effectuer les calculs des propri#t~s thermodynamiques pour les m#langes non az~otropiques sous-refroidis, en r~gime biphasique, sur- chauff#s et permet d'analyser un cycle frigorifique simple.

The growing interest in non-azeotropic re- frigerant mixtures for heat pumps and ap- pliances has led to the need for a method of analysing the refrigeration cycle properties of such mixtures. Thermodynamic property tables are not sufficient, because they give no infor- mation about the two-phase region, where

temperature and vapour and liquid compo- sitions are not constant at constant pressure. This paper describes the computer program CYCLE, which can perform thermodynamic pro- perty calculations for subcooled, two-phase, and superheated non-azeotropic mixtures and can analyse a simple refrigerating cycle.

In a recent paper Connon and Drew 2 reported the successful calculation of thermodynamuc property tables for a non-azeotropic misture of R 13B1 and R 152a using the Redlich-Kwong-Soave (RKS) equa- tion of state. 11,12 The equations and constants are reproduced for the reader's convenience in the appen- dix of this paper. The saturation tables for the com- position of each mixture were presented in even increments of pressure, and the bubble-point and dew- point temperatures were given at each pressure. Unlike pure refrigerants and azeotropic mixtures, the tempera- ture and the vapour and liquid compositions of non- azeotropic mixtures do not remain constant at constant pressure as the refrigerant evaporates or condenses. Thus, the saturation tables gave only bubble-point and dew-point information but no information about the two-phase region. An additional computer routine was written to provide the information necessary to fill in the temperature and quality lines in the two-phase region of a pressure-enthalpy diagram. That routine calculated pressure, enthalpy, and vapour and liquid compositions from the temperature, refrigerant quality, and total composition.

The completed p-/-/ diagrams were useful for

The author is a Research Chemist with Freon Products Laboratory, E.I. du Pont de Nemours and Company, Wilmington, DE, USA. Reprinted by permission of the American Society of Heating, Refrigerating and Air-Conditioning Engineers, Inc.. from ASHRAE Transactions 89 2A (1983) (Atlanta: American Society of Heating. Refrigerating and Air-Conditioning Engineers, Inc. 1983),

analysing simple refrigerating cycles, but it was felt that a computer program capable of such an analysis would be faster and more accurate than reading values from the diagram. This led to the development of the program CYCLE.

Program descr ip t ion

Fig.1 represents a simplified pressure-enthalpy diagram for a non-azeotropic refrigerant mixture. Con- stant temperature lines are included in the two-phase region to illustrate the temperature change at constant pressure. Although the total refrigerant composition (vapour+ liquid) is fixed, the vapour and liquid com- positions are always different in the two-phase region and change continuously at constant pressure be- tween the bubble point and the dew point. The CYCLE program performs the thermodynamic property calcu- lations that are necessary to analyse the refrigerating cycle of a mixture like the one drawn on Fig.1. The calculations include condensation (1-3), liquid sub- cooling (3-4), isenthalpic expansion into the two- phase region (4-5), evaporation (5-7), superheating (7-8), and isentropic compression (8-9).

CYCLE is designed as a modular program in which each of the subroutines performs a discrete thermody- namic property calculation. The program has been generalized to allow analysis of binary mixtures of fourteen different refrigerants. All the pure-component properties are available in the form needed by the program. The user need only supply at least one

Volume 7 Numero 3 Mai 1984 0140-7007/84/0031 67-06S3.00 © 1984 Butterworth S" Co (Publishers) Ltd and IIR 167

7OO~-

65155 Weight percent Freon 15BI/152a pressure-enfhalpy diagram

r~

n

1u u wu ~u o v -~J ou ou ru ou ~u = w ~U

Enthalpy, Btu Ib "~ above saturated l iquid a t - 4 0 = F

Fig. 1 A s impl i f ied pressure-entha lpy d iagram for non -azeo t rop i c refr igerant mix ture

Fig. 1 Diagramme pression-enthalpie simplifi~ pour un m~lange de frigorig~nes non az~otropique

pressure-temperature-composition point for the mix- ture to calculate the interaction constant. This in- teraction constant is an empirical parameter charac- terizing the binary mixture. If the data are not available, the interaction constant can be set at 0; however, this does not mean the mixture is being treated as ideal. It has been shown ln4 that for ideal mixtures, the in- teraction constant is not, in general, equal to 0, since the mixing rules do not always result in an ideal solution. Subroutines RFLASH, PFLASH, TFLASH, HSCALC, AND CaM PRS are general routines that can be used as they are to incorporate the RKS mixture equations into equipment-modelling programs (for heat pumps, air conditioners, appliances, etc.) that currently contain pure-component equations of state. CONDNS. EVAP, and EVAP2 were designed specifi- cally for analysis of the simple refrigerating cycle shown in Fig.l. However. modifications of these routines to suit other applications would not be dif- ficult. The function of each subroutine is described below. The f low chart of Fig.2 is a schematic repre- sentation of the program logic.

In addition to the refrigerant property information described in the appendix, the simple refrigerating- cycle calculation requires knowledge of the circulating composition and only four pieces of input data: the condenser mean temperature, TCON (2): the expansion-device inlet temperature, TSUB (4) (or the degrees of subcooling); the evaporator mean tempera-

ture, TE (6) or the evaporator pressure, PE; and the compressor suction temperature, TSUP (8) (or the degrees of superheat).

Circu/ating composition. After the refrigerant data is read into the program, CYCLE asks for the composition of the circulating refrigerant. In most cases, the refrigerant composition will be known and can be entered by the user. For the special case where refrigerant may have collected in an accumulator and the actual circulating composition is unknown, it can be calculated from the refrigerant charge composition by either of three routines: RFLASH, PFLASH, or TFLASH.

RFLASH calculates refrigerant quality and va- pour and liquid compositions from the known tem- perature, pressure, and total composition. If there is liquid in an accumulator, RFLASH can use the tem- perature and pressure above the accumulator to cal- culate the circulating refrigerant composition, the ac- cumulator liquid composition, and the weight fraction of the original charge that is in the accumulator.

PFLASH and TFLASH are complements to RF- LASH and are more useful for design work. PFLASH calculates the pressure and the vapour and liquid compositions from a known temperature, refrigerant quality, and total refrigerant composition. Similarly, TFLASH calculates the temperature and the vapour and liquid compositions at a specified pressure, re- frigerant quality, and the total composition. Thus, the circulating composition can be calculated for a desired weight fraction of accumulation at a desired tempera- ture or pressure. Once the circulating composition has been established, it is used in all subsequent calcu- lations in CYCLE.

Ca/cu/ation at specific points. Before the pro- gram proceeds to the analysis of the refrigeration cycle, the user is given the opportunity to use RFLASH, PFLASH, or TFLASH to perform specific calculations anywhere in the subcooled, two-phase, or superheated regions. Each call to RFLASH, PFLASH, or TFLASH is followed by a call to sub-routine HSCALC, which calculates the liquid and vapour enthalpies and en- tropies from known temperature, pressure, and liquid and vapour compositions. Thus, a call to RFLASH can be used simply to check the refrigerant quality at any point in the two-phase region or to calculate the refrigerant enthalpy and entropy for a given pressure and temperature in either the subcooled, two-phase, or superheated regions. PFLASH orTFLASH can be used, respectively, to calculate the pressure or temperature, along with the enthalpies and entropies, at the dew point, bubble point, or any other point where a specific refrigerant quality is desired. PFLASH and TFLASH are also useful for predicting refrigerant composition cha- nges resulting from vapour leaks.

Refrigerating-cyc/eana/ysis. A fair-performance comparison between a non-azeotropic mixture and a pure refrigerant requires the proper choice of conde- nser and evaporator conditions for each. There appears to be no generally accepted convention at this time. Jakobs and Kruse 7 have compared pure refrigerants to mixtures by using the same condenser dew-point and evaporator-inlet temperatures for each. Stoecker and Walukas (1 981 ) performed their comparison by requir-

168 International Journal of Refrigeration

N° I I Circulating composition must be calculated I

Enter original charge composition CALL CHOICE

,..Choose RFLASH, PFLASH, OR TFLASH

I

i "°a° °°re c°m°°°en' °ar'm=°r" } Read mixture interaction coefficient k12

I I Is the composition of the circulating refrigerant

known? Is it the same as the charge composition? I

i

'Dew point or bubble point calulation?' 'Or calculation at some other point?'

~ Yes

CALL CHOICE Choose RFLASH, PFLASH, or TFLASH Perform T-P-X calculation in two-phase, superheated or subcooled regions

CALL RSCALC Calculate enthalpy and entropy

I Pressure from condenser mean temperature I CALL CONDNS I Input desired condenser mean temperature

I Yes

Enter weight fraction of component one in the circulating composition

i

No

Calculate the pressure at which the dew point I and bubble point temperatures average to give I the desired temperature

I 'Enter expansion valve inlet temperature' CALL RSCALC Use condenser pressure and specified expansion valve inlet temperature to calculate enthalpy of the liquid entering the expansion valve

I Analyse evaporator 'Choose type of calculation' 'Evaporator inlet temperature from known pressure

'Pressure from mean evaporator temperature - 2'

CALL EVAP Enter evaporator pressure Calculate inlet temperature and dew point temperature

I Isentropic compression Enter compressor inlet pressure and temperature CALL DEWPIT Check to see that the lefrigerant is superheated CALL RSCALC Calculate enthalpy and entropy before compression Enter compressor outlet pressure CALL COMPRS Calculate compressor outlet temperature for isentropic compression to outlet pressure

Fig. 2 Program logic flow chart for CYCLE

Fig.2 Orgonigrarome Iogiquo du programme "CYCLE"

ing the same condenser and evaporator UA values for pure refrigerants and mixtures. The initial version of the CYCLE program has been set up to compare the pure- refrigerant condenser and evaporator temperatures with the mean condenser and evaporator temperatures for a non-azeotropic mixture. The mean condenser temperature refers to the average of the dew-point and bubble-point temperatures at the condenser pressure, and the mean evaporator temperature is the average of the dew-point and evaporator-inlet temperatures at the

2

CALL EVAP 2 Enter desired mean evaporator temperature Calculate the pressure at which the evaporator inlet temperature and dew point temperature average to give the desired temperature

evaporator pressure. This convention is an arbitrary one and can be changed when the CYCLE program is incorporated into an equipment-modelling program.

Condenser pressure, The user must enter the desired condenser mean temperature, TCON (point 2 on Fig.l). Subroutine CONDNS wil l calculate the pressure, PCON, at which the dew-point (1) and bubble-point (3) temperatures average to give the desired mean temperature.

Isenthalpic expansion. The evaporator analysis must be preceded by calculation of the enthalpy of the

Volume 7 Number 3 May 1984 169

Table 1. Comparison of calculated and experimental refrigeration-cycle results for 65/35 weight % R 13B1/R 152a Tableau 1. Comparaison des r~su/tats calculus et exp#rimentaux du cycle frigorifique pour un m~/ange de 65/35% en masse de R 13B1/R 152a

Calculated Experimental

1. Condenser

Enter: mean temperature, TCON (2)=40°C (104°F) - - 40°C

Calculate: pressure, PCON 1.61 x 1 0 3 kPa (234 psia) 1.60x 103 kPa (232 psia)

Dew point temp (1) 44.2°C (111.6°F) - - Bubble point temp (3) 35.8°C (96.4°F) - -

2. Expansion device

Enter: inlet temperature TSUB (4)=17°C (63~F) - - 17°C

Calculate: enthalpy (4-5) 56.0 kJ kg -1 (24.1 Btu Ib -1) - -

3. Evaporator

Enter: pressure, PE=446 kPa (64.7 psia) - - 446 kPa

Calculate: inlet temperature (5) - 11.40C (11.5°F) - 1 t .70C (1 lOF) % vapour at inlet (5) 23.4% - - inlet composition, weight % R 13B1/R 152a

liquid 59.2/40.8 - - vapour 84.0/16.0 - -

dew point temp (7) ~0.1°C (31.9°F) - -

4. Compressor

Enter: suction temp, TSUP (8)=24°C (76"F) 24°C

Calculate: discharge temp (9) 81,60C (178.8°F) 108°C (226°F)

subcooled liquid (4). The expansion-device inlet tem- perature, TSUB (4), must be entered. If an exact amount of subcooling is desired, TSUB can be ob- tained by subtracting the degrees of subcooling from the bubble-point temperature of the condenser (3) calculated by subroutine CONDNS. TSUB and conde- nser pressure PCON are then used by HSCALC to obtain the expansion enthalpy.

Evaporator analysis. The user may choose be- tween two evaporator subroutines, EVAP and EVAP2.

Evaporator subroutine EVAP uses the expansion enthalpy calculated above and a desired evaporator pressure, PE, to calculate the evaporator-inlet tempera- ture (5), the refrigerant quality and the vapour and liquid composit ions at the evaporator inlet, and the evaporator dew-point temperature (7).

Subroutine EVAP2 is similar to CONDNS in that the enthalpy of the subcooled l iquid and a desired mean evaporator temperature, TE (6), are used to calculate the pressure at which the evaporator-inlet temperature (5) and the dew-point temperature (7) average to give the desired mean evaporator tempera- ture (6). The refrigerant quality and vapour and liquid composit ions at the evaporator inlet are also calculated.

Isentropic compression. Compressor suction temperature TSUP (8) must be supplied by the user. If desired, TSUP can be selected by adding the desired degrees of superheat to the dew-point temperature (7) calculated in the evaporator routine. Temperature TSUP and evaporator pressure, PE, are used by HS- CALC to compute the entropy at point 8. Subroutine COMPRS then uses the entropy and the condenser pressure, PCON, to calculate the isentropic compres- sion temperature (9).

Results

Table 1 shows a typical comparison between calculated CYCLE results and experimental calorimeter data for a heat pump run with R 13B1 and 152a at 65/35 weight percent. The table shows very good agreement between calculated and experimental con- denser pressures. The expansion enthalpy, calculated from the high-side pressure, PCON, and temperature TSUB (4), is used wi th the low-side pressure, PE, to give a calculated evaporator-inlet temperature (5) that agrees wi th the experimental value to well wi th in experimental error. The calculated compressor dis- charge temperature is for an isentropic compression

170 Revue Internationale du Froid

and is expected to be lower than the experimental value.

Conclusions The CYCLE computer program provides a con-

venient method for performing a variety of specific thermodynamic property calculations for two-phase, subcooled, and superheated non-azeotropic refri- gerant mixtures. The program has been generalized to allow analysis of binary mixtures of fourteen re- frigerants, and pure-refrigerant properties can also be calculated. For binary mixtures, at least one experimen- tal pressure-temperature-composition data point must be specified to calculate the interaction constant; however, data covering the ranges of composition, temperature, and pressure that are of interest will give better results.

So far, the agreement between the Redlich- Kwong-Soave calculated values and experimental ca- lorimeter data has been very good. Future plans include measurement of pressure-temperature-composition data for more mixtures (in order to generate the interaction constants) and incorporation of the CYCLE routines into equipment-modelling computer programs.

Program listings are available upon written request.

References

1 Asselineau, L., Bogdanic, G.. Vidal, J. Calculation of thermodynamic properties and vapor-liquid equilibria of re- frigerants Chem Engng Sci 33 (1978) 787-791

2 Connon, H. A., Draw, D. W. Estimation and application of thermodynamic properties for a non-azeotropic refrigerant mixture, IIR Meeting of Commissions BI, B2, El, and E2, Essen, West Germany (September 1981 )

3 Downing, R. C. Refrigerant equations ASHRAE Transac- tions 80, Part 2 (1974)

4 Edmister, W. C. Applied hydrocarbon thermodynamics. Part 32 - Compressibility factors and fugacity coefficients from the Redlich-Kwong equation of state, Hydrocarbon Processing 47 9 (1968) 239-244

5 Edminster, W. C. Applied hydrocarbon thermodynamics. Part 33- Isothermal changes in enthalpy and entropy from the Redlich-Kwong Equation, Hydrocarbon Processing 47 10 (1968) 145-149

6 Graboski, M. S., Daubert, T. E. A modified Soave equation of state for phase equilibrium calculations. 2. Systems con- taining CO 2, H2S, N 2 and CO IndEng Chem Process Des Dev 17 4 (1978) 448-454

7 Jakobs, R., Kruse, H. The use of non-azeotropic refrigerant mixtures in heat pumps for energy saving. IIR Meeting of Commission B2. Delft, The Netherlands (September 1978)

8 Mart in, J. J., Welshans, L, M., Chou, C. H.. Gyrka, G. E. R 13B1, University of Michigan Engineering Research Institute Project M777 (1953)

9 Mears, W. H., Stahl, R. F., Orfeo, S. R., Shair, R. C., Kells, L. F., Thompson, W., McCann, H. R 152a, Thermo- dynamic properties of halogenated ethanes and ethylenes Ind Eng Chem 47 7 (1955) 1449-1454

10 Redlich, O., Kwong, J. N. S. On the thermodynamics of solutions. V, An equation of state. Fugacities of gaseous solutions Chem Rev44 (1949) 233-244

11 Scare, G. Equilibrium constants from a modified Redlich- Kwong equation of state Chem Engng Sci 27 (1972) 1197- 1203

12 Scare, G. Rigorous and simplified procedures for determin- ing the pure-component parameters in the Redlich-Kwong- Soave equation of state Chem Engng Sci 35 (1980) 1725- 1729

13 Stoecker, W. F., Walukes, D. J. Conserving energy in domestic refrigerators through the use of refrigerant mixtures ASHRAE Transactions 87 1 (1981 ) 279-291

14 Vidal, J. Mixing rules and excess properties in cubic equa- tions of state Chem Engng Sci 33 (1978) 1269-1276

Appendix

Computations The Redlich-Kwong-Soave 11.12 equation of state

was chosen for calculation of the thermodynamic properties of R 1 3B1/R 1 52a because it allows calcu- lation of a variety of thermodynamic properties (en- thalpy, entropy, fugacity) of a mixture from a very limited set of experimental pure-component and mix- ture data. For each pure component, one needs the critical temperature and pressure, molecular weight, ideal-gas heat-capacity coefficients, and some pressure-temperature (p-T) data. Obviously, more ex- tensive p-T data will give better results. The information listed above is readily available for a wide variety of pure refrigerants. 3 For binary mixtures, at least one experimental pressure versus composition (p-x) data point must be specified. Ideally, however, data should be available to cover the ranges of composition, temperature, and pressure that are of interest.

The original Redlich-Kwong 10 equation of state was:

RT a P=v-b v(v+b) (1)

where p= pressure, T= temperature, v= molar volume, and a,b=composition and temperature-dependent parameters.

Soave's 1980 version, known as the Redlich- Kwong-Soave equation of state 12 has expanded the a parameter for pure components to include a temperature-dependent function with two constants, m and n. For each pure component, i:

ai=oliaci (2)

where

e,=l+ 1 ~ - kin,+.,7)

aci = 0.42748 R 2 Toi2/pci (4)

bi=bc~ =0.08664R Tc~/P~i (5)

The constants m and n must be found for each pure component by minimizing the rms deviation between experimental and calculated vapour pressures. The pure-component constants, m and n, were found to be 0.6209 and 0.1684, respectively, for R 13B1 and 0.9331 and 0.677 for R 152a. These constants were verified by using the RKS equation of state to generate saturation tables for each of the pure components.

For determining mixture properties, Soave's me- thod 11 also includes an experimentally determined interaction constant, kij, in the mixing rule for a. Thus, for a mixture of components i and j:

a,i= (1 -ki i ) (aiai)l/2 (6)

a= Z~,x~P~i (7) i j

Volume 7 Num6ro 3 Mai 1984 171

/'bi + bi'~ b = Z :Z :x~ , [ -~ - - } (8)

i i ", /

The interaction constant, k,i, is determined by minimiz- ing the sum of the errors between calculated and experimental bubble-point pressures. The R 13B1/152a interaction parameter, k~i=0.833, mini- mized the average absolute error over 22 data points to 1.16%.

Vapour density. The Redlich-Kwong-Soave equation of state 1~ has been shown to be a cubic polynomial in compressibility factor, Z:

Z3- Z2 + (A - B - B2)Z- AB =O (9)

where

z=p v RT'

Edmister 4

a P b P A - R2T2 and B = ~ (10)

has shown that in some regions of temperature and pressure, there are three real roots. The largest is vapour density (compressibility), the smallest is liquid density, and the third has no physical signifi- cance. Comparison of experimental and calculated vapour densities in a temperature range from 64.4 to 150.1"F (18.0 to 65.6°C) gave an average absolute error of 0.90% for a 70/30 weight percent R 13B1/152a mixture.

Liquid density. It is generally recognized that the equation of state does not have the desired accuracy for calculating liquid density. The liquid density was calculated using a simple mixing equation and the known densities 8.9 of the pure components (p = I bs ft-3 T=°R):

p ( 1 381 ) = 46.50 + 0.039808 T+ 0.82202 T 1/2 + 9.30627T~J3+ 9.22 x 10-7T 2 (11

#(1 52a) = 22.7862 + 36.8350T+ 33.9508T 2 -8 ,2530T (12)

Volume=weight fraction (1 3B1 )/#(1 3B1 ) +weight fraction (1 52a)/p(1 52a) (1 3)

Density = 1/volume (14)

Experimental and calculated liquid densities over a temperature range from - 60.0 to + 130.0°F ( - 51.1 to + 54.4°C) agreed to within 0.68% average absolute error.

Fugacity. The fugacity coefficient is derived from the general thermodynamic relationship:

P

f ( -P) l In ¢i = ~-~- (1 5)

0

This has been solved by Graboski and Daubert (1978) for the Redlich-Kwong-Soave equation:

bl tn ~ ,=~-(Z- 1 ) - I n (Z-B)

A- ~xja,j B - ~ [ 2 J a bi ] ln(1 + ~ ) (16)

Entha/py. Following the method given by Ed- mister. 5 isothermal change in enthalpy is derived from this general equation:

v~./ /#P'~ p'~dv+ A / (17)

v 1

Combination with the Redlich-Kwong-Soave equa- tion of state leads to the integrated form:

RT bRT In + Z - 1 (18)

where

. xixi(1 _k,i) " \1/2 ( a c , a ~ j ~ , , a = ~ Z -2- - - (~,~j +~j~, ) (19)

and

n i f c i T 2

Vapour or liquid computed by

/-

rn I T~, (20)

enthalpy. H. at temperature. T. is

AH H=ICpdT+(-R-~)TRT-fAH~ RT \ R / ] D D (21) T D

where Cp= mixture ideal-gas heat capacity at constant pressure; TD=reference datum temperature, -40°F ( -40°C); and subscript D=evaluate at datum tem- perature. Each table is referenced to 0.0 Btu Ib (0.0 kJ kg -1) liquid enthalpy at -40°F (-40°C).

Entropy. The derivation for entropy parallels the derivation for enthalpy. The thermodynamic equation

V2

$2-S~ = ~ dv v

v]

can be transformed to

/ Z - B , a' ( B ) (23) &S- In t~P- - )+b -R InR 1 + ~

Finally, vapour or liquid entropy at temperature, T, and pressure, p, is computed by

T R AS AS S=f~dT+ ( ( R - ) T - - ( ~ ) D ) ( 2 4 )

TD

where PD is the bubble-point pressure at -40°F ( -40°C).

172 International Journal of Refrigeration