Modeling, simulation and optimization of a vapor compression refrigeration system dynamic and steady state response

Applied Energy 158 (2015) 540–555

Contents lists available at ScienceDirect

Applied Energy journal homepage: www.elsevier.com/locate/apenergy

Modeling, simulation and optimization of a vapor compression refrigeration system dynamic and steady state response T.K. Nunes a, J.V.C. Vargas a,⇑, J.C. Ordonez b, D. Shah b, L.C.S. Martinho a a b

Department of Mechanical Engineering, Graduate Program in Mechanical Engineering, PGMEC, Federal University of Paraná, UFPR, CP 19011, CEP: 81531-980 Curitiba, PR, Brazil Department of Mechanical Engineering, Energy Sustainability Center and Center for Advanced Power Systems, Florida State University Tallahassee, FL, USA

h i g h l i g h t s
We introduce a transient vapor compression refrigeration (VCR) system model.
Dimensionless groups are identified, and results presented in normalized charts.
Model adjustment and experimental validation are conducted.
A system performance comparison is done for refrigerants R12, R134a, and R1234yf.
Optimal heat transfer area allocation is found for maximum system performance.

a r t i c l e

i n f o

Article history: Received 5 February 2015 Received in revised form 13 August 2015 Accepted 21 August 2015

Keywords: Refrigeration Optimal area allocation Heat exchanger Total heat transfer area

a b s t r a c t This paper introduces a dimensionless simplified mathematical model of a vapor compression refrigeration (VCR) system, in order to optimize the system dynamic response. The model combines principles of thermodynamics, heat and mass transfer applied to the system components with empirical correlations, assigning thermodynamic control volumes to each component, which yield a system of ordinary differential equations with respect to time that is integrated explicitly and accurately with low computational time. Appropriate dimensionless groups are identified, and the results are presented in the form of normalized charts for general application to similar systems. Experiments were conducted using an industrial chiller that was instrumented for real time data acquisition in order to determine model adjustment parameters through the solution of the inverse problem of parameters estimation (IPPE) for a 300 W thermal load. The model was then experimentally validated using the adjusted parameters by direct comparison of simulation results to the system measured thermal response for a 1200 W thermal load with good agreement. After experimental validation, the study addressed the optimization of the heat exchangers heat transfer area inventory for minimum pull down time. A system performance comparison is conducted for three refrigerant fluids: (i) a banned refrigerant (R12); (ii) its original ozone depletion harmless substitute, but with a high global warming potential, GWP (R134a), and (iii) one of the current substitutes of R134a, i.e., R1234yf. The dynamic results show that, for a system originally designed for R12, substitution with R1234yf depicts a closer performance to R12 than R134a. The optimal configuration that leads to steady state maximum system first (or coefficient of performance, COP), and second law efficiencies is also pursued. The normalized results for refrigerants R12, R134a, and R1234yf show that an optimal heat transfer area distribution in both evaporator and condenser, represented by an evaporator to total system heat exchanger area ratio x4;opt ffi 0:55, leads the system to minimum pull down time and maximum system 2nd law efficiency, whereas x4;opt ffi 0:4 to maximum COP, when the heat exchangers global heat transfer coefficients are of the same magnitude. The difference in the optimum location when the objective function is the pull-down time (and second law efficiency) and COP is due to the fact that the first law analysis does not fully capture the thermodynamic losses due to changes in heat exchangers’ areas, which stresses the importance of a system second law assessment for realistic results. However, changes in the optima location are observed when the ratio of heat exchangers global heat transfer coefficients departs from 1. The obtained maxima and minima are
 sharp f min  100, in the searching interval ð0:1 6 x4 6 0:9Þ, showing up to a 189.6% and 93.37% variation, f max f max=min

⇑ Corresponding author. Tel.: +55 41 3361 3307; fax: +55 41 3361 3129. E-mail address: [email protected] (J.V.C. Vargas). http://dx.doi.org/10.1016/j.apenergy.2015.08.098 0306-2619/Ó 2015 Elsevier Ltd. All rights reserved.

T.K. Nunes et al. / Applied Energy 158 (2015) 540–555

541

Nomenclature area, m2 dimensionless set point of expansion valve opening area corresponding to DT dsh a constant for a gas in the van der Waals equation, kPa m6 kmol2 a1 ; . . . ; a5 polynomial coefficients, Eqs (14) and (15) b constant for a gas in the van der Waals equation, m3 kmol1 b1 ; b2 polynomial coefficients, Eq. (39) c1 ; c2 ; c3 polynomial coefficients, Eq. (25) c specific heat, J kg1 K1 c0 dead volume ratio in the compressor Cv compressor valve coefficient CV control volume E energy, J G thermostatic valve adjustment constant, m2 K1 h specific enthalpy, J kg1 M molecular weight, kg kmol1 m mass, kg _ m mass flow rate, kg s1 n compressor polytropic coefficient p pressure, N m2 Q_ heat transfer rate, W R reliability coefficient in polynomial fitting R universal gas constant, 8.314 kJ kmol1 K1 rps compressor revolutions per second s evaporator to condenser global heat transfer ratio, Eq. (3) t time, s T temperature, K U global heat transfer coefficient, W m2 K1 v specific volume, m3 kg1 v0 specific volume of the refrigerant at ðp0 ; T 0 Þ Vc compressor volumetric displacement, m3 _ W power, W x heat exchanger area fraction, Eq. (4) ~ x vector of dimensionless system variables y space-averaged quality z dimensionless wall thermal conductance, Eq. (7) A e set A

Greek symbols D variation, difference g efficiency gv compressor volumetric efficiency l dimensionless conversion factor, Eq (44) n dimensionless conversion factor, Eq. (41) s dimensionless temperature, Eq. (5) / dimensionless compressor speed w dimensionless heat capacity rate, Eq. (5) Subscripts a air a, b, c, d points in the VCR system cycle, Fig. 1 cp compressor crit critical point C Carnot dsh desired degree of superheat GS global system heat heating mode l liquid p constant pressure r refrigerant sat saturation sc subcooled liquid zone in the condenser sh superheated vapor set set point v constant volume; vapor w wall 0 exterior ambient; reference conditions, initial conditions 1, 2, . . . , 7 control volume number 2p two phase zone (refrigerant liquid and vapor mixture) I first law of thermodynamics II second law of thermodynamics Superscript  dimensionless variable

in which f is either the calculated system pull down time or second law efficiency, respectively, which was observed with refrigerant R1234yf, which points out their importance for actual HVAC–R (heating, ventilation, air conditioning and refrigeration) vapor compression systems design. Ó 2015 Elsevier Ltd. All rights reserved.

1. Introduction The statistics on the total electrical energy consumption in the United States of America, USA, in 2009 [1] show that the residential, commercial and industrial sectors were responsible for 38%, 36%, and 26% of that total, respectively. The heating, ventilation, and air conditioning (HVAC) and refrigeration (R) sectors are responsible for 42% and 7% of the electrical energy consumption of the residential sector in 2008, respectively [2]. For the commercial sector, HVAC and R were responsible for 38% and 7% of the sector’s electrical energy consumption in 2008, respectively [3]. Facility heating, ventilating, and air conditioning (HVAC), and refrigeration (R) systems account for about 9% of all industrial electricity use, which somehow varies according to the type of industry [4]. Using such data, it is possible to calculate that the HVAC–R

systems electrical energy consumption, in the residential, commercial and industrial sectors, correspond to 18.62%, 16.20%, and 2.34% of the total electrical energy consumption in the USA, respectively, which adds up to 37.16% of the total electrical energy consumption in the country. Hence, it is demonstrated that HVAC–R systems affect considerably the USA energetic matrix, which could be reasonably extrapolated to the rest of the world, accounting for regional particular differences. As a result, any technical or scientific action that aims at the energy consumption reduction of HVAC–R systems will greatly contribute to the search of solutions to the current growing world energetic demand. Additionally, HVAC–R systems are widely used for the thermal management of large manmade devices, which comprise the so called field of systems engineering. Modeling and simulation are tools that are currently used for such complex design [5]. Systems

542

T.K. Nunes et al. / Applied Energy 158 (2015) 540–555

