Abstract

A thermosyphon-based modular cooling approach offers an energy efficient cooling solution with an increased potential for waste heat recovery. Central to the cooling system is an air-refrigerant finned tube heat exchanger (HX), where air is cooled by evaporating refrigerant. This work builds on a previously published two-dimensional (2D) model for the finned-tube HX by updating and validating the model using in-house experimental data collected from the proposed system using R1233zd(E) as the working fluid. The results show that key system variables such as refrigerant outlet quality, air and refrigerant outlet temperatures, and exchanger duty agree within 20% of their experimental counterparts. The validated model is then used to predict the mean heat transfer coefficient on the refrigerant side for each tube in the direction of airflow, indicating a maximum heat transfer coefficient of nearly 1200 W/(m2 K) for a HX duty of 5.3 kW among the tested cases. The validated model therefore enables accurate predictions of HX performance and provides insights into improving the heat exchange efficiency and the corresponding system performance.

1 Introduction

Cooling and other secondary equipment in a data center (DC) consume an average of about 43% of total DC power [1]. These equipment, combined with information technology (IT) equipment, generate a significant amount of low-grade waste heat with little economic value. Moreover, no major improvement in DC non-IT power consumption has occurred since 2013 [1]. Therefore, novel approaches for DC thermal management and waste heat recovery are needed.

Data centers have been traditionally cooled by supplying cold air through an underfloor plenum [2,3] with ASHRAE recommended IT inlet air temperatures in the range of 18–27 °C [4]. The exhaust air temperatures are, in turn, below 50 °C, as recommended by major original equipment manufacturers. This scheme presents disadvantages in terms of the electric energy consumed by the cooling system and the discharge of heat to the ambient without any intermediate use. The waste heat is typically high in quantity but low in quality. Therefore, the potential to both (1) reduce the cooling power requirement of legacy DCs and (2) capture, transport, and reuse this heat for various purposes remains.

Furthermore, an increasing DC average rack power density over the past three years [5] represents a growing opportunity to switch to liquid cooling. Liquid cooling allows thermal engineers to operate with condensing temperatures in the 40–60 °C range, depending on the choice of a single or multiphase liquid cooling system [6,7] while maintaining chip case temperatures within reliability limits. The use of vapor-compression driven liquid cooling allows for condensing temperatures up to 90 °C [8]. A higher condensing temperature leads to effective waste heat harvesting, monetary benefits gained by a reduction in cooling system energy costs and further capital gain via waste heat recovery, and reduced carbon footprint via lesser CO2 emissions. These reductions in emissions further merit possible financial gain if a local carbon tax is in place.

Figure 1 depicts a novel system that not only efficiently cools and captures heat from hardware components in DCs and telecom central offices but also allows for heat recovery via a vapor recompression (i.e., heat pump) system. The equipment to be cooled (computer servers in this case) are contained in an air-cooled rack. The IT equipment is cooled by air that is continuously circulated within the cabinet, while the air itself is cooled via an air-to-refrigerant heat exchanger (HX) that is the evaporator of a thermosyphon loop. In particular, the evaporator is a finned tube HX that is widely used in refrigeration and air-conditioning applications. The refrigerant flows in a thermosyphon loop between a source and sink (evaporator and condenser, respectively). A reservoir tank (LV-separator) connects the two HXs, which physically separates the refrigerant liquid and vapor phases.

Fig. 1
Proposed thermal management and heat recovery system
Fig. 1
Proposed thermal management and heat recovery system
Close modal

The vapor refrigerant either can reject heat directly through a water-cooled condenser or be driven through a heat pump, thereby boosting its pressure (and temperature) for enhanced heat recovery. The high temperature condenser cooling water can then be used for various applications such as district heating, as service fluid within process industries, or as potable water in the same or colocated buildings. Alternatively, the condenser of Fig. 1 can be replaced or augmented with an absorption chiller to produce cold water, which can be supplied to an heating, ventilating, and air conditioning system or used as potable water in the summer months, after further treatment. The hot and/or cold water thus produced can be a source of revenue for the DC owner. If priced correctly, the incoming revenue can more than offset the cost of running the compressor, thereby generating a yearlong profit for the DC owner. The reader is referred to Ref. [9] for case studies detailing the applications and efficacy of the heat recovery system. This study, however, focuses on the evaporator cooling loop only and therefore assumes that the heat pump loop is not active since enhanced heat recovery adds complexity to the system and does not affect evaporator performance.

Figure 1 indicates that the evaporator as a finned-tube HX is a key component in the proposed cooling and waste heat recovery system. System performance in the form of thermal management of the rack-mounted IT equipment and waste heat recovery ability is strongly dependent upon evaporator performance. Figure 2 shows the pilot scale experimental system with key components labeled. The evaporator is located inside the cabinet directly behind the fan array (upstream of the cabinet fans) and cannot be seen without opening the cabinet. Our previous work [10] introduced an evaporator model for the system but did not include experimental validation nor a detailed analysis of heat transfer coefficient (HTC) variation within the HX. This study therefore extends the previous work by determining the appropriate HTC correlations to yield agreement with experimental measurements on a pilot-scale rig using R1233zd(E) as the working fluid (refrigerant), and then exploring the HTC variation within the evaporator for different HX duty values. HX model validation against data from a pilot-scale experiment is vital since model accuracy is extremely important for predicting overall system thermal performance.

Fig. 2
Proposed thermal management and heat recovery system. Evaporator is located directly behind the fan array to the left of the servers.
Fig. 2
Proposed thermal management and heat recovery system. Evaporator is located directly behind the fan array to the left of the servers.
Close modal

Finned-tube HXs have been widely studied over the past 50 years. Greater interest was initially placed in studying the airside due to its higher thermal resistance, thereby influencing the overall system performance. Eminent among these early studies are the works of McQuiston [11,12] and Rich [13] in quantifying the air-side performance of these HXs, in particular the heat transfer coefficient using Colburn j-factors and the pressure drop, evaluated separately for the tubes and fins. The work of Schmidt [14] and subsequently Hong and Webb [15] in modeling the continuous fin geometry as discrete hexagons and further approximating the geometry as circular/annular fins is of paramount importance. This allowed the determination of the fin efficiency with ease using trigonometric functions as opposed to more complex Bessel functions.

On the refrigerant side, flow boiling of refrigerants continues to be actively explored to date, with the works of Kandlikar [16,17], Shah [18], and Kays and London [19] standing out for their contribution to Nusselt number correlation development and for their comprehensive work on characterizing the performance of compact HXs.

The current model is based on the work of Qiao et al. [20] to predict the performance of a finned-tube HX. The HX in this study functions as an evaporator, designed to operate as part of a thermosyphon for a novel thermal management solution. Our model builds upon the Qiao et al. model by incorporating three new features lacking in the earlier model, which are used to predict the performance of the evaporator for a given set of inlet stream conditions. The model highlighted in this study is presented in detail in Ref. [10] and features the following novel aspects:

  • The new model uses the geometry of the heat exchanger to predict the refrigerant pressure drop in the inlet and outlet headers of the exchanger, primarily due to wall friction and sudden expansion and contraction of the refrigerant.

  • Treatment of two-phase pressure drop is different in the new model. First, the model artificially incorporates the momentum pressure drop on the refrigerant side by computing the difference in kinetic energy of the liquid and vapor phase, as opposed to using the discretized momentum equation. Second, the model determines the frictional pressure drop using the widely known and validated Müller-Steinhagen and Heck [21] correlation.

  • The new model incorporates the pressure drop on the airside as a means to judge whether the cabinet fans can maintain flow circulation inside the enclosed cabinet.

The advantages of creating a custom, detailed model lies in its ability to capture the flow physics and problem domain with a high level of accuracy. For example, in this case, the pressure drop through the headers can be accurately computed by assigning the exact geometry. Moreover, such user-created models also allow the flexibility of tailoring the correlations used to better match the flow physics, an ability that commercial software generally do not allow.

2 Evaporator Geometry

The evaporator is a finned-tube HX incorporated in a thermosyphon loop. The evaporator geometry as fabricated divides the HX into two identical finned tube sections, with three rows of tubes in each section, where each row is seven tubes wide. Thus, each section consists of 21 vertical, single-pass tubes with the entire evaporator containing 42 tubes. A computer-aided design model of the overall finned-tube HX (evaporator) is shown in Fig. 3, while Fig. 4 shows a close-up of the inlet manifold and tubing.

Fig. 3
Evaporator geometry
Fig. 3
Evaporator geometry
Close modal
Fig. 4
Evaporator inlet manifold with the rows and tubes in each row labeled
Fig. 4
Evaporator inlet manifold with the rows and tubes in each row labeled
Close modal

3 Modeling Strategy

The evaporator geometry is modeled using several assumptions for the air and refrigerant sides. On the refrigerant side:

  • No flow maldistribution between tube bundles, tube rows, and individual tubes

  • Incoming flow splits equally among all tubes

  • Modeling takes place in the yz-plane

  • Only one tube modeled from each row

  • Each tube in a given row faces the same conditions inside and outside

  • Viscous heating is neglected through the inlet and outlet headers. Thus, the refrigerant temperature and enthalpy entering the tubes are the same as that entering the HX

