- Email: [email protected]
A computer simulation program for mixed-refrigerant air conditioning
A computer simulation program for mixed-refrigerant air conditioning Geoffrey G. Haselden and J. Chen Department of Chemical Engineering, University of Leeds, Leeds LS2 9JT, UK Received I December 1992; revised 13 December 1993
A steady-state design simulation program for air-conditioning systems using binary refrigerant mixtures has been developed and tested. The Redlich-Kwong-Soave (RKS) equation of state is used to model refrigerant vapour mixtures and equilibrium behaviour. Specific heat data are used directly in the liquid phase. The pinch method is introduced to facilitate the refrigerant cycle simulation. The program is able to handle all mixtures in which liquid solution behaviour is near ideal. All the main models and subroutines are provided. For a given air-conditioning duty and specified mixed refrigerants, the program will predict system COP, the size of compressor, and the required areas of the evaporator and condenser, together with relevant temperatures, pressures and flows. The program can also be adapted for heat pump design and simulation. (Keywords:air conditioning;air conditioner;refrigerantmixtures;simulation;calculation)
Programme de simulation par ordinateur des conditionneurs d'air utilisant des mrlanges de frigorigSnes On a mis au point et essayb un programme de simulation en rbgime stable des installations de conditionnement d'air utilisant des mblanges binaires de frigorig~nes. On utilise l'bquation d'btat de Redlich-Kwong-Soave (RKS) pour modbliser les m~langes gazeux de frigorig~nes et leur ~quilibre. Les donn~es sur la chaleur massique sont utilisbes directement dans la phase liquide. On prbconise l'utilisation de la mbthode des kcarts pour faciliter la simulation des cycles de frigorigOnes. Le programme peut traiter tous les mblanges pour lesquels la solution liquide se comporte presque idbalement. On donne tousles modOles prineipaux et les sous-programmes. Pour une installation de conditionnement d'air donnbe et pour des mblanges de frioorig~nes spbeifiques, le programme sera en mesure de prbvoir le COP, la dimension du eompresseur et les emplacements appropribs de l~vaporateur et du condenseur, ainsi que les tempbratures, les pressions et les beoulements correspondants. Le programme peut bgalement ~tre adapt~ en rue de la conception et de la simulation des pompes ~ ehaleur.
(Mots cl6s: conditionnementd'air; conditionneurd'air; frigorig6ne;m61ange;simulation;calcul) In recent years there has been a growing interest in mixed refrigerants, either as near-azeotropic mixtures to replace CFCs, or as wide-boiling non-azeotropic refrigerant mixtures for power saving. A prerequisite of any such study is a simulation of the proposed system to indicate its probable power consumption and capital cost. Referring to a given air-conditioning application, the required outputs are: 1. the best refrigerant components and compositions; 2. the optimized system COP; 3. the size of the component parts. Although there are some papers related to mixedrefrigerant modelling I and simulation 2, no one has yet published a complete mixed-refrigerant simulation and design program. This paper presents an effective model for mixed-refrigerant air-conditioning systems.
Mixed-refrigerant cycle and programming The following assumptions were made. There is no refrigerant accumulation in the system, therefore the 0140- 7007/94/050343 -08 © 1994 Butterworth-Heinemann Ltd and IIR
composition of the refrigerant circulating in the system is constant and equal to z. There is no heat transfer and pressure drop in the connecting pipelines and fittings. The pressure drop in the condenser and evaporator is dominated by the two-phase region and pressure drops in the superheated vapour and subcooled liquid flow are neglected. With these assumptions a mixed-refrigerant air-conditioning cycle can be simplified as shown in Figure 1, with four main processes. The evaporation process 7-8-1 starts with two-phase refrigerant mixture at 7 and finishes at 8 or 1 with or without superheating. The compressor sucks refrigerant vapour at 1 and compresses it to 2. The condensation process 2--6 consists of desuperheating 2-4, condensing 4--5 and subcooling 5-6. The expansion process 6-7 is assumed to be isenthalpic. Referring to the heat transfer processes, the air and refrigerant temperature profiles in both evaporator and condenser can be expressed in a T - H diagram (Figure 2). In contrast to a pure refrigerant system, a mixed refrigerant has gliding temperatures during evaporation and condensation. To develop a robust computer program for a
Rev. Int. Froid 1994 Volume 17 Num6ro 5
343
Computer simulation of mixed-refrigerant air conditioning: G. G. Haselden and J. Chen Greek letters
Nomenclature
Specific heat (kJ kg-1 K-1) Condenser air flow volume (m 3 s-1) Condenser fin area (m 2) Tabulated constants in liquid density equation Refrigerant pressure drop (Pa) Temperature differences (K) Evaporator air flow volume (m 3 s-1) Evaporator fin area (m 2) Convergency limits Specific enthalpy (kJ kg-1) Mixed refrigerant interactive parameter Refrigerant mass flow rate (kg s-1) Refrigerant pressure (Pa) Condenser duty (kW) Total cooling duty (kW) Specific entropy (kJ kg-1 K-1) Compressor swept volume (m 3 h-1) Air temperature (K) Refrigerant temperature (K) Overall heat transfer coefficient (J m -2 K - l ) Refrigerant or air specific volume (m 3 kg- 1) Compressor work (kW) Humidity ratio Moisture condensation rate (kg s-1) Liquid refrigerant composition (mol mol- 1) Vapour refrigerant composition (mol mol- 1) Refrigerant input composition (mol mol- 1)
Cp
DAV DFA dl-d 4 AP At EAV EFA e
h kq P Qd
Qo s
SW t T U
W w
WD
x(i) y(i) z(i)
P
A 6 ........... ~ ~
.....
~
4
3 2
/
~ / /
/1
Refrigerant vapour weight fraction (kg kg- 1) Compressor isentropic efficiency Relative humidity Compressor volumetric efficiency Liquid density (kg m-3)
/']ise
2 P Subscripts
1,2,3,4,5,6,7,8 Cycle state shown in Figure 1 a Air c Critical properties ca Calculated value d Condenser dp Condenser pinch e Evaporator ep Evaporator pinch ex Experimental value i Refrigerant component i m Mixed refrigerant or molar property p Pinch point sb Subcooling sp Superheating tp Two phase 0 Reference state Superscripts
id v 1
Ideal term Vapour Liquid
In most cases there will be a pinch point between the air and refrigerant temperature profiles in both the evaporator and condenser, at which the temperature differences ITp-T,[ reach minimum values. For mixed refrigerants, with gliding temperatures in the two-phase region less than 10 °C, the pinch points are probably at T4 and t,4 in the condenser and at T7 and tat in the evaporator. To simulate the heat transfer process correctly the model should guarantee that the refrigerant and air temperature profiles never cross: that is, Atdp=T4-ta4 and Atcp=ta7-T 7 should always be positive. The proposed method is to use Atdp and Atep as control parameters and assign them positive values. The refrigerant temperature profiles are then dependent o n ta4 a n d taT.
r
h
The flowchart for the system is shown in Figure 3. To start the design program the required input data are:
Figure 1 Simplified cycle diagram for an air-conditioning system Figure 1 Diagramme simplifib du cycle pour une installation de conditionnement d'air
mixed-refrigerant system it is necessary to analyse the temperature profile diagrams give in Figure 2 and to define a characteristic that specifies heat exchange performance.
1. refrigerant critical data, reference values, RKS coefficients and mixture composition; 2. air flow and temperature specifications, compressor performance coefficients and heat transfer coefficients; 3. specifications for pinch-point values, subcooling in the condenser, superheating in the evaporator, and pressure-drop assumptions.
Computer simulation of mixed-refrigerant air conditioning." G. G. Haselden and J. Chen
, ~ - , - , air flow
I'
refrigerant flow
ta6~~erarum
T2
~/~ T4
,6 ~ ~ . T ' ;
outlet variables, P7, f17 are then available and T7 from the evaporator pinch. At the evaporator outlet end P8 =P~, and P~ can be evaluated from the specified pressure drop AP1, as P t = P 7 - A P 1 . The dew-point temperature Ta can be calculated by DEWP. The evaporator outlet temperature is evaluated by adding the specified superheat Atsp to T8 to give Ti = T8 + Atsp. The outlet enthalpy h~ and entropy sl are calculated by HSV. The refrigerant mass flow rate is then calculated from the total cooling duty Qe and the refrigerant enthalpy change:
~
,4. dn've in condenser . , . tal
"' "' ~ ' ~ * ~ "
"......
m=
"1"8
Q~ h 1- h 6
(4)
With the inlet state (T1, Pt) and outlet pressure P2 available, the subroutine COMP is called to calculate
H
Figure 2 Simplified air and refrigerant temperature profiles for a mixed-refrigerant air-conditioning system Figure 2 Profilsde tempbrature simplifibs de Fair et dufrigorig~ne pour
un conditionneur d'air utilisant des mblanges de frigorio~nes [Refrigerant property data input[
4,
[System duty and efficiency input I
Although it is possible to calculate refrigerant pressure drops, this requires tube and fin arrangement details of each heat exchanger as a precondition. However, it is sensible and practical to assume overall pressure drop values, say 0.5 bar, for both evaporator and condenser, at the initial design stage. Based on the required air flow in the evaporator the program calculates the required cooling duty Q~. A subroutine called AIRHV was developed based on ASHRAE moist-air models 3 to calculate ha, w,, v,, h, from known t, and ~ba. The total cooling Qe is calculated from the air enthalpy changes less the enthalpy of the condensed water:
Q,= EAV[hal-ha7--(Wal-- WaT)hw7 ]
4' [ Control parameters input l
4, Cooling duty: Call AIRHV(tal,¢l,hal,val,wl,hwl) Call AIRHV(ta7,c:7,ha7,va7,w7,hw7) Qe=(hal-ha7-[wl-w7].hw7).EAV/val
4' Estimate Cond.pinch:ta4=f(ta6,Oe) Con& refrig, pineh:T4=ta4+dtdp Call DEWP(z,T4.P2) P5 =P2- AP~ Call BUBT(z,P5,T5) T6=T5-At,~ Call HSL(z,T6,h6)
(1)
I Call EXPA(z,h6,T7,p7,:37)
Va1
To initiate the cycle calculation, the condenser-side air flow pinch value was estimated as a function of ta6 and Q~: t~4 = t~6 +
2.4Qe
Cp.DAV/1)a6
Evap.: PI=P7-AP1 Call DEWT(z,P1,T8) TI=T8+ At,p Call HSV(z,T1,PI,hl.sl m=Oe/(hl-h6)
(2)
4,
The refrigerant pinch-point temperature in the condenser is given by T4 = ta4 + Atdp
Comp: Call COMP(z,T1,P1,P2,~,~, T2, h2)J Call AIRHV(ta6,¢6,ha6,va6,wa6,hw6) 1 ta4p= f(ta6,DAV,m(h4-h6)) I
(3)
By calling subroutine DEWP the refrigerant dew-point pressure P4 in the condenser can be calculated. The condenser inlet pressure P2 is assigned the same value as P4. Similarly, at the condenser outlet end P6 = P5 and P5 can be evaluated from the specified pressure drop APz: P5 = P2 - AP2. The bubble point Ts of the refrigerant mixture is then calculated in subroutine BUBT. Adding subcooling At,b, the condenser outlet temperature T6 is given by T6=Ts-At,b. The outlet enthalpy h 6 is calculated in HSL. Through subroutine EXPA, the refrigerant is isenthalpically expanded from state 6 to 7, and expansion
Succesive iteration ta4=ta4p
/" N Output Calculation: W,COP,Pr,SW,EFA,WD,Qd,DFA,ta2
Figure 3 Flowchart for a mixed-refrigerant air-conditioning cycle Figure 3 Organigramme d'un cycle de conditionneur d'air utilisant des
mJlanges de frigorigknes
Computer simulation of mixed-refrigerant air conditioning: G. G. Haselden and J. Chen the refrigerant outlet temperature T2 and enthalpy h 2 due to compression. By now the whole refrigerant cycle calculation has been completed. The condenser air pinch-point temperature ta4p can be evaluated more accurately based on the calculated refrigerant mass flow m, refrigerant enthalpy change and condenser inlet air state: rh(h4 - - h6)
ta4 p = ta6 -k
CpaDA VIVa6
Compressor W = th(h2
2.
C O P - Qe w
3. P r -
-hi)
T"~0.6667
Pi=tOc[1"t-dl(l-Tc/# -t-d2tl-Tc) +d3(ii~e)+d4(i__~} ZX) 1"333] J
(9)
and the densities of the liquid mixtures were formed from 1
2 -
Pm
1
~
(10)
x,-
i= 1 Pi
as appropriate mixtures will have near-ideal behaviour. The enthalpy and entropy models for a mixed refrigerant are 2 hm= Z Yihld-Ah' (11) i=i
2 Sm= 2 i=1
Pz P1
2
Y,Sld - A s ' - R ~ ylx In Yi
(12)
i=1
3600rhv
hii d and siid can be integrated from ideal gas heat capacity
Qd=rh(h2-h6)
Cpd for each component 5. For each refrigerant the reference states were chosen as saturation conditions at one atmospheric pressure. Hence the values ofho and So can be taken from ASHRAE data 3. The departure functions Ah' and As' for vapour mixtures were integrated from the RKS equation. Liquid-mixture enthalpies and entropies were calculated directly from pure component values at the same temperature, assuming no heat or entropy of mixing:
Condenser 2. DFA = Asb + Atp + Asp Qd 3. ta2 = ta6 -~ CpaDAV/v,6 Evaporator
[Lo}
(13)
sl= Lxi[slo+ IrC~dTl
(14)
h'= 1.
//
and
4. S W - - -
1.
r ~ 0.355
(5)
The new condenser air pinch temperature t.4 p is then compared with ta4, the previous air pinch temperature. If Ita4p-ta41 > e another iteration will be required using ta4 = ta4p. When the pinch-point convergence is satisfied, the calculation is complete and the required parameters can be evaluated:
1.
be calculated from the smallest root of the cubic equation, this can lead to large errors. Instead, the densities of the pure liquid components were calculated as empirical functions of temperature 6 by
2. WD=
DAV
(w1-w7)
Va6
i=l
The simulation can then be repeated for different refrigerant blends or pinch values to optimize the chosen parameters. Modelling
To predict thermodynamic properties for mixed refrigerants the RKS equation 4'5 was chosen. The mixture equilibrium calculations can be carried out by solving the following simultaneous equations:
Tt=T v
(6)
p,=pv
(7)
Xitpli=yiq9 v
xi hi+ i=l
EFA = Atp+ Asp
i=
1,2
(8)
The fugacity coefficients are calculated using the RKS equation 5. Although the liquid density of a mixed refrigerant can
Cl dt
,Jror
C~ is a function of temperature, and can be obtained by curve-fitting ASHRAE data. For heat exchanger surface design a complete heat exchanger can be subdivided into subcooled, two-phase and superheated regions based on refrigerant state. Therefore the required heat transfer surface is the summation of these regions. For the condenser the required fin area is calculated by
V
DFA = ! mAh56" + rhAh45 q LU,6LMTDs6 U45LMTD,,s
mAhz4
]1
U24L---M~D24J(a" (15)
assuming that there are sufficient tube passes for the heat transfer to be essentially counter-current. Similarly, the evaporator fin area is calculated by f rhAhvs mAhsx 1 1 EFA=L uTs-~---M-~D78 + Us,LMT DslA
(16)
Computer simulation of mixed.refrigerant air conditioning." G. G. Haselden and d. Chen data input: P,y(i)] i
initialvalue: T=280K calculate: Pv(i)=f(T,Tc(i),Pc(i)) x(2)=Pv(1)y(2)/fPv(1)y(2)+Pv(2)y(1)) x(1)=l-x(2) Ps=x(1)Pv(1)+x(2)Pv(2)
,,@
new T from Secant method
call FUGAV(T,P,y(i),q~) [ x(i)=y(i)~i/ tpI ]
except that the second loop is controlled by the convergence of the calculated vapour composition ly(1) + y(2)- 11< e.
Subroutine BU BP BUBP (shown in Figure 5) is used to calculate bubble-point pressure P and vapour composition y(i) at a given T and x(i). It employs Raoult's law to obtain an initial estimate of P and y(/). Then accurate values of P and y(i) are achieved by fugacity balance. Subroutine D E W P DEWP is used to calculate the dew-point pressure P and liquid composition x(i) at known temperature T and vapour composition y(/), and is similar to BUBP.
new T from N-R method
Subroutine F L A S H
Figure 4 Flowchart for subroutine DEWT
Figure 4 Organigramme du sOUSoprogran~neDEWT
Here rh is the refrigerant mass flow rate, and Ahi~ is the enthalpy difference of the refrigerant mixture for the given region. (d and ~e are the factors for cross-flow. LMTDij is the local log mean temperature difference of a specific region. U 0 is the local overall heat transfer coefficient based on fin area, and can be calculated by summing the reciprocals of the local coefficients, in the normal way. By keeping the two-phase refrigerant velocity in the annular flow regime the effect of mass transfer resistance could be ignored.
Both phases of a non-azeotropic refrigerant mixture will change composition during equilibrium condensation or evaporation. This requires a subroutine to predict vapour composition y(i), liquid composition x(i) and vapour weight fraction for a given overall composition z(i) at pressure P and varying temperature T. The inside loop iterates the vapour fraction, and the outside loop calculates fugacity coefficients g0~', ~oI and mixture compositions. The adopted vapour fraction algorithms 7 are:
~2= 1 (K i - 1)zi
1 +/~m(K,- 1)
tim=tim+ X~=l[1
(17)
K,-1 ]2z, + ,Bin(K,- 1)A
where Main subroutines
Subroutine D E W T Subroutine DEWT (shown in Figure 4) is used to calculate dew-point temperature T and liquid composition x(i) from known pressure P and vapour composition y(i). Two iteration loops are used, the first to obtain reasonable starting values of T and x(i) and the second to refine these values. At the beginning of the first loop an initial value of 280 K is assigned to T. The vapour pressure of each component Pv(i) is evaluated at T. Then Raoult's law is applied to estimate the liquid composition and the corresponding total pressure Pc,. The first loop is controlled by the deviation in total pressure [Po, - P[ > e. In the second loop the liquid dew composition x(i) is calculated from the fugacity balance. The NewtonRaphson method is used to iterate the dew-point temperature T on the basis of the liquid composition summation Ix(I) + x(2)- I I < e.
Subroutine B U B T The function of the BUBT subroutine is to calculate the temperature T and vapour composition y(i) of the equilibrium bubble at a known P and liquid composition x(i). The flowchart for BUBT is similar to that for DEWT
K,-
(18)
data input: T,x(i)I
, initial estfmateP & y(i) : Pv(i)=f(T,Tc(i),Pc(i)) P=x(1)Pv(1) + x(2)Pv(2) y(i)=x(i)Pv(i)/P
¢
l
]call FUGAV(T,P,y(i),cpi )[
new P from N-R method
Figure 5 Flowchart for subroutine BUBP
Figure 5 Organigramme du sous-progran~'neBUBP
Computer simulation of mixed-refrigerant air conditioning." G. G. Haselden and J. Chen Data: z , h 6 , T 7 /
+ Guess P7:T7b=T7-1.5 I Call BUBP(z,T7b,P7)
q, •q Call FLASH(z,T7,P7 y,x,~)
(====~Call HSL(x,hca,sca)
I---] Call HSV(y,hca,sca) Call H S V ( y , T , P7,hv,sv) Call I--ISL(x,T, hl,sl) hca=fl hv+(l-/3) hl
q, ,-[ dh=hca-h6 J New P from NR Method
÷v
/ utput. P7,77 Figure 6 Flowchartfor subroutine EXPA Figure 6 Organigramme du sous-programme EXPA xl-
2i 1 + flm(Ki- 1)
Yi = Kixi
(19)
the mixed refrigerant R22/R142b at a mass fraction of 0.55 R22. The best fitted value of the interactive parameter for the RKS equation was k12 = -0.0036. A comparison of experimental and calculated dewpoint pressures showed an absolute average pressure deviation of 0.35%, with the largest deviation being less than 1%0. The dew-point temperature calculations were also checked. The average deviation was 0.05 °C over the pressure range of 1.6-21.6 bar, and the largest deviation was 0.2 °C. For the enthalpy of saturated vapour the absolute average deviation was 1.0 kJ kg-1, with the largest deviation being less than 1%0 over the dew-point pressure range from 1.6 bar to 21.6 bar. For the entropy of saturated vapour, the absolute average deviation was 0.0035 kJ kg-1 K - ~ and the largest deviation was less than 0.6% over the same dew-point pressure range. Similar comparisons were made for saturated liquid enthalpy and entropy. In the temperature range of - 20.6 to 71.1 °C the absolute enthalpy deviation was 0.22 kJ kg-1 with the largest deviation 0.5%. For the entropy of the saturated liquid, the absolute average deviation was 0.002 kJ kg- ~ K-~ with the largest deviation less than 1.5% at the lowest temperature. Comparisons were then made in the superheated vapour region. At a pressure of 2.07 bar, where the temperature rises from a dew-point of - 5 . 2 °C to superheated vapour of 18.3 °C, the calculated enthalpy showed an absolute average deviation of 0.2 kJ kg- ~ with a largest deviation of 0.12%; for entropy the deviations were 0.005 kJ kg -~ K -~ and 0.56% in the same superheated region. Another test was made in the pressure range of 2.1-20.7 bar, with between 10 and 35 °C
(20)
//~am: z,TI,PI,P2,Nis~ /
Subroutine E X P A
Subroutine EXPA (shown in Figure 6) is used to simulate the isenthalpic expansion process 6-7 (referring to Figure 1). It calculates expansion outlet pressure P and vapour fraction fl from known refrigerant composition z, inlet enthalpy h 6 and required outlet temperature T7. The calculation is initialized by guessing an evaporator inlet pressure P7 from subroutine BUBP. Subroutine C O M P
The subroutine C O M P (shown in Figure 7) is used to calculate compression outlet temperature T2 and enthalpy h2 for specified inlet state (P1, Tt) and outlet pressure P2. It contains two main iteration loops, with the first one being the isentropic process and the second one using isentropic efficiency r/is=to calculate the actual outlet state. The actual compression outlet enthalpy is then evaluated from the given isentropic efficiency value ~/ise. An initial outlet temperature T2 is estimated from rhs=. The Newton-Raphson method is used to modify T2 to reach convergency.
Call HSV(z,TI,Pl,hl,sl)
J
I •
~Call HSV(z,T3,P2,h3,s3) ds=sT-s3
1
I New T3 from ] N-R method I
h2=hl+(h3-hl)/rhj. Guess "i'9-: T2='I'3+('l'3-T1)/rh.=
"f Call ItSV(z,T2,P2,h2c,s2e) 1
dh=h2-h2c INew "I"2from .d [N-R method "
Program test
The calculation accuracy of the thermodynamic property program was checked with experimental results* using *Supplied by the Pennwalt Corporation of USA.
I
guess T3=f(k,Pr) . I
Figure 7 Flowchartfor subroutine COMP Figure 7 Oroanigramme du sous-prooramme COMP
I
Computer simulation of mixed-refrigerant air conditioning." G. G. Haselden and J. Chen
;
5.0
.
i
and evaporator fin areas of about 80 and 50 m 2 respectively, would have a C O P of approximately 3.1. F i g u r e 8 shows how the calculated system C O P varies with the overall composition of R22, with the condenser pinch value as a parameter• C O P values reach a maximum at composition z(R22)=0.5 at each condenser pinch. Using a condenser pinch of 6.5 °C it is seen that a maximum C O P of 4.2 is attainable, representing a power saving of 35% compared with the standard unit. Increasing the condenser pinch to 10 °C reduces the C O P to 3.77 and the power saving to 22%. The corresponding condenser fin areas required are presented in F i g u r e 9. For the 50% mixture and a condenser pinch of 6.5 °C the required fin area is 82 m 2, while for 10 °C the area is 59 m 2.
i
4.5
Jl
00
I
i \
(") 4•0
3.5
j
~
J
3.0 0.2
..~ I~
0.3
~to=6.,5'*C
i
At~'=lo*c
'1
{ I
i
0.6
0.7
0.8
0.5
0.4
Overall R22 c o m p o s i ~ : i o n ( m a s s ) Figure
g COP versus condenser temperature pinch: At,o=5.0°C;
At,b =O°C; At,p = 2.0°C
Figure 8 Ecart entre le COP et la tempbrature du condenseur: At,~=0°C; At,v= 2.0°C
4.6
~
•
;
4.2
i
!
!
n 0
II
~-al / "
250
i >-/ / " I !
3,t / •~
200
3•4
o_ h m rtO ~
-~-
/xtn~=10*C ~t~=l*C
I
'
', ~
i
150
3.01
I
100
,so o.1
0.2
o.4
o.s
o.6
!
I
~
i';
i;
' r-;~,.,_-,0oo ~ti "If" At.~ =5°C
i i \i i
!
!
-*- Et;,-z*c
!
i
"f
i
t
\
x ~ \
i
!
i i : *
I
t
!
0.2
0.3
0.4
t
i
I
0.5
0•6
0.7
i
!
T'~
!
,
i; 0.8
Overall R22 c o m p o s i t i o n ( m a s s )
' o.3
',
!
'
'
V--I
321-J
i
o.7
o.a
o.s
Figure lg COP versus evaporator temperature pinch: Atdp=6.5°C; At,b=2°C; At,p=0°C Figure 10 Ecart entre le COP et la tempbrature de l'bvaporateur: Atdp=6.5°C; At,b=2°C; At,p=0°C
Overall R22 c o m p o s i t i o n ( m a s s ) Figure 9 Required condenser fin area versus condenser pinch: At.p = 5.0°C; At,b = 2.0°C; At., = 0°C .'120
Figure 9 Ecart entre la zone des ailettes du condenseur et le condenseur : At,p = 5.0°C; At,b = 2.0°C; At,p = 0 ° C
200
i
180
E "~
superheat for each point. The superheated vapour enthalpy gave an absolute average deviation of 0.8 kJ k g - z and the largest deviation was 0.8 %. The superheated vapour entropy gave the absolute average and the largest deviation of 0•003 kJ k g - 1 K - ~ and 0.5%. On the basis of these tests it is expected that the calculated properties of this mixture will be within 1% over the range required. The overall model was then used to design a unit for the required duty of cooling 0.95 m a s - z of air from 26.7 °C to 13.3 °C in the evaporator, while the condenser takes in 1.25 m a s - z of air at a temperature of 35 °C. The superheat at compressor suction and the subcooling at condenser discharge are kept at zero and 2 °C. For this duty a standard unit operating on R22, with condenser
180
J
P I
t.= ,,:10c I
O 14.0 i-"
,5.top--5* C At.~=I*C
i I
LL 120 "~ 0 t-~
! ;
100
i i
!
8(
44:
2O
0.2
0.3 0•4 0.5 0.6 0.7 Overall R22 c o m p o s i t i o n ( m a s s )
o.e
Figure ll Evaporator fin area versus evaporator pinch: Atdp= 6.5°C; .Atsb= 2.0°C; Atsp= 0oc Figure 11 Ecart entre la zone des ailettes de l~vaporateur et l~vaporateur: Atdp= 6.5°C; At,b= 2.0°C; Atsp= 0°C
Computer simulation of mixed-refrigerant air conditioning: G. G. Haselden and J. Chen The corresponding results when the evaporator pinch is changed while the condenser pinch is held steady at 6.5 are given in Figures 10 and 11. An evaporator pinch of 5 °C again gives the optimum energy performance (COP = 4.2) and the required fin area is 80 m 2. Increasing the pinch to 10 °C reduces the COP to 3.57 and decreases the fin area to 46 m2. The optimum design must take into account both the cost of power and the appropriate capital charges, but the above figures suggest that significant overall savings are likely when mixed refrigerants are properly applied. The run time of the full simulation program on a typical PC (Dell 200) is about 5 s. Conclusion
It is possible to predict mixed-refrigerant thermodynamic properties with an accuracy of better than 1% by using the RKS equation for the vapour phase and direct combination methods for the liquid phase. The pinch-point method, together with other programming approaches presented in this paper, provides a robust and reliable simulation program for air-conditioning systems using binary refrigerant mixtures. The simulation program has been used in an airconditioning system optimization study a, and shows that
power savings of the order of 30% should be achieved by using similar increases in heat transfer area. The program can be readily adapted for use with heat pumps, and can handle other mixtures provided that the liquid solution behaviour is nearly ideal and condenser pressure does not come within 40% of the critical pressure of either component.
References Domanski, P., Didion, D.A. Simulationofa heat pump operating
with a non-azeotropic mixture ASHRAE Trans (1985) 91(2B) 250-259 Connan, H.A. A generated computer program for analysis of mixed refrigerantcyclelnt J Refrig (1984) 7(3) 167-172 ASHRAE Fundamentals Handbook (1989) 17 Kruse, H., Kuver, M., Quast, U., Schroeder, M. and Upmeier, B.
Theoretical and experimental investigations of advantageous refrigerant mixture applicationsASHRAE Trans (1985) 91(2B) Reid, R.C.,Prausaitz,J.M., Poling,B.E. The Properties of Gases and Liquids 4th edn, McGraw-Hill,New York (1987) MeLiaden,M.O.Thermodynamicpropertiesof CFC alternatives: a surveyof the availabledata Int J Refrio (1990) 13 149-162 Wains, S.M. Phase Equilibria in Chemical Engineering Butterworth (1985) Chen, J. Optimization of a vapour compression air conditioning/ heat pump system using mixed refrigerants PhD Thesis, Leeds University, UK (1992)
A steady-state design simulation program for air-conditioning systems using binary refrigerant mixtures has been developed and tested. The Redlich-Kwong-Soave (RKS) equation of state is used to model refrigerant vapour mixtures and equilibrium behaviour. Specific heat data are used directly in the liquid phase. The pinch method is introduced to facilitate the refrigerant cycle simulation. The program is able to handle all mixtures in which liquid solution behaviour is near ideal. All the main models and subroutines are provided. For a given air-conditioning duty and specified mixed refrigerants, the program will predict system COP, the size of compressor, and the required areas of the evaporator and condenser, together with relevant temperatures, pressures and flows. The program can also be adapted for heat pump design and simulation. (Keywords:air conditioning;air conditioner;refrigerantmixtures;simulation;calculation)
Programme de simulation par ordinateur des conditionneurs d'air utilisant des mrlanges de frigorigSnes On a mis au point et essayb un programme de simulation en rbgime stable des installations de conditionnement d'air utilisant des mblanges binaires de frigorig~nes. On utilise l'bquation d'btat de Redlich-Kwong-Soave (RKS) pour modbliser les m~langes gazeux de frigorig~nes et leur ~quilibre. Les donn~es sur la chaleur massique sont utilisbes directement dans la phase liquide. On prbconise l'utilisation de la mbthode des kcarts pour faciliter la simulation des cycles de frigorigOnes. Le programme peut traiter tous les mblanges pour lesquels la solution liquide se comporte presque idbalement. On donne tousles modOles prineipaux et les sous-programmes. Pour une installation de conditionnement d'air donnbe et pour des mblanges de frioorig~nes spbeifiques, le programme sera en mesure de prbvoir le COP, la dimension du eompresseur et les emplacements appropribs de l~vaporateur et du condenseur, ainsi que les tempbratures, les pressions et les beoulements correspondants. Le programme peut bgalement ~tre adapt~ en rue de la conception et de la simulation des pompes ~ ehaleur.
(Mots cl6s: conditionnementd'air; conditionneurd'air; frigorig6ne;m61ange;simulation;calcul) In recent years there has been a growing interest in mixed refrigerants, either as near-azeotropic mixtures to replace CFCs, or as wide-boiling non-azeotropic refrigerant mixtures for power saving. A prerequisite of any such study is a simulation of the proposed system to indicate its probable power consumption and capital cost. Referring to a given air-conditioning application, the required outputs are: 1. the best refrigerant components and compositions; 2. the optimized system COP; 3. the size of the component parts. Although there are some papers related to mixedrefrigerant modelling I and simulation 2, no one has yet published a complete mixed-refrigerant simulation and design program. This paper presents an effective model for mixed-refrigerant air-conditioning systems.
Mixed-refrigerant cycle and programming The following assumptions were made. There is no refrigerant accumulation in the system, therefore the 0140- 7007/94/050343 -08 © 1994 Butterworth-Heinemann Ltd and IIR
composition of the refrigerant circulating in the system is constant and equal to z. There is no heat transfer and pressure drop in the connecting pipelines and fittings. The pressure drop in the condenser and evaporator is dominated by the two-phase region and pressure drops in the superheated vapour and subcooled liquid flow are neglected. With these assumptions a mixed-refrigerant air-conditioning cycle can be simplified as shown in Figure 1, with four main processes. The evaporation process 7-8-1 starts with two-phase refrigerant mixture at 7 and finishes at 8 or 1 with or without superheating. The compressor sucks refrigerant vapour at 1 and compresses it to 2. The condensation process 2--6 consists of desuperheating 2-4, condensing 4--5 and subcooling 5-6. The expansion process 6-7 is assumed to be isenthalpic. Referring to the heat transfer processes, the air and refrigerant temperature profiles in both evaporator and condenser can be expressed in a T - H diagram (Figure 2). In contrast to a pure refrigerant system, a mixed refrigerant has gliding temperatures during evaporation and condensation. To develop a robust computer program for a
Rev. Int. Froid 1994 Volume 17 Num6ro 5
343
Computer simulation of mixed-refrigerant air conditioning: G. G. Haselden and J. Chen Greek letters
Nomenclature
Specific heat (kJ kg-1 K-1) Condenser air flow volume (m 3 s-1) Condenser fin area (m 2) Tabulated constants in liquid density equation Refrigerant pressure drop (Pa) Temperature differences (K) Evaporator air flow volume (m 3 s-1) Evaporator fin area (m 2) Convergency limits Specific enthalpy (kJ kg-1) Mixed refrigerant interactive parameter Refrigerant mass flow rate (kg s-1) Refrigerant pressure (Pa) Condenser duty (kW) Total cooling duty (kW) Specific entropy (kJ kg-1 K-1) Compressor swept volume (m 3 h-1) Air temperature (K) Refrigerant temperature (K) Overall heat transfer coefficient (J m -2 K - l ) Refrigerant or air specific volume (m 3 kg- 1) Compressor work (kW) Humidity ratio Moisture condensation rate (kg s-1) Liquid refrigerant composition (mol mol- 1) Vapour refrigerant composition (mol mol- 1) Refrigerant input composition (mol mol- 1)
Cp
DAV DFA dl-d 4 AP At EAV EFA e
h kq P Qd
Qo s
SW t T U
W w
WD
x(i) y(i) z(i)
P
A 6 ........... ~ ~
.....
~
4
3 2
/
~ / /
/1
Refrigerant vapour weight fraction (kg kg- 1) Compressor isentropic efficiency Relative humidity Compressor volumetric efficiency Liquid density (kg m-3)
/']ise
2 P Subscripts
1,2,3,4,5,6,7,8 Cycle state shown in Figure 1 a Air c Critical properties ca Calculated value d Condenser dp Condenser pinch e Evaporator ep Evaporator pinch ex Experimental value i Refrigerant component i m Mixed refrigerant or molar property p Pinch point sb Subcooling sp Superheating tp Two phase 0 Reference state Superscripts
id v 1
Ideal term Vapour Liquid
In most cases there will be a pinch point between the air and refrigerant temperature profiles in both the evaporator and condenser, at which the temperature differences ITp-T,[ reach minimum values. For mixed refrigerants, with gliding temperatures in the two-phase region less than 10 °C, the pinch points are probably at T4 and t,4 in the condenser and at T7 and tat in the evaporator. To simulate the heat transfer process correctly the model should guarantee that the refrigerant and air temperature profiles never cross: that is, Atdp=T4-ta4 and Atcp=ta7-T 7 should always be positive. The proposed method is to use Atdp and Atep as control parameters and assign them positive values. The refrigerant temperature profiles are then dependent o n ta4 a n d taT.
r
h
The flowchart for the system is shown in Figure 3. To start the design program the required input data are:
Figure 1 Simplified cycle diagram for an air-conditioning system Figure 1 Diagramme simplifib du cycle pour une installation de conditionnement d'air
mixed-refrigerant system it is necessary to analyse the temperature profile diagrams give in Figure 2 and to define a characteristic that specifies heat exchange performance.
1. refrigerant critical data, reference values, RKS coefficients and mixture composition; 2. air flow and temperature specifications, compressor performance coefficients and heat transfer coefficients; 3. specifications for pinch-point values, subcooling in the condenser, superheating in the evaporator, and pressure-drop assumptions.
Computer simulation of mixed-refrigerant air conditioning." G. G. Haselden and J. Chen
, ~ - , - , air flow
I'
refrigerant flow
ta6~~erarum
T2
~/~ T4
,6 ~ ~ . T ' ;
outlet variables, P7, f17 are then available and T7 from the evaporator pinch. At the evaporator outlet end P8 =P~, and P~ can be evaluated from the specified pressure drop AP1, as P t = P 7 - A P 1 . The dew-point temperature Ta can be calculated by DEWP. The evaporator outlet temperature is evaluated by adding the specified superheat Atsp to T8 to give Ti = T8 + Atsp. The outlet enthalpy h~ and entropy sl are calculated by HSV. The refrigerant mass flow rate is then calculated from the total cooling duty Qe and the refrigerant enthalpy change:
~
,4. dn've in condenser . , . tal
"' "' ~ ' ~ * ~ "
"......
m=
"1"8
Q~ h 1- h 6
(4)
With the inlet state (T1, Pt) and outlet pressure P2 available, the subroutine COMP is called to calculate
H
Figure 2 Simplified air and refrigerant temperature profiles for a mixed-refrigerant air-conditioning system Figure 2 Profilsde tempbrature simplifibs de Fair et dufrigorig~ne pour
un conditionneur d'air utilisant des mblanges de frigorio~nes [Refrigerant property data input[
4,
[System duty and efficiency input I
Although it is possible to calculate refrigerant pressure drops, this requires tube and fin arrangement details of each heat exchanger as a precondition. However, it is sensible and practical to assume overall pressure drop values, say 0.5 bar, for both evaporator and condenser, at the initial design stage. Based on the required air flow in the evaporator the program calculates the required cooling duty Q~. A subroutine called AIRHV was developed based on ASHRAE moist-air models 3 to calculate ha, w,, v,, h, from known t, and ~ba. The total cooling Qe is calculated from the air enthalpy changes less the enthalpy of the condensed water:
Q,= EAV[hal-ha7--(Wal-- WaT)hw7 ]
4' [ Control parameters input l
4, Cooling duty: Call AIRHV(tal,¢l,hal,val,wl,hwl) Call AIRHV(ta7,c:7,ha7,va7,w7,hw7) Qe=(hal-ha7-[wl-w7].hw7).EAV/val
4' Estimate Cond.pinch:ta4=f(ta6,Oe) Con& refrig, pineh:T4=ta4+dtdp Call DEWP(z,T4.P2) P5 =P2- AP~ Call BUBT(z,P5,T5) T6=T5-At,~ Call HSL(z,T6,h6)
(1)
I Call EXPA(z,h6,T7,p7,:37)
Va1
To initiate the cycle calculation, the condenser-side air flow pinch value was estimated as a function of ta6 and Q~: t~4 = t~6 +
2.4Qe
Cp.DAV/1)a6
Evap.: PI=P7-AP1 Call DEWT(z,P1,T8) TI=T8+ At,p Call HSV(z,T1,PI,hl.sl m=Oe/(hl-h6)
(2)
4,
The refrigerant pinch-point temperature in the condenser is given by T4 = ta4 + Atdp
Comp: Call COMP(z,T1,P1,P2,~,~, T2, h2)J Call AIRHV(ta6,¢6,ha6,va6,wa6,hw6) 1 ta4p= f(ta6,DAV,m(h4-h6)) I
(3)
By calling subroutine DEWP the refrigerant dew-point pressure P4 in the condenser can be calculated. The condenser inlet pressure P2 is assigned the same value as P4. Similarly, at the condenser outlet end P6 = P5 and P5 can be evaluated from the specified pressure drop APz: P5 = P2 - AP2. The bubble point Ts of the refrigerant mixture is then calculated in subroutine BUBT. Adding subcooling At,b, the condenser outlet temperature T6 is given by T6=Ts-At,b. The outlet enthalpy h 6 is calculated in HSL. Through subroutine EXPA, the refrigerant is isenthalpically expanded from state 6 to 7, and expansion
Succesive iteration ta4=ta4p
/" N Output Calculation: W,COP,Pr,SW,EFA,WD,Qd,DFA,ta2
Figure 3 Flowchart for a mixed-refrigerant air-conditioning cycle Figure 3 Organigramme d'un cycle de conditionneur d'air utilisant des
mJlanges de frigorigknes
Computer simulation of mixed-refrigerant air conditioning: G. G. Haselden and J. Chen the refrigerant outlet temperature T2 and enthalpy h 2 due to compression. By now the whole refrigerant cycle calculation has been completed. The condenser air pinch-point temperature ta4p can be evaluated more accurately based on the calculated refrigerant mass flow m, refrigerant enthalpy change and condenser inlet air state: rh(h4 - - h6)
ta4 p = ta6 -k
CpaDA VIVa6
Compressor W = th(h2
2.
C O P - Qe w
3. P r -
-hi)
T"~0.6667
Pi=tOc[1"t-dl(l-Tc/# -t-d2tl-Tc) +d3(ii~e)+d4(i__~} ZX) 1"333] J
(9)
and the densities of the liquid mixtures were formed from 1
2 -
Pm
1
~
(10)
x,-
i= 1 Pi
as appropriate mixtures will have near-ideal behaviour. The enthalpy and entropy models for a mixed refrigerant are 2 hm= Z Yihld-Ah' (11) i=i
2 Sm= 2 i=1
Pz P1
2
Y,Sld - A s ' - R ~ ylx In Yi
(12)
i=1
3600rhv
hii d and siid can be integrated from ideal gas heat capacity
Qd=rh(h2-h6)
Cpd for each component 5. For each refrigerant the reference states were chosen as saturation conditions at one atmospheric pressure. Hence the values ofho and So can be taken from ASHRAE data 3. The departure functions Ah' and As' for vapour mixtures were integrated from the RKS equation. Liquid-mixture enthalpies and entropies were calculated directly from pure component values at the same temperature, assuming no heat or entropy of mixing:
Condenser 2. DFA = Asb + Atp + Asp Qd 3. ta2 = ta6 -~ CpaDAV/v,6 Evaporator
[Lo}
(13)
sl= Lxi[slo+ IrC~dTl
(14)
h'= 1.
//
and
4. S W - - -
1.
r ~ 0.355
(5)
The new condenser air pinch temperature t.4 p is then compared with ta4, the previous air pinch temperature. If Ita4p-ta41 > e another iteration will be required using ta4 = ta4p. When the pinch-point convergence is satisfied, the calculation is complete and the required parameters can be evaluated:
1.
be calculated from the smallest root of the cubic equation, this can lead to large errors. Instead, the densities of the pure liquid components were calculated as empirical functions of temperature 6 by
2. WD=
DAV
(w1-w7)
Va6
i=l
The simulation can then be repeated for different refrigerant blends or pinch values to optimize the chosen parameters. Modelling
To predict thermodynamic properties for mixed refrigerants the RKS equation 4'5 was chosen. The mixture equilibrium calculations can be carried out by solving the following simultaneous equations:
Tt=T v
(6)
p,=pv
(7)
Xitpli=yiq9 v
xi hi+ i=l
EFA = Atp+ Asp
i=
1,2
(8)
The fugacity coefficients are calculated using the RKS equation 5. Although the liquid density of a mixed refrigerant can
Cl dt
,Jror
C~ is a function of temperature, and can be obtained by curve-fitting ASHRAE data. For heat exchanger surface design a complete heat exchanger can be subdivided into subcooled, two-phase and superheated regions based on refrigerant state. Therefore the required heat transfer surface is the summation of these regions. For the condenser the required fin area is calculated by
V
DFA = ! mAh56" + rhAh45 q LU,6LMTDs6 U45LMTD,,s
mAhz4
]1
U24L---M~D24J(a" (15)
assuming that there are sufficient tube passes for the heat transfer to be essentially counter-current. Similarly, the evaporator fin area is calculated by f rhAhvs mAhsx 1 1 EFA=L uTs-~---M-~D78 + Us,LMT DslA
(16)
Computer simulation of mixed.refrigerant air conditioning." G. G. Haselden and d. Chen data input: P,y(i)] i
initialvalue: T=280K calculate: Pv(i)=f(T,Tc(i),Pc(i)) x(2)=Pv(1)y(2)/fPv(1)y(2)+Pv(2)y(1)) x(1)=l-x(2) Ps=x(1)Pv(1)+x(2)Pv(2)
,,@
new T from Secant method
call FUGAV(T,P,y(i),q~) [ x(i)=y(i)~i/ tpI ]
except that the second loop is controlled by the convergence of the calculated vapour composition ly(1) + y(2)- 11< e.
Subroutine BU BP BUBP (shown in Figure 5) is used to calculate bubble-point pressure P and vapour composition y(i) at a given T and x(i). It employs Raoult's law to obtain an initial estimate of P and y(/). Then accurate values of P and y(i) are achieved by fugacity balance. Subroutine D E W P DEWP is used to calculate the dew-point pressure P and liquid composition x(i) at known temperature T and vapour composition y(/), and is similar to BUBP.
new T from N-R method
Subroutine F L A S H
Figure 4 Flowchart for subroutine DEWT
Figure 4 Organigramme du sOUSoprogran~neDEWT
Here rh is the refrigerant mass flow rate, and Ahi~ is the enthalpy difference of the refrigerant mixture for the given region. (d and ~e are the factors for cross-flow. LMTDij is the local log mean temperature difference of a specific region. U 0 is the local overall heat transfer coefficient based on fin area, and can be calculated by summing the reciprocals of the local coefficients, in the normal way. By keeping the two-phase refrigerant velocity in the annular flow regime the effect of mass transfer resistance could be ignored.
Both phases of a non-azeotropic refrigerant mixture will change composition during equilibrium condensation or evaporation. This requires a subroutine to predict vapour composition y(i), liquid composition x(i) and vapour weight fraction for a given overall composition z(i) at pressure P and varying temperature T. The inside loop iterates the vapour fraction, and the outside loop calculates fugacity coefficients g0~', ~oI and mixture compositions. The adopted vapour fraction algorithms 7 are:
~2= 1 (K i - 1)zi
1 +/~m(K,- 1)
tim=tim+ X~=l[1
(17)
K,-1 ]2z, + ,Bin(K,- 1)A
where Main subroutines
Subroutine D E W T Subroutine DEWT (shown in Figure 4) is used to calculate dew-point temperature T and liquid composition x(i) from known pressure P and vapour composition y(i). Two iteration loops are used, the first to obtain reasonable starting values of T and x(i) and the second to refine these values. At the beginning of the first loop an initial value of 280 K is assigned to T. The vapour pressure of each component Pv(i) is evaluated at T. Then Raoult's law is applied to estimate the liquid composition and the corresponding total pressure Pc,. The first loop is controlled by the deviation in total pressure [Po, - P[ > e. In the second loop the liquid dew composition x(i) is calculated from the fugacity balance. The NewtonRaphson method is used to iterate the dew-point temperature T on the basis of the liquid composition summation Ix(I) + x(2)- I I < e.
Subroutine B U B T The function of the BUBT subroutine is to calculate the temperature T and vapour composition y(i) of the equilibrium bubble at a known P and liquid composition x(i). The flowchart for BUBT is similar to that for DEWT
K,-
(18)
data input: T,x(i)I
, initial estfmateP & y(i) : Pv(i)=f(T,Tc(i),Pc(i)) P=x(1)Pv(1) + x(2)Pv(2) y(i)=x(i)Pv(i)/P
¢
l
]call FUGAV(T,P,y(i),cpi )[
new P from N-R method
Figure 5 Flowchart for subroutine BUBP
Figure 5 Organigramme du sous-progran~'neBUBP
Computer simulation of mixed-refrigerant air conditioning." G. G. Haselden and J. Chen Data: z , h 6 , T 7 /
+ Guess P7:T7b=T7-1.5 I Call BUBP(z,T7b,P7)
q, •q Call FLASH(z,T7,P7 y,x,~)
(====~Call HSL(x,hca,sca)
I---] Call HSV(y,hca,sca) Call H S V ( y , T , P7,hv,sv) Call I--ISL(x,T, hl,sl) hca=fl hv+(l-/3) hl
q, ,-[ dh=hca-h6 J New P from NR Method
÷v
/ utput. P7,77 Figure 6 Flowchartfor subroutine EXPA Figure 6 Organigramme du sous-programme EXPA xl-
2i 1 + flm(Ki- 1)
Yi = Kixi
(19)
the mixed refrigerant R22/R142b at a mass fraction of 0.55 R22. The best fitted value of the interactive parameter for the RKS equation was k12 = -0.0036. A comparison of experimental and calculated dewpoint pressures showed an absolute average pressure deviation of 0.35%, with the largest deviation being less than 1%0. The dew-point temperature calculations were also checked. The average deviation was 0.05 °C over the pressure range of 1.6-21.6 bar, and the largest deviation was 0.2 °C. For the enthalpy of saturated vapour the absolute average deviation was 1.0 kJ kg-1, with the largest deviation being less than 1%0 over the dew-point pressure range from 1.6 bar to 21.6 bar. For the entropy of saturated vapour, the absolute average deviation was 0.0035 kJ kg-1 K - ~ and the largest deviation was less than 0.6% over the same dew-point pressure range. Similar comparisons were made for saturated liquid enthalpy and entropy. In the temperature range of - 20.6 to 71.1 °C the absolute enthalpy deviation was 0.22 kJ kg-1 with the largest deviation 0.5%. For the entropy of the saturated liquid, the absolute average deviation was 0.002 kJ kg- ~ K-~ with the largest deviation less than 1.5% at the lowest temperature. Comparisons were then made in the superheated vapour region. At a pressure of 2.07 bar, where the temperature rises from a dew-point of - 5 . 2 °C to superheated vapour of 18.3 °C, the calculated enthalpy showed an absolute average deviation of 0.2 kJ kg- ~ with a largest deviation of 0.12%; for entropy the deviations were 0.005 kJ kg -~ K -~ and 0.56% in the same superheated region. Another test was made in the pressure range of 2.1-20.7 bar, with between 10 and 35 °C
(20)
//~am: z,TI,PI,P2,Nis~ /
Subroutine E X P A
Subroutine EXPA (shown in Figure 6) is used to simulate the isenthalpic expansion process 6-7 (referring to Figure 1). It calculates expansion outlet pressure P and vapour fraction fl from known refrigerant composition z, inlet enthalpy h 6 and required outlet temperature T7. The calculation is initialized by guessing an evaporator inlet pressure P7 from subroutine BUBP. Subroutine C O M P
The subroutine C O M P (shown in Figure 7) is used to calculate compression outlet temperature T2 and enthalpy h2 for specified inlet state (P1, Tt) and outlet pressure P2. It contains two main iteration loops, with the first one being the isentropic process and the second one using isentropic efficiency r/is=to calculate the actual outlet state. The actual compression outlet enthalpy is then evaluated from the given isentropic efficiency value ~/ise. An initial outlet temperature T2 is estimated from rhs=. The Newton-Raphson method is used to modify T2 to reach convergency.
Call HSV(z,TI,Pl,hl,sl)
J
I •
~Call HSV(z,T3,P2,h3,s3) ds=sT-s3
1
I New T3 from ] N-R method I
h2=hl+(h3-hl)/rhj. Guess "i'9-: T2='I'3+('l'3-T1)/rh.=
"f Call ItSV(z,T2,P2,h2c,s2e) 1
dh=h2-h2c INew "I"2from .d [N-R method "
Program test
The calculation accuracy of the thermodynamic property program was checked with experimental results* using *Supplied by the Pennwalt Corporation of USA.
I
guess T3=f(k,Pr) . I
Figure 7 Flowchartfor subroutine COMP Figure 7 Oroanigramme du sous-prooramme COMP
I
Computer simulation of mixed-refrigerant air conditioning." G. G. Haselden and J. Chen
;
5.0
.
i
and evaporator fin areas of about 80 and 50 m 2 respectively, would have a C O P of approximately 3.1. F i g u r e 8 shows how the calculated system C O P varies with the overall composition of R22, with the condenser pinch value as a parameter• C O P values reach a maximum at composition z(R22)=0.5 at each condenser pinch. Using a condenser pinch of 6.5 °C it is seen that a maximum C O P of 4.2 is attainable, representing a power saving of 35% compared with the standard unit. Increasing the condenser pinch to 10 °C reduces the C O P to 3.77 and the power saving to 22%. The corresponding condenser fin areas required are presented in F i g u r e 9. For the 50% mixture and a condenser pinch of 6.5 °C the required fin area is 82 m 2, while for 10 °C the area is 59 m 2.
i
4.5
Jl
00
I
i \
(") 4•0
3.5
j
~
J
3.0 0.2
..~ I~
0.3
~to=6.,5'*C
i
At~'=lo*c
'1
{ I
i
0.6
0.7
0.8
0.5
0.4
Overall R22 c o m p o s i ~ : i o n ( m a s s ) Figure
g COP versus condenser temperature pinch: At,o=5.0°C;
At,b =O°C; At,p = 2.0°C
Figure 8 Ecart entre le COP et la tempbrature du condenseur: At,~=0°C; At,v= 2.0°C
4.6
~
•
;
4.2
i
!
!
n 0
II
~-al / "
250
i >-/ / " I !
3,t / •~
200
3•4
o_ h m rtO ~
-~-
/xtn~=10*C ~t~=l*C
I
'
', ~
i
150
3.01
I
100
,so o.1
0.2
o.4
o.s
o.6
!
I
~
i';
i;
' r-;~,.,_-,0oo ~ti "If" At.~ =5°C
i i \i i
!
!
-*- Et;,-z*c
!
i
"f
i
t
\
x ~ \
i
!
i i : *
I
t
!
0.2
0.3
0.4
t
i
I
0.5
0•6
0.7
i
!
T'~
!
,
i; 0.8
Overall R22 c o m p o s i t i o n ( m a s s )
' o.3
',
!
'
'
V--I
321-J
i
o.7
o.a
o.s
Figure lg COP versus evaporator temperature pinch: Atdp=6.5°C; At,b=2°C; At,p=0°C Figure 10 Ecart entre le COP et la tempbrature de l'bvaporateur: Atdp=6.5°C; At,b=2°C; At,p=0°C
Overall R22 c o m p o s i t i o n ( m a s s ) Figure 9 Required condenser fin area versus condenser pinch: At.p = 5.0°C; At,b = 2.0°C; At., = 0°C .'120
Figure 9 Ecart entre la zone des ailettes du condenseur et le condenseur : At,p = 5.0°C; At,b = 2.0°C; At,p = 0 ° C
200
i
180
E "~
superheat for each point. The superheated vapour enthalpy gave an absolute average deviation of 0.8 kJ k g - z and the largest deviation was 0.8 %. The superheated vapour entropy gave the absolute average and the largest deviation of 0•003 kJ k g - 1 K - ~ and 0.5%. On the basis of these tests it is expected that the calculated properties of this mixture will be within 1% over the range required. The overall model was then used to design a unit for the required duty of cooling 0.95 m a s - z of air from 26.7 °C to 13.3 °C in the evaporator, while the condenser takes in 1.25 m a s - z of air at a temperature of 35 °C. The superheat at compressor suction and the subcooling at condenser discharge are kept at zero and 2 °C. For this duty a standard unit operating on R22, with condenser
180
J
P I
t.= ,,:10c I
O 14.0 i-"
,5.top--5* C At.~=I*C
i I
LL 120 "~ 0 t-~
! ;
100
i i
!
8(
44:
2O
0.2
0.3 0•4 0.5 0.6 0.7 Overall R22 c o m p o s i t i o n ( m a s s )
o.e
Figure ll Evaporator fin area versus evaporator pinch: Atdp= 6.5°C; .Atsb= 2.0°C; Atsp= 0oc Figure 11 Ecart entre la zone des ailettes de l~vaporateur et l~vaporateur: Atdp= 6.5°C; At,b= 2.0°C; Atsp= 0°C
Computer simulation of mixed-refrigerant air conditioning: G. G. Haselden and J. Chen The corresponding results when the evaporator pinch is changed while the condenser pinch is held steady at 6.5 are given in Figures 10 and 11. An evaporator pinch of 5 °C again gives the optimum energy performance (COP = 4.2) and the required fin area is 80 m 2. Increasing the pinch to 10 °C reduces the COP to 3.57 and decreases the fin area to 46 m2. The optimum design must take into account both the cost of power and the appropriate capital charges, but the above figures suggest that significant overall savings are likely when mixed refrigerants are properly applied. The run time of the full simulation program on a typical PC (Dell 200) is about 5 s. Conclusion
It is possible to predict mixed-refrigerant thermodynamic properties with an accuracy of better than 1% by using the RKS equation for the vapour phase and direct combination methods for the liquid phase. The pinch-point method, together with other programming approaches presented in this paper, provides a robust and reliable simulation program for air-conditioning systems using binary refrigerant mixtures. The simulation program has been used in an airconditioning system optimization study a, and shows that
power savings of the order of 30% should be achieved by using similar increases in heat transfer area. The program can be readily adapted for use with heat pumps, and can handle other mixtures provided that the liquid solution behaviour is nearly ideal and condenser pressure does not come within 40% of the critical pressure of either component.
References Domanski, P., Didion, D.A. Simulationofa heat pump operating
with a non-azeotropic mixture ASHRAE Trans (1985) 91(2B) 250-259 Connan, H.A. A generated computer program for analysis of mixed refrigerantcyclelnt J Refrig (1984) 7(3) 167-172 ASHRAE Fundamentals Handbook (1989) 17 Kruse, H., Kuver, M., Quast, U., Schroeder, M. and Upmeier, B.
Theoretical and experimental investigations of advantageous refrigerant mixture applicationsASHRAE Trans (1985) 91(2B) Reid, R.C.,Prausaitz,J.M., Poling,B.E. The Properties of Gases and Liquids 4th edn, McGraw-Hill,New York (1987) MeLiaden,M.O.Thermodynamicpropertiesof CFC alternatives: a surveyof the availabledata Int J Refrio (1990) 13 149-162 Wains, S.M. Phase Equilibria in Chemical Engineering Butterworth (1985) Chen, J. Optimization of a vapour compression air conditioning/ heat pump system using mixed refrigerants PhD Thesis, Leeds University, UK (1992)