engineering includes commercial or military ships, airplanes, oil refinery, power plants, and similar processes. For example, allelectric ships present a series of thermal management challenges imposed by the energy levels and energy dissipation of advanced power systems together with space and weight limitations. The ships are divided in zones that need a particular HVAC–R system for each of them, besides dedicated cooling systems for the high heat generating components that they contain (e.g., electromagnetic gun, radar, electric motors, motor drives) that need to be driven by a finite fuel energy source [6]. Such systems need specific design strategies based on imposed high heat flux (with possible abrupt dynamic changes) at the evaporators differently from traditional vapor compression cycle systems design that is typically based on a fixed thermal load (or temperature differences). Therefore, the HVAC–R design challenge is to minimize energy consumption utilized for thermal management, and to develop VCR dynamic models that deal realistically with such boundary conditions for control and optimization [7]. In order to reduce energy consumption, directions that have been pursued by the industry and academia are related to system control, exergy analyses, and optimization. For example, Buzelin et al. [8] developed experimentally a novel power law continuous control system, and tested it in a vapor compression refrigeration system, obtaining up to 35% energy consumption reduction as compared to the same system operating with the traditional on– off control, under the same thermal load and environmental conditions. Kanoglu [9] developed a steady state exergy analysis of the multistage cascade refrigeration cycle used for natural gas liquefaction, and using a typical actual work input value found the system exergetic efficiency to be 38.5%, which indicated great potential for improvements. With respect to HVAC–R optimization, numerous studies were found in the literature, from which a few are herein discussed. Shiba and Bejan [10] showed how the internal geometric configuration of a component could be deduced by optimizing the global performance of the installation that uses the component, using the counterflow heat exchanger that serves as condenser in a vapor-compression-cycle refrigeration system for environmental control of aircraft. The optimization of global performance was achieved by minimizing the total power requirement or the total entropy generation rate at steady state operation. Sanaye and Malekmohammadi [11] minimized the total cost per unit cooling load of a vapor compression system including capital investment for components as well as the required electricity cost, using the Lagrange multipliers method to find optimal design parameters for steady state operation. The effects of changing the cooling load on optimal design parameters were studied. Hovsapian et al. [12] using a solar concentrator, and Vargas et al. [13] using a solar collector, performed the thermodynamic optimization of hybrid solar and combustion chamber driven absorption refrigerators and found optimal coupling thermal fluid flow rates to minimize water heating pull-up and refrigeration pull-down times, as well as maximum second law efficiency at steady state. However, only the cold chamber was treated dynamically since an endoreversible model was used for the absorption refrigerator by neglecting its thermal inertia with respect to the other system components’ thermal inertias. Sharp optima were found, stressing their importance for actual system design. Tracy et al. [14] addressed the fundamental question of the splitting of a hot exhaust into two heat recovery heat exchangers that are part of a trigeneration system (electricity, refrigeration, and hot water by recovering waste energy from an internal combustion engine). First and second law analyses were performed to find the optimal splitting of the available hot exhaust stream between the refrigerator and the hot water heat exchangers. Ordonez and Bejan [15] conducted a thermodynamic optimization study to find the cabin inlet temperature that

minimizes power consumption in an aircraft environmental control system. The problem of optimal area allocation has been studied [16–22] in various types of refrigeration and power systems. Radcenco et al. [16] optimized (1) the frequency of on/off operation, and (2) heat exchanger allocation between evaporator and condenser through the entropy generation minimization methodology. Khan and Zubair [17] studied various steady state HVAC–R vapor compression models to determine the optimal heat transfer area allocation for maximum performance. Analogously, Chen et al. [18] developed analytical expressions for the optimal allocation of heat-exchanger area for the endoreversible and irreversible Carnot cycles with constant- and variable-temperature heat-reservoirs, and found optimal values for real vaporcompression cycles numerically. Area allocation was optimized in solar driven refrigeration systems for maximum heat extraction [19] and in a hot-gases driven power plant (pre-heating, boiling and superheating sections) for maximum power output [20]. Ordóñez et al. [21] optimized the area allocated to counter flow heat exchangers in an installation for maximum combined heat and power production. Luo et al. [22] gave guidelines to allocate the heat transfer surface area for maximum cooling load in thermoelectric refrigerators. The thermodynamic optimization (or entropy generation minimization, finite-time thermodynamics) methodology is a common point in these studies for steady state operation. Finally, Dingec and Ìleri [23] developed a thermoeconomic optimization study for refrigerators based on local and total irreversibilities to find optimal area allocation for minimum cost and thermodynamic losses. Several comprehensive mathematical models have been published to simulate the dynamic response of vapor compression refrigeration systems. Some well known benchmark models were developed by Dhar and Soedel [24], Chi and Didion [25], MacArthur [26], Rajendran and Pate [27], and Vargas and Parise [28], which allowed for the investigation of startup transients, as well as continuous feedback control, and global system energy consumption reduction. In order to evaluate the performance of refrigeration plants under all aspects of operation, such as pressures, temperatures and refrigerant mass distribution, dynamic simulation models are more appropriate than steady state models [7,29,30]. The design of the refrigeration equipment depends strongly on the properties of the selected refrigerant [31]. Accurate knowledge of thermodynamic properties is of great importance in the efficiency of the refrigeration equipment [32]. Pressures in the system are temperature-dependent and are different for each particular refrigerant. Refrigerant transport properties, such as liquid and vapor density, define heat transfer coefficients and consequently temperature differences in heat exchangers thus directly influencing pressures in the system as well as necessary heat transfer surface of heat exchangers [31]. The calculation methods of thermodynamic properties of refrigerants can be divided into two types [33]: (1) accurate methods, and (2) fast methods. In the first type, the equations of state are for the accurate prediction of thermodynamic properties in a wide range. Using appropriate refrigerants’ properties data banks, it is possible to evaluate either pure or mixture refrigerants as possible substitutes for banned substances, such as CFCs (chlorofluorocarbons) [32]. However, the calculation speed is limited by unavoidable iterations in calculation when the equations of state are used in refrigeration system simulation [33]. When the fast calculation methods are used instead of the equations of state method, the time duration for a single simulation can be reduced considerably. However, the majority of the existing fast methods is not accurate to all of the refrigerants and the physical states of the working fluid [33].

543

T.K. Nunes et al. / Applied Energy 158 (2015) 540–555

The literature review shows that numerous vapor compression HVAC–R mathematical models, and several related optimization studies have been published. However, no dynamic response optimization study was found in the literature for such plants that addresses in general terms, i.e., non-dimensionally, the variation of design and operating parameters, and their effect on performance, and that has been utilized for general system thermodynamic optimization (minimum system pull-down time, and maximum first and second law efficiencies). The objectives of this paper therefore are: (i) to introduce a general (dimensionless) transient mathematical model for vapor compression HVAC–R plants, and (ii) to minimize system pull-down time and extract maximum exergy input rate from the energy source. The approach is crossdisciplinary and pursues simultaneously: (i) the local optimization of components and processes with (ii) the optimal global integration and configuration of the system. The model is represented by a system of ordinary differential equations with respect to time, the solution of which consists of the temperatures of each control volume, from which, the resulting cooling rate, first and second law combined efficiencies are computed as functions of time, to be able to analyze the system dynamic response. The model is simple enough to ensure small computational time requirements, so that it is possible to simulate the system steady state and transient response in a large number of competing geometric and operating configurations. 2. Mathematical model Fig. 1 shows a schematic diagram of the vapor compression refrigeration system, and the thermodynamic control volumes assigned to each component. On the air side, three control volumes are defined as follows: the refrigerated space (CV1), the air side of the evaporator (CV2), and the air side of the condenser (CV3). The refrigerant side consists of the compressor (CV5), the expansion valve (CV6), and the refrigerant side of evaporator (CV4) and condenser (CV7). The HVAC–R system optimization problem is formulated as follows: 1. One objective function is defined to evaluate the system dynamic response: the pull-down time to achieve a prescribed setpoint temperature T 1;set , i.e., for operation as a cooling system, which needs to be minimized;

2. For designing purposes, two other objective functions are defined to be maximized at steady state operation, i.e.: (i) the coefficient of performance (COP) or 1st law efficiency, and (ii) the second law efficiency, and 3. Since the HVAC–R hardware, i.e., heat exchangers heat transfer areas, shown in Fig. 1, are commodities in short supply, it makes sense to recognize the total area inventory as a physical constraint for the optimization problem. The total HVAC–R area inventory, A, is defined as follows:

A ¼ A4 þ A7

Assuming that the heat exchangers global heat transfer coefficients are of the same order of magnitude, it is recognized that:

U4  U7  U

& Q w

& Q 1

ð2Þ

in which U is a reference global heat transfer coefficient. However, in order to evaluate the possibility of distinct global heat transfer coefficients in the evaporator and condenser, the following evaporator to condenser global heat transfer ratio is introduced:

U 7 ¼ U;

s¼

U4 U7

A4 A

ð4Þ

Therefore, the condenser heat transfer area fraction is given by 1  x4 . Similarly, the dimensionless heat transfer and work rate, temperature, time and heat capacity rate are defined by:

_ _ e i; W f i Þ ¼ ðQ i ; W i Þ ; ðQ UAT0 Air flow

CV2

& W cp

CV5 Compressor

a

CV7

c

CV3

~i ¼ m CV6 Expansion valve

Condenser

& Q 7

Ti ; T0

~t ¼

t ; m1 cv;a =ðUAÞ

wi;j ¼

_ i cj m UA

The dimensionless mass, enthalpy, specific heat, pressure ratio and area are given by

Evaporator

Refrigerant flow

si ¼

ð5Þ

& Q 4

CV4 d

ð3Þ

therefore U 4 ¼ sU 7 ¼ sU. There is a need to simulate the system steady state and dynamic response in a large number of competing geometric and operating configurations. Hence, the model should be simple enough to ensure small computational time requirements. The dimensionless model aims to achieve two main objectives. The first is numerical stability, since the nondimensionalization of the variables is done based on appropriate physical scales for the original ones, i.e., by using scale analysis. In this way, the calculated dimensionless variables are as close as possible to the unity, therefore avoiding the divergence of the algorithm, which is possible to occur with the original dimensional variables. The other objective of the dimensionless model is the generalization of the results because with dimensionless variables, the results are normalized, i.e., the graphs and numerical tables are valid to any geometric configuration (or architecture) with functional and physical characteristics similar to the system analyzed by the model. Thus, with the fixed total heat transfer area inventory in mind, the dimensionless evaporator heat transfer area fraction is defined by:

x4 ¼ CV1 – Controlled ambient

ð1Þ

mi ; m1

~ ¼ h i

hi ; cv;a T 0

~ci ¼

ci ; cv;a

~i ¼ p

pi ; p0

e i ¼ Ai A A

ð6Þ

In Eqs. (5) and (6), the subscripts i and j refer to a particular substance, CV or location as shown in Fig. 1.

b

2.1. Controlled ambient (CV1) Air flow

Fig. 1. Schematic diagram of the vapor compression refrigeration (VCR) system.

The model adopts the following assumptions for CV1: uniform e 1 , is properties, dry air, and constant mass. The thermal load, Q

544

T.K. Nunes et al. / Applied Energy 158 (2015) 540–555

estimated based on the internal surfaces, roof and fenestration, heat gain due to ventilation, infiltration and internal heat sources (e. g., people, equipment). A balance of energy (first law of thermodynamics) combined with mass conservation in CV1 states that:

ds1 e 1 þ w ðs2  s1 Þ ew þ Q ¼Q 2;a d~t

ð7Þ

e w ¼ zð1  s1 Þ, and z ¼ ðUAÞw . in which Q UA

configuration. However, for system level analyses, control, and optimization, simpler and low computational time demanding approaches suffice [7,24–28]. Therefore, in this model, similarly to previously published system oriented HVAC–R models [7,25,28], an homogeneous model is adopted, i.e., the variation of the refrigerant quality, y4 , across CV4 is assumed to be linear between the refrigerant inlet and outlet conditions resulting that:

y4 ¼

2.2. Air side of the evaporator (CV2)

y4;in þ 1 2

ð11Þ

The CV4 refrigerant liquid–vapor mixture enthalpy is given by:

For simplicity, in this model, the possibility of the use of outside air through the air side of the evaporator is not considered. However, the model could be easily amended to include such aspect. Heat transfer through the evaporator walls to CV1 is neglected in presence of the heat transfer with the refrigerant side of the evaporator (CV4). Therefore, the balance of energy in CV2 states that

~4 ¼ y h ~ ~ h 4 4v þ ð1  y4 Þh4l

e4 ds2 w2;a ðs1  s2 Þ  Q ¼ ~2 m d~t

in which the left-hand side term reflects the assumption of du=dt  dh=dt since dh ¼ du þ pdv þ v dp, under negligible pressure drop, i.e., v dp ! 0, and pdv  du in a phase changing process. This assumption was also previously used by MacArthur [26]. The dimensionless evaporator saturation temperature, s4 , is obtained from Eq. (13). For that, initially, the properties of saturated refrigerant are approximated through suitable polynomial interpolations, as follows:

ð8Þ

2.3. Air side of the condenser (CV3) In the air side of the condenser (CV3), with similar assumptions, the balance of energy states that:

e7 ds3 w3;a ð1  s3 Þ þ Q ¼ ~ ~ m3 dt

ð9Þ

Along the evaporator, the flow pressure drop is assumed negligible. A balance of energy applied to CV4 states that:


 ~ w h ~ ~ dm~ 4 e ~ ~4 w5;r h b 6;r d =c p;r þ Q 4  h4 d~t dh ¼ ~4 m d~t

ð14Þ

h4l ¼ a4 T 4 þ a5

ð15Þ

Therefore

With the system in operation, the refrigerant thermodynamic state varies along the evaporator. The refrigerant enters the evaporator (CV4) as a liquid–vapor mixture and exits as superheated vapor. The space-averaged and time varying values of the thermodynamic properties of the refrigerant are approximated as the properties of a liquid–vapor mixture with space-averaged and time varying quality y4 , i.e., the vapor mass fraction in two-phase flow (liquid and vapor). Since the model is homogeneous (uniform properties in CV4), y4 is also assumed equal to the void fraction, i.e., the vapor volume fraction in two-phase flow (liquid and vapor). At the expansion valve outlet, i.e., point c in the diagram shown in Fig. 1, the refrigerant is a 2-phase mixture of liquid and vapor. For that condition, the refrigerant dimensionless specific enthalpy ~ ¼y h ~ . Hence, the refrigerant quality at the ~ 4v  ð1  y Þh is h 4;in

4;in

4l

evaporator inlet is calculated as follows:

y4;in

~ h ~ h b 4l ¼ ~ ~ h4v  h

ð13Þ

h4v ¼ a1 T 24 þ a2 T 4 þ a3

2.4. Refrigerant side of evaporator (CV4)

b

ð12Þ

dh4v dT 4 dT 4 ¼ 2a1 T 4 þ a2 dt dt dt

ð16Þ

dh4l dT 4 ¼ a4 dt dt

ð17Þ

dy4 dt

Next deriving Eq. (12) with respect to time, and assuming
  dhdt4v ; dhdt4l , it follows that:

dh4 dh4v dh4l ffi y4 þ ð1  y4 Þ dt dt dt Using Eqs. (16) and (17) into Eq. (18), and solving for result is given by:

dT 4 dh4 =dt ¼ c4 dt

ð18Þ dT 4 , dt

the

ð19Þ

where

ð10Þ

4l

in which the refrigerant dimensionless specific enthalpy at the con~ ¼h ~ c since the expansion valve underdenser outlet is taken as h b goes an isoenthalpic process. At this point, it is important to note that void fraction detailed models for HVAC–R heat exchangers are abundant in the literature [34], but even the more complex correlations have been reported to be inaccurate in the past [35]. The various models include different assumptions about velocity ‘‘slip” between liquid and vapor phases and the role of total refrigerant mass flux. The prediction of refrigerant charge level has been shown to be strongly influenced by the choice of void fraction models [36]. Such models are particularly important for component development purposes, and can demand large computational time even for a single system

c4 ¼ ð1  y4 Þa4 þ 2a1 y4 T 4 þ a2 y4

ð20Þ

which has specific heat dimensions, i.e., J kg1 K1. Finally, the dimensionless evaporator saturation temperature, s4 , is obtained as follows:

~4 ds4 1 dh ¼ ~c4 d~t d~t

ð21Þ

~ , The dimensionless specific enthalpy at the compressor inlet, h d is calculated by assuming a linear variation of the enthalpy through the evaporator, as follows:

~ ¼ 2h ~ h ~ h 4 d b

ð22Þ

In order to obtain the refrigerant dimensionless temperature at the evaporator outlet, ideal gas behavior is assumed for the refrigerant in the evaporator superheated section, as follows:

545

T.K. Nunes et al. / Applied Energy 158 (2015) 540–555

~4v dsd 1 dh ¼ ~cp;r d~t d~t

ð23Þ

~ =d~t  dh ~4v =d~t, which is calculated with Eq. (23) assumes that dh d Eq. (16). The dimensionless specific volume at the evaporator outlet, i.e., ~ d ¼ vv d , is estimated at point d in the diagram shown in Fig. 1, v 0 based on empirical correlations for the selected refrigerant, if available, or through the van der Waals equation, so that:

v d ¼ f ðT d ; p4 Þ or

(

pþ

ð24aÞ

)

a ðv MÞ

2

ðv M  bÞ ¼ RT;

a¼

27R2 T 2crit 64pcrit

and b ¼

 crit RT 8pcrit ð24aÞ

in which, T crit should be in K, and pcrit in kPa, for the selected refrigerant. In Eq. (24), whenever the specific volume appears implicitly (e. g., van der Waals equation), the equation has to be solved numerically to obtain v d . ~4 ; is The refrigerant dimensionless pressure at the evaporator, p estimated through a suitable polynomial interpolation of the saturated refrigerant thermodynamic properties table, as follows:

p4 ¼ c1 T 24 þ c2 T 4 þ c3 ;

~4 ¼ and p

refrigerant average temperature in the condenser, similarly to CV4, as follows:

~7 ¼ y h ~ ~ h 7 7v þ ð1  y7 Þh7l

ð25Þ

Eq. (23) neglects the pressure effect on the variation of enthalpy in the superheated region, within a narrow range of temperatures that constitutes the degree of superheat at the evaporator outlet. The degree of superheat at the evaporator outlet, Dssh , is calculated by:

Dssh ¼ sd  s4

ð26Þ

The dimensionless heat interaction with the air side of the evaporator is given by:

e 4 ¼ sx4 ðs2  s4 Þ Q

ð27Þ

Conservation of mass applied to CV4 dictates that:

~4 dm 1 ¼ ðw  w5;r Þ ~cp;r 6;r d~t

ð28Þ

2.5. Refrigerant side of condenser (CV7) In the refrigerant side of the condenser, i.e., CV7 in the diagram shown in Fig. 1, the fluid refrigerant enters as superheated vapor and exits it as subcooled liquid. The space-averaged and time varying values of the refrigerant thermodynamic properties are approximated as the properties of a liquid–vapor mixture with a constant time- and space-averaged quality y7 ¼ 0:5. This is a result of assuming an homogeneous model for the refrigerant liquid–vapor mixture, and a linear variation of the quality, and refrigerant entering CV7 superheated ðy7 ¼ 1Þ and leaving CV7 always subcooled ðy7 ¼ 0Þ. Since the model is homogeneous (uniform properties in CV7), y7 is also assumed equal to the void fraction, i.e., the vapor volume fraction in two-phase flow (liquid and vapor). Along the condenser, the flow pressure drop is assumed negligible. The resulting balance of energy applied to CV7 states that:

~a  w h ~ ~ ~ dm7 e ~7 ðw5;r h dh 6;r b Þ=cp;r  Q 7  h7 d~t ¼ ~7 m d~t ~

ð29Þ

Note that the left-hand side term of the Eq. (29) reflects the assumption of du=dt  dh=dt, also used in Eq. (13) for CV4. The

ð30Þ

where the same polynomials of Eqs. (14) and (15) are used to evaluate h7v and h7l replacing T 4 by T 7 . Next, based on the same reasoning as Eqs. (16)–(20), replacing T 4 by T 7 , and y4 by y7 , the dimensionless condenser saturation temperature is calculated by:

~7 ds7 1 dh ¼ ~c7 d~t d~t

ð31Þ

Assuming negligible refrigerant mass variation with respect to time in the expansion valve and in the compressor, in presence of the refrigerant mass variation in the heat exchangers (condenser and evaporator), conservation of mass applied to the global system states that:

~7 ~4 dm dm ¼ ~ dt d~t

ð32Þ

e 7 ; is calThe total dimensionless condenser heat transfer rate, Q culated with the global energy balance on the HVAC–R system. First, it is recognized that the expansion valve is adiabatic, there
 fore dE6 ¼ 0; dE5  dE4 ; dE7 , since Q_ 5 is negligible in presence of dt

p4 p0

s7 , is calculated

dt

dt

dt

7 < 0, therefore assuming Q_ 4 and Q_ 7 . Next, note that dEdt4 > 0 and dE dt     dE4  dE7  dEGS dE4 dE7 that  dt  ffi  dt  in time, dt ffi dt þ dt ffi 0, and as a result:

e4 þ W e7 ¼ Q f cp Q

ð33Þ

The refrigerant undergoes a change of phase along the condenser, i.e., superheated vapor, two-phase (liquid + vapor), and subcooled liquid. The dimensionless heat transfer rate between the refrigerant and air side of the condenser is divided in two contributions in this model: (i) superheated vapor and two-phase e sh2p , and (ii) subcooled liquid zone, Q e sc . The result is given zones, Q by:

e sh2p þ Q e sc e7 ¼ Q Q

ð34Þ

in which

e sh2p ¼ ð1  x4 Þðs7  s3 Þ Q

ð35Þ

e sh2p is dominated by the two-phase where it is assumed that Q zone. Using Eqs. (33) and (34), the dimensionless heat transfer rate in e sc , is obtained as follows: the subcooled liquid zone, Q

e7  Q e sh2p e sc ¼ Q Q

ð36Þ

Applying the first law of thermodynamics to the subcooled liquid zone, the degree of subcooling is calculated as follows:

Dssc ¼

e sc Q w5;r

ð37Þ

Eq. (37) also assumes quasi steady conditions in the condenser subcooled zone, i.e., dE7;sc =dt ! 0. As a result, the refrigerant dimensionless temperature at the condenser outlet is calculated by:

sb ¼ s7  Dssc

ð38Þ

The refrigerant subcooled liquid dimensionless specific volume and enthalpy at the condenser outlet, i.e., point b in the diagram shown in Fig. 1, are estimated based on a suitable linear

546

T.K. Nunes et al. / Applied Energy 158 (2015) 540–555

interpolation of the saturated liquid refrigerant specific volume thermodynamic properties table evaluated at the same temperature, T b , as follows:

v b ¼ b1 T b þ b2 ;

hb ¼ a4 T b þ a5 ;

and v~ b ¼

vb v0

ð39Þ

The refrigerant dimensionless saturation pressure at the con~7 , is estimated through Eq. (25) replacing T 4 by T 7 . denser, p 2.6. Expansion valve (CV6)

1. Pull-down time minimization in transient operation: for each dimensionless evaporator heat transfer area fraction tested, x4 , the model is used to determine the dimensionless pull-down time, ~t pd , to reach the controlled ambient temperature setpoint, s1;set . 2. For designing purposes, the 1st and 2nd law efficiencies maximization are pursued: for each dimensionless evaporator heat transfer area fraction tested, x4 , the model is used to determine the steady state system first law, gI , and second law, gII , efficiencies as follows:

COP ¼ gI ¼

In the expansion valve, the refrigerant undergoes an isoenthalpic process, and enters as subcooled liquid and exits as a liquid–vapor mixture. A thermostatic valve is assumed by the model, therefore the expansion orifice area varies according to the desired degree of superheat at the evaporator outlet, as follows:

e set þ G e6 ¼ A e ðDssh  Dsdsh Þ A

constant. The dimensionless mass flow rate of the refrigerant through the expansion valve is calculated by

 1=2 e 6 2 ðp ~7  p ~4 Þ w6;r ¼ A n v~ b in which n ¼

cp;r U


1=2 p0

v0

ð41Þ

.

2.7. Compressor (CV5)

 1=n ~ p gv ¼ 1 þ c0  c0 ~7 p4

( )  n1 ~7 n n p g p~4 / 1  ~ 1n v p4

v~ d

in which

ð43Þ

ð44Þ

l ¼ cpv;a0 vT00 .

The first law of thermodynamics applied to CV5 states that:

f ~a ¼ h ~ þ W cp ~cp;r h d w5;r

where gC ¼ 1s1s1 is the system Carnot efficiency. 2.9. Numerical method The numerical problem consisted of integrating Eqs. (7), (8), (11), (13), (17) and (23), with all auxiliary algebraic equations in time, and solving a nonlinear algebraic equation with the secant method, Eq. (24), at each time step, using an adaptive time step, 4th–5th order Runge–Kutta method [37]. In order to solve Eq. (24) with the secant method, for fast convergence, an initial guess for the specific volume of the superheated refrigerant at the evaporator outlet was estimated with the ideal gas equation of state, i. RT v d;0 ¼ Mp , at each time step. The time step was adjusted autod

d

matically according to the local truncation error, which was kept below a specified tolerance of 106. Steady state conditions were ~4 ; h ~ 7 ; s4 ; s7 ; m ~ 4; m ~ 7 ; s1 ; x ¼ ðh verified when @@~x~t 6 104 , where ~

s2 ; s3 Þ is the vector of dimensionless system variables defined by

2.10. Experiments

The dimensionless heat capacity rate in the compressor is given by:

gv C v ~ /cp;r l

ð47Þ

ð42Þ

p0 V c where / ¼ rps UAT is the dimensionless compressor speed. 0

w5;r ¼

gI gC

the model.

in which c0 is the dead volume ratio in the compressor, and n is the compressor polytropic coefficient. The dimensionless compressor power is given by:

f cp ¼ W

gII ¼

e.,

In the compressor, CV5, quasi steady conditions, negligible heat exchange across the cylinder walls, and a polytropic compression process are assumed. The compressor volumetric efficiency, gv , is calculated as follows:

for operation as a cooling system ð46Þ

ð40Þ

e ¼ G is the dimensionless thermostatic valve adjustment where G A=T 0

e4 Q ; e4 e Q7  Q

ð45Þ

2.8. HVAC–R system optimization In order to solve the optimization problem formulated in the beginning of Section 2, the search for the HVAC–R system optimization is conducted as follows:

A vapor compression refrigeration system was set up and instrumented in the Center for Advanced Power Systems (CAPS) at Florida State University (FSU), USA. The objective was to produce transient experimental data for the mathematical model adjustment and experimental validation. The system consists of a chiller, an external tank, water heaters, thermally insulated flexible tubing, and valves, as shown in Fig. 2. The chiller was manufactured by MTA, Tribano, Italy (model TAE-TWE) in 2003, as Fig. 3 shows with connections to the external tank and a data acquisition system. The system is loaded with 2.65 kg of refrigerant fluid R22. Since the chiller has a sealed water storage tank, an 115-L external water tank was added to the system, as shown in Figs. 2 and 3. The objective was to allow for testing the system with different thermal loads. For that, water heaters were installed inside the tank. The system was instrumented with 10 high precision thermistors, manufactured by Measurement Specialties, Model YSI 44004RC (2.252 X @ 25 °C), with a recommended operating range from 80 °C to 120 °C. The thermistors were calibrated in the laboratory with 24 temperature measurements in a controlled temperature thermal bath in the range from 9.5 °C to 86.5 °C, in order to determine constants for each of them. A National Instruments data acquisition system together with the software LabView 8.2 was used to monitor and record the experimental data. The system used the boards NI PXI 1010, NI SCXI-1140 and NI SCXI-1300. The refrigeration system high and

547

T.K. Nunes et al. / Applied Energy 158 (2015) 540–555

Valve Water input Valve

Electric heaters Water tank

Air output

Air input

Air input

Chiller

Insulated flexible tubing

Fig. 2. Schematic diagram of the experimental rig.