On the airside:

  • Incoming air mass flow rate splits equally among all tubes in a given row

  • No mixing on the airside takes place in the z-direction due to the presence of fins

  • Mixing in the y-direction due to back flow does not impact results—same stream packet impacts the next refrigerant cell

  • No significant temperature gradients exist in the x-direction

  • The output from one row of tubes becomes the input for the next row of tubes

The above assumptions simplify the geometry from three-dimensional to two-dimensional (2D). The resulting geometry can then be depicted as shown in Fig. 5. The simplified geometry is used for further analysis of the evaporator. It is advantageous to use this geometry since it makes the resulting computation easier and faster while still preserving the flow physics, and therefore the model predictive capabilities, as shown later in this paper.

Fig. 5
Simplified evaporator geometry and applied grid
Fig. 5
Simplified evaporator geometry and applied grid
Close modal

A segment-by-segment approach is utilized to solve the simplified 2D geometry. A grid can then be overlaid on the geometry as shown in Fig. 5. Each hatched square is a segment and forms the basic unit on which the transport equations for the air and refrigerant side are evaluated. The height of each rectangular cell (z-direction) is chosen to ensure an integer number of fins are captured in each grid cell. The cells occupying the region between the two tube bundles do not take part in the simulation and the air temperature is taken to be constant as it flows across those cells.

The modeling strategy is presented in detail in an earlier work by the authors [10] and is summarized below for completeness. An exception is the calculation method for the refrigerant void fraction, which is now determined using the well-known Rouhani–Axelsson correlation [22] for vertical flow [23], which is based on the drift flux model. It is understood to be accurate in predicting the gas volume fraction (void fraction) [24]
(1a)
where
(1b)
and
(1c)
The homogenous void fraction,γhom, is determined by setting the slip ratio equal to 1 in the standard definition of void fraction, thereby yielding
(1d)

The void fraction is used in determining the density of the two-phase mixture, which in-turn is used to determine the momentum pressure drop during phase-change, which contributes to the total refrigerant-side pressure drop.

The general procedure utilizes the equations in Table 1 and assumes that the system is at steady-state. The refrigerant-side equations are discretized using either a second-order or the upwind differencing scheme, depending on the term, and solved on a back-staggered grid to decouple velocity from pressure. The continuity equation acts as a check for the mass flowrate of the refrigerant, while the momentum and energy equations yield the pressure and enthalpy (p and h¯) at the centroid of each cell, which are the primary state variables. Other state variables such as temperature can then be readily determined.

Table 1

Summary of key equations in evaporator model

DescriptionEquationNotes
Air-side energy balancema˙cp,adTadyΔy=αa(Ao,tube+ηfinAo,fin)(TwallTa)Fin efficiency is based on annular representation of hexagonal fins
Air-side mass balancem˙adωadyΔy=αm(Ao,tube+ηfinAo,fin)min(0,ωw,satωa)αm=αacp,aLe2/3 is the mass transfer coefficient, where Le2/30.90 is the Lewis number
ωw,sat is the humidity ratio of saturated air, evaluated at Twall
Tube wall temperaturem˙a[cp,a(Ta,inTa,out)+(ωa,inωa,out)Δhfg]αrAh(TwTr)=0Tube wall assumed thin
Equation based on a surface energy balance on the tube wall
Ah is surface area based on heated perimeter, i.e., tube outer diameter
Refrigerant mass balanceddz(GA)=0G is mass flux
Refrigerant momentum balanceddz(G2Aρ¯m)=Adpdzτ¯wPρ¯gAsinθThe frictional pressure drop is determined using correlation [21] and results from the wall shear stress,τ¯w
Refrigerant energy balanceddz(Gh¯A)=Phq¯+Aq¯GgGAsinθThe volumetric heat generation term, q¯G, is zero for this application
Frictional pressure dropΔPfric=12ρfum2LDiAcum=4m˙ρπDi2Single-phase (Moody friction factor using Colebrook formula)
Equations (1)(7) in Ref. [21]
ΔPfric=(dPdL)LAc
Two-phase
Momentum pressure dropΔPmom=G2{[(1x̂)2ρf(1γ)+x̂2ρgγ]out[(1x̂)2ρf(1γ)+x̂2ρgγ]in}For two-phase (zero for single-phase)
Average momentum density1ρ¯m=1ρgx2γ+1ρf(1x2)(1γ)x is flow quality (gas weight-rate fraction)
Average densityρ¯=ρgγ+ρf(1γ)Found using property routines for both single and two-phase flow
Average density-weighted enthalpyh¯ρ=hgx̂+hf(1x̂)x̂ is static quality (gas mass-fraction)
Average flow-weighted enthalpyh¯=hgx+hf(1x)This is definition but h¯ is determined by solving the energy balance
Refrigerant pressure dropΔpfriction1ϕ=0.5fρum2LAc/Dhum=4m˙/(πρDh2)Δpfriction2ϕ Müller-Steinhagen and Heck [21]
ΔPmomentum1ϕ=0ΔPmomentum2ϕ=G2Ac(AB)A=xi2ρvγi+(1xi)2ρl(1γi)B=xi12ρvγi1+(1xi1)2ρl(1γi1)(γi10)B=1ρl(γi1=0)
Darcy friction factor, f, using Colebrook Formula
For two-phase frictional pressure drop, flow length is equal to cell height; refrigerant quality is taken at the exit of each cell
Difference in kinetic energy between liquid and vapor phase gives rise to momentum pressure drop
Static pressure drop obtained from discretized momentum equation
i represents cell number within any tube row
Δpstatic=ρ¯igAcsinθL(i=1)
0.5(ρ¯i+ρ¯i1)gAcsinθL(i>1)
Airside pressure dropRich [13] method for pressure drop over fins
Zukauskas and Ulinskas [25] method for pressure drop over tubes
Total pressure drop is the sum of tube and fin pressure drop
Manifold pressure dropSudden expansion and contraction K-factors from Ref. [26], Fig. 6.3(a)
Horizontal section same as core section (as outlined in refrigerant pressure drop row)
DescriptionEquationNotes
Air-side energy balancema˙cp,adTadyΔy=αa(Ao,tube+ηfinAo,fin)(TwallTa)Fin efficiency is based on annular representation of hexagonal fins
Air-side mass balancem˙adωadyΔy=αm(Ao,tube+ηfinAo,fin)min(0,ωw,satωa)αm=αacp,aLe2/3 is the mass transfer coefficient, where Le2/30.90 is the Lewis number
ωw,sat is the humidity ratio of saturated air, evaluated at Twall
Tube wall temperaturem˙a[cp,a(Ta,inTa,out)+(ωa,inωa,out)Δhfg]αrAh(TwTr)=0Tube wall assumed thin
Equation based on a surface energy balance on the tube wall
Ah is surface area based on heated perimeter, i.e., tube outer diameter
Refrigerant mass balanceddz(GA)=0G is mass flux
Refrigerant momentum balanceddz(G2Aρ¯m)=Adpdzτ¯wPρ¯gAsinθThe frictional pressure drop is determined using correlation [21] and results from the wall shear stress,τ¯w
Refrigerant energy balanceddz(Gh¯A)=Phq¯+Aq¯GgGAsinθThe volumetric heat generation term, q¯G, is zero for this application
Frictional pressure dropΔPfric=12ρfum2LDiAcum=4m˙ρπDi2Single-phase (Moody friction factor using Colebrook formula)
Equations (1)(7) in Ref. [21]
ΔPfric=(dPdL)LAc
Two-phase
Momentum pressure dropΔPmom=G2{[(1x̂)2ρf(1γ)+x̂2ρgγ]out[(1x̂)2ρf(1γ)+x̂2ρgγ]in}For two-phase (zero for single-phase)
Average momentum density1ρ¯m=1ρgx2γ+1ρf(1x2)(1γ)x is flow quality (gas weight-rate fraction)
Average densityρ¯=ρgγ+ρf(1γ)Found using property routines for both single and two-phase flow
Average density-weighted enthalpyh¯ρ=hgx̂+hf(1x̂)x̂ is static quality (gas mass-fraction)
Average flow-weighted enthalpyh¯=hgx+hf(1x)This is definition but h¯ is determined by solving the energy balance
Refrigerant pressure dropΔpfriction1ϕ=0.5fρum2LAc/Dhum=4m˙/(πρDh2)Δpfriction2ϕ Müller-Steinhagen and Heck [21]
ΔPmomentum1ϕ=0ΔPmomentum2ϕ=G2Ac(AB)A=xi2ρvγi+(1xi)2ρl(1γi)B=xi12ρvγi1+(1xi1)2ρl(1γi1)(γi10)B=1ρl(γi1=0)
Darcy friction factor, f, using Colebrook Formula
For two-phase frictional pressure drop, flow length is equal to cell height; refrigerant quality is taken at the exit of each cell
Difference in kinetic energy between liquid and vapor phase gives rise to momentum pressure drop
Static pressure drop obtained from discretized momentum equation
i represents cell number within any tube row
Δpstatic=ρ¯igAcsinθL(i=1)
0.5(ρ¯i+ρ¯i1)gAcsinθL(i>1)
Airside pressure dropRich [13] method for pressure drop over fins
Zukauskas and Ulinskas [25] method for pressure drop over tubes
Total pressure drop is the sum of tube and fin pressure drop
Manifold pressure dropSudden expansion and contraction K-factors from Ref. [26], Fig. 6.3(a)
Horizontal section same as core section (as outlined in refrigerant pressure drop row)

