Experimental and Numerical Testing of Heat Pump Evaporator

Abstract

When designing a heat pump evaporator, it is necessary to use correlations that ensure small deviations of the designed and realized process parameters for specific input data. The aim of the work is to propose a suitable mathematical model for the physical process in the tubular evaporator of the heat pump. The applicability of the proposed mathematical model in the design of the heat pump was evaluated by comparing the results obtained from the experimental tests of the tubular evaporator of the heat pump with the numerical results obtained from the application of the proposed mathematical model. For the experimental tests, a tubular evaporator was made and 10 measuring points were set up, where the process parameters were measured (temperature and pressure drop of the working media R134a and water). Theoretical results were obtained by dividing the evaporator into control volumes and solving the corresponding system of equations of the proposed mathematical model using the Runge-Kutta and Adams Moulton predictor-corrector method. As an independent parameter, the water temperature at the inlet to the evaporator was varied in the range of 10 °C to 18 °C. The test results show that the largest deviation of the calculated and measured water temperature is +0.41 °C to −0.58 °C, while the refrigerant temperature is +0.43 °C to + 0.52 °C. The largest deviation of the evaporator thermal capacity based on the calculations and experimental tests is +9.39% to −6.31%. Based on the obtained results, it is possible to recommend the use of the proposed mathematical model for the design of the tubular evaporator of a heat pump.

1. Introduction

The evaporator plays an important role in the design and operation of the heat pump, as it has a great impact on the energy efficiency of the system. The structure of the evaporator determines the type of flow both on the side of the refrigerant and the cooled medium. The type of flow, in turn, affects the pressure drop of the media and the intensity of heat transfer. In addition to the design of the evaporator, the throttling rate, compression, and combination of well water condition indicators determine the rate of evaporation of the refrigerant flowing through the evaporator. The evaporator environment consists of the thermal expansion valve, compressor, and wells (heat sinks).

A variety of steady-state and transient mathematical models have been developed in the literature to describe heat pump systems under various conditions and with various working fluids that are more or less neglected.

MacArthur [1], MacArthur and Grald [2], Nyers and Stojan [3], Mithraratne et al. [4], and Jong Won Choi et al. [5] have developed a transient mathematical model with split parameters. The description of the steady state was given by Koury et al. [6], Belman et al. [7], Kassai [8], and Nyers [9]. Zhao Lei et al. [10], Chi and Didion [11], Welsby et al. [12], Nyers [13], Santa [14], and Santa [15] have developed a mathematical model with lumped parameters to study the behavior of the evaporator.

Esbri et al. [16] present a model of a shell and tube evaporator using R1234yf and R134a as working fluids. The model uses the effectiveness NTU method to predict the evaporating pressure and the refrigerant and secondary fluid temperatures at the evaporator outlet. The evaporator geometry, the refrigerant mass flow rate, the secondary fluid volume flow rate, and the evaporator inlet temperature are inputs.

Kaddoura et al. [17] present mathematical modeling of a “Convection-Enhanced Evaporator (CEE) system”. Comparison between evaporation rates predicted by the model and experimental results shows good agreement; for 11 of 16 operating conditions, the error was within 10%, for 4 conditions within 20%, and for 1 condition within 35%. The effects of the different operating conditions on evaporation performance are also investigated.

Abadi and Bahrami [18] investigate the feasibility of a combined evaporator and condenser. Tests were conducted with a water vapor pressure of 0.61 −5.63 kPa and a heat transfer fluid (HTF) inlet temperature of 0–35 °C. Comparison of tubes with different fins a fin height of 1.42 mm and 40 parallel fins per inch (FPI), have a higher overall heat transfer coefficient as evaporators (40 W/K), and the tubes with 0.96 mm fin height (40 FPI) slightly outperform the other tubes as condensers (47 W/K).In the paper by Yadav et al. [19], the approximate solution of the nonlinear dynamic energy model of a multiple effect evaporator (MEE) using Fourier series and metaheuristics is presented. Before solving the dynamic model, the nonlinear steady-state model is solved to obtain the optimal steady-state process parameters. The optimization goal was to find the best estimates of unknown coefficients in the Fourier series expansion using two preeminent metaheuristic approaches: Particle swarm optimization and harmony search. The optimization results show that both approaches provide minimal constraint violation, proving their efficiency in solving such complex nonlinear energy models

Lin et al. [20] investigated the counterflow dew point evaporative cooler, which provides a much better approach to air cooling instead of the conventional vapor compression cooler. The model was able to accurately predict the product air temperature, cooling effectiveness, cooling capacity, and COP with a maximum deviation of ±5.0%. The main results of this study showed that the transient responses of the duct plate and cooler were in good agreement with an exponential decay function.