whereas for the second set 1200 W, which was used for model experimental validation. The two tests were conducted three times so that uncertainties could be calculated, and the sampling rate was 100 ms. Model adjustment was conducted to determine appropriately system geometric and operating parameters used by the mathematical model, that were not accurately known either by direct measurements, or by the manufacturer’s specs. The procedure consisted of the solution of the inverse problem of parameters estimation (IPPE) [40] using the measured first data set (300 W), and the proposed mathematical model. In this way, measured temperatures were imposed as input parameters, and the not accurately known system geometric and operating parameters calculated by the model as unknowns. Next, the adjusted parameters were used by the model to simulate the system response operating in the conditions of the second data set (1200 W). The simulation results were then compared to the measured temperatures, in order to experimentally validate the model. 3. Results and discussion In this section, initially, a set of experimental data is used to adjust the mathematical model. Then, a new set of experimental data is used to experimentally validate the adjusted model. Finally, a thermodynamic optimization of a vapor compression refrigeration system is conducted for three different refrigerant fluids. 3.1. Uncertainty in experimental data The temperature uncertainties in 10 measuring points were estimated according to the methodology described in Section 2.10. Additionally, for higher accuracy, the largest calculated uncertainty in each of the measured temperatures in the two experiments during transient operation was considered for all points. Accordingly, error bars were drawn in the experimental points in the graphs with the calculated values shown in Table 1. 3.2. Mathematical model adjustment Three parameters were selected to adjust the mathematical model: (i) the overall heat transfer coefficient in the condenser; (ii) the fraction of the total heat transfer area allocated to the

Table 1 Dimensionless uncertainties calculated for the measured temperatures. Fig. 3. Photograph of the experimental rig and data acquisition system.

low pressures were recorded with analog manometers for the refrigerant R22 installed in the chiller’s panel. The average of the measurements taken in all runs was utilized to estimate the temperature at each point in the refrigeration system. The precision limit of the measurements was calculated as twice the standard deviation of the experimental runs assuming that the population follows a symmetric unimodal normal distribution, within a 95% confidence interval [38,39]. The maximum calculated bias limit observed during the thermistors calibration was ±0.005 °C, which was negligible in presence of the temperature measurements precision limit calculated during transient operation, which was therefore assumed as the measurements uncertainty. The setpoint temperature was the same for two experimental conditions. The first data set was used for model adjustment, in which the water heaters dissipated 300 W inside the external tank,

Symbol

Temperature

Dimensionless uncertainty

sa

Dimensionless refrigerant temperature at compressor’s exit Dimensionless refrigerant temperature at condenser’s inlet Dimensionless refrigerant saturation temperature at condenser Dimensionless refrigerant temperature at condenser’s exit. Dimensionless refrigerant temperature – expansion valve Dimensionless refrigerant temperature at evaporator’s inlet Dimensionless refrigerant temperature at evaporator’s exit Dimensionless water temperature at the expansion tank Dimensionless water temperature at evaporator’s exit Dimensionless air temperature at condenser

±0.11

sa;7 s7 sliq;7 sliq;7 sc sd s1 s2 s3

±0.12 ±0.03 ±0.01 ±0.01 ±0.04 ±0.01 ±0.02 ±0.03 ±0.03

548

T.K. Nunes et al. / Applied Energy 158 (2015) 540–555

Table 2 Values of adjusted parameters – using the 300 W data set. Symbol

Parameter

Value

U7 x4

Overall heat transfer coefficient in the condenser Fraction of total heat transfer area allocated to evaporator Refrigerant quality at condenser

0.01 W m2 K1 0.38

y7

0.55

Table 3 Geometric and operating parameters used for model adjustment and experimental validation. Q_ 1 ¼ 0:3 kW T 0 ¼ 295:15 K

A ¼ 63:83 m2 Aset ¼ 106 m2 Aw ¼ 1:32 m2

U ¼ 0:01 kW m2 K1 1

K1

1

1

cp;a ¼ 1:005 kJ kg cv;a ¼ 0:716 kJ kg cv ¼ 1:008 c0 ¼ 0:03743

K

G ¼ 0:5  107 m2 K1 m1 ¼ 113 kg

U w ¼ 1:123  102 kW m2 K1 DT dsh ¼ 16 K rpm = 2900 System initial conditions: ~ ¼ hðx4;0 ;T 4;0 Þ ; with y ¼ 0:60 h 4;0 4;0 cv;a T 0

m2 ¼ 110 kg

~ ¼ hðx7;0 ;T 7;0 Þ ; with y ¼ 0:55 h 7;0 7;0 cv;a T 0

m3 ¼ 0:269 kg

~ 4;0 ¼ 2:12  102 m

_ 1 ¼ 0:636 kg s1 m _ 3 n ¼ 0:02 kg s1 m _ 1 ¼ 2:30 kg s1 m V c ¼ 6:789  105 m3 n ¼ 1:132 p0 ¼ 0:1 MPa

Fig. 4. Dimensionless temperatures in refrigerated space (CV1) with 300 W thermal load (model adjustment).

~ 7;0 ¼ 6:02  103 m

s1;0 ¼ 1:011 s2;0 ¼ 0:997 s3;0 ¼ 1:005 s4;0 ¼ 0:9922 s7;0 ¼ 1:003

evaporator; and (iii) the quality of the refrigerant in the condenser. Table 2 shows the adjusted values for these parameters according to the methodology described in Section 2.10, i.e., the inverse problem of parameters estimation (IPPE). Table 3 shows the geometric and operating parameters of the experimental refrigeration system, and the initial conditions for the model adjustment that was conducted with the mathematical model. Fig. 4 illustrates the temporal evolution of the experimental dimensionless temperature for the refrigerated space (CV1) and the corresponding temperature obtained with the adjusted mathematical model. Good agreement is observed both in the trend (cooling) and quantitative values (simulated values fall within the uncertainty in the experimental data). Fig. 5 shows the dimensionless saturation temperature for the refrigerant in CV7. Again, the adjusted model follows the experimental trends and for the most part stays within the uncertainty limits.

Fig. 5. Dimensionless refrigerant temperatures in the condenser (CV7) with 300 W thermal load (model adjustment).

3.3. Mathematical model experimental validation After the adjustment of the mathematical model a second data set for a thermal load of 1200 W was used for the model experimental validation. Figs. 6 and 7 illustrate the experimental data and the simulation results. The experimental and simulation results for the dimensionless saturation temperature in the condenser (CV7) are shown in Fig. 7. The comparison indicates that the model predicts the system behavior taking into account the uncertainties in the experimental data, therefore the model is considered experimentally validated according to the methodology described in Section 2.10. 3.4. System thermodynamic optimization After experimental validation, the model was considered reliable to perform a thermodynamic optimization study. For that,

Fig. 6. Dimensionless temperatures in refrigerated space (CV1) with 1200 W thermal load (experimental validation).

549

T.K. Nunes et al. / Applied Energy 158 (2015) 540–555 Table 5 Refrigerants’ coefficients, constants and properties. Variable

R12

R134a

R1234yf

Obs

a1 a2 a3 R_hv a4 a5 R_hl b1 b2 R_vl c1 c2 c3 R_psat pcrit ðkPaÞ T crit ðKÞ M (kg kmol1)

0.00079898 0.42074 187.58 0.99996 0.94094 36.772 0.99975 1.7979e6 0.00072218 0.99442 0.00013694 0.011086 0.30459 0.99949 4115 385.15 120.93 0.612

0.0019937 0.56359 247.85 0.99919 1.5408 54.332 0.99451 2.3877e6 0.00078627 0.98632 0.00025045 0.011438 0.20888 0.99801 4060.3 374.23 102.03 0.9

0.0022087 0.62917 364.15 0.99876 1.4527 203.66 0.99755 2.8292e6 0.00086781 0.98426 0.00020706 0.011836 0.26954 0.99891 3382.2 367.85 114.04 0.92618

233.15– 353.15 K

0.935

1.34

1.2893

1

cp;r ðkJ kg Fig. 7. Dimensionless refrigerant temperatures in the condenser (CV7) with 1200 W thermal load (experimental validation).

three refrigerant fluids were considered: R12 (Dichlorodifluoromethane), R134a (1,1,1,2-Tetrafluoroethane), and R1234yf (2,3,3,3-Tetrafluoropropene). Although R12 (also known as Freon12) has been banned due to atmospheric ozone layer depletion effects, it is well known to lead to excellent HVAC–R performance. Therefore, R12 is herein used as a reference, in order to investigate the impact on performance of a system originally designed for R12, operating with R134a, and R1234yf. The system parameters and constant values were selected, ~4 ; h ~7 ; s4 ; s7 ; m ~ 4; according to Table 4, and initial conditions ~ x0 ¼ ðh ~ 7 ; s1 ; s2 ; s3 Þ0 were established to complete the initial value probm lem formulation. For clarity, the dimensional values used in this study are presented in Table 4, from which the dimensionless input parameters were calculated. However, the simulation results are expected to be valid for any set of dimensional parameters that lead to the dimensionless parameters used in the simulations performed in this study, as explained previously in the text. Aiming at conducting an optimization study, a fast computational methodology was selected to calculate the refrigerants’ properties. As described by the model, polynomial curve fitting was used in the saturation region, and for the superheated vapor

K1 Þ

cr , liquid 1

ðkJ kg

K1 Þ

233.15– 353.15 K 233.15– 333.15 K 233.15– 363.15 K