Solving the refrigerant-side equations requires determining the algebraic variables, comprising of the refrigerant mass flowrate,m˙, densities,ρ¯,ρ¯m, and the average density-weighted enthalpy, h¯ρ. The refrigerant mass flowrate is determined by solving the refrigerant mass balance.

Determination of the densities and enthalpies is phase-dependent. For the single-phase region, both densities and enthalpies collapse into one and can easily be determined using property routines. All fluid properties are determined using REFPROP [27]. For the two-phase region, the values of the two state variables upstream of the current node are used for determining ρ¯ and in addition the saturation temperature,Tsat. The saturation temperature is then used to determine the saturated densities and enthalpies ρl,ρg,hl,andhg. The saturated enthalpies can be compared to the current enthalpy,h¯, which then yields the relative enthalpy (static quality). The static quality and densities at saturated state are then used to solve for the void fraction γ using equation set (1). The value of the average density-weighted enthalpy, h¯ρ, at the current node can then be computed based on the static quality through the relevant relation outlined in Table 1.

Finally, the average momentum density,ρ¯m, contributes to the acceleration term of the pressure drop and is not computed since the discretized momentum equation requires knowing ρ¯m at the next cell [10]. Rather, it is determined from first principles as the difference in kinetic energy between the two-phases per the expression provided in Table 1. Note that the momentum pressure drop has no contribution under single-phase flow. The frictional component of pressure drop results from the wall shear stress, which is difficult to compute directly and is generally determined via empirical correlations. Hence, the widely known correlation from Müller-Steinhagen and Heck [21] is used to determine the contribution of the frictional pressure drop under two-phase flow and is presented in Table 1.

The heat transfer coefficients on both the air and refrigerant-side are determined using known correlations, which are discussed in detail in Secs. 6 and 7. Similarly, the pressure drop for each fluid is determined as per Table 1. With the values of the algebraic variables, HTCs and pressure drop established for each cell, the two state variables could be determined using the momentum and energy equation for the refrigerant-side.

For the airside, the mass balance determines the amount of moisture lost, via the humidity ratio,ω, while the energy balance determines the temperature rise of the air as it flows across the evaporator. The reader is encouraged to refer to Ref. [10] for specific details of the calculation scheme.

4 Airside Sensor Calibration

The air temperature and velocity sensors were calibrated in a small wind tunnel against a known velocity before installation and testing. The blower motor was connected to a variable frequency drive and set to a predetermined frequency to achieve a known air velocity. In addition, the air temperatures were also measured to provide a reference for comparison with sensor temperature measurements. The calibration study found the sensors to read lower than the expected velocity, except at the lowest tested speed of 0.5 m/s, where the sensors agreed with the set velocity. As a result, individual calibration curves for each sensor were plotted, yielding the corrected/actual velocity as a function of the sensor velocity. A linear trend fit the data well (R2 > 0.99 for all ten sensors). The calibration curve for one of the sensors on the evaporator airside is provided in Fig. 6.

Fig. 6
Calibration curve for sensor at point 7a
Fig. 6
Calibration curve for sensor at point 7a
Close modal

5 Performance Measurement

Air temperature and velocity measurements were made at four distinct locations (state points 7a, 8a, 9a, and 10a per Fig. 7) across the evaporator using air temperature and velocity sensors. These sensors consist of a hot-wire anemometer and a thermocouple in a low-profile shape to measure the air velocity and temperature at the inlet and exit of the HX (evaporator) under study. Air temperature and velocity measurements were collected simultaneously for all sensors, sampled at a frequency of 0.2 Hz for all 15 data points.

Fig. 7
Relative position of airside sensors
Fig. 7
Relative position of airside sensors
Close modal

Measurements of the air pressure at various points inside the cabinet were made using a handheld barometric pressure recorder before any tests were ran. The air pressure inside the cabinet was found to be at most 0.2 mm Hg higher than the room ambient. Thus, the cabinet pressure for any test was taken to be the ambient pressure plus the small difference noted above.

The airside relative humidity was measured using a handheld hygrometer at the outlet of the load banks. The location was chosen based on practicality so that the display could be placed outside the enclosed cabinet while the measuring wand remained inside. Repeated measurements of the air's relative humidity at the load bank exhaust were manually recorded during each test, and at steady-state, were found to be nearly constant.

The three measured psychrometric properties (T, p, and RH) were used to determine the humidity ratio of the air at the load bank outlet. Moreover, the lowest air temperature inside the cabinet is measured at the evaporator outlet, and since that did not drop below the dew point, no moisture should have been removed in the evaporator. This was also confirmed while running the system, where no moisture was seen collecting below the evaporator. Hence, the humidity ratio stayed constant throughout the cabinet and could be used as the third psychrometric parameter instead of the relative humidity for processing of the airside data.

Refrigerant temperature and pressure measurements were taken at the inlet of the evaporator, while only temperature was measured at the evaporator outlet in the vertical section of the riser. Refrigerant flowrate in the evaporator was measured on the downcomer using an ultrasonic flowmeter. The relative locations of the instrumentation used on the refrigerant and waterside are depicted in Fig. 8, while instrumentation specifications are detailed in Table 2. For the evaporator, the state points of interest are 1r and 2r only, specifying the evaporator inlet and outlet, respectively.

Fig. 8
Cooling system flow loop
Fig. 8
Cooling system flow loop
Close modal
Table 2

Instrumentation specifications

InstrumentLocation per Fig. 8 Measurement principleMeasuring rangeMeasurement accuracy
PDURight of cabinet wall 3c (inside cabinet)Power meterUp to 8.6 kW±3% of reading +1 for last digit
Thermocouple (K-type)1r–5r, 1w–2w, 1c–3c, 1tnVoltage difference between measured and reference junction0–1038 °C±0.375% of reading
Pressure transmitter1r, 3r, 5r (1w: pressure gauge)Strain gauge (piezoresistive ceramic diaphragm)0–150 psig0.5% of best fit straight line, i.e., end value
Air velocity sensorInlet and exit of load banks (three each) and evaporator (two each)Hot-wire anemometer0–5 m/s• Between 15–35 °C
±0.015 m/s or ±3% of reading (greater value)
• Beyond 35 °C
increases by ±0.25% per deg and ±0.005 m/s
Air temp. sensorSame as air velocity sensor (each air sensor measures temp. and velocity)Voltage difference0–70 °C±1 °C (>0.5 m/s)
±1.5 °C (<0.5 m/s)
Ultrasonic flowmeter1rLiquid velocity using ultrasoundUp to 300 L/min
1.25 DIN pipe
±2% of reading [28]
Magnetic flowmeter1wFaraday's law (induced voltage fluid velocity in a magnetic field)0.03–6.6 gpm±(2% MV + 0.5% EV)
InstrumentLocation per Fig. 8 Measurement principleMeasuring rangeMeasurement accuracy
PDURight of cabinet wall 3c (inside cabinet)Power meterUp to 8.6 kW±3% of reading +1 for last digit
Thermocouple (K-type)1r–5r, 1w–2w, 1c–3c, 1tnVoltage difference between measured and reference junction0–1038 °C±0.375% of reading
Pressure transmitter1r, 3r, 5r (1w: pressure gauge)Strain gauge (piezoresistive ceramic diaphragm)0–150 psig0.5% of best fit straight line, i.e., end value
Air velocity sensorInlet and exit of load banks (three each) and evaporator (two each)Hot-wire anemometer0–5 m/s• Between 15–35 °C
±0.015 m/s or ±3% of reading (greater value)
• Beyond 35 °C
increases by ±0.25% per deg and ±0.005 m/s
Air temp. sensorSame as air velocity sensor (each air sensor measures temp. and velocity)Voltage difference0–70 °C±1 °C (>0.5 m/s)
±1.5 °C (<0.5 m/s)
Ultrasonic flowmeter1rLiquid velocity using ultrasoundUp to 300 L/min
1.25 DIN pipe
±2% of reading [28]
Magnetic flowmeter1wFaraday's law (induced voltage fluid velocity in a magnetic field)0.03–6.6 gpm±(2% MV + 0.5% EV)

The voltage and current signals generated by the instrumentation employed on the refrigerant and waterside were acquired using a National Instrument data acquisition system. Based on test conditions, it took the overall system about 40–45 min to reach statistically steady-state conditions, following which all data channels were scanned simultaneously at a frequency of 1 Hz for a duration of 15 min. A total of 15 data points were collected over a 45 day period, spanning a cabinet heat load range between 2 and 7.5 kW in increments of 2 kW. An arithmetic average of the acquired data samples for each heat load tested was used to process the data.

It should be noted that temperature measurements were obtained indirectly as the difference in voltage (generated by a thermocouple) between the state point of interest and an ice bath, taken as the reference voltage, and converted to a temperature reading via a ninth-order polynomial [29]. This is the reason for specifying the thermocouple uncertainty as a percentage rather than a value in Table 2.

6 Data Reduction

6.1 Air Temperature and Flowrate.