Salman et al. [21] investigated a dynamic model of a domestic refrigerator developed in this study and used it to analyze both steady-state and hot-load cooling operations. The results show that the power consumption value is well captured, with the average power consumption deviating by less than 10% in all tests. Power consumption is within ±5% for all tests. For load processing, the water temperature deviation did not exceed 1.2 °C throughout the cooling period.

Santos et al. [22] propose a rigorous model to represent the operation of an air-cooled chiller of a large office building and some simple correlations to predict the global heat exchange coefficient and the pressure drop of the evaporator. The rigorous model shows that large variations on the order of 10% are observed in the prediction of the global heat exchange coefficient of the evaporator.

We have found that the structures of mathematical models in the literature, the results for the properties of the processes, and the formulas have an extraordinary degree of dispersion, the limits of their usability are usually not disclosed, and their comparison is difficult. The developed mathematical models describe the operating processes that take place in the system elements of the heat pump only partially or with certain cutbacks.

Within the scope of this work, the goal is to formulate a complete mathematical model of the shell-and-tube evaporator that describes the steady-state behavior of the system with sufficient accuracy.

In the vast majority of models, researchers assume that the heat transfer coefficients and tube friction factors are constant. They do not take into account that the change in the process is a function of temperature, quantity, and vapor quality, or they use inaccurate, outdated equations created in previous decades for other refrigerants. The mathematical model presented also takes into account the pressure drop of the single-phase and two-phase refrigerant and the change in the heat transfer coefficient. The equation of state and thermodynamic properties of other refrigerants can be integrated into the presented model. However, the model is limited to shell-and-tube heat exchangers and does not support plate heat exchangers, for example.

The refrigerant R134a is selected because of its very favorable thermodynamic properties. In view of future restrictions (e.g., EU F-Gas Regulation), it can be relatively easily replaced by non-flammable (A1) HFO/HFC alternatives.

2. Materials and Methods

The aim of this study is to define a mathematical model of the evaporator suitable for operation with the refrigerant R134a and chilled water in a steady state. The proposed mathematical models of the heat exchanger are described by coupled differential equations with boundary conditions. The mathematical model of the evaporator cannot be solved analytically, so numerical simulation must be applied.

2.1. Mathematical Model

The mathematical model of the evaporator (Figure 1) is a coupled system of partial differential equations under certain boundary conditions. For solving the mathematical model the numerical recursive procedure Runge-Kutta and iterative Mac Adams were applied. The Runge-Kutta procedure is used only for getting initial solutions at discrete points, while the Mac Adams procedure converges the approximate initial solutions in an iterative manner to the accurate values. The following neglections were applied in the mathematical model with distributed parameters in a steady state:

-

the refrigerant flow in the evaporator is one-dimensional and steady. Only the radial heat transfer is taken into account, while the axial heat transfer by conduction through the tube walls and baffle is ignored,

-

in the two-phase flow of a refrigerant, the liquid, and the vapor are in thermodynamic equilibrium: their pressures and temperatures are the same,

-

in heat exchanger, only axial flow is considered,

-

the diameter of the tubes is the constant with length,

-

the effect of gravity on heat transfer is neglected,

-

evaporator’s heat gain from the environment is neglected.

The governing laws of fluid motion is derived using a control volume approach. In this case, the control volume is fixed for simplicity. When analyzing a control volume problem, there are three laws that are always valid:

• Conservation of Mass
• Conservation of Momentum
• Conservation of Energy

In addition to the above laws, the equations for heat transfer and heat dissipation are applied.

The mass conservation equation for a one-dimensional compressible fluid flow:

$$\frac{\partial \rho }{\partial t}+\frac{\partial \left(\rho ·w\right)}{\partial z}=0.$$

(1)

The steady state operation mode:

$$\frac{\partial \left(\rho ·w\right)}{\partial z}=0.$$

(2)

$$\rho ·w=G \to w=G·v.$$

(3)

$$\frac{\partial w}{\partial z}=G·\frac{\partial v}{\partial z},$$

(4)

$$\frac{\partial v}{\partial z}=\frac{\partial v}{\partial p}·\frac{\partial p}{\partial z}+\frac{\partial v}{\partial h}·\frac{\partial h}{\partial z},$$

(5)

$$\frac{\partial v}{\partial z}=v_p·\frac{\partial p}{\partial z}+v_h·\frac{\partial h}{\partial z},$$

(6)

By substituting Equations (4) and (6):

$$\frac{\partial w}{\partial z}=G·\left(v_p·\frac{\partial p}{\partial z}+v_h·\frac{\partial h}{\partial z}\right).$$