Saturation @ 273.15 K Saturation @ 273.15 K

region, an equation of state reported by Gatecliff and Lady [41] for R12, and the van der Waals equation for R134a and R1234yf were utilized. For that, the necessary coefficients, constants and properties are listed in Table 5. In order to determine the polynomial coefficients of Table 5, all refrigerant thermodynamic properties in the saturation region were obtained from the NIST REFPROP Database [42]. The following analysis is conducted for the system defined in Fig. 1, with the parameters listed in Table 4. The transient evolution of the dimensionless temperatures that characterize the system response and computed in the simulation are shown in Fig. 8 for a system with equal heat transfer areas in the evaporator and condenser, x4 ¼ 0:5, and same global heat transfer coefficients in both heat exchangers as well, s = 1, which are parameters also used in Figs. 9–13. At ~t ¼ 0, all dimensionless temperatures are equal to 1, since the system departs from thermal equilibrium with the external ambient. The system is allowed to run until steady state conditions are achieved. The response of the same system operating with the three different refrigerants is depicted in the graph. On the condenser side, heat is rejected to the external ambient, and since the input air flow is always at the ambient temperature, a similar response is observed with the three refrigerants. On the other hand, on the evaporator side, due to the confined thermal

Table 4 System parameters, constant values, and ODE’s initial values for the optimization. A ¼ 12 m

1.2

T 0 ¼ 298:15 K

e set ¼ 8:33  10 Aset ¼ 10 m ð A Aw ¼ 54 m2 6

8

Þ

2

U ¼ 0:1 kW m

K

U w ¼ 1:472  103 kW m2 K1

1

K1 ð~cp;a ¼ 1:4Þ

V c ¼ 7  10

1

K1

z ¼ 1:227  103 DT dsh ¼ 15 K ðDsdsh ¼ 0:0503Þ / ¼ 0:0039 (for 1200 rpm) ODE’s initial values:

cp;a ¼ 1:005 kJ kg cv;a ¼ 0:716 kJ kg cv ¼ 0:808 c0 ¼ 0:05

G ¼ 0:5  107 m2 K1 e ¼ 1:24  106 Þ ðG

4

m3

m1 ¼ 27 kg

~ ¼ hðx4;0 ;T 4;0 Þ ; with x ¼ 0:6 h 4;0 4;0 cv;a T 0

~ 2 ¼ 0:0011Þ m2 ¼ 0:03 kg ðm

~ ¼ hðx7;0 ;T 7;0 Þ ; with x ¼ 0:5 h 7;0 7;0 cv;a T 0

~ 3 ¼ 0:0011Þ m3 ¼ 0:03 kgðm _ 1 ¼ 0:5 kg s1 m ðw2;a ¼ w3;a ¼ 0:419Þ n ¼ 1:1 p0 ¼ 0:1 MPa e 1 ¼ 0:0112Þ Q_ 1 ¼ 4 kW ð Q

~ 4;0 ¼ 0:009815 m ~ 7;0 ¼ 0:001852 m

s1;0 ¼ s2;0 ¼ s3;0 ¼ s4;0 ¼ s7;0 ¼ sd;0 ¼ 1

R134a R12 R1234yf

x 4 = 0.5

1

s= 1 1.1

1

0.9

0.8 0

4

8

12

16

20

24

28

32

~ t Fig. 8. The VCR dimensionless temperatures dynamic response in CV1, CV2, CV3, CV4 and CV7 with R134a, R12, and R1234yf.

550

T.K. Nunes et al. / Applied Energy 158 (2015) 540–555 9 10 -8

20

R134a R12 R1234yf

x 4 = 0.5

~ p7

s=1

15 8 10 -8

~ p4 , ~ p7 10

~ A6

R134a R12 R1234yf

x 4 = 0.5

s=1

7 10 -8

5

~ p4 0 0

4

8

12

16

20

24

28

32

6 10 -8

~ t

0

4

8

12

16

20

24

28

32

~ t

Fig. 9. The VCR dimensionless saturation pressures dynamic response in CV4 and CV7 with R134a, R12, and R1234yf.

Fig. 12. The VCR dimensionless thermostatic expansion valve orifice area dynamic response with R134a, R12, and R1234yf.

0.18 3

R134a R12 R1234yf

x 4 = 0.5

0.16

s=1

0.14

0.04

R134a R12 R1234yf

x 4 = 0.5

s=1

0.03 0.12

2

0.1

COP

~

0.06

~

0.02 Wcp , Q 4

COP

0.08

~ Q4

1

0.01

0.04

~ Wcp

0.02 0

0

0

0

4

8

12

16

20

24

28

32

~ t

0.01

~ m 7

0.007

R134a R12 R1234yf

x 4 = 0.5

s=1 0.004

~ m 4 0.001 0

4

8

12

16

4

8

12

16

20

24

28

32

~ t

Fig. 10. The VCR dimensionless degrees of superheating and subcooling dynamic response with R134a, R12, and R1234yf.

~ ,m ~ m 4 7

0

20

24

28

32

~ t Fig. 11. The VCR dimensionless refrigerant mass dynamic accommodation in the evaporator and condenser with R134a, R12, and R1234yf.

load in the controlled ambient, the response is quite different depending on the refrigerant. The lowest steady state refrigerant temperatures, s4 , are observed with R12, followed by R1234yf and R134a. The same trend is observed for the evaporator air side temperature, s2 , and the controlled ambient temperature, s1 ,

Fig. 13. The VCR dimensionless compressor power input, refrigeration rate and COP dynamic response with R134a, R12, and R1234yf.

which indicates that the highest refrigeration rate is obtained with R12. As reported earlier in the text, the system was originally sized for refrigerant R12, and it is well known that a vapor compression HVAC–R system should be sized according to the selected refrigerant, which will dictate appropriate heat transfer areas and other design parameters. Therefore, the results should not be understood in a way that R12 is a refrigerant that performs better than R134a and R1234yf. However, in case of existing systems designed for R12, the simulations anticipate how the system would respond to the substitution of a particular refrigerant for an alternative one. In the current analysis, it is seen that R1234yf depicts a closer response to R12 than R134a. The transient evolution of the high and low pressure lines in the system, p4 and p7 , is shown in Fig. 9. The pressure difference is directly responsible for the resulting refrigeration effect since the refrigerant temperatures increase and decrease in a similar manner, as it was shown in Fig. 8, i.e., through s4 and s7 . The pressure difference results from the compressor power input, which in the current model results from the compressor speed, geometry, and the properties of the selected refrigerant. Therefore, the smaller the pressure difference, the lower the compressor power input and energy consumption rate. For the three tested refrigerants, R12 shows the smallest pressure ratio, followed by R1234yf and R134a. Fig. 10 depicts the transient behavior of the refrigerant degrees of superheating at the evaporator inlet and subcooling at the

551

T.K. Nunes et al. / Applied Energy 158 (2015) 540–555

The larger the orifice area, the larger the refrigerant mass flow rate, and a larger refrigeration rate is expected, which is also confirmed by the lowest evaporation temperature, s4 , observed with R12 in Fig. 8. The transient evolution of the system performance with the three tested refrigerants is displayed in Fig. 13. The R12 and R1234yf required compressor power input is smaller than what is required by R134a. However, R12 shows the highest refrigeration rate, followed by R1234yf and R134a. As a result, R12 leads to the highest COP for the tested system configuration, followed closely by R1234yf, and more distantly by R134a. The analysis of the curves show that R1234yf could substitute R12 for the tested system without a large performance drop, but for R134a the system should have the operating and design parameters reevaluated. For that, the simulation model could be readily used. All curves in Figs. 8–13 captured the expected physical trends for the tested system, and therefore this gives reliability to the model to evaluate the system dynamic response until steady state is achieved, and to perform the system thermodynamic optimization. A transient model is also mandatory if one is to evaluate system performance operating under fluctuations in thermal load and operating parameters. The search for system thermodynamic optimization opportunities started by monitoring the behavior of s1 in time, for three different evaporator area fractions, i.e., x4 ¼ 0:3; 0:5; and 0:85, for the three tested refrigerants. Fig. 14 shows that there is an

condenser outlet. The thermostatic expansion valve has its orifice regulated by a desired degree of superheating at the evaporator outlet, so that no liquid enters the compressor, therefore, for the three tested refrigerants, the transient response at the evaporator outlet is similar. However, the obtained degrees of subcooling depend on the refrigerants’ properties, and are distinct for each refrigerant, with R134a showing the lowest degree of subcooling (less heat is rejected by the condenser than with R12 and R1234yf), which result in the lowest refrigeration rate for the tested system. The transient refrigerant mass distribution in the system (evaporator and condenser according to the model) is assessed through Fig. 11. The initial conditions were set for a system at rest, in which more refrigerant mass is found in the evaporator since the compressor is off, and the refrigerant migrated from the high pressure region (condenser) to the low pressure region (evaporator) until equilibrium is reached. When the compressor is turned on, the pressure builds up in the condenser, and drops in the evaporator, so that more refrigerant mass accommodates in the condenser, which is captured by the model. The largest refrigerant mass gradient is observed with R12, followed by R1234yf and R134a, which is consistent with the lowest evaporation temperature, s4 , observed with R12 in Fig. 8. The dynamic variation of the expansion valve orifice is shown in Fig. 12, which is adjusted according to the desired degree of superheating at the evaporator outlet. The curves show that for R12 the orifice area stabilizes at a higher value than R1234yf and R134a.