The airside volumetric flowrate is obtained at both cross section of the evaporator (inlet and outlet) by taking the product of the area-weighted average velocity, Eq. (2), and total cross-sectional area, as given by Eq. (3). The average air temperature at each cross section is determined by the arithmetic mean of the individual temperature values at that cross section, as per Eq. (4). This is reasonable given that each cross section is divided into segments of equal area and the velocities (and by extension mass flowrate) through each section are similar. Temperature and velocity readings at state points 7a and 8a are used for the evaporator inlet. Similarly, readings at state points 9a and 10a are used for the evaporator outlet
(2)
(3)
(4)
Finally, the mass flowrate of the air inside the cabinet is determined as per Eq. (5), where the average density at a given cross section is determined using the average temperature at that cross section, and the air pressure and the humidity ratio inside the cabinet
(5)

6.2 Evaporator Duty.

The electric power supplied to the load banks (used in place of the computer servers in Fig. 2) is dissipated as heat to the air moving inside the cabinet, due to the resistive nature of the load banks. Assuming negligible losses through the evaporator's shrouding, the heat gained by the air from the load banks is completely transferred to the refrigerant as the air flows across the evaporator. The refrigerant in turn transfers that heat to the cooling water as it flows through the condenser. Since the intermediate tubing, separator tank, and condenser are well insulated, there is negligible heat loss to the environment through these components. Therefore, the heat gained by the waterside is taken as a measure of the heat lost by the air, and by extension, the evaporator duty, determined as per the following equation:
(6)

Using the waterside has the advantage that it is easy to instrument and measure and, therefore, provides the least uncertainty among the three fluids (as compared to instrumenting the airflow or refrigerant flow inside the enclosed cabinet). In comparison, the large uncertainty in air temperature measurements compared to its temperature change across the evaporator leads to a large uncertainty in evaporator duty predictions.

6.3 Evaporator Outlet Quality.

The refrigerant quality at the evaporator outlet cannot be measured physically and therefore has to be calculated indirectly using measured variables. Given the size of the evaporator and the physical limitations of the enclosed cabinet, the differential pressure drop across the evaporator was not measured and only a temperature reading was obtained at the evaporator outlet, as shown in Fig. 8.

Assuming the temperature to be saturated, the corresponding saturation enthalpies can be readily determined via REFPROP [27]. The actual enthalpy of the refrigerant at the evaporator outlet is determined as the sum of the refrigerant enthalpy at the evaporator inlet and the heat lost by the air, which is equal to the waterside heat gain per the argument presented in Sec. 6.2. A comparison of the calculated actual enthalpy at the evaporator outlet versus the two saturation enthalpies yields the relative enthalpy, more commonly known as the quality (mass-based). This is determined as per equation set (7)
(7a)
(7b)

A value between 0 and 1 mathematically confirms that the fluid is a two-phase mixture and validates the saturation temperature assumption. In addition, an observation of the refrigerant state via a sight glass installed in the riser visually confirms that the fluid at the evaporator outlet is indeed a two-phase mixture, as discussed in the section on model validation. Finally, the quality (mass fraction) cannot be visually gauged, but the void fraction (volume fraction) can be calculated via equation set (1) and directly judged via visual observations of the vapor space in the riser.

7 Results and Uncertainty Analysis

The processed data for the measured and calculated variables are presented in Tables 3 and 4, representing evaporator inlet and outlet conditions, respectively. Note that all values in Tables 3 and 4 are quoted to three significant figures. The last two columns of Table 4 are provided for reference for validating the arguments presented in Sec. 6.2.

Table 3

Evaporator input performance data

AirRefrigerant
Data pointTemp. (°C)Pressure (bar)Rel. hum (%)Flow (kg/s)Temp. (°C)Pressure (bar)Flow (kg/s)
137.51.0133.50.51921.61.400.596
243.01.0125.40.58227.31.650.674
350.91.0118.00.52833.32.000.750
441.41.0127.60.57427.21.660.747
535.51.0141.60.51122.31.440.716
636.10.99832.80.48422.31.380.731
736.50.99333.40.53021.91.360.703
832.71.0043.10.47224.21.430.555
934.81.0040.20.47920.91.310.673
1041.11.0032.80.55725.21.520.678
1149.71.0022.90.53431.91.870.704
1250.20.99923.20.53031.91.870.660
1341.60.99832.70.60325.61.570.650
1434.80.99842.30.52118.71.280.572
1531.80.99942.30.35220.61.360.487
AirRefrigerant
Data pointTemp. (°C)Pressure (bar)Rel. hum (%)Flow (kg/s)Temp. (°C)Pressure (bar)Flow (kg/s)
137.51.0133.50.51921.61.400.596
243.01.0125.40.58227.31.650.674
350.91.0118.00.52833.32.000.750
441.41.0127.60.57427.21.660.747
535.51.0141.60.51122.31.440.716
636.10.99832.80.48422.31.380.731
736.50.99333.40.53021.91.360.703
832.71.0043.10.47224.21.430.555
934.81.0040.20.47920.91.310.673
1041.11.0032.80.55725.21.520.678
1149.71.0022.90.53431.91.870.704
1250.20.99923.20.53031.91.870.660
1341.60.99832.70.60325.61.570.650
1434.80.99842.30.52118.71.280.572
1531.80.99942.30.35220.61.360.487
Table 4

Evaporator output performance data

Refrigerant
Data pointAirTemp. (°C)Temp. (°C)Quality (%)HXDuty (kW)Load bankPower input (kW)Percentage differenceLoad bank versus HX
130.221.83.263.903.961.55
234.127.54.085.335.757.37
338.832.35.366.577.4812.1
432.427.43.625.285.666.65
528.022.62.723.973.872.49
628.022.62.684.033.932.77
729.422.12.693.843.901.69
828.725.11.241.931.920.80
927.122.02.233.823.820.06
1032.026.03.585.275.738.10
1138.032.24.686.437.4113.2
1238.632.05.056.357.4214.4
1333.225.94.015.235.687.90
1427.619.33.123.853.850.10
1526.521.11.681.881.922.0
Refrigerant
Data pointAirTemp. (°C)Temp. (°C)Quality (%)HXDuty (kW)Load bankPower input (kW)Percentage differenceLoad bank versus HX
130.221.83.263.903.961.55
234.127.54.085.335.757.37
338.832.35.366.577.4812.1
432.427.43.625.285.666.65
528.022.62.723.973.872.49
628.022.62.684.033.932.77
729.422.12.693.843.901.69
828.725.11.241.931.920.80
927.122.02.233.823.820.06
1032.026.03.585.275.738.10
1138.032.24.686.437.4113.2
1238.632.05.056.357.4214.4
1333.225.94.015.235.687.90
1427.619.33.123.853.850.10
1526.521.11.681.881.922.0

A first-order uncertainty analysis is then applied to characterize the propagation of errors from the measured values to the final results, as per the approach outlined by ASME [30] (based primarily on the work of Ref. [31]). The expanded uncertainty in the measured variables are based on instrument uncertainties specified in Table 2. Further, the uncertainty values for various thermodynamic and transport properties used for data reduction are stated in Table 5. The expanded uncertainty range was determined for each value for the high-impact variables based on a 95% confidence interval, and are provided in Table 6. Note that since the evaporator duty is taken as the condenser waterside heat gain, the uncertainties associated with the waterside measurements and heat gain are also included as part of Table 6.

Table 5

Uncertainty in thermodynamic and transport properties of R1233zd(E) obtained via REFPROP

PropertySourceUncertainty
Temperature (T)New data obtained to determine Helmholtz equation of state [32]±6 mK
Density (ρ)New Helmholtz equation of state fitted to empirical data [32]0.020% of value
Vapor pressure (Pv)0.223% of value
Speed of sound (c)0.131% of value
Thermal conductivity (k)New correlation for thermal conductivity based on experimental data [33]Liquid: 1% of value
Vapor: 3% of value (P < 1 MPa)
Viscosity (μ)Unpublished data from personal communications of authors in Ref. [27]Liquid: 4% of value
243–433 K up to 40 MPa
Vapor: 4% of value (P < 1 MPa)
Specific heat (cp)REFPROP database [27]
Considering low GWP fluids, e.g., R1234yf, R1234ze(E), and traditional R134a
Assume 1%
Specific enthalpy (h)Assume same as specific heat
Δh=cpavgΔT
Assume same as cp or cpavg
Surface tension (σ)New experimental data for low GWP refrigerants [34]0.14 mN/m
PropertySourceUncertainty
Temperature (T)New data obtained to determine Helmholtz equation of state [32]±6 mK
Density (ρ)New Helmholtz equation of state fitted to empirical data [32]0.020% of value
Vapor pressure (Pv)0.223% of value
Speed of sound (c)0.131% of value
Thermal conductivity (k)New correlation for thermal conductivity based on experimental data [33]Liquid: 1% of value
Vapor: 3% of value (P < 1 MPa)
Viscosity (μ)Unpublished data from personal communications of authors in Ref. [27]Liquid: 4% of value
243–433 K up to 40 MPa
Vapor: 4% of value (P < 1 MPa)
Specific heat (cp)REFPROP database [27]
Considering low GWP fluids, e.g., R1234yf, R1234ze(E), and traditional R134a
Assume 1%
Specific enthalpy (h)Assume same as specific heat
Δh=cpavgΔT
Assume same as cp or cpavg
Surface tension (σ)New experimental data for low GWP refrigerants [34]0.14 mN/m
Table 6

Uncertainty range for key variables (95% confidence interval)