(7)

The momentum conservation equation:

$$\frac{\partial \left(\rho ·w\right)}{\partial t}+\frac{\partial \left(\rho ·w·w+p\right)}{\partial z}+f_x=0.$$

(8)

where f_x is the effect of shear stress on the fluid and is actually the pressure drop of the refrigerant and can be denoted as ∆p/dz. It is calculated according to Equations (52)–(57).

The steady state operation mode:

$$\frac{\partial \left(\rho ·w^2\right)}{\partial z}+\frac{\partial p}{\partial z}+f_x=0,$$

(9)

$$\frac{\partial \left(G·w\right)}{\partial z}+\frac{\partial p}{\partial z}+f_x=0,$$

(10)

$$G·\frac{\partial w}{\partial z}+\frac{\partial p}{\partial z}+f_x=0.$$

(11)

By combining Equations (7) and (11):

$$G·\left(G·\left(v_p·\frac{\partial p}{\partial z}+v_h·\frac{\partial h}{\partial z}\right)\right)+\frac{\partial p}{\partial z}+f_x=0,$$

(12)

$$G^2·\left(v_p·\frac{\partial p}{\partial z}+v_h·\frac{\partial h}{\partial z}\right)+\frac{\partial p}{\partial z}+f_x=0,$$

(13)

$$G^2·v_p·\frac{\partial p}{\partial z}+G^2·v_h·\frac{\partial h}{\partial z}+\frac{\partial p}{\partial z}+f_x=0,$$

(14)

$$\left(G^2·v_p+1\right)·\frac{\partial p}{\partial z}+G^2·v_h·\frac{\partial h}{\partial z}+f_x=0.$$

(15)

The energy conservation equations of the refrigerant:

$$\frac{\partial \left(\rho ·h_0\right)}{\partial t}+\frac{\partial \left(\rho ·w·h_0\right)}{\partial z}-\frac{\partial p}{\partial t}-\dot{q_i}\frac{O_i}{A_i}=0.$$

(16)

Stagnation enthalpy is defined by Equation (17).

$$h_o=h+\frac{1}{2}{w}^2.$$

(17)

The steady state operation mode:

$$\frac{\partial \left(G·\left(\frac{w^2}{2}+h\right)\right)}{\partial z}-\frac{\partial {\dot{q}}_i}{\partial z}=0,$$

(18)

$$G·w·\frac{\partial w}{\partial z}+G·\frac{\partial h}{\partial z}=\frac{\partial {\dot{q}}_i}{\partial z},$$

(19)

$$G^2·v·\frac{\partial w}{\partial z}+G·\frac{\partial h}{\partial z}=\frac{\partial {\dot{q}}_i}{\partial z}.$$

(20)

By combining Equations (7) and (20):

$$G^3·v·v_p·\frac{\partial p}{\partial z}+\left(G^3·v·v_h+G\right)·\frac{\partial h}{\partial z}=\frac{\partial {\dot{q}}_i}{\partial z}.$$

(21)

The heat transfer conservation equation between the water and the refrigerant:

The wall has negligible thickness: d_o = d_i, so it is T_wall,I = T_wall,o = T_wall

$${\rho }_{wall}·c{p}_{wall}·A_{wall}·n·\frac{\partial T_{wall}}{\partial t}={\dot{q}}_o-{\dot{q}}_i .$$

(22)

$$A_{wall}=\left(d_o^2-d_i^2\right)·\pi /4$$

(23)

A_wall is the tube’s cross-sectional area.

$${\dot{q}}_o={\alpha }_w·d_o·\pi ·n·\left(T_w-T_{wall}\right)$$

(24)

$${\dot{q}}_i={\alpha }_r·d_i·\pi ·n·\left(T_{wall}-T_r\right)$$

(25)

Mean temperatures of water and refrigerant is calculated acc. Equations (26) and (27).

$$T_w=\frac{T_{w,in}+T_{w,out}}{2},$$

(26)

$$T_r=\frac{T_{r,in}+T_{r,out}}{2}$$

(27)

The steady state operation mode:

$${\dot{q}}_o={\dot{q}}_i$$

(28)

$${\alpha }_w·d_o·\pi ·n·\left(T_w-T_{wall}\right)={\alpha }_r·d_i·\pi ·n·\left(T_{wall}-T_r\right).$$

(29)

$${\alpha }_w·\left(T_w-T_{wall}\right)={\alpha }_r·\left(T_{wall}-T_r\right),$$

(30)

$$T_{wall}=\frac{{\alpha }_w·T_w+{\alpha }_r·T_r}{{\alpha }_w+{\alpha }_r}.$$

(31)