1 1

τ1,set = 0.97 ; s = 1; R12

τ 1,set = 0.97 ; s = 1; R134a

τ1

τ1

0.95

0.3 = x 4

0.3 = x 4

0.85

0.85 0.5

0.5 0.9

0.95 0

4

1 424 3 8

12

16

20

24

28

0

32

~ t

~ tpd

4 { ~

8

12

16

20

24

28

32

~ t

tpd

(a)

(b) 1

τ1,set = 0.97; s =1; R1234yf

τ 1 0.95

0.3 = x 4

0.85 0.5

0.9 0

4

{ ~

8

12

16

tpd

(c)

20

24

28

32

~ t

Fig. 14. The effect of the variation of the VCR system evaporator area fraction, x4 , on the system pull-down time to achieve s1;set ¼ 0:97 with R134a (a), R12 (b), and R1234yf (c).

552

T.K. Nunes et al. / Applied Energy 158 (2015) 540–555

(3) For each evaporator heat exchanger area fraction, x4 , use the model to determine the pull-down time, ~t pd , to reach the refrigerated space temperature setpoint, s1;set . (4) From all tested evaporator heat exchanger area fractions, x4 , select x4;opt that leads to the one-way minimized pulldown time, ~tpd;min . (b) System first (COP) and second law efficiency maximization: (1) Fix the system total heat exchangers area, A, and the evaporator to condenser global heat transfer ratio, s ¼ UU47 ¼ UU4 , for a given U, at any desired values.

intermediate value of the evaporator area fraction, between 0.3 and 0.85, such that the temporal temperature gradient is maximum, minimizing the pull-down time to achieve a prescribed setpoint temperature (s1;set ¼ 0:97) for the three tested refrigerants, i.e., for R134a in Fig. 14a, R12 in Fig. 14b, and R1234yf in Fig. 14c. The system thermodynamic optimization for optimal heat exchanger area allocation was carried out in this study for minimum pull-down time, and maximum system performance. The range of variation of the evaporator area fraction to be optimized was: 0:1 6 x4 6 0:9. The chosen discretization in that range for the pull-down minimization, and performance maximization was the coarsest set for which the optimal value of the parameter did not change as the sets became finer, while the relative error was kept below 1% in all cases. The VCR system optimization procedure utilized to obtain the results shown in Figs. 15–17 is summarized by the following brute force algorithms: (a) Pull-down time minimization: (1) Fix the system total heat exchangers area, A, and the evaporator to condenser global heat transfer ratio, s ¼ UU47 ¼ UU4 , for a given U, at any desired values. (2) Vary the evaporator heat exchanger area fraction, x4 : 0:1 6 x4 6 0:9.

(2) Vary the evaporator heat exchanger area fraction, x4 : 0:1 6 x4 6 0:9. (3) For each evaporator heat exchanger area fraction, x4 , use the model to determine the steady state system first and second law efficiencies, gI (COP) e gII , Eqs. (46) and (47), respectively. (4) From all tested evaporator heat exchanger area fractions, x4 , select x4;opt that leads to the one-way maximized system first and second law efficiencies, gI;max ðor COPmax Þ and gII;max . Fig. 15 shows the VCR system optimization results obtained for the three tested refrigerants. The results for R134a, R12 and R1234yf are shown in Fig. 15a–c, respectively.

(a)

(b)

(c) Fig. 15. The minimization of the VCR system pull-down time, and maximization of second law efficiency, COP and dimensionless refrigeration rate with respect to the VCR system evaporator area fraction, x4 , with R134a (a), R12 (b), and R1234yf (c).

553

T.K. Nunes et al. / Applied Energy 158 (2015) 540–555 6 2

x ~ 4 ,opt tpd , min

5

x 4,opt COPmax

1.5

4

~ tpd ,min

COPmax

3

1

2 0.5

1

x 4 , opt

x 4, opt

0 1 R134a

2 R12

0

3 R1234yf

1 R134a

(a)

2 R12

3 R1234yf

(b)

0.6

0.5

x 4,opt 0.4

x 4 ,opt

η II , max

0.3

0.2

0.1

η II , max 0

1 R134a

2 R12

3 R1234yf

(c) Fig. 16. The optimal evaporator area fraction, x4;opt , for minimum pull-down time, ~tpd;min (a), COPmax (b), and maximum second law efficiency, gII;max (c), with R134a, R12, and R1234yf.

It is clear that there is a sharp minimum pull-down time for an optimal evaporator area fraction allocation for all refrigerants. The optimal evaporator heat exchanger area fraction, x4 , is ‘‘robust” with respect to the variation of refrigerant, being the same for all of them. The one way minimized pull-down time happens at the optimal value x4;opt ffi 0:55, for the three refrigerants, according to the results of Fig. 15. The steady state system first law efficiency gI ðor COPÞ also shows a maximum with respect to evaporator heat exchanger area allocation, but not as sharp as for the pull-down time. Also in Fig. 15 it is noted that the evaporator refrigeration rate shows a maximum for all tested refrigerants. Therefore, the question to be answered is: why does the system show a sharp pull-down minimum? The answer is because the first law analysis does not fully capture the losses experienced by the system as the system area allocation is varied. As the heat transfer area is reduced, thermal contact either with the cold or hot reservoirs through the condenser and evaporator worsens, which shows there must be an optimal heat transfer area allocation for maximum thermal contact in both sides of the system, so that the system gets as close as possible to the ideal Carnot behavior. This discussion brings to light the need to evaluate the performance of the VCR system (or thermal systems in general) on a more concrete basis, which is provided by exergy analysis, i.e., by evaluating the resulting second law efficiency of the VCR system, according to Eq. (47). It is then clear that, by maximizing the second law efficiency, the potential for use of the system cooling rate (at s1;set ¼ 0:97) will be maximized. The optimization results are summarized in Fig. 16, which shows that minimum pull-down time and second law efficiency

occur for the three refrigerants at approximately the same evaporator heat transfer area allocation, x4;opt ffi 0:55, as seen in Fig. 16a and c. However, maximum COP occurs at x4;opt ffi 0:4 for the three refrigerants, which is slightly shifted from the actual optimal area allocation for minimum pull-down time, as it was physically explained in the previous paragraph. Next, the design issue to be addressed concerns the possibility of distinct global heat transfer coefficients in the evaporator and condenser. For that, simulations were conducted with the model for the system operating with refrigerant R1234yf and the results are shown in Fig. 17. The variation of the evaporator to condenser global heat transfer coefficient ratio (s = 0.67, 1, and 2) is addressed in Fig. 17a showing the system performance through the pulldown time, COP, second law efficiency and dimensionless evaporator heat transfer rate with respect to the evaporator heat transfer area fraction, x4 . There are optimal values x4;opt for minimum ~t pd , e 4 . However, the optima locations and maximum COP, g , and Q II

vary with respect to the evaporator to condenser global heat transfer coefficient ratio, s. This issue is further investigated in Fig. 17b, in which only the optimal evaporator area fractions, x4;opt , and the obtained extrema are shown. The lower the values of s, the higher x4;opt will be, since the evaporator global heat transfer coefficient is smaller than the condenser one. Conversely, as s increases above unity, the U 4 > U 7 , therefore less area allocation is required by the evaporator for maximum performance, which is characterized by ~tpd;min and gII;max . Note that, as s decreases, more area is allocated to the evaporator to compensate for low global heat transfer coefficient, and the condenser rejects heat accordingly even with smaller area than the evaporator, due to high global heat transfer

554

T.K. Nunes et al. / Applied Energy 158 (2015) 540–555 10

s = 0.67 s=1 s=2

8

0.2

τ1,set = 0.97; R1234yf

~ t pd

(2)

0.15

η II

6

~

0.1 η II , Q 4

~ , t pd COP 4

(3) 0.05

COP

2

~ Q4

0 0

0.2

0.4

0.6

0.8

0

(4)

1

(5)

1

x4

(a) 4

τ1,set = 0.97; R1234yf

0.9 0.8 0.7

x 4,opt

0.6

~ t pd,min 3

0.5

~ t pd,min

x 4,opt

0.4

plant, and therefore minimum pull-down time and maximum second law efficiency, no matter how complicated the actual design may be; Since this optimization principle is general, the model could be used as a preliminary project tool, to locate x4;opt for any set of system design parameters that are different from Table 4, for maximum second law efficiency, and therefore optimal operation; It was shown that x4;opt is robust with respect to the variation of the type of refrigerant. This is an important finding for the purpose of scalability, and designing systems using different refrigerants; Changes in the optima location, x4;opt , are observed when the ratio of heat exchangers global heat transfer coefficients departs from 1. Therefore, an accurate assessment of the heat exchangers global heat transfer coefficients is mandatory for finding the appropriate optimal design, and The minimum pull-down time, and maximum second law efficiencies found with respect to the optimized system area allocation are sharp, stressing their importance for practical design, and therefore must be identified accurately in the quest for achieving VCR systems higher net efficiencies, and size reduction in order to make these systems more compact for applications in systems with high heat generating components, and less energy consuming to make them more commercially competitive.

0.3 0.2

η II, max

0.1 η II,max

Acknowledgements

0

2 0.5

1

1.5

2

s

(b) Fig. 17. The effect of the variation of evaporator to condenser global heat transfer rate ratio, s, on the minimization of the VCR system pull-down time, and maximization of second law efficiency, COP and dimensionless refrigeration rate with respect to evaporator area fraction, x4 , with R1234yf (a), and the obtained minimum pull-down time, and maximum second law efficiency as functions of s.

coefficient, therefore ~t pd;min drops and gII;max increases. On the other hand, as s increases, the evaporator global heat transfer increases, less area is allocated to the evaporator, but still keeping a large portion of the total system constrained area, and when s > 1 the condenser does not reject heat as effectively due to insufficient heat transfer area and low global heat transfer coefficient, therefore ~t pd;min increases and gII;max drops. 4. Conclusions The basic thermodynamics problem of how to obtain maximum exergy output rate from a VCR plant through optimal heat exchangers’ area allocation has been considered. A dynamic VCR system mathematical model was developed to obtain the system response in time and to calculate the system second law efficiency, as functions of operating and design parameters. Appropriate dimensionless groups were identified and the generalized results reported in dimensionless charts. Based on the results of Figs. 8–17, obtained for three different refrigerants, R12, R134a and R1234yf, the main conclusions of this study are summarized as follows: (1) There exists a fundamental optimal system heat exchangers area allocation, x4;opt , that characterizes the system such that maximum exergy output rate is obtained in the refrigeration

The Office of Naval Research, Contract N00014-14-1-0198. This work was also supported through the National Council for Scientific and Technological Development, CNPq, Brazil, Grant – 483467/2009-0. References [1] U.S. Energy Information Administration (EIA). Electric Power Monthly. Table 5.1; March 11, 2011. [2] U.S. Department of Energy, DOE. Buildings Energy Data Book, Section 2.1.5. 2011; 2010. [3] U.S. Department of Energy, DOE. Buildings Energy Data Book, Section 3.1.4, 2011; 2010. [4] Energy Information Administration (EIA). International Energy Outlook; 2010. [5] Sage AP. Systems engineering. New York: Wiley; 1992. [6] Dilay E, Vargas JVC, Souza JA, Ordonez JC, Yang S, Mariano AB. A volume element model (VEM) for energy systems engineering. Int J Energy Res 2015;39(1):46–74. http://dx.doi.org/10.1002/er.3209. [7] Catano J, Zhang T, Wenc JT, Jensen MK, Peles Y. Vapor compression refrigeration cycle for electronics cooling – Part I: Dynamic modeling and experimental validation. Int J Heat Mass Transf 2013;66:911–21. [8] Buzelin LOS, Amico SC, Vargas JVC, Parise JAR. Experimental development of an intelligent refrigeration system. Int J Refrig 2005;28:165–75. [9] Kanoglu M. Exergy analysis of multistage cascade refrigeration cycle used for natural gas liquefaction. Int J Energy Res 2002;26(8):763–74. [10] Shiba T, Bejan A. Thermodynamic optimization of geometric structure in the counterflow heat exchanger for an environmental control system. Energy 2001;26(5):493–512. [11] Sanaye S, Malekmohammadi HR. Thermal and economical optimization of air conditioning units with vapor compression refrigeration system. Appl Therm Eng 2004;24(13):1807–25. [12] Hovsapian R, Vargas JVC, Ordonez JC, Krothapalli A, Parise JAR, Berndsen JC. Thermodynamic optimization of a solar system for cogeneration of water heating and absorption cooling. Int J Energy Res 2008;32(13):1210–27. [13] Vargas JVC, Ordonez JC, Dilay E, Parise JAR. Modeling, simulation and optimization of a solar collector driven water heating and absorption cooling plant. Sol Energy 2009;83(8):1232–44. [14] Tracy T, Ordonez JC, Vargas JVC. First and second law thermodynamic analysis of a domestic scale trigeneration system. ASME 2007 energy sustainability conference; 2007. [15] Ordonez JC, Bejan A. Minimum power requirement for environmental control of aircraft. Energy 2003;28(12):1183–202. [16] Radcenco V, Vargas JVC, Bejan A, Lim JS. Two design aspects of defrosting refrigerators. Int J Refrig 1995;8(2):76–86.

T.K. Nunes et al. / Applied Energy 158 (2015) 540–555 [17] Khan JR, Zubair SM. Thermodynamic optimization of finite time vapor compression refrigeration systems. Energy Convers Manage 2001;42 (12):1457–75. [18] Chen L, Sun F, Wu C. Optimal allocation of heat-exchanger area for refrigeration and air-conditioning plants. Appl Energy 2004;77:339–54. [19] Vargas JVC, Sokolov M, Bejan A. Thermodynamic optimization of solar-driven refrigerators. J Solar Energy Eng – Trans ASME 1996;118(2):130–5. [20] Vargas JVC, Ordonez JC, Bejan A. Power extraction from a hot stream in the presence of phase change. Int J Heat Mass Transfer 2000;43(2):191–201. [21] Ordóñez JC, Vargas JVC, Bejan A. Combined power and refrigeration from a hot stream. Int J Thermodyn 1999;2(2):49–57. [22] Luo J, Chen L, Sun F, Wu C. Optimum allocation of heat transfer surface area for cooling load and COP optimization of a thermoelectric refrigerator. Energy Convers Manage 2003;44(20):3197–206. [23] Dingec H, Ìleri A. Thermoeconomic optimization of simple refrigerators. Int J Energy Res 1999;23:949–62. [24] Dhar M, Soedel W. Transient analysis of a vapor compression refrigeration system, Part I: The mathematical model. Venice, Italy: XV Int Congr Refrig; 1979. p. 1035–48. [25] Chi J, Didion D. A simulation model of the transient performance of a heat pump. Int J Refrig 1982;5(3):176–84. [26] MacArthur JW. Transient heat pump behavior: a theoretical investigation. Int J Refrig 1984;7(2):123–32. [27] Rajendran H, Pate M. A computer model of the start-up transients in a vapor compression refrigeration system. In: Preprints of the Int Inst Refrig Meeting, Purdue University, West Lafayette, IN, USA; 1986. p. 129–40. [28] Vargas JVC, Parise JAR. Simulation in transient regime of a heat pump with closed loop and on–off control. Int J Refrig 1995;18(4):235–43. [29] Estrada-Flores S, Cleland DJ, Cleland AC, James RW. Simulation of transient behaviour in refrigeration plant pressure vessels: mathematical models and experimental validation. Int J Refrig 2003;26:170–9. [30] Kim M, Kim MS, Chung JD. Transient thermal behavior of a water heater system driven by a heat pump. Int J Refrig 2004;27:415–21. [31] Sarbu I. A review on substitution strategy of non-ecological refrigerants from vapour compression-based refrigeration, air-conditioning and heat pump systems. Int J Refrig 2014;46:123–41.

555

[32] Arcaklioglu E, Çavusoglu A, Erisen A. An algorithmic approach towards finding better refrigerant substitutes of CFCs in terms of the second law of thermodynamics. Energy Convers Manage 2005;46:1595–611. [33] Zhao D, Ding G, Wu Z. Extension of the implicit curve-fitting method for fast calculation of thermodynamic properties of refrigerants in supercritical region. Int J Refrig 2009;32:1615–25. [34] Milián V, Navarro-Esbrí J, Ginestar D, Molés F, Peris B. Dynamic model of a shell-and-tube condenser. Analysis of the mean void fraction correlation influence on the model performance. Energy 2013;59:521–33. [35] Janssen MJ, de Witt JA, Kuijpers LJM. Cycling losses in domestic appliances: an experimental and theoretical analysis. Int J Refrig 1992;15(3):152–8. [36] Assawamartbunlue K, Brandemuehl MJ. The effect of void fraction models and heat flux assumption on predicting refrigerant charge level in receivers. In: Proceedings of international refrigeration and air conditioning conference. Paper 519. Purdue University, West Lafayette, IN, USA; 2000. [37] Kincaid D, Cheney W. Numerical analysis mathematics of scientific computing. 1st ed. Belmont, CA: Wadsworth; 1991. [38] Kim JH, Simon TW, Viskanta R. Journal of heat transfer policy on reporting uncertainties in experimental measurements and results [editorial]. J Heat Transfer 1993;115:5–6. [39] Lipschutz S, Lipson ML. Theory and problems of probability. 2nd ed. New York: McGraw-Hill; 2000. [40] Minkowycz WJ, Sparrow EM, Schneider GE, Pletcher RH. Handbook of numerical heat transfer. 2nd ed. Chapter 17. John Wiley & Sons Inc, New York; 2006. [41] Gatecliff GW, Lady ER. Explicit representation of the thermodynamic properties of refrigerants 12 and 22 in the superheat region. Purdue Compressor Tech Conf 1974:187–90. [42] McLinden MO, Klein SA, Lemmon EW, Peskin AP. NIST Standard Reference Database 23, NIST thermodynamic and transport properties of refrigerants and refrigerant mixtures – REFPROP version 7.0. Standard Reference Data Program, National Institute of Standards and Technology; 2002.