VariableUncertainty
Air temperature±1.4 °C
Air mass flowrate±(0.04–0.07) kg/s
Refrigerant temperature±(0.1–0.2) °C
Refrigerant pressure±0.1 bar
Refrigerant mass flowrate±(0.02–0.03) kg/s
Refrigerant outlet quality±(1.7–1.9)%
Evaporator duty±(0.3–0.4) kW
Load bank power input±(0.1–0.6) kW
Water temperature±0.3 °C
Water mass flowrate±(0.04–0.07) kg/s
Condenser duty (waterside)±(0.1–0.9) kW
VariableUncertainty
Air temperature±1.4 °C
Air mass flowrate±(0.04–0.07) kg/s
Refrigerant temperature±(0.1–0.2) °C
Refrigerant pressure±0.1 bar
Refrigerant mass flowrate±(0.02–0.03) kg/s
Refrigerant outlet quality±(1.7–1.9)%
Evaporator duty±(0.3–0.4) kW
Load bank power input±(0.1–0.6) kW
Water temperature±0.3 °C
Water mass flowrate±(0.04–0.07) kg/s
Condenser duty (waterside)±(0.1–0.9) kW

The uncertainty in airside pressure readings was not determined since there is little to no variation in this parameter. Similarly, uncertainty in the airside relative humidity was not determined since it is only used to fix the air state (determine humidity ratio) and since no condensation occurs inside the evaporator, the humidity ratio remains constant throughout the cabinet for all data points. Hence, these variables are deemed to be low-impact and specific uncertainty values not required which may impact further results.

Finally, the reason for the relatively high uncertainty in the quality values is due to the propagation of error resulting from the outlet quality being determined from the inlet enthalpy, refrigerant mass flowrate through the evaporator and evaporator duty, as opposed to the outlet enthalpy being determined directly based on measured temperature and pressure at the evaporator outlet.

8 Prevalidation

8.1 Pressure Correction.

The pressure at the physical outlet of the evaporator was determined prior to validating the model with the experimental data collected from the test system. The calculated outlet pressure was based on the difference between the saturation pressure at the measured state point 2r (Fig. 8) and the end of the outlet header of the evaporator, 2r. The piping upstream of state point 2r involves a vertical copper pipe 12.5 in (31.75 cm) in length terminating in a 90 deg bend, followed by a 2-in (5.08 cm) horizontal pipe stub, as shown in Fig. 9. The nominal tube diameter for this section of the riser is 1.5 inch (3.81 cm), as indicated in Fig. 8.

Fig. 9
Relative location of state point 2r to evaporator outlet
Fig. 9
Relative location of state point 2r to evaporator outlet
Close modal
Fig. 10
Evaporator model grid resolution study
Fig. 10
Evaporator model grid resolution study
Close modal

The pressure difference was computed using the Chisholm correlation of the Lockhart–Martinelli approach [35]. To determine the friction factor in the straight sections of the copper tubing, a value of 1.5 μm [36] was used for the surface roughness. A value of 0.36 was considered for the single-phase loss factor for the 90 deg bend [37]. The pressure drop between points 2r and 2r came out to be about 0.2 psi (1.1 kPa), which is small (about 5%) compared to the actual pressure drop in the evaporator computed as the difference in pressure between state points 1r and 2r.

8.2 Grid Independence Study.

A grid resolution study for the evaporator was performed based on the first three data points collected, as shown in Fig. 10. The figure shows that 2254 grid cells are sufficient for model validation against the collected experimental data. Results for the predicted refrigerant and air outlet temperatures, evaporator outlet quality and exchanger duty are presented in Sec. 9, Validation and Results. It should be noted that since the airside results are not used for further analysis, the exchanger duty value for the evaporator is taken as the heat exchange value for the waterside of the low-pressure condenser.

Fig. 11
Predicted and experimental refrigerant outlet temperature using Kandlikar's correlation
Fig. 11
Predicted and experimental refrigerant outlet temperature using Kandlikar's correlation
Close modal
Fig. 12
Predicted and experimental air outlet temperature using Kandlikar's correlation
Fig. 12
Predicted and experimental air outlet temperature using Kandlikar's correlation
Close modal

9 Validation and Results

The model is validated using existing flow boiling correlations found in the literature. In particular, the correlations from Refs. [17] and [3840] are used. These correlations are selected based on the following reasons:

  • The correlations capture the flow physics accurately. The first three are flow boiling correlations and include the combined effect of nucleate and forced convection boiling, while the last one (Cooper's correlation) is based only on pool boiling

  • The first three correlations have been validated to a large set of experimental data (>5000 data points)

  • Shah's correlation is recent (year 2022) and includes R1233zd(e) in its experimental database, which is the working fluid for this study

  • Kandlikar's correlation includes a fluid dependent parameter, Ffl, to extend the correlation to fluids other than the base fluid, water, which was assigned a fluid factor of 1.0

  • Cooper's correlation has been validated in several pool boiling studies and captures the physics of the pool boiling process by inclusion of the surface roughness

Results for the refrigerant-side using Kandlikar's correlation with a Ffl value of 2.2 are presented in Figs. 1114. The statistics outlined in these figures are explained as follows:

  • n is the number of data points used for model comparison

  • |δ| is the mean absolute error, computed as
    (8)
  • λ is the number of points that fall within the ±20% band, expressed as a percentage

For each of the four plots, 12 out of the 15 data points collected experimentally are used for purposes of model validation. This is done intentionally to keep the same number of data points for validating the evaporator and condenser models, since the condenser model could only be simulated for 12 out of the 15 data points2. Data points 8, 14, and 15 of Tables 3 and 4 were not considered for the evaporator model validation.

As seen from the plots in Figs. 1114, the model can predict the refrigerant and air outlet temperatures within a 20% band for all data points. This indicates proper operation of the core solver for the evaporator model. The evaporator outlet quality and exchanger duty, however, are based on the refrigerant (and airside) HTCs, which are computed using empirical correlations obtained from published literature. Note that the heat exchanger duty is captured well, albeit for two data points at the lower end of the duty values, which are slightly off. Further, the HX duty follows the trend of the air outlet temperature, i.e., a less than ideal prediction for the air temperature results in a higher than ideal prediction for the duty and vice versa, since the HX duty is proportional to the airside temperature difference, with the inlet air temperature being constant and coming from experiment. This leaves the refrigerant outlet quality, which is predicted well for 2/3 data points using Kandlikar's correlation with a Ffl value of 2.2; however, it shows the largest scatter. An initial guess of 1.5 for the fluid factor is obtained by determining the heat transfer coefficient using Cooper's correlation, as per Kandlikar's [16] suggestion, and then using a hit-and-trial approach in its vicinity till the mean absolute error is reduced to a minimum. This yields the fluid factor value given above.

In an attempt to improve the HX duty and refrigerant quality prediction, the correlations by Shah and Liu and Winterton were subsequently applied, with both yielding the same result in terms of prediction accuracy and mean absolute error. Results for the same four variables using Shah's correlation are presented in Figs. 1518.

As can be seen from Figs. 1518, air and refrigerant outlet temperature are again predicted with 100% accuracy, while the HX duty is again predicted accurately for 5/6 data points. However, the mean absolute error improves marginally for the air temperature and exchanger duty. The HX duty again follows the same trend as the air outlet temperature, as explained above. The refrigerant outlet quality shows the same scatter as before; however, it is predicted with lower accuracy, as compared to Kandlikar's correlation with Ffl of 2.2 (Fig. 13).

Fig. 13
Predicted and experimental refrigerant outlet quality using Kandlikar's correlation
Fig. 13
Predicted and experimental refrigerant outlet quality using Kandlikar's correlation
Close modal
Fig. 14
Predicted and experimental heat exchanger duty using Kandlikar's correlation
Fig. 14
Predicted and experimental heat exchanger duty using Kandlikar's correlation
Close modal
Fig. 15
Predicted and experimental refrigerant outlet temperature using Shah's correlation
Fig. 15
Predicted and experimental refrigerant outlet temperature using Shah's correlation
Close modal
Fig. 16
Predicted and experimental air outlet temperature using Shah's correlation
Fig. 16
Predicted and experimental air outlet temperature using Shah's correlation
Close modal
Fig. 17
Predicted and experimental refrigerant outlet quality using Shah's correlation
Fig. 17
Predicted and experimental refrigerant outlet quality using Shah's correlation
Close modal

The low quality obtained in the experiment (<10%) and the high HTCs predicted by the three flow boiling correlations used above suggested the lack of flow boiling and the dominance of pool boiling in the experimental setup of Fig. 2. As a result, the widely used pool boiling correlation from Ref. [40] was next attempted to predict the experimental data, in particular, to better predict the refrigerant quality. Cooper's correlation incorporates the surface roughness of the tubing through which the fluid flows, thus better capturing the physics of the pool boiling process. The surface roughness for copper tubing generally lies between 1 and 2 μm, and an average value of 1.5 μm [36] was employed for analysis. Results using Cooper's correlation are presented in Figs. 1922.

Fig. 18
Predicted and experimental heat exchanger duty using Shah's correlation
Fig. 18
Predicted and experimental heat exchanger duty using Shah's correlation
Close modal
Fig. 19
Predicted and experimental refrigerant outlet temperature using Cooper's correlation
Fig. 19
Predicted and experimental refrigerant outlet temperature using Cooper's correlation
Close modal
Fig. 20
Predicted and experimental air outlet temperature using Cooper's correlation
Fig. 20
Predicted and experimental air outlet temperature using Cooper's correlation
Close modal
Fig. 21
Predicted and experimental refrigerant outlet quality using Cooper's correlation
Fig. 21
Predicted and experimental refrigerant outlet quality using Cooper's correlation
Close modal
Fig. 22
Predicted and experimental heat exchanger duty using Cooper's correlation
Fig. 22
Predicted and experimental heat exchanger duty using Cooper's correlation
Close modal

The plots in Figs. 1922 show that the refrigerant and air outlet temperature are still predicted with 100% accuracy within a ±20% band and the same absolute error as with Kandlikar's correlation. However, the HX duty is now predicted with 100% accuracy (has the same error as with Kandlikar's correlation) but comes at a loss of accuracy in predicting the refrigerant quality, which is heavily under-predicted (values lie below the solid line).