The heat transfer equation between the water and the tube wall:

$${\rho }_w·c{p}_w·A_w·\frac{\partial T_w}{\partial t}=-{\dot{V}}_w·{\rho }_w·c{p}_w·\frac{\partial T_w}{\partial z}-{\alpha }_w·d_o·\pi ·n·\left(T_w-T_{wall}\right)$$

(32)

The steady state operation mode:

$${\dot{V}}_w·{\rho }_w·c{p}_w·\frac{\partial T_w}{\partial z}=-{\alpha }_w·d_o·\pi ·n·\left(T_w-T_{wall}\right),$$

(33)

$$\frac{\partial T_w}{\partial z}=-\frac{{\alpha }_w·d_o·\pi ·n·\left(T_w-T_{wall}\right)}{{\dot{V}}_w·{\rho }_w·c{p}_w}$$

(34)

The mathematical model of the evaporator is obtained by organizing the system of equations consisting of Equations (15) and (21) and derivatives:

$$\frac{dp}{dz}=-\frac{\left(f_x·G^2·v+\frac{\partial {\dot{q}}_i}{\partial z}·G\right)·v_h+f_x}{G^2·v_p+G^2·v·v_h+1}.$$

(35)

$$\frac{dh}{dz}=\frac{\left(f_x·G^3·v+\frac{\partial {\dot{q}}_i}{\partial z}·G\right)·v_p+\frac{1}{G}\frac{\partial {\dot{q}}_i}{\partial z}}{G^2·v_p+G^2·v·v_h+1}.$$

(36)

The proposed equations, supplemented by auxiliary equations, together with the equations of state for refrigerant and water and the physical properties of the evaporator (length, diameter, number of evaporator tubes, etc.) form the mathematical model of the evaporator for the steady state. Equations (29) and (34)–(36) are suitable for Runge-Kutta and Adam-Moulton solutions. To solve the basic equations, auxiliary equations and initial conditions are also needed. The auxiliary equations serve as correlations for the determination of the heat transfer coefficients (

Table 1. Correlations for heat transfer coefficient.

Region	Tube Side Heat Transfer Correlation	Equation
Single phase flow: Dittus-Boelter [24]	${\alpha }_l=0.023·R{e}^{0.8}·P{r}^{0.3}·\frac{\lambda }{d_i}.$ The Reynolds number (Re) and Prandtl number (Pr) for the flow in a tubes can be expressed as: $Re=\frac{w\cdot d_{\mathrm{i}}}{{\nu }_{\mathrm{r}}}$$, Pr=\frac{cp\cdot {\mu }_{\mathrm{r}}}{{\lambda }_{\mathrm{r}}}$	(37) (38)
Two phase flow: Kandlikar [25]	${\alpha }_r=max\left({\alpha }_{nb, }{\alpha }_{cb, }\right).$	(39)
	${\alpha }_{nb}=\left(0.6683·C{o}^{-0.2}f\left(F{r}_{lo}\right)+1058·B{o}^{0.7}·F_f\right)·{\alpha }_l$	(40)
	${\alpha }_{cb}=\left(1.136·C{o}^{-0.9}f\left(F{r}_{lo}\right)+667.2·B{o}^{0.7}·F_f\right)·{\alpha }_l$	(41)
	$f\left(F{r}_{lo}\right)={\left(25·F{r}_{lo}\right)}^{0.3}$ for horizontal tubes with $F{r}_{lo}\le 0.04$$\mathrm{where} {\alpha }_l$ is calculated by Equation (38). $\mathrm{For} \mathrm{R}134\mathrm{a} F_f$ = 1.63.	(42)
	$Bo=\frac{\dot{q}}{G{i}_{lg}} Co={\left(\frac{1-x}{x}\right)}^{0.8}{\left(\frac{\dot{q}}{G{i}_{lg}}\right)}^{0.5} F{r}_{lo}=\frac{G^2}{{\rho }_l^2·g·d_i}$	(43)
Single phase flow: Bell-Delaware [26]	Shell side heat transfer correlation
	$\alpha ={\alpha }_{id}·J_c·J_L·J_B·J_S·J_{\mu }$	(44)
	$\mathrm{Idealized} \mathrm{heat} \mathrm{transfer} \mathrm{coefficient} {\alpha }_{id}:$
	${\alpha }_{id}=j·c_p·m·P{r}^{\raisebox{1ex}{$-2$}\!\left/ \!\raisebox{-1ex}{$3 $}\right.}$Sm−1	(45)
	Baffle cut correction factor: $J_c=0.55+0.72·F_c$,	(46)
	Baffle Leakage correction factor: $J_L=0.44·\left(1-r_s\right)+\left. [1-0.44·\left(1-r_s\right)\right]exp\left(-2.2·r_m\right)$	(47)
	Bundle bypass correcrtion factor: $J_B=\exp \left(-C_{bh}\cdot F_{sbp}\cdot \left(1-\sqrt[3]{2\cdot r_{ss}}\right)\right)$ $\mathrm{The} \mathrm{empirical} \mathrm{factor} C_{bh}=1.35 \mathrm{for} \mathrm{laminar} \mathrm{flow} \mathrm{and} C_{bh}=1.25$ for transition and turbulent flows.	(48)
	Unequal baffle spacing correction factor: $J_S=\frac{\left(N_b-1\right)+{\left(L_{bi}/L_{bc}\right)}^{1-n}+{\left(L_{bo}/L_{bc}\right)}^{1-n}}{\left(N_b-1\right)+\left(L_{bi}/L_{bc}\right)+\left(L_{bo}/L_{bc}\right)}$	(49)
	Sm in Equation (45) is the cross-flow area.
	$S_m=L_b·\left. [D_s-D_{otl}+\frac{\left(D_{otl}-D_0\right)·\left(P_t-D_0\right)}{P_t}\right].$	(50)
	The viscosity correction factor: $J_{\mu }={\left(\frac{T+273}{T_{wall}+273}\right)}^{0.14}$	(51)

) and the pressure losses (

Table 2. Pressure drops correlations.

Region	Tube Side Pressure Drop Correlation	Equation
Single and two-phase flow: Grönnerud [27]	$\frac{\mathsf{\Delta }p}{dz}=\mathsf{\Phi }·{\left(\frac{\mathsf{\Delta }p}{dz}\right)}_l.$	(52)
	$\mathsf{\Phi }=1+{\left. [\frac{dp}{dz}\right]}_{Fr}·\left. [\frac{{\rho }_l/{\rho }_v}{{\left({\rho }_l/{\rho }_v\right)}^{0.25}}-1\right],$	(53)
	${\left. [\frac{dp}{dz}\right]}_{Fr}=f_{Fr}·\left. [x+4·\left(x^{1.8}-x^{10}·f_{Fr}{}^{0.5}\right)\right],$	(54)
	$f_{Fr}=F{r}^{0.3}+0.0055·\left. [ln\frac{1}{F{r}_l}\right].$	(55)
	${\left(\frac{dp}{dz}\right)}_l=\frac{2·{\lambda }_{wall}·G^2}{d·{\rho }_l},$	(56)
	${\lambda }_{wall}=\frac{0.316}{R{e}^{0.25}}$	(57)
	Shell side pressure drop correlation
Single phase flow: Bell-Delaware [25]	$$\begin{gathered} \mathsf{\Delta }{p}_s=\left. [\left(N_B-1\right)·\mathsf{\Delta }{p}_c·R_B+N_B·\mathsf{\Delta }{p}_w\right]·R_L+2·\mathsf{\Delta }{p}_c·R_B \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}·\left(1+\frac{N_{cw}}{N_C}\right). \end{gathered}$$ The ideal cross-flow pressure drop through one baffle space is obtained with the use of Equation (60)	(58)
	$\mathsf{\Delta }{p}_c=\left(K_a+N_c·K_f\right)·\left(\frac{{\rho }_l·{\dot{V}}^2}{2}\right),$	(59)
	The window zone pressure drop for Re ˃ 100 $\mathsf{\Delta }{p}_w=\frac{\left(2+0.6·N_{cw}\right)·{\dot{m}}_l^2}{2·S_m·S_w·{\rho }_l},$	(60)
	The number of effective cross-flow rows in the window zone: $N_{cw}=\frac{0.8·L_c}{p_p}, N_c=\frac{d_o·\left(1-2·L_c/D_s\right)}{p_p},$	(61)
	Equation (62) define the net cross-flow area through one baffle. $S_w=\frac{D_i{}^2}{4}·\left. [co{s}^{-1}{D}_b-D_b·{\left(1-D_b{}^2\right)}^{1/2}\right]-\frac{n}{8}· \left(1-F_t\right)·\pi ·D_s^2$ $D_b=\frac{D_s-2·L_c}{D_s}$	(62) (63)

). The equations of state of the thermodynamic properties of R134a are taken according to Huber and Ely [23].

The following table summarizes the correlations for calculating the pressure drop of a single-phase and two-phase refrigerant and the pressure drop of water in the evaporator shell side.

2.2. Numerical Analysis