The result of Fig. 21 nullifies the assumption that nucleate boiling alone is contributing to the flow boiling process in the evaporator; rather, there is a vapor core in the copper tube where flow boiling is taking place at the liquid–vapor interface, and the solid–liquid interface at the tube's inner surface aids the nucleate boiling process. This hypothesis is confirmed by the images presented in Fig. 23, which are taken at point 2r in Fig. 8. The figure shows an increase in vapor fraction with increasing evaporator heat exchange as expected.

Fig. 23
Vapor core generated at different heat loads. Vapor region is designated by dashed line.
Fig. 23
Vapor core generated at different heat loads. Vapor region is designated by dashed line.
Close modal

The surface roughness of the tube is the primary catalyst for this nucleate boiling to take place, providing nucleation sites for bubble growth and departure. Given that Cooper's correlation only considers the pool boiling (nucleate boiling) process, the low quality prediction of Fig. 21 reveals that an additional boiling phenomenon is taking place, which is convective boiling (due to forced convection of the flow). The lower quality results from a lower prediction of the HTC using Cooper's correlation, which results in lower HX duty and higher airside outlet temperatures. Due to this lower duty prediction, Cooper's correlation is therefore able to capture the two data points for the HX duty at the cost of under-predicting the refrigerant outlet quality.

In summary, the Shah, Kandlikar, and Liu and Winterton correlations capture the physics of the boiling process in the evaporator, with Kandlikar's correlation with Ffl of 2.2 yielding the best prediction for all four variables, followed equally well by the Shah and Liu and Winterton correlations. This is attributed to the fact that Kandlikar's correlation assigns a factor to the nucleate boiling term, which in this case, enhances the impact of nucleate boiling (Ffl>1), while the Liu and Winterton correlation assigns a suppression factor, S, to the nucleate boiling term (S < 1) and an enhancement factor, E, to the convective boiling term (E > 1). Moreover, although Cooper's correlation predicts the HX duty and fluid outlet temperatures with 100% accuracy, it severely under-predicts the refrigerant outlet quality and hence cannot be reliably used for future design or prediction purposes.

Further, Kandlikar's correlation best meets the simulation convergence criteria for all data points (converges with least error), closely followed by Shah's correlation and then by Liu and Winterton's correlation. The simulation is deemed to converge if the total heat lost by the airside is within 0.20% of the total heat gained by the refrigerant side. The heat gain and loss are computed using two methods: aggregated on a per cell basis (cells shown in Fig. 5) and overall for the HX using each fluid's inlet and outlet temperature and flowrate (Q˙=m˙cpΔT). The higher of the two differences is compared against the converge criteria for a conservative comparison.

Finally, Kandlikar's correlation with Ffl of 2.2 may not predict 100% of the data points accurately, but the mean absolute error is 0.7%, meaning that on average, the exit quality is in fact predicted with an error of 0.007. Given that the maximum exit quality is close to 5% (very low), the prediction is remarkably accurate, with 2 out of every 3 points (67%) lying within ±20% of their experimental counterpart, and one data point being borderline accurate. If the exit quality were to be within the 30–60% range, the margin for error would be significantly larger, with a higher probability that a greater number of points would fall within a ±20% band (possibly closer to 100%). In effect, the exit quality predictions are reasonable enough for future use given the low experimental exit quality values and the high uncertainty associated with them. This conclusion is important since the model is accurate yet neglects any effects of refrigerant maldistribution, and accounting for refrigerant maldistribution greatly enhances the complexity and computational cost of the model.

10 Evaporator Performance Analysis

Using the validated model with Kandlikar's correlation, the performance of the evaporator on each fluid's side is analyzed. On the refrigerant-side, the fluid flows from the bottom of the separator tank through the downcomer into a horizontal section of tubing, as shown in Fig. 8, eventually entering the evaporator as a subcooled liquid. The liquid splits through the six rows of tubing in the evaporator as it flows upward, while picking up heat from the airside and changing phase along the way. The heat transfer coefficients achieved on a per tube basis are plotted in Figs. 2427 for select HX duty values. In addition, the boiling and convection number on a per tube basis for the maximum heat load (6.35 kW, data point 12) are plotted in Figs. 28 and 29, respectively. These figures reveal the following about the refrigerant flow through the evaporator:

Fig. 24
Predicted HTC by row for 1.9 kW HX duty
Fig. 24
Predicted HTC by row for 1.9 kW HX duty
Close modal
Fig. 25
Predicted HTC by row for 3.9 kW HX duty
Fig. 25
Predicted HTC by row for 3.9 kW HX duty
Close modal
Fig. 26
Predicted HTC by row for 5.3 kW HX duty
Fig. 26
Predicted HTC by row for 5.3 kW HX duty
Close modal
Fig. 27
Predicted HTC by row for 6.4 kW HX duty
Fig. 27
Predicted HTC by row for 6.4 kW HX duty
Close modal
Fig. 28
Predicted boiling number by row for 6.4 kW HX duty
Fig. 28
Predicted boiling number by row for 6.4 kW HX duty
Close modal
Fig. 29
Predicted convection number by row for 6.4 kW HX duty
Fig. 29
Predicted convection number by row for 6.4 kW HX duty
Close modal
  • Single-phase HTCs are steady, while two-phase HTCs are increasing as the fluid progresses downstream.

  • Single-phase HTCs are lower than their two-phase counterparts for each tube. This allows the system designer to take advantage of the higher HTCs afforded by the phase-change process.

  • The minimum length of the single-phase region is close to 30% of the overall tube length for the flowrate achieved (Table 3). This fraction of overall tube length occurs for the maximum heat load case and increases with decreasing HX duty.

  • For a given heat load, refrigerant stays in single-phase the longest for tubes in row 6. This is because row 6 has the lowest heat flux due to precooled air from upstream. Row 1 has the shortest length of a single-phase region because the row contains the highest heat flux.

  • Both convective and nucleate boiling contribute to the boiling mechanism (Figs. 28 and 29). Nucleate boiling dominates, as evident by the steady increase in boiling number along the length of the tube. Convective boiling dominates at the onset of saturation, but its contribution rapidly diminishes to a minimum yet steady amount, as evident by the plot of Fig. 29. The boiling number is on the order of 104 to 105 for the range of heat loads tested. The convection number is larger than 0.65 for all heat loads (chosen as threshold for predominantly nucleate boiling region by Kandlikar for correlation development [17])

  • A potential for cooling higher heat loads exists by improving the evaporator design, leading to a lower refrigerant flowrate to achieve a saturated state achieved sooner (i.e., a shorter single-phase length). This result would take greater advantage of the higher HTCs afforded by phase change.

On the airside, the overall heat transfer coefficient, which varies between 55 and 60 W/(m2 K) for the heat loads studied, is calculated using the Colburn j-factor approach. McQuiston and Parker's correlation [41] is used as the default for a four row HX. The fin efficiency is calculated using the Bessel function approach from Ref. [36]. The effective airside HTC is lower than the overall HTC by about 8% due to a less-than-ideal fin efficiency for the given fin geometry. The effective airside HTC,h¯eff, is related to the overall HTC,h¯, as per the following equation:
(9)

where Afseg denotes the surface area of fins per segment in m2 and Awseg is the surface area of tubes between fins per segment in m2.

The overall evaporator design heat load was 10 kW with air inlet and exit temperatures of 50 °C and 40 °C, respectively. The relatively high temperatures were chosen for enhancing heat recovery. The cell-level heat flux is low (below 2 W/cm2) for the highest heat load (6.4 kW). The highest minimum HTC is achieved for an evaporator duty of 5.3 kW per Fig. 30 and decreases slightly for increasing heat loads.

Fig. 30
Maximum HTC achieved versus HX duty
Fig. 30
Maximum HTC achieved versus HX duty
Close modal

11 Conclusions and Future Work

The advantage of creating an in-depth, custom model lies in its ability to capture the flow physics and problem domain with a high level of accuracy. For example, in this case, the heat transfer coefficient for the air and refrigerant side can be determined as the fluid flows through and across the heat exchanger, respectively. In addition, the pressure drop through the headers can be accurately computed by assigning the exact geometry. Moreover, such user-created models also allow the flexibility of tailoring the correlations used to better match the flow physics. The validated model is used to predict the heat transfer coefficients achieved in the evaporator.

It was determined that the evaporator outlet quality is low enough to dissipate significantly higher heat loads before reaching dry-out, provided the evaporator (and condenser) are adequately sized with a sufficient flow of air and water through them. Furthermore, it was found that nucleate boiling dominates over convective boiling of the working fluid in this study. The maximum HTC appears to occur at a duty of 5.3 kW of the data points studied, or approximately half of the HX design capacity.

Finally, areas of further improvement include better instrumentation of the airside to accurately capture temperature and velocity measurements (e.g., Refs. [42] and [43]) and an accurate measurement of the evaporator refrigerant side pressure drop needs to be made using a differential pressure sensor (e.g., Ref. [44]).

Acknowledgment

Special thanks to Mr. S. G. Schon of QuantaCool Corp. for designing the evaporator and for permission to use Figs. 24 upon modification. The authors would also like to thank Mr. Ken Allen of Super Radiator Coils for fabricating the evaporator and providing design parameters, which was instrumental in validating the model.

Funding Data

  • NSF IUCRC Award No. IIP-1738782. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of the National Science Foundation. Directorate for Engineering (Funder ID: 10.13039/100000084).

Data Availability Statement

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

Nomenclature

A =

area (m2)

A =

air (used to denote state point on airside)

Ac =

cross-sectional area based on wetted perimeter (m2)

Ah =

cross-sectional area based on heated perimeter (m2)

c =

cabinet (enclosure to mount information technology equipment)

c =

speed of sound (m/s)

cp =

specific heat capacity (J/kg K)

Co =

constant

D =

diameter

Dh =

hydraulic diameter (m)

E =

enhancement factor

f =

Darcy (Moody) friction factor (dimensionless)

Ffl =

fluid dependent parameter

g =

acceleration due to gravity (m/s2)

G =

mass flux (kg/(m2s))

h =

specific enthalpy (J/kg)

h¯ =

average flow-weighted enthalpy (J/kg)

h¯ρ =

average density-weighted enthalpy (J/kg)

Δhfg =

enthalpy of vaporization (J/kg)

k =

thermal conductivity (W/m K)

L =

length (denotes flow length) (m)

Le =

Lewis number

m˙ =

mass flow rate (kg/s)

n, N =

number of data points

p =

pressure (Pa)

Ph =

heated perimeter (m)

Q =

flowrate (m3/s)

Q˙ =

heat flow (W)

q¯G =

average rate of volumetric heat generation (W/m3)

q¯ =

average heat flux (W/m2)

r =

refrigerant (used to denote state point on refrigerant side)

Ref =

refrigerant (system working fluid)

RH =

relative humidity

S =

suppression factor

Sat =

saturated

T =

temperature

tn =

tank (liquid–vapor separator tank)

ugi =

gas/vapor velocity at inlet (m/s)

um =

mean velocity (m/s)

V =

velocity (m/s)

w =

water (used to denote state point on waterside)

x =

flow quality (vapor weight-rate fraction)

x̂ =

static quality (vapor mass fraction)

Greek Symbols
α =

thermal diffusivity (m2/s)

αa =

heat transfer coefficient (W/m2 K)

αm =

mass transfer coefficient (K s)

γ =

void fraction

|δ| =

mean absolute error

ηfin =

fin efficiency

θ =

angle of flow or tube to horizontal (deg)

λ =

number of data points within a specified range

μ =

dynamic viscosity, Pa s

ρ =

density (kg/m3)

ρ¯ =

average density (kg/m3)

ρl =

density of liquid phase in a two-phase mixture (kg/m3)

ρv =

density of vapor phase in a two-phase mixture (kg/m3)

ρ¯m = =

average momentum density (kg/m3)

σ =

surface tension (N/m)

τw =

wall shear stress (Pa)

ω =

humidity ratio

Subscripts
a =

air

avg =

average

eff =

effective

evap =

evaporator

f =

saturated liquid

fin =

fin-side

g =

saturated vapor

h =

heated

hom =

homogenous

i =

index of current cell, inner (refers to diameter; used in Table 1, single-phase pressure drop)

in =

inlet

l =

liquid

m =

mean

o,out =

outside/outlet

R, r =

refrigerant

sat =

saturation

seg =

segment

tube =

tube-side

v =

vapor

w, wtr =

water

wall =

wall

1ϕ =

single-phase (fluid is either all liquid or all vapor)

2ϕ =

two-phase (fluid exists as a liquid–vapor mixture)

Superscript
X¯ (overbar) =

average value for parameter X

Abbreviations
ASHRAE =

American Society for Heating, Refrigerating, and Air Conditioning Engineers

DC =

data center

EV =

end value

HTC =

heat transfer coefficient

HX =

heat exchanger

IT =

information technology

LV =

liquid–vapor

MV =

measured value

PDU =

power distribution unit

Footnotes

2

The lowest heat loads were challenging to simulate since the simulation would either not converge or crash due to the property routine package (REFPROP) unable to determine properties at a given state point. This is due to the higher error in measured values at the lower loads; measurement accuracy and repeatability increases at higher loads due to the values of the measured variables falling in the middle range of the sensor.

References