The Runge-Kutta and Adams-Moulton predictor-corrector methods were used for the numerical solution of the differential equations. The developed simulation algorithm is coded in C++ and uses the REFPROP 10: Thermodynamic and Transport Properties Database developed by the National Institute of Standards and Technology (NIST). This database describes the properties of refrigerants such as density, enthalpy, liquid and vapor-specific heat, entropy, saturation temperature and pressure, thermal conductivity, dynamic viscosity, and surface tension. The numerical simulation is able to determine the thermodynamic parameters and the state parameters of the refrigerant at any point of the evaporator and, of course, at the endpoint of the evaporator with sufficient accuracy. The parameters such as temperature, pressure, vapor quality, and tube wall temperature of the refrigerant are studied as a function of the chilled water temperature. The heat transfer process is also analyzed. Both sections of the thermodynamic process in the horizontal tubes of the evaporator, evaporation, and superheating, are considered. The developed mathematical model is tested by a series of numerical simulations to obtain the values for the refrigerant and the water. The results of the numerical simulation were validated by experimental tests.

The input data of the mathematical model describe the steady-state operation of the evaporator:

• Refrigerant R134a
• Mass flow rate of water: m_w = 0.28 kg/s.
• Water inlet temperature: 13–19 °C
• Inner diameter of the shell: D_s = 28 mm.
• Number and diameter of tubes in the bundle: 5 × Ø8 × 1 mm.
• Tubes material: cooper
• Hydraulically smooth tubes
• Number of baffles: N_b = 15.
• The window of the baffles: 25%.

The heat exchanger shell is insulated from the outside, so that the heat exchange with the environment is neglected.

3. Experimental Investigation

In the experiment, the shell and tube evaporator is used for heat transfer between the water and the refrigerant. In the present case, the refrigerant R134a flows in the tubes while the cooling medium, i.e., the water, flows on the shell side. The baffles support the tube bundle and regulate the flow in the shell.

The measurement setup is shown in Figure 2. Measurement points were set up at ten discrete points on the evaporator to measure the temperature of the refrigerant and water. The measuring devices and their accuracy are listed in

Table 3. The measuring instruments’ characteristics.

Measuring Instruments	Designation	Accuracy
Thermometer sensor	DS18S20.	$\pm 0.3 \mathrm{K}$
Pressure gauge and pressure sensor	Dial indicator, Transducers.	1%, 0.5%.
Volume flow meter	B10 water flow sensor	0.6% and 0.2%
Coriolis effect mass flow meter	Krohne Optimas 6400	0.1%.

.

The initial conditions are presented in

Table 4. The initial conditions.

Measuring Variable	Designation	Value
The temperature of the cooled water	Tw	10–18 °C
The mass flow of cooled water	${\dot{m}}_w$	0.28 (kg/s)
Type of refrigerant	R134a	-
The mass flow of the refrigerant	${\dot{m}}_r$	0.023 (kg/s)

.

3.1. Experimental Procedure

The system was turned on and operated until steady-state heat exchange was achieved. The steady-state procedure shall include a period, usually twenty minutes, with a sampling time of one minute, during which the deviation of the vapor pressure must not exceed ±2.5 kPa. The deviation of the measured temperature must not exceed ±0.5 K, and the deviation of the flow rate of water must not exceed ±0.0005 kg/s.

Then, the inlet and outlet parameters of the evaporator were also read on the gauges, such as the temperatures and pressures of the refrigerant and the water, as well as the mass flow rate of the refrigerant and the volume flow rate of the water. At the beginning of the measurements, the inlet temperature of the water in the evaporator was 10 °C and the mass flow rate was 0.28 kg/s, a constant value during the series of measurements. Then the temperature of the cooled water was changed to 11 °C. The system was released again to establish a steady state so that the measurement results could be read again. The procedure was repeated for all other cooled water temperatures up to 18 °C.

3.2. Measurement Uncertainty

The experiment was performed under repeatable conditions, but the results were affected by small factors that were difficult to take into account, so the measured values were subject to some errors. Therefore, the value of the standard deviation σ was determined, which indicates how much the measured values y_k deviate from the arithmetic mean $\overline{y}$ of the data set.

$$\sigma =\sqrt{\frac{{{\displaystyle \sum }}_1^{60}{\left(y_k-\overline{y}\right)}^2}{60}}$$

(64)

The calculated uncertainty for cooling capacity is σ = 0.94%, for the temperature of refrigerant σ = 0.87%, for the temperature of water wall σ = 0.84%, and for wall temperature σ = 0.83%.

4. Result and Discussion