1.
Bizo
,
D.
,
Ascierto
,
R.
,
Lawrence
,
A.
, and
Davis
,
J.
,
2021
, “
Uptime Institute Global Data Center Survey 2021
,” UII-511, Uptime Institute, New York.
2.
Tradat
,
M.
,
Mohsenian
,
G.
,
Manaserh
,
Y.
,
Sammakia
,
B.
,
Mendo
,
D.
, and
Alissa
,
H. A.
,
2020
, “
Experimental Analysis of Different Measurement Techniques of Server-Rack Airflow Predictions Towards Proper DC Airflow Management
,”
Proceedings of the 2020 19th IEEE Intersociety Conference on Thermal and Thermomechanical Phenomena in Electronic Systems (ITherm)
,
Orlando, FL
, July 21–23, pp.
366
373
.10.1109/ITherm45881.2020.9190584
3.
Tradat
,
M. I.
,
Mohammad
,
Y.
,
Manaserh
,
A.
,
Sammakia
,
B. G.
,
Hoang
,
C. H.
, and
Alissa
,
H. A.
,
2021
, “
An Experimental and Numerical Investigation of Novel Solution for Energy Management Enhancement in Data Centers Using Underfloor Plenum Porous Obstructions
,”
Appl. Energy
,
289
, p.
116663
.10.1016/j.apenergy.2021.116663
4.
ASHRAE
,
2015
,
Thermal Guidelines for Data Processing Environments
, 4th ed.,
ASHRAE
,
Atlanta, GA
.
5.
Miller
,
R.
,
2020
, “
Rack Density Keeps Rising at Enterprise Data Centers
,” Data Center Frontier, accessed June 4, 2020, https://datacenterfrontier.com/rack-density-keeps-rising-at-enterprise-data-centers/
6.
Fujitsu Limited
,
2015
, “
Fujitsu Cool-Central Liquid Cooling Technology
,”
Presented at the ISC15
, Frankfurt, Germany, July 12–16.https://www.fujitsu.com/dk/imagesgig5/fujitsu-coolcentral-liquid-cooling-technology.pdf
7.
Wu
,
D.
,
Marcinichen
,
J. B.
, and
Thome
,
J. R.
,
2013
, “
Experimental Evaluation of a Controlled Hybrid Two-Phase Multi-Microchannel Cooling and Heat Recovery System Driven by Liquid Pump and Vapor Compressor
,”
Int. J. Refrig.
,
36
(
2
), pp.
375
389
.10.1016/j.ijrefrig.2012.11.011
8.
Marcinichen
,
J. B.
,
Olivier
,
J. A.
, and
Thome
,
J. R.
,
2012
, “
On-Chip Two-Phase Cooling of Datacenters: Cooling System and Energy Recovery Evaluation
,”
Appl. Therm. Eng.
,
41
, pp.
36
51
.10.1016/j.applthermaleng.2011.12.008
9.
Khalid
,
R.
,
Schon
,
S. G.
,
Ortega
,
A.
, and
Wemhoff
,
A. P.
,
2019
, “
Waste Heat Recovery Using Coupled 2-Phase Cooling Heat-Pump Driven Absorption Refrigeration
,”
Proceedings of the 2019 18th IEEE Intersociety Conference on Thermal and Thermomechanical Phenomena in Electronic Systems (ITherm)
, Las Vegas, NV, May 28–31, pp.
684
692
.10.1109/ITHERM.2019.8757465
10.
Khalid
,
R.
,
Amalfi
,
R. L.
, and
Wemhoff
,
A. P.
,
2021
, “
Rack-Level Thermosyphon Cooling and Vapor-Compression Driven Heat Recovery: Evaporator Model
,” ASME Paper No. IPACK2021-73269. 10.1115/IPACK2021-73269
11.
McQuiston
,
F. C.
,
1980
, “
Finned Tube Heat Exchangers: State of the Art for the Air Side
,”
Presented at the 5th Annual Heat Pump Technology Conference
,
Oklahoma State University
,
Stillwater, OK
, Apr. 14–15.https://ui.adsabs.harvard.edu/abs/1980ahpt.confR..14M/abstract
12.
McQuiston
,
F. C.
,
1978
, “
Heat, Mass and Momentum Transfer Data for Five Plate-Fin-Tube Heat Transfer Surfaces
,”
ASHRAE Trans.
,
84
(
1
), pp.
266
293
.https://www.techstreet.com/standards/at-2486-rp-155-heat-mass-and-momentum-transfer-data-for-five-plate-fin-tube-heat-transfer-surfaces?product_id=1854047
13.
Rich
,
D. G.
,
1973
, “
The Effect of Fin Spacing on the Heat Transfer and Friction Performance of Multi-Row, Smooth Plate-Fin-and-Tube Heat Exchangers
,”
ASHRAE Trans.
,
79
(
2
), pp.
137
145
.https://books.google.co.in/books/about/The_Effect_of_Fin_Spacing_on_the_Heat_Tr.html?id=XcKErgEACAAJ&redir_esc=y
14.
Schmidt
,
T. E.
,
1945
, “
La Production Calorifique Des Surfaces Munies D'ailettes
,”
Annexe Du Bulletin De L'Institut International Du Froid, Annexe G-5
,
International Institute of Refrigeration (IIR)
,
Paris, France
.
15.
Hong
,
K. T.
, and
Webb
,
R.
,
1996
, “
Calculation of Fin Efficiency for Wet and Dry Fins
,”
HVACR Res.
,
2
(
1
), pp.
27
41
.10.1080/10789669.1996.10391331
16.
Kandlikar
,
S. G.
,
1983
, “
An Improved Correlation for Predicting Two-Phase Flow Boiling Heat Transfer Coefficient in Horizontal and Vertical Tubes
,”
Heat Exchangers for Two-Phase Applications
,
Seattle, WA, July 24
.https://www.semanticscholar.org/paper/An-improvedcorrelation-for-predicting-two-phase-in-Kandlikar/9a000f116abba6014b3c00dbba565306d12109a0
17.
Kandlikar
,
S. G.
,
1990
, “
A General Correlation for Saturated Two-Phase Flow Boiling Heat Transfer Inside Horizontal and Vertical Tubes
,”
ASME J. Heat Mass Transfer-Trans. ASME
,
112
(
1
), pp.
219
228
.10.1115/1.2910348
18.
Shah
,
M. M.
,
1979
, “
A General Correlation for Heat Transfer During Film Condensation Inside Pipes
,”
Int. J. Heat Mass Transfer
,
22
(
4
), pp.
547
556
.10.1016/0017-9310(79)90058-9
19.
Kays
,
W. M.
, and
London
,
A. L.
,
1984
,
Compact Heat Exchangers
, 3rd ed.,
McGraw-Hill
,
New York
.
20.
Qiao
,
H.
,
Aute
,
V.
, and
Radermacher
,
R.
,
2015
, “
Transient Modeling of a Flash Tank Vapor Injection Heat Pump System—Part I: Model Development
,”
Int. J. Refrig.
,
49
, pp.
169
182
.10.1016/j.ijrefrig.2014.06.019
21.
Müller-Steinhagen
,
H.
, and
Heck
,
K.
,
1986
, “
A Simple Friction Pressure Drop Correlation for Two-Phase Flow in Pipes
,”
Chem. Eng. Process.: Process Intensif.
,
20
(
6
), pp.
297
308
.10.1016/0255-2701(86)80008-3
22.
Rouhani
,
S. Z.
, and
Axelsson
,
E.
,
1970
, “
Calculation of Void Volume Fraction in the Subcooled and Quality Boiling Regions
,”
Int. J. Heat Mass Transfer
,
13
(
2
), pp.
383
393
.10.1016/0017-9310(70)90114-6
23.
Albertsen
,
B.
, and
Schmitz
,
G.
,
2021
, “
Experimental Parameter Studies on a Two-Phase Loop Thermosyphon Cooling System With R1233zd(E) and R1224 yd(Z)
,”
Int. J. Refrig.
,
131
, pp.
146
156
.10.1016/j.ijrefrig.2021.07.036
24.
Woldesemayat
,
M. A.
, and
Ghajar
,
A. J.
,
2007
, “
Comparison of Void Fraction Correlations for Different Flow Patterns in Horizontal and Upward Inclined Pipes
,”
Int. J. Multiphase Flow
,
33
(
4
), pp.
347
370
.10.1016/j.ijmultiphaseflow.2006.09.004
25.
Zukauskus
,
A.
, and
Ulinskas
,
R.
,
1998
, “
Banks of Plain and Finned Tubes,” Heat Exchanger Design Handbook
,
G. F.
Hewitt
, ed.,
Begell House
,
New York
.
26.
Shah
,
R. K.
, and
Sekulić
,
D. P.
,
2003
,
Fundamentals of Heat Exchanger Design
,
Wiley
,
Hoboken, NJ
.
27.
Lemmon
,
E. W.
,
Bell
,
I. H.
,
Huber
,
M. L.
, and
McLinden
,
M. O.
,
2018
, “
NIST Standard Reference Database 23: Reference Fluid Thermodynamic and Transport Properties-REFPROP
,” National Institute of Standards and Technology, Standard Reference Data Program, Gaithersburg, MD.
28.
Frei
,
M.
,
Hischier
,
I.
,
Deb
,
C.
,
Sigrist
,
D.
, and
Schlueter
,
A.
,
2021
, “
Impact of Measurement Uncertainty on Building Modeling and Retrofitting Decisions
,”
Front. Built Environ.
,
7
.10.3389/fbuil.2021.675913
29.
Omega Engineering
,
Omega Temperature Measurement Handbook
, 6th ed.
30.
ASME
,
2006
, “
Test Uncertainty
,”
ASME
,
New York
, Report No. PTC-19.1.
31.
Coleman
,
H. W.
, and
Steele
,
W. G.
,
2009
,
Experimentation, Validation, and Uncertainty Analysis for Engineers
, 3rd ed.,
Wiley
,
Hoboken, NJ
.
32.
Mondéjar
,
M. E.
,
McLinden
,
M. O.
, and
Lemmon
,
E. W.
,
2015
, “
Thermodynamic Properties of Trans -1-Chloro-3,3,3-Trifluoropropene (R1233zd(E)): Vapor Pressure, (p, ρ, T) Behavior, and Speed of Sound Measurements, and Equation of State
,”
J. Chem. Eng. Data
,
60
(
8
), pp.
2477
2489
.10.1021/acs.jced.5b00348
33.
Perkins
,
R. A.
,
Huber
,
M. L.
, and
Assael
,
M. J.
,
2017
, “
Measurement and Correlation of the Thermal Conductivity of Trans -1-Chloro-3,3,3-Trifluoropropene (R1233zd(E))
,”
J. Chem. Eng. Data
,
62
(
9
), pp.
2659
2665
.10.1021/acs.jced.7b00106
34.
Kondou
,
C.
,
Nagata
,
R.
,
Nii
,
N.
,
Koyama
,
S.
, and
Higashi
,
Y.
,
2015
, “
Surface Tension of Low GWP Refrigerants R1243zf, R1234ze(Z), and R1233zd(E)
,”
Int. J. Refrig.
,
53
, pp.
80
89
.10.1016/j.ijrefrig.2015.01.005
35.
Ghiaasiaan
,
S. M.
,
2007
, “
Two-Phase Flow, Boiling and Condensation
,”
Conventional and Miniature Systems
,
Cambridge University Press
,
Cambridge, UK
.
36.
Incropera
,
F. P.
, and
DeWitt
,
D. P.
, eds.,
2007
,
Fundamentals of Heat and Mass Transfer
, 6th ed.,
Wiley
,
Hoboken, NJ
.
37.
Munson
,
B. R.
,
Young
,
D. F.
, and
Okiishi
,
T. H.
,
2002
,
Fundamentals of Fluid Mechanics
, 4th ed.,
Wiley
,
New York
.
38.
Shah
,
M. M.
,
2022
, “
New General Correlation for Heat Transfer During Saturated Boiling in Mini and Macro Channels
,”
Int. J. Refrig.
,
137
, pp.
103
116
.10.1016/j.ijrefrig.2022.02.019
39.
Liu
,
Z.
, and
Winterton
,
R.
,
1991
, “
A General Correlation for Saturated and Subcooled Flow Boiling in Tubes and Annuli, Based on a Nucleate Pool Boiling Equation
,”
Int. J. Heat Mass Transfer
,
34
(
11
), pp.
2759
2766
.10.1016/0017-9310(91)90234-6
40.
Cooper
,
M. G.
,
1984
, “
Saturation Nucleate Pool Boiling—A Simple Correlation
,”
First U.K. National Conference on Heat Transfer
, Vol.
86
,
Elsevier
,
Leeds, UK
, pp.
785
793
.
41.
McQuiston
,
F. C.
, and
Parker
,
J. P.
,
1994
,
Heating Ventilating and Air-Conditioning-Analysis and Design
,
Wiley
,
New York
.
42.
Degree Controls
,
2022
, “
°C Grid Multi-Point Air Measurement System
,” Degree Controls, Nashua, NH, accessed Dec. 31, 2022, https://www.degreec.com/products/airflow-tools-instruments/airflow-measurement-instrumentation/c-grid/
43.
Degree Controls
,
2022
, “
°C Grate Air Velocity and Temperature Measurement
,” Degree Controls, Nashua, NH, accessed Dec. 31, 2022, https://www.degreec.com/products/airflow-tools-instruments/airflow-tools/c-grate/
44.
Setra Systems
,
2022
, “
Model 230 Wet-to-Wet Differential Water Pressure Transducer
,” Setra Systems, Boxborough, MA, accessed Dec. 31, 2022, https://www.setra.com/product/pressure/model-230