In this study, the change in the parameters of the refrigerant R134a flowing in the tubes of a counterflow evaporator is investigated as a function of the water temperature flowing on the shell side. The mass flow rate of the cooled water and the mass flow rate of the refrigerant are constant and are 0.28 kg/s and 0.023 kg/s, respectively. The simulation results are shown graphically in Figure 3, Figure 4, Figure 5, Figure 6, Figure 7, Figure 8, Figure 9 and Figure 10.

The decrease in water temperature was continuous along the evaporator (Figure 3). The change in water temperature ranged from 2.92 °C to 3.47 °C flowing into the shell space of the evaporator along its entire length.

Figure 4 shows the change in temperature of the tube’s outer wall along the length of the evaporator. At a water inlet temperature of 13 °C, the change in wall temperature along the length of the evaporator is 2.36 °C, while at a water temperature of 19 °C at the inlet of the evaporator, the change in tube wall temperature is 2.29 °C.

The change in the quality (content) of the refrigerant vapor is shown in Figure 5.

Vapor quality can be determined from Equation (65).

$$x=\frac{h-h_l}{h_v-h_l}$$

(65)

where:

h—specific enthalpy of mixture

h_l—specific enthalpy of saturated liquid

h_v—specific enthalpy of saturated vapor

Enthalpy h_l and h_v can be found in the thermodynamic properties tables. Enthalpy h can be determined on the basis of temperature and pressure. At the entrance of the evaporator, part of the refrigerant is already in the vapor phase. The evaporation process starts in the expansion valve due to the pressure drop.

The change in vapor content showed almost the same trend at all water inlet temperatures. At a length of the evaporator of z = 1.83 m, the vapor content takes the value x = 0.9, while in the last 25 cm of the evaporator length, overheating of the refrigerant (x = 1) already occurs.

The applied correlation of the heat transfer coefficient for water consists of two parts: The first part contains correction factors related to the geometry of the heat exchanger, which affects the turbulence of the flow. These values are constant along the tube, so the turbulence along the tube changes a little, and therefore the part of convective heat transfer changes a little as well. The second part of the correlation considers the physical properties of water as a working fluid. The values of thermal conductivity, heat capacity, viscosity, and density depend on temperature and pressure. As the temperature of the water changes along the pipe, these values also change and so do the dimensionless Reynolds number and Prandtl number. As the water temperature decreases, the density and dynamic viscosity values increase. As a result, the Reynolds number decreases, and the Prandtl number increases, decreasing the heat transfer coefficient. Since the influence of this part of the correlation is greater than the influence of the first part of the correlation, the general trend is a decrease in the value of the heat transfer coefficient along the tube (Figure 6).

Cooling capacity is the measure of a system’s ability to absorbing heat from water (brine). It is equivalent to the heat supplied to the refrigerant.

$${\dot{Q}}_w={\dot{m}}_w·C{p}_w·\left(T_{w,in}-T_{w,out}\right)$$

(66)

The cooling capacity increases as the temperature of the water at the inlet of the evaporator increases (Figure 7). The cooling capacity at a water inlet temperature of 13 °C was 3420 W, while the cooling capacity at a water inlet temperature of 19 °C was 4103 W. In general, the cooling capacity increases with the increase of the surface (length of evaporator) through which heat is exchanged.

At the entrance of the evaporator, the tubes are mostly filled with the liquid phase of the refrigerant. As a result of the heat flow through the tube wall, the internal energy of the liquid increases, which leads to an increase in turbulence and the heat transfer coefficient in the inlet area. It is characteristic of the Kandlikar model that the maximum value of the heat transfer coefficient occurs at a vapor quality of x = 0.5–0.6 for bubbly flow and nucleate boiling. The dominant region of nucleate boiling (bubbly region) is characterized by high heat transfer coefficients and low dependence on steam quality, since the wall temperature is nearly constant. The region where weak nucleate boiling predominates (slug region) is characterized by a decrease in heat transfer coefficients with increasing steam quality. In this region cyclic drying and rewetting of the tube wall occurs at the top and lateral ends of the tube. In the flow boiling region, convective boiling becomes more important, with a smaller contribution from nucleate boiling, where the liquid phase is usually present only in the lower half of the tube. As the wall temperature increases, the vapor density also increases, while the viscosity of the vapor decreases in contrast. These two facts together result in a significant increase in Reynolds number, which leads to increased forced convection heat transfer through the steam at the top wall of the tube. The bottom half of the tube continues to contain the stratified liquid in which nucleate boiling takes place. When the quality of the steam increases, the contribution of nucleate boiling decreases, and the heat transfer by forced convection of the steam is dominant. In this region, the heat transfer coefficient continues to decrease. This is followed by a drying region of high steam quality, characterized by a small decrease in the heat transfer coefficient. In the region of steam superheating, the heat transfer coefficient is low and almost constant, while the temperature of the tube wall increases faster. It can be seen that the refrigerant is superheated in the last 45 cm of the evaporator. In the superheated region, the values of the single-phase heat transfer coefficient are low and almost the same for all temperatures of inlet water.

The pressure drop of the refrigerant is continuous along the evaporator (Figure 9). The pressure drop of the refrigerant flowing in the evaporator tubes was Δp = 0.41 bar at a water inlet temperature of 13 °C, while at a water inlet temperature of 19 °C it was Δp = 0.297 bar.

Considering the equation of the state of an ideal gas and the p-h diagram of a heat pump, the explanation of Figure 10 is clear. When the temperature of the water at the inlet of the evaporator rises, the temperature of the refrigerant also rises. This increases the pressure in the system. When the pressure increases, the condensing temperature also increases. Consequently, the temperature and vapor content of the refrigerant is also higher after the expansion valve (inlet to the evaporator). The temperature drop of the evaporating refrigerant along the evaporator length is 3 °C at a water temperature of 13 °C, while at 19 °C the temperature drop is only 1.6 °C. The values for the temperature of the two-phase refrigerant show a decreasing tendency because the viscous friction is present that cause pressure drop. But the temperature of the refrigerant in the vapor phase in the superheating region increases spectacularly (2.76 m to 3 m evaporator length).

In the following Figure 11, Figure 12, Figure 13 and Figure 14, measurement results are compared with the values obtained from the mathematical model presented.

It is noted that the deviation between the experimental and predicted values of evaporating temperature ranged from −0.43 °C to +0.52 °C, while the experimental and predicted values of wall temperature were very similar, ranging from −0.48 to +0.51 °C.

Belman, et al. [7] compare predicted and measured values for the refrigerant R134a. The predicted temperature of the water-glycol mixture at the evaporator outlet has a maximum error of ±0.5 K. Relative error predictions for cooling capacity are approximately ±5%. Prediction of refrigerant mass flow rate ±5%. For evaporating pressure, the prediction error is within ±5%.

Navarro-Esbría et al. [16] tested the shell-and-tube evaporator experimentally and compared the results with correlations from the literature. They also chose R134a as the working fluid. The thermal analysis of the heat exchange was based on the ε-NTU method. The evaporator was modeled as a counter-flow heat exchanger. Among others, they also use the Kandlikar correlation, which achieves the highest accuracy and the lowest absolute mean error. The deviation between calculated and experimental values of cooling capacity was ±5%. The output parameter with the largest deviations between the predicted and experimental data is the evaporating pressure with 93.94% of the data points below 10%. The deviation for vapor pressure was 4.87%, for refrigerant outlet temperature ±1%, and for secondary liquid outlet temperature 0.03%.

Santos et al. [22] proposed a simplified physical model that assumes a constant heat exchange coefficient and pressure drop of the evaporator. This model is used to simulate and optimize the operation of the chiller. The rigorous model shows that large variations of the order of 10% are observed in the prediction of the global heat exchange coefficient of the evaporator.

These experimental test results are similar to the test results presented in this article but, according to the authors, show a slightly smaller deviation from the results obtained by numerical calculations.

5. Conclusions

The paper analyzes the thermodynamic properties of the counterflow tubular evaporator as a component of a heat pump with the working fluid R134a. The paper proposes a mathematical model that defines the physical processes in the evaporator. The mathematical model consists of basic conservation laws and auxiliary equations, as well as boundary and initial conditions. The basic laws are the laws of conservation of momentum, energy, and mass. The auxiliary equations are correlations of heat conduction and heat transfer, correlations of pressure drop, and correlations of thermodynamic and thermophysical properties of water and coolant. The mathematical model was verified by numerical simulations and experiments. The numerical simulation was performed using software originally developed in C++. For this purpose, the evaporator was divided into control volumes and the thermodynamic properties along the heat exchanger for the steady state of heat exchange were calculated.

Numerical simulation was used to study the variation of temperature, pressure, vapor quality of the refrigerant R134a, and temperature of the tube wall as a function of water temperature from 13 °C to 19 °C over the length of the evaporator, which was 3 m. The validation of simulation results was verified by experimental tests in a heat evaporator designed for this purpose and by the use of appropriate measuring equipment. The results of the experimental test confirm that the selected mathematical model is applicable to the selected type of heat evaporator. The deviation between the calculated and the measured value of the cooling capacity has a maximum error of 9.39%, while the average error is 3.05%. The deviations between the experimental data and the calculated values of refrigerant temperature along the evaporator range from −0.43 °C to +0.52 °C, while the deviation between the experimental data and the calculated values of water temperature range from −0.41 °C to +0.58 °C.