**Feasibility Study of a Centralised Electrically Driven Air Source Heat Pump Water Heater to Face Energy Poverty in Block Dwellings in Madrid (Spain)**

### **Roberto Barrella, Irene Priego, José Ignacio Linares \*, Eva Arenas, José Carlos Romero and Efraim Centeno**

ICAI School of Engineering, Comillas Pontifical University, 28015 Madrid, Spain; roberto.barrella@iit.comillas.edu (R.B.); 201807511@alu.comillas.edu (I.P.); earenas@icai.comillas.edu (E.A.); jose.romero@iit.comillas.edu (J.C.R.); efraim.centeno@iit.comillas.edu (E.C.) **\*** Correspondence: linares@comillas.edu

Received: 16 April 2020; Accepted: 25 May 2020; Published: 28 May 2020

**Abstract:** Energy poverty can be defined as the inability to pay the bills that are required for maintaining the comfort conditions (usually in winter) in dwellings. The use of energy efficient systems is one way forward to mitigate this problem, with one option being the electrically driven air source heat pump water heater. This paper assesses the performance of a centralised heat pump (200 kW of heating capacity) to meet the space heating demand of block dwellings in Madrid (tier four out of five in winter severity in Spain). Two models have been developed to obtain the following variables: the hourly thermal energy demand and the off-design heat pump performance. The proposed heat pump is driven by a motor with variable rotational speed to modulate the heating capacity in an efficient way. A back-up system is also considered to meet the peak demand. A levelised cost of heating of 92.22 €/MWh is obtained for a middle-level energy efficiency in housing (class E, close to D). Moreover, the following energy-environmental parameters have been achieved: more than 74% share of renewable energy in primary energy and 131.7 g CO2 avoided per kWh met. A reduction of 60% in the heating cost per dwelling is obtained if an energy retrofit is carried out, improving the energy performance class from E to C. These results prove that the proposed technology is among the most promising measures for addressing energy poverty in vulnerable households.

**Keywords:** energy poverty; centralised heat pump; hourly heating demand; off-design heat pump model

### **1. Introduction**

Although there is no agreed definition of energy poverty, there is some consensus in identifying it as the situation suffered by "individuals or households that are unable to adequately heat, cool or provide other necessary energy services in their homes at an affordable cost" [1]. Approximately 40 million people in the Union (8% of the total population) suffer from this situation, according to data from the latest European Living Conditions Survey [2,3]. In Spain, the incidence ranges from 7.2% to 16.9% of the population [4], depending on the indicator used, according to the latest update of the National Strategy against Energy Poverty.

In this context, there is a wide consensus in recognising that energy poverty is a key social challenge that must be addressed by Member States through their public policies. The private sector can also contribute to tackle this issue, by means of social entrepreneurship and the involvement of large corporations through their Corporate Social Responsibility [5].

Energy poverty, indeed a manifestation of general poverty, is directly linked to low income, and many low-income households are energy poor. Nevertheless, energy poverty does not fully

overlap with income poverty. Two energy-related variables, namely housing energy efficiency and energy prices, are also key features characterising the problem [6].

In terms of consequences, focusing only on health impacts, it is well documented that living in households with inadequate heating or cooling has detrimental consequences for respiratory, circulatory and cardiovascular systems, as well as for mental health and well-being [7]. Additionally, energy poverty impacts on other aspects of people life, namely economic, social, and educational ones, among others [8]. In the thick of energy needs that remain totally or partially unmet in vulnerable households, heating needs are particularly noteworthy. A potential improvement in this respect could be the installation of clean and efficient heating systems, such as the centralised heat pump analysed in this paper.

The air source heat pumps are a promising alternative to replace the old heating installations of block dwellings suffering energy poverty. These dwellings usually have oversized radiators, which enables them to be fed with low temperature water, then having a temperature range that is suitable for air source heat pumps water heaters (ASHPWH) [9]. European Union considers that the energy taken from the air in the evaporator can be accounted as renewable energy if the heating seasonal performance factor exceeds certain values [10]. Due to this fact, the replacement of centralised gas boilers by heat pumps makes it possible to introduce renewable energies in heating, so achieving an efficient heating system more de-carbonised than the former boiler. Such replacement should take the impact of the electricity taken from the grid in terms of CO2 emissions, which depends on the electricity generation mix of the country, into account. Furthermore, heat pumps must face two specific environmental issues, namely, the ozone deployment potential (solved some time ago with the hydrofluorocarbons, HFC) and the global warming potential (GWP). In fact, the GWP of the HFCs used as refrigerants is usually too high. European Union has regulated the use of fluorinated greenhouse gases (F-gases) [11], establishing an agenda to move from HFC to natural refrigerants (propane, butane, ammonia, CO2, etc.), passing through hydrofluoroolefins (HFO) as intermediate fluids. Spanish regulation [12] recently included heat pumps as a renewable option to supply thermal services to buildings, especially domestic hot water. Despite this, the use of heat pumps for heating purposes is far away from common practice in Spain. Conversely, space heating from heat pumps is usually a by-product of cooling services, using the so-called reversible heat pumps. This situation is even considered by the European Union [10], which penalises reversible heat pumps when counting renewable energy due to the fact that they are usually designed to operate in cooling mode.

Curve fitting from actual machines is a common methodology for modeling the behaviour of heat pumps. Accordingly, Underwood et al. [13] develop a model that is based on refrigerant-side variables, which makes it suitable for the analysis of the performance of heat pumps in service. So, this type of model is used when the focus is on the representation of the entire system building-heat pump. For instance, Lohani et al. [14] integrate correlation curves of coefficient of performance (COP) for ground and air source heat pumps in a building modelling software that lacks heat pump models. In some cases [15], the curve fitting is based on the theoretical behaviour of Carnot efficiency, as compared with actual performance. These models are used to identify key performance parameters in a site, as carried out by Vieira et al. [16].

Physical models of heat pumps are based on energy and heat transfer equations of the different components. The key components are the heat exchangers, which are modelled using the effectiveness-number of transfer units method or the equivalent logarithmic mean temperature difference method [17]. These types of methods usually consider the single phase and phase change zones. For example, Fardoun et al. [18] propose a quasi-steady state model that is based on an iterative procedure for the heat exchangers. In this sense, Patnode [19] develops a model of heat exchangers that is based on the Dittus–Boelter correlation, which can obtain the overall heat transfer coefficient as a function of the mass flow rate. This type of models is the base of rules and standards that determine the COP of ASHPWH as a function of water temperature and air dry and wet temperature [20].

Other studies use dynamic models, which are usually based on commercial software as TRNSYS. This type of analyses might pursue the improvement of the control strategies in the operation of heat pumps [21] or enhance the coupling with the thermal inertia of the building [22]. In short, dynamic models are useful when the seasonal performance is sought [23].

The thermal load, which is required as an input data to simulate the behavior of the heat pump, can be calculated using hourly average data or load predictions. Xu et al. [24] carry out a case study while using data that were recorded along two years of operation of a data center, evaluating the performance of the combined cooling, heat, and power (CCHP) system, also performing transient modeling. Gadd et al. performed an analysis of the heat meter readings that were obtained on an hourly basis along one year [25], aiming to provide heat load patterns, whereas Bacher et al. [26] attempted to develop a method for predicting the space heating thermal load in a dwelling. The latter model is based on measured data from actual houses in combination with local climate measurements and weather forecasts. Noussan et al. [27] propose a thermal demand model built up while using heat measurements taken every six minutes along several years of operation. Energy Plus software from the U.S. Department of Energy [28] was used by many researchers to determine the thermal load. In this sense, Wood et al. [29] and Michopoulos et al. [30] employed this software to analyse the use of biomass in space heating. Other authors used regressive models to obtain thermal energy demand prediction methods [31] or autoregression analysis with exogenous time and temperature indexes [32]. In this sense, Powell et al. [33] selected nonlinear autoregressive models with exogenous inputs as the best methodology based on artificial neural networks.

Several scholars used the degree-day method to obtain the hourly thermal energy demand profile. For example, Büyükalaca et al. [34] estimate the energy needs in a building in Turkey by calculating the heating and cooling degree-days using variable-base temperatures. Furthermore, Martinaitis et al. [35] perform an exergy analysis of buildings based on the degree-days method. Moreover, Layberry [36] analysed the errors in the degree-day method that may affect the building energy demand analysis. Carlos et al. [37] combined solar radiation data with the degree-day method and compared several different simplified methodologies for building energy performance assessment in winter. In Spain, the winter climatic severity index is defined from the winter degree-days and solar radiation measurements [38]. This index is used to determine the thermal energy demand and it makes it possible to characterise the country's climatic zones in order to assess the energy requirements, according to the building energy performance regulations [39]. This thermal energy demand modelling is also suitable for assessing the performance of other thermal devices for the evaluated scenario.

This paper analyses the efficiency of an electrically driven a heat pump as a realistic alternative to achieve winter thermal comfort in vulnerable households' dwellings. Thermally driven heat pumps (driven by absorption or internal combustion engines) also constitute a feasible option; nevertheless, this research is focused on electrically driven heat pumps, due to its higher commercial deployment in Spain. In order to do so, the baseline case is defined for a block of dwellings with a middle-level of energy efficiency and, additionally, retrofitting alternatives are considered for under average situations. The performance of the heat pump in terms of cost and CO2 emissions is compared with other alternatives (centralized and decentralized boilers). The cost breakdown of the heat pump is detailed, allowing for the evaluation of subsidy schemes in order to fund this kind of devices.

The novelty of this paper lies on the assessment of the use of centralised electrically driven heat pumps to meet the heating demand in Madrid. Such solutions are usually employed in northern Europe areas, especially with ground source heat pumps. However, the use of air source heat pumps is not common in Spain, as it can be followed from the support to these devices in the last release of the Spanish Technical Building Code [12]. In this sense, the developed heat pump model is not sophisticated, although it is accurate enough, as can be derived from the comparison with a commercial machine. Regarding the heating demand forecasting model, a simpler version has been proposed by the authors [40], but, as a novelty, the version that is used in this manuscript is able to retain information regarding thermal insulation of the building. This information made it possible to assess the effect of energy building retrofitting. Therefore, the aim of this paper is to analyse the feasibility of ASHPWH as active measure to fight energy poverty in dwelling blocks. The environmental and economic assessment performed in this paper can eventually provide insightful information to policy makers for implementing clean and efficient measures to tackle energy poverty.

### **2. Methodology**

### *2.1. Hourly Heating Demand*

Building Technical Code (BTC) in Spain establishes a procedure for assessing the energy demand in both winter and summer as a function of the climate severity index (CSI) [41]. In this work, only heating demand is assessed, so focusing on winter energy needs, because energy poverty studies have been mainly focused on this season [42].

The reference specific demand (*RD*, kWh/m2) in winter is given by Equation (1), where *WSI* stands for the winter severity index and Table 1 provides the coefficients α and β. The winter severity index is defined in Equation (2), where *RAD* is the average accumulated global radiation over horizontal surface during January, February, and December (Equation (3)), and *DD* is the average degree-days (at base temperature Tb = 20 ◦C) for the same months (Equation (4a)) [43]. Table 2 provides the coefficients for Equation (2). The calculation of *RAD* requires the global hourly radiation over horizontal surface (*rk*), whereas the calculation of *DD* requires the hourly temperature difference (Δ*Tk*), as defined by Equation (4b). Hourly values of *Tk* are available for each climatic zone in the web site of the BTC [44].

$$RD = \alpha + \beta \cdot \text{WSI} \tag{1}$$

$$\text{WSI} = a \cdot \text{RAD} + b \cdot \text{DD} + c \cdot \text{RAD} \cdot \text{DD} + d \cdot \text{RAD}^2 + e \cdot \text{DD}^2 + f \tag{2}$$

$$RAD = \frac{\sum\_{k=1}^{24 \times 90} r\_k}{3} \tag{3}$$

$$DD = \frac{\sum\_{k=1}^{24 \times 90} \Delta T\_k}{24 \times 3} \tag{4a}$$

$$
\Delta T\_k = \begin{cases} T\_b - T\_k & \text{if } T\_b > T\_k \\ 0 & \text{otherwise} \end{cases} \tag{4b}
$$

**Table 1.** Coefficients required to obtain the reference specific demand in winter [41].


**Table 2.** Coefficients required to obtain the winter severity index (*WSI*) [43].


*RD* and *WSI* are overall values, that is, they are calculated for the complete winter season, as it is observed in Equations (1) and (2). Based on these equations, a Taylor series expansion of first order around (*RAD*, *DD*) has been carried out over *WSI*. This procedure leads to an hourly expression of the specific reference demand (Equation (5a)), where: (1) the index "*j*" is extended from 1 to 4368 (number of hours from January to March and from October to December), (2) *Nd* is equal to 182 days (number of days in the same period), (3) *Nm* is equal to 6 (number of months), and (4) the star denotes that this specific reference demand needs to be corrected. This correction has to be done to take into account two different issues. Firstly, the effect of the radiation and second, the fact that three additional months have been included with respect to the original correlation.

The correction of the radiation is performed because, sometimes, its value is high enough to result in a cooling demand (negative heating demand). In this case, the heating demand is set to zero. On the other hand, to consider the inclusion of additional months in the formulation, a reduction coefficient (*Cr*) is defined as the ratio of the actual specific seasonal demand (*RDa*, given in regulations [45] and equal to 53 kWh/m2 for Madrid) to the summation of *RD*<sup>∗</sup> *<sup>j</sup>* over the 4368 h. Accordingly, the corrected hourly specific reference demand *RDj* is given in Equation (5d).

$$RD\_{\dot{j}}^{\*} = \frac{\alpha + \beta \cdot (\mathcal{WSI} - \rho \cdot \text{RAD} - \delta \cdot DD)}{24 \cdot N\_d} + \left(\frac{\beta \cdot \rho}{N\_m}\right) \cdot r\_{\dot{j}} + \left(\frac{\beta \cdot \delta}{24 \cdot N\_m}\right) \cdot \Delta T\_{\dot{j}} \tag{5a}$$

$$
\rho = a + 2 \cdot d \cdot RAD + c \cdot DD \tag{5b}
$$

$$\delta = b + \mathbf{c} \cdot \mathbf{R}AD + \mathbf{2} \cdot \mathbf{c} \cdot DD \tag{5c}$$

$$RD\_{\dot{j}} = RD\_{\dot{j}}^{\*} \underbrace{\begin{pmatrix} RD^{a} \\ \sum\_{j=1}^{4368} RD\_{\dot{j}}^{\*} \end{pmatrix}}\_{\mathcal{C}\_{\mathcal{I}}} \begin{cases} 1 & \text{if } RD\_{\dot{j}}^{\*} > 0 \\ 0 & \text{otherwise} \end{cases} \tag{5d}$$

Once the hourly specific reference demand is obtained, it should be corrected according to the energy performance index (*EPI*) and the ratio of the reference demand of the whole stock of reference buildings to the 10-th percentile of this stock (*R*) [41]. This correction leads to the calculation of the hourly absolute demand (*Dj*, Equation (6)), where the heated area (*A*) has been included. In Equation (6), the *EPI* is obtained from the energy performance certificate of the building and *R* is given at Table 3, where the climatic zone ranges from mild winter (A) to severe winter (E). For the current research, the *EPI* values have been calculated by cross-correlating the CENSUS 2011 data [46] and the buildings-energy-certification data [47]. The average *EPI* values for D zone (Madrid) are found to be 3.53 for block dwellings built before 1981, 2.18 if they were built between 1981 and 2007 and 0.92 for the ones built after 2008. Vulnerable households typically live in old buildings and use inefficient heating installations, typically electric radiators, as shown in the literature [48–50]. For this reason, the first two building age categories are those in which this collective of households is commonly located. In this study, as explained in Section 3.2, the block dwellings built between 1981 and 2007 are chosen as the baseline case.

$$D\_{\dot{j}} = A \cdot RD\_{\dot{j}} \cdot \left(\frac{1 + (EPI - 0.6) \cdot 2 \cdot (R - 1)}{R}\right) \tag{6}$$



Figure 1 displays the hourly demand for each ambient wet bulb temperature (Wet bulb temperature is used as a measurement of the enthalpy of humid air). It shows a cloud of points following a linear regression curve (enclosed in a red dashed line), along with a set of disperse data with lower heating demand. The data density increases with the temperature, in agreement with the radiation issue previously explained.

**Figure 1.** Hourly heating demand profile for each ambient wet bulb temperature for a set of block dwellings in Madrid built between 1981 to 2007 with 6000 m2 of total heated surface.

Finally, the hourly demand is sorted from maximum to minimum, obtaining the annual cumulative heating demand profile, as shown in Figure 2.

**Figure 2.** Annual cumulative heating demand profile for a set of block dwellings in Madrid built between 1981 to 2007 with 6000 m2 of total heated surface.

#### *2.2. Heat Pump Model*

Two models have been developed for obtaining the performance of the heat pump: one for the best efficiency point (BEP) and another for the off-design operation. The former is used to size the main components and the latter to obtain the performance map. The heat pump uses the air as thermal source. Thus, a rotational speed control (inverter) driving the compressor motor is assumed to avoid the loss of heating capacity when the ambient temperature decreases. The nominal rotational speed is taken as 1490 rpm, with a range of variation from 745 to 2235 rpm (±50%). Out of these limits, a back-up system is necessary, assumed as a condensing boiler with modulation, fuelled by natural gas. The efficiency of this boiler has been taken as 95% (based on higher heating value, HHV), assumed constant when considering the usual seasonal performance factor values that were reported by manufacturers [51]. The rotational speed of the evaporator fan is also controlled to keep constant the temperature drop in the air to further improve the heat pump efficiency. The heat pump is an air/water system, so water of the existing radiators heating loop is heated in the condenser. The heat pump takes advantage of the usual oversizing employed in the radiators heat transfer area to make them behave as low temperature radiators, thus enabling the use of heat pump as heater. Domestic hot water demand is covered by solar thermal energy that is supported by natural gas, with this aspect being out of the scope of this analysis. Cooling demand is not considered, according to the common trend in energy poverty studies [42].

Figure 3 shows a scheme of the heat pump. An adiabatic reciprocating compressor is chosen, in accordance with actual commercial trends [52]. The isentropic efficiency (Equation (7)) is used to model the compressor and the polytropic exponent (*n*) can be obtained, as shown in Equation (8). The volumetric efficiency (η*v*, Equation (9)) is used to set the refrigerant mass flow rate . *m* . In Equations (7)–(9), *h* stands for enthalpy, *s* entropy, *p* pressure, *v* specific volume, vs. swept volume, *r* pressure ratio, and α relative clearance volume.

**Figure 3.** Layout of the heat pump.

The nominal parameters to solve the best efficiency point are:

	- Heating capacity . *Qc* : 200 kWth
	- Water conditions: 45 ◦C at condenser inlet (*Twi*) and 55 ◦C at condenser outlet (*Two*). A counterflow heat exchanger is assumed, feeding the water to an existing low temperature radiator system.
	- Air conditions: 5 ◦C as ambient wet bulb temperature (taken as evaporator inlet, *Tai*)) and −5 ◦C as wet bulb temperature at evaporator outlet (*Tao*). An ambient pressure of 95 kPa is assumed (Madrid exhibits an elevation over the sea level of 655 m).
	- Outlet working conditions: superheating (Δ*Tv*) of 5 ◦C
	- Inlet working conditions: approach temperature (Δ*Te*) of 10 ◦C
	- Speed (*N*): 1490 rpm
	- Isentropic efficiency (η*s*): 75% (includes motor and inverter efficiencies)
	- Relative clearance volume (α): 3%.

$$\eta\_s = \frac{h(s\_1; p\_2) - h\_1}{h\_2 - h\_1} \tag{7}$$

$$r = \left(\frac{p\_2}{p\_1}\right) = \left(\frac{v\_1}{v\_2}\right)^n \tag{8}$$

$$\eta\_{\mathcal{V}} = \frac{\dot{m}}{\left(\frac{V\_r}{v\_1}\right) \cdot \left(\frac{N}{60}\right)} = 1 - \alpha \cdot \left(r^{1/n} - 1\right) \tag{9}$$

The condenser is a refrigerant/water counter-flow heat exchanger. Equation (10) shows the energy balance, where . *mw* stands for mass flow rate of water. The approach temperature results in a relationship between both fluids (Equation (11)), where refrigerant outlet temperature is the condensation temperature (refrigerant is saturated liquid at state 3), linked to the condensation pressure (Equation (12)). Enthalpy for water is assessed as enthalpy of saturated liquid at the water temperature. No pressure drop is considered (Equation (13)).

$$
\dot{Q}\_{\mathbb{C}} = \dot{m} \cdot (h\_2 - h\_3) = \dot{m}\_{w^\cdot} (h\_{wo} - h\_{wi}) \tag{10}
$$

$$
\Delta T\_{\mathcal{E}} = T\_{\mathcal{B}} - T\_{\text{wi}} \tag{11}
$$

$$T\mathfrak{g} = T\_{\text{sat}}(p\mathfrak{g})\tag{12}$$

$$p\_2 = p\_3 \tag{13}$$

The valve is considered to be adiabatic, so, as kinetic and potential energy are neglected, it is modelled as iso-enthalpic (Equation (14)). Besides, the valve is assumed to be a thermostatic expansion valve, so it keeps the superheating at the compressor suction (Equation (15)) constant by acting on the refrigerant mass flow rate.

$$h\mathfrak{z} = h\_{\mathfrak{4}}\tag{14}$$

$$
\Delta T\_{\overline{v}} = T\_1 - T\_4 = \text{constant} \tag{15}
$$

The evaporator is an air/refrigerant cross-flow heat exchanger. Equation (16) gives the energy balance, where . *ma* stands for air mass flow rate and . *QE* stands for the rate of heat transfer. As in the condenser, the approach temperature establishes a relationship between the fluids (Equation (17)), where refrigerant inlet temperature is the evaporation temperature (refrigerant is a liquid-vapor

mixture at state 4), linked to the evaporation pressure (Equation (18)). No pressure drop is considered (Equation (19)). .

$$
\dot{Q}\_{\rm E} = \dot{m} \cdot (h\_1 - h\_4) = \dot{m}\_{\rm a} \cdot (h\_{\rm ai} - h\_{\rm ao}) \tag{16}
$$

$$
\Delta T\_{\varepsilon} = T\_{\text{ao}} - T\_4 \tag{17}
$$

$$T\_4 = T\_{sat}(p\_4) \tag{18}$$

$$p\_4 = p\_1 \tag{19}$$

The main performance parameters of the heat pump in its nominal point (BEP) are the heating capacity (useful heat released in the condenser), previously defined, the compressor consumption ( . *W*, Equation (20)), and the COP (Equation (21)). In the design point, all of them are instantaneous values; later on, in the off-design operation, they will be redefined as seasonal parameters, i.e., both the power rates ( . *QC* and . *W*) and the instantaneous COP will be time-integrated.

$$
\dot{\mathcal{W}} = \dot{m} \cdot (\hbar\_2 - h\_1) \tag{20}
$$

$$\text{COP} = \frac{\dot{Q}\_{\text{C}}}{\dot{W}} \tag{21}$$

The refrigerant must comply with regulations about the global warming potential (GWP) and ozone depletion potential (ODP) [11] of the European Union. Accordingly, R-290 (propane) is chosen as fluid valid for long term, due to the fact that it is a natural refrigerant with null ODP and a GWP value of 3. On the other hand, it is included in the A3 class of the flammability safety classification (ASHRAE), thus considered highly flammable. Figure 4 shows the pressure-enthalpy diagram for R-290, where the assumptions previously stated for the best efficiency point are represented. This diagram shows that the compressor discharge temperature is not too high and the de-superheating interval is small (around 20% of the overall heating capacity). Both of the values are in accordance with the fact that domestic hot water is heated by using other procedures.

**Figure 4.** P-h diagram for the design point.

The main parameters of the heat pump are fixed once the design point has been solved. Table 4 summarises these values. At this point (inlet air temperature of 5 ◦C and outlet water temperature of 55 ◦C, usually known as A5/W55), the heating capacity is set as 200 kWth, the compressor consumption is found to be 76.15 kWe, and then the COP results in 2.626.

**Table 4.** Main parameters of the heat pump calculated in the design point and set as constant in off-design operation.


The modelling of the BEP and the off-design operation have been carried out for different purposes. On one hand, the model of the BEP is developed to define the heat pump parameters that determine its size (listed in Table 4). On the other hand, the off-design model aims to obtain the performance of the heat pump as a function of both the ambient temperature and the heating load, in order to work out the seasonal performance of the device. The input variables are the ambient wet bulb temperature and the rotational speed of the compressor. Equations (8)–(21) are now solved using the parameters that are given in Table 4. It should be noted that the isentropic efficiency of the compressor is replaced by the polytropic relationship (Equation (8)) and the mass flow rates of refrigerant, air, and water are now unknown. Approach temperatures in the heat exchangers have been fixed, while assuming that the number of transfer units are high enough to work in the asymptotic range of effectiveness. The main functions of the off-design model are listed in Equations (22) and (23), which lead to Equation (24). Another limit should be imposed: the maximum driving power of the motor, being assumed as 1.5 times the power consumed at BEP. .

$$
\dot{Q}\_{\mathbb{C}} = \dot{Q}\_{\mathbb{C}}(T\_{ai}, N) \tag{22}
$$

$$\text{COP} = \text{COP}(T\_{\text{ai}}, N) \tag{23}$$

$$
\dot{\mathcal{W}} = \frac{\dot{\mathcal{Q}}\_{\mathbb{C}}(T\_{\dot{a}i\prime} \text{ N})}{\text{COP}(T\_{\dot{a}i\prime} \text{ N})} = \dot{\mathcal{W}}(T\_{\dot{a}i\prime} \text{ N})\tag{24}
$$

### *2.3. Coupling of the Demand and Heat Pump Models*

Once the hourly heating demand (Equation (6)) and the heat pump performance (Equations (22) and (24)) are determined, the map of the device coupled to the demand can be obtained, producing a diagram similar to the one that is shown in Figure 5. In this chart, the heating demand only considers the linear regression curve for the sake of clarity. At the nominal rotational speed, a positive slope line determines the heating capacity of the heat pump for each temperature. This line cuts to the heating load line in one point. To modulate the response of the heat pump, the rotational speed is varied, therefore sweeping the operation zone and cutting to the load curve in a large range. For temperatures out of that range, the back-up system operates (if the load is higher than the heat pump capacity) or the machine is operated in on/off mode, in order to adapt the excess of capacity to the low demand.

**Figure 5.** Coupling between demand and performance of the heat pump.

Thus, the consumption of the heat pump ( . *Wj*) and the back-up system ( . *F bkp <sup>j</sup>* ) (if any) to meet the demand are calculated, for each operation hour, while using the functions that are given in Equation (25a–d). In Equation (25c), *bkp* stands for the back-up system, assumed as a natural gas condensing boiler with a constant efficiency (η*bkp*) value of 95% on higher heating value basis. In the same equation, . *F* stands for the natural gas consumption (again based on HHV) of the back-up system. Equation (25d) determines the power of the back-up boiler ( . *Qbkp*). In the on/off operation range, the demand is covered by the heat pump working at the minimum rotational speed, as it is derived from the algorithm described in Equation (25). Defrosting cycles are neglected, due to the usually low air moisture content in Madrid [53].

$$\dot{Q}\_{\text{C,j}}^{\text{max}} = \max \left\{ \dot{Q}\_{\text{C}} \left( T\_{\text{ai,j}}, N\_{\text{max}} \right), \text{COP} \left( T\_{\text{ai,j}}, N\_{\text{max}} \right) \cdot 1.5 \cdot \dot{W} \left( T\_{\text{ai,j}}^{\text{design}}, N\_{\text{design}} \right) \right\} \tag{25a}$$

$$
\dot{\mathcal{W}}\_{\dot{j}} = \begin{cases}
\dot{Q}\_{\mathcal{C},\dot{j}} / \text{COP}\{T\_{\dot{a}\dot{i},\dot{j}\prime} \ N\_{\dot{j}}\} & \text{if } D\_{\dot{j}} \le \dot{Q}\_{\mathcal{C},\dot{j}}^{\text{max}} \\
\dot{Q}\_{\mathcal{C},\dot{j}}^{\text{max}} / \text{COP}\{T\_{\dot{a}\dot{i},\dot{j}\prime} \ N\_{\dot{j}}\} & \text{otherwise}
\end{cases} \tag{25b}
$$

. *F bkp <sup>j</sup>* = ⎧ ⎪⎪⎪⎨ ⎪⎪⎪⎩ <sup>0</sup> *if Dj* <sup>≤</sup> . *Q max C*,*j Dj* <sup>−</sup> . *Q max C*,*j* /η*bkp otherwise* (25c)

$$
\dot{Q}\_{bkp} = \max \left\{ \dot{F}\_j^{bkp} \right\} \tag{25d}
$$

Some seasonal performance indexes have been defined: heating seasonal performance factor (*HSPF*, Equation (26)), CO2 avoided ratio (*AVCO2*, Equation (27)) and renewable input to heating demand ratio (*R2H*, Equation (28)). Numerical coefficients that are employed in Equations (27) and (28a) make it possible to consider the environmental impact of the electricity coming from the grid to drive the heat pump. They have been taken from [54], according to the Spanish energy sector. Equation (28b) comes from the current EU regulation regarding the support to heat pumps [10].

$$HSPF = \frac{\sum\_{j=1}^{4368} \dot{W}\_j \cdot COP\left(T\_{ai,j\text{\textquotedblleft N}j}\right)}{\sum\_{j=1}^{4368} \dot{W}\_j} \tag{26}$$

$$AVCO2 = \frac{\left(\frac{\sum\_{j=1}^{4368} D\_j}{\eta\_{k\bar{p}}}\right) \cdot 0.252 - \left[0.331 \cdot \sum\_{j=1}^{4368} \dot{W}\_j + 0.252 \cdot \sum\_{j=1}^{4368} \dot{F}\_j^{k\bar{p}}\right]}{\sum\_{j=1}^{4368} D\_j} \tag{27}$$

$$R2H = \frac{\sum\_{j=1}^{4368} RES\_j + 0.414 \cdot \sum\_{j=1}^{4368} \dot{W}\_j}{\sum\_{j=1}^{4368} D\_j} \tag{28a}$$

$$\begin{array}{c} \text{RES}\_{j} = \begin{cases} \dot{Q}\_{\text{C},j} \cdot \left(1 - \frac{1}{\text{HSPF}}\right) \text{if } HSPF \ge 2.5275\\ 0 & \text{otherwise} \end{cases} \tag{28b} \end{array} \tag{28b}$$

### *2.4. Economic Model*

The main indicator to assess the economic feasibility of this kind of device is the levelised cost of heating (LCOH), which integrates both the investment and operating costs [55]. In this research, two LCOH have been calculated: one referred to the whole heating demand (*LCOHDB*[/*MWh*], (Equation (29a)) and another referred to the whole heated area (*LCOHAB* /*m*<sup>2</sup> (Equation (29b)). In Equation (29a), *INV* stands for investment, *C* for annual cost, superscript *M* refers to maintenance, *W* to power consumption, *F* to fuel consumption, subscript *0* to costs in year zero, *CRF* to the capital recovery factor (Equation (29c)), *CELF* stands for the constant escalation levelisation factor (Equation (29d)). In Equation (29c–e), *rx* represents the nominal escalation rate of the item *x*, *wacc* the weighted average capital cost, and *Ny* is the life span of the project.

$$LCOH\_{DB} = \frac{\left(INV\_{HP} + INV\_{bkp}\right) \cdot \text{CRF} + \text{C}\_0^W \cdot \text{CELF}^W + \text{C}\_0^F \cdot \text{CELF}^F + \text{C}\_0^M \cdot \text{CELF}^M}{\sum\_{j=1}^{4368} D\_j} \tag{29a}$$

$$LCOOH\_{AB} = LCOOH\_{DB} \cdot \left[ \frac{\sum\_{j=1}^{4368} D\_j}{A} \right] \tag{29b}$$

$$CRF = \frac{wacc \cdot (1 + wacc)^{Ny}}{(1 + wacc)^{Ny} - 1} \tag{29c}$$

$$\text{CELF}\_{\text{x}} = \left[ \frac{k\_{\text{x}} \cdot \left( 1 - k\_{\text{x}}^{Ny} \right)}{1 - k\_{\text{x}}} \right] \cdot \text{CRF} \tag{29d}$$

$$k\_x = \frac{1 + r\_x}{1 + uacc} \tag{29e}$$

The investment for the heat pump (inverter model with 200 kWth of heating capacity) has been taken as 24,753 € [52], and a scale law has been fit for the investment of the back-up boiler (Equation (30)).

$$[INV\_{bkp}[]] = 1087.6 \cdot \dot{Q}\_{bkp}^{0.506} [k\mathcal{W}\_{th}] \tag{30}$$

For consumptions with installed power higher than 15 kWe, the cost of electricity includes a term for maximum consumed power in a year and a term for annual consumed energy. Moreover, this consumption considers hourly discrimination in three periods (P1, P2, and P3), according to Table 5.


**Table 5.** Electrical tariff [56].

<sup>1</sup> The cost associated to the power term is derived from multiplying the tariff by the maximum power consumed along the year at each period.

For comparison purposes, two additional scenarios have been considered: a decentralised system and a centralised one, both using only condensing natural gas boilers with the same efficiency as the back-up boiler (95% based on HHV). In the de-centralised boiler case, 60 single dwellings of 100 m<sup>2</sup> are considered, each employing a natural gas boiler with an investment of 800 € and maintenance costs of 150 €/year. In the centralised boiler scenario, a larger boiler is used to meet the overall demand, being given its investment by Equation (30) and a maintenance cost assumed of 2000 €/year.

All of the costs include taxes. Regarding the natural gas costs, two tariffs have been selected, depending on the annual consumption. Accordingly, for large consumptions (centralised cases in both heat pump with back-up and boiler alone) the tariff is 971.64 €/year plus 40.66 €/MWh (based on HHV), whereas for small consumptions (decentralised boilers) the tariff is 115.98 €/year plus 42.4 €/MWh (again based on HHV) [57].

Maintenance cost has been set to 2000 €/year for the proposed system (heat pump supported by a boiler).

Table 6 shows the economic parameters used to calculate the levelised costs.


**Table 6.** Assumed economic parameters.

### **3. Results**

### *3.1. Heat Pump Model*

Once the model is completed, Equations (22)–(24) are solved, obtaining, respectively, Equations (31) to (33). For a given rotational speed, the relation between heat capacity and air temperature is linear, and the slope increases with the speed, as it can be seen in Equation (31). Regarding the COP, Equation (32) is obtained, but it does not depend on rotational speed. This result is foreseeable, because the COP is the ratio between heat capacity and compressor consumption, both proportional to the mass flow rate, which is directly dependent on the rotational speed. Combining both Equations (31) and (32), the expression for the compressor consumption is obtained (Equation (33)).

$$\dot{Q}\_{\text{C}} = (0.1106 \cdot \text{N} + 0.0002) + (0.0053 \cdot \text{N} + 0.0003) \cdot T\_{\text{ai}} \tag{31}$$

$$\text{COP} = 2.49499 + 0.02981 \cdot T\_{\text{ai}} \tag{32}$$

$$\dot{W} = \frac{(0.1106 \cdot N + 0.0002) + (0.0053 \cdot N + 0.0003) \cdot T\_{ai}}{2.49499 + 0.02981 \cdot T\_{ai}} \tag{33}$$

Figure 6 shows the iso-lines of heating capacity as a function of the rotational speed and the ambient wet bulb temperature. The limitation of maximum electrical power (Equation (25a)) has been taken into account.

**Figure 6.** Iso-lines of heating capacity of the ASHPWH.

A comparison with an actual machine has been performed in order to validate the results of the heat pump model. A literature search has been carried out to find an existing device that fits with the characteristics of the modelled one. The selected equipment has been HERA 190-2-2 [58], an air to water heat pump from EUROKLIMAT. A reciprocating compressor using an inverter drives it and it has a design heating capacity of 190 kWth. The refrigerant is R290 as the one selected in this paper. The water outlet temperature of the existing heat pump is 55 ◦C, as in the model. Figure 7 plots the comparison of the COP, showing a good match, especially at low temperatures. This behavior is due to the fact that the HERA 190-2-2 data are based on dry bulb air temperature, whereas the data from the model are based on the wet bulb temperature, with both temperatures being very similar at low dry bulb temperature.

**Figure 7.** Comparison between the coefficient of performance (COP) obtained with the model and the COP of HERA 190-2-2 from EUROKLIMAT.

### *3.2. Baseline Case*

A block dwelling with an overall heated area of 6000 m<sup>2</sup> and EPI of 2.18, i.e., built between 1981 and 2007, located in Madrid, has been simulated. It is important to point out that the selected baseline-building-age represents a considerable percentage of vulnerable households in Spain and the 41% of the overall main houses with heating [46]. Figure 2 shows its annual cumulative heating demand profile. Applying the algorithm given by Equation (25) to each winter hour, the results that are shown in Table 7 are obtained. From these values, the performance indexes that are given by Equations (26)–(28) are calculated, and Table 8 summarises the results. Furthermore, Figure 8 shows the contribution to the heating demand of both the heat pump (96%) and the back-up boiler (4%).


**Table 7.** Energy results in the baseline case.



**Figure 8.** Contribution of heat pump (ASHPWH) and back-up boiler (BOILER) to meet the annual cumulative heating demand profile in the baseline case (block dwellings in Madrid built between 1981 to 2007 with 6000 m2 of total heated surface).

At low thermal demands, the heat pump is controlled in on/off mode due to the minimum rotational speed limit. This type of control causes important energy losses when an inverter is not available, but this driving control is expected to reduce them. Accordingly, Figure 9 shows the cloud of points for the thermal demand under the minimum heat capacity line. Once these points are sorted and moved into the ASHPWH contribution to the annual cumulative heating demand profile, they take up the end tale of the distribution (Figure 10), accounting for 8% of the overall demand that is met by the heat pump. Figure 9 also shows that the motor consumption ranges between 31.2 to 49.9 kWe (41 to 65.5% as compared to the consumption at design point). This allows for us to assume that low current peaks will take place in the startups and, moreover, they also will be smoothed by the use of the speed modulation. On the other hand, the ratio of the thermal demand to the minimum heating capacity of the heat pump at each hour determines the time fraction when the heat pump has been working over that hour. Figure 11 shows this ratio, after being sorted. It can be observed that in 58% of the time when the heat pump was working in on/off mode was "on" more than 30 min (hour fraction 0.5), therefore minimizing the transitory effect on the components of the heat pump. In summary, the availability of the inverter is expected to significantly contribute to reducing the energy losses of the on/off operation mode and, consequently, it make sense to neglect it for the scope of the current analysis.

Regarding the costs, Table 9 summarises the results.

**Figure 9.** Thermal demand points in on/off operation mode versus the ambient wet bulb temperature. Heating capacity and compressor consumption at minimum rotational speed are also plotted.

**Figure 10.** Thermal demand met by the heat pump at inverter control mode and on/off one.

**Figure 11.** Hour fraction (time in hours) when the heat pump was "on", while working in on/off mode.


**Table 9.** Levelised costs in the baseline case for heat pump, de-centralised, and centralised boilers.

Figure 12 shows the breakdown of the levelised costs that are given in Table 9. As the demand and the heated area are established, the percentage breakdown is the same for both levelised costs (see Equation (29a,b)). Moreover, Figure 12 points out the predominance of operating costs (especially the energy term, being the maintenance less significant) over investment (heat pump and boiler).

The proposed active measure should be supplemented by the application of a social tariff based on the cost breakdown. This would make it possible to take advantage of the environmental benefits of the proposed system at the same levelised cost of the most economical system, i.e., the centralised boiler. In this case, a discount of 23.75% in the electricity cost (overall cost, including power and energy terms) would reduce the LCOH to 73.54 €/MWh, matching the cost of the centralised boiler. This discount is in accordance with the current social electricity tariff in Spain [59], which ranges from 25% for vulnerable consumers, 40% for severely vulnerable consumers, and up to 100% for consumers at risk of social exclusion. The average discount was 32% when considering the distribution of each cluster in 2019 (respectively, 648,826 and 630,086 households in the first two categories, and 4545 in the third one [60]). Therefore, the required discount to spread the proposed system would be even lower than the one currently applied to vulnerable-households' electricity bill.

**Figure 12.** Percentage breakdown of the levelised costs.

### *3.3. Influence of Energy Retrofitting*

The baseline case has been selected to be representative of the middle-energy efficiency level pool of dwellings, with an energy performance assessment of E, close to D (EPI = 2.18). However, in the context of energy poverty, it is usual to find dwellings built before 1980 with worse insulation condition than this middle energy-efficiency-level. In these cases, an energy retrofitting of the dwelling is recommended, which makes it possible to reduce the heating cost per dwelling. To assess this fact, Figure 13 has been obtained, varying EPI and maintaining the levelised cost of heating in demand base with the same value as the baseline case. This condition requires varying the heated area, obtaining a potential fitting curve (Equation (34)). In Figure 13 the average EPI obtained for dwelling blocks in Spain have been represented over such fitting curve, so obtaining their heated area to maintain the same levelised cost than in the baseline case.

$$A \;= 11,139 \cdot EPI^{-0.789} \tag{34}$$

### **4. Discussion**

The proposed demand model makes it possible to obtain the hourly demand profile of a building, which is essential for working out the instantaneous consumptions of the heating system. The shape of the annual cumulative heating demand curve is consistent with the profiles that were obtained by other simulation tools [40]. This model only requires the location of the building, its useful area and the energy performance certification (EPC). A detailed simulation of the building thermal behaviour is necessary to assess the EPC, which leads to the EPI that is integrated in the proposed model. The fact that no specific additional information about the building is required by the model makes it very easy-to-use and helpful for planning purposes. The flexibility of the model has been used in order to solve a case study representative of vulnerable dwellings.

**Figure 13.** Required heated area (A) for each energy performance index (EPI) to maintain the same *LCOHDB* than in the baseline case (92.22 €/MWh) and *LCOHAB* obtained. Energy performance classes (A to G) limits are indicated in dashed lines.

Regarding the heat pump, R290 (propane) has revealed as a suitable refrigerant for this application, and it is also supported by many manufacturers. It is a natural fluid, with zero ODP and very low GWP. The low compressor discharge temperature and the low de-superheating zone in the condenser advise taking this fluid into consideration for the future. Its high flammability is not relevant in the current case, since it is a centralised unit, with roof-top allocation and maintenance performed by professional workers. The speed control in the compressor gives to the heat pump the capability to meet nearly all of the heating demand (96% in the baseline case) with good efficiency. In the baseline case, the heating seasonal performance factor achieved is 2.585, exceeding the minimum that is required to consider the thermal energy taken from the air as renewable energy. Furthermore, the overall renewable energy to meet the heating demand is obtained, taking both the energy supplied by the air and the renewable share in the electricity mix into account (applied to the electricity consumption of the heat pump). Therefore, each kWh-th met with the proposed system avoids the emission of 131.7 g CO2, whereas 265.3 g CO2 are emitted if a boiler is used.

Regarding the cost for the baseline case, the proposed technology (heat pump plus boiler back-up) reduces the levelised cost of heating by 15% with respect to the de-centralised boiler. However, the proposed system increases the cost by 25% with respect to the centralised boiler. On the other hand, with centralised gas neither renewables sources are employed nor carbon dioxide emissions are avoided. The cost breakdown reveals that the energy cost, which facilitates the integration of subsidy policies to the operation, is the most important contribution to the overall cost. Therefore, a discount of 23.75% in electricity bill (comparable with the current social-tariff average-discount) would be enough to equalise the levelised cost of heating of the proposed system with the centralised boiler (the most economical scenario). Furthermore, the innovative and social aspect of this system makes it possible to integrate funding by non-profit organisations, energy companies, Public Administration, or even the actual consumers (depending on their vulnerability level). These subsidies to the proposed system would be endorsed due to its excellent environmental performance, balancing this way de-carbonisation with energy affordability.

Finally, the effect of energy retrofitting has been investigated. The results show that, if the energy performance assessment is improved from E to C due to a retrofit, the same heat pump would be able to meet the demand of 160 dwellings of 100 m<sup>2</sup> instead of 60 of the baseline case. Regarding the

heating cost per square meter, the cost drops 4 €/m<sup>2</sup> for every EPI unit reduction. Accordingly, in the same example, the levelised cost of heating decreases 500 € for a dwelling of 100 m2 (with a baseline cost of 848 €).

### **5. Conclusions**

The feasibility of heating a block dwelling in Madrid by an air-source heat-pump water heater has been investigated in the framework of energy poverty research. The implementation of the heat pump is planned as a retrofit of the heating system. Therefore, an existing system that is based on radiators is assumed, and the operation temperatures of the heat pump are adapted to this configuration. A global methodology has been used to forecast the hourly heating demand while using an expansion method that is based on the Spanish regulation. The performance map has been obtained after selecting an eco-friendly refrigerant (R290), and considering a speed variation control over the compressor. This control makes it possible to meet nearly the whole demand with a high efficiency. Finally, the demand has been coupled to the performance of the heat pump.

When considering environmental and efficiency indicators, the obtained results show an excellent behaviour of the heat-pump as compared with the classical solution of the natural gas centralised boiler system. However, the levelised cost of heating is 25% higher than the centralised boiler, due to the low prices of gas for high consumption volumes. In this sense, a discount in electricity bill for vulnerable households might be taken into account, considering both social and environmental benefits, to equalise the energy costs of this technology with the one of a centralised gas boiler. The situation is reversed in the case of comparing the heat pump with the de-centralised boiler solution, with the heating cost of the heat pump being 15% lower than the boiler one. In all the scenarios, the main contribution to the overall cost is the operation, especially the energy cost.

Therefore, the heat pump has revealed as an efficient and sustainable system to tackle energy poverty in a typical case of vulnerable households, such as the building blocks analysed, especially when compared with a de-centralised boiler heating system. In any case, a new electricity tariff frame that reduces the operation costs of heat pumps with respect to the centralised boiler solution would be highly recommended.

In future works, the model will be expanded to summer season (cooling demand), taking advantage of the reversible ability of heat pumps. In the context of energy poverty, this is a new trend, especially in countries in Southern Europe. Other technologies, such as thermally driven heat pumps (both absorption and internal combustion engine), might be also considered to be alternatives. Moreover, the application to the different climatic zones of Spain will be carried out and policy implications will be pointed out.

**Author Contributions:** Conceptualization, J.I.L. and R.B.; methodology, E.A.; software, I.P.; validation, J.C.R. and E.C.; formal analysis, J.I.L. and R.B.; investigation, I.P.; resources, E.A.; writing—original draft preparation, J.I.L.; writing—review and editing, R.B. and E.A.; supervision, J.C.R.; project administration, J.C.R.; funding acquisition, E.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Chair of Energy and Poverty of Comillas Pontifical University.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Applicability of Swaging as an Alternative for the Fabrication of Accident-Tolerant Fuel Cladding**

### **Dae Yun Kim 1, You Na Lee 1, Joon Han Kim 2, Yonghee Kim 3,\* and Young Soo Yoon 1,\***


Received: 25 May 2020; Accepted: 17 June 2020; Published: 19 June 2020

**Abstract:** We suggest an alternative to conventional coating methods for accident-tolerant fuel (ATF) cladding. A Zircaloy-4 tube was inserted into metal tubes of different materials and the inserted tubes were subjected to physical force at room temperature. The manufactured tube exhibited a pseudo-single tube (PST) structure and had higher thermal stability than a Zircaloy-4 tube. Optical microscopy and scanning electron microscopy images showed that the PST had a uniform and well-bonded interface structure, i.e., no gaps or voids were found at the interface between the inner and outer tubes. Energy-dispersive X-ray spectroscopy analysis confirmed that the metal components did not interdiffuse at the interface of the PST, even after being kept at 600 and 900 ◦C for 1 h and rapidly cooled to room temperature. Unlike pure Zircaloy-4 tubes, Zircaloy-4/stainless use steel (SUS) 316 PST did not show significant structural collapse, even after being stored at 1200 ◦C for 1 h. Based on these results, if a PST was fabricated using a Zircaloy-4 tube thinner than the Zircaloy-4 tube used in this study and an outer tube of micron-scale thickness, swaging may be a feasible alternative to Zircaloy-4-based ATF cladding.

**Keywords:** alternative process; non-coating method; room-temperature swaging; pseudo-single tube (PST); accident-tolerant fuel (ATF) cladding

### **1. Introduction**

Nuclear fuel claddings are crucial core materials in nuclear power plants. They effectively transfer the heat generated by the fission reaction to the coolant, while preventing the fuel and fission products from leaking into the coolant. Therefore, selecting cladding materials with high corrosion resistance, suitable mechanical properties, and the ability to withstand high pressures, temperatures, and irradiation is vital. Zirconium alloys have shown high potential as cladding materials, owing to their low neutron absorption rate, high corrosion resistance, and stable mechanical properties under the operating conditions of the reactor [1–4]. However, zirconium alloys were found to lack sufficient physical and chemical stability in the Fukushima nuclear power plant accident, highlighting the necessity for research on the development of accident-tolerant fuel (ATF) claddings with improved physicochemical stability. Several alternative ATF claddings have been developed in and outside of Korea, including the M5 physical vapor deposition (PVD)-coated cladding by AREVA [5], the Fe-Cr-Al alloy by ORNL [6], the SiC/SiCf composite by MIT [7], and the Cr laser coating by KAERI [8,9]. However, the ATF cladding fabrication techniques described above have the following disadvantages: it is difficult to achieve uniform surface modifications, such as coatings, on long fuel rods, and issues such as peeling can occur due to thermal shock [10]; FeCrAl alloys experience neutron absorption and

tritium emission problems [11]; SiC/SiCf composite claddings consist of ceramics and are therefore exceedingly unstable in terms of fatigue behavior; and, finally, the abovementioned claddings are difficult to mass-produce, and also incur high production costs [12,13]. To overcome these shortcomings, we propose swaging as an alternative method for the production of ATF claddings. Swaging is a simple process and can be carried out at room temperature. In general, the term "swaging" refers to the process of creating a tube with an outer diameter of a desired size by applying a physical force toward the center of a tube with a larger outer diameter. In this study, swaging indicates the process whereby a tube is inserted into another tube that has an inner diameter greater than the outer diameter of the inserted tube; subsequently, these tubes are joined by applying compressive and tensile stresses toward the center and along the lengths to form a pseudo-single tube (PST) without gaps. If a metal or alloy material with physicochemical properties suitable for nuclear power plant operations is produced in the form of a tube into which Zircaloy-4 can be inserted, there is a possibility that a PST of a length of several meters (more than 1 m) can be easily manufactured using swaging. This process requiresa4t/cm2 load rolling press, and the finished PST has an increased length and a reduced outer diameter compared with the starting tube. The entire swaging process is performed at room temperature; hence, there are no heat-induced formations of microstructures or phase deformations during the process. A heat treatment test confirmed that the interface of the PST remained thermally stable at high temperatures. Furthermore, with swaging, surface modification is possible without length constraints and the results of heat treatment suggest that the PST is suitable for use as an ATF cladding [14–16]. The biggest advantage of swaging compared to the coating method is that uniform cladding can be manufactured in a short time at low cost. This advantage can be increased as the size of cladding increases, which leads to the expectation of mass productivity. The coating method may be more advantageous than swaging to set the optimum point of the cladding thickness but there may be many shortcomings in terms of practical use [17]. When cladding is formed over the entire area of a metric cylindrical tube, it is necessary to verify whether the optimal thickness can be formed uniformly. This is important because it is directly related to the safety of nuclear energy. The produced PSTs can also be applied in other industries in the form of metal tubes, which can be manufactured with the desired physical and chemical properties, depending on the choice of materials and the thickness of the outer and inner tubes.

Although in this study, a sample with an ATF geometry (cladding length and thickness) that is directly applicable to nuclear operations has not been used, a 90 cm-long PST was constructed of stainless use steel (SUS) 316 (outside) and Zircaloy-4 (inside) and its high-temperature stability was confirmed. Based on this, we suggest that swaging can be an alternative process for the fabrication of ATF claddings. For this purpose, the interface between the two metals of a PST tube was analyzed before and after high-temperature heat treatment by optical microscopy (OM) and scanning electron microscopy (SEM). Furthermore, Zircaloy-4 cladding and Zircaloy-4 with SUS 316 on the surface, i.e., PST, were exposed to high temperatures for the same durations and their physicochemical stabilities were compared. The core of this study is to suggest a simple and practical process method that can compensate for the disadvantages of the existing coating method. The most frequently mentioned cladding materials are Cr or Cr-based alloys and studies to optimize them are still incomplete. For the swaging process, a cladding tube manufactured in the form of a pipe is required and the difficulty of its manufacture serves as a barrier to the progress of this study [18]. In this situation, we have conducted prior studies on the applicability of the swaging process for ATF cladding by focusing on the process method using SUS 316, which has excellent oxidation resistance, heat resistance, and workability and is readily available. Although SUS 316 is not the optimal material for ATF cladding, we think it is persuasive to use it to confirm the possibility of the swaging process.

### **2. Experimental**

### *2.1. Swaging Method*

The inner tube comprised Zircaloy-4 (outer diameter: 9.57 mm, inner diameter: 8.30 mm, length: 1 m) and the outer tube comprised SUS 316 (outer diameter: 11.90 mm; inner diameter: 9.82 mm; length: 1 m). The surfaces of Zircaloy-4 and SUS 316 tubes were cleaned by wiping with ethanol. A double tube was prepared by inserting the Zircaloy-4 tube into the SUS 316 tube. The inner tube was completely filled with water-soluble KNO3 filler powder, so that the tube would not experience distortion during compression in the direction perpendicular to the tube surface and also so that the stress could be uniformly applied in all directions. The PST was synthesized by applying a pressure of 4 t/cm<sup>2</sup> in the direction of the central axis of the double tube. Meanwhile, it was difficult to quantify the magnitude of the tensile stress in the transverse direction, which arose from the 4 t/cm<sup>2</sup> compressive stress applied in the longitudinal direction. Each time the tube was passed through the shaft jig, the speed of the tube was manually controlled, so that the outer diameter of the outer tube was reduced by 0.5 mm of the total. Figure 1 shows a schematic of the swaging method for the double tube and the stress type applied to the tube during the process. Figure 2 shows that the inner surface of the outer tube and the outer surface of the inner tube adhere closely to each other. The compressive stress applied during swaging and the tensile force applied in the longitudinal direction result in a strong adhesive force. In addition, the thickness, inner and outer diameters, and length of the final tube can be adjusted, depending on the number of repetitions. The KNO3 left inside the PST was removed by immersing the PST in water at room temperature.

**Figure 1.** Schematic view of the swaging method for a double tube and the stress type applied to the tube during the process.

**Figure 2.** Mechanism of swaging.

### *2.2. Heat Treatment*

Cut pieces of Zircaloy-4 cladding tubes and PSTs with a thickness of 2 cm were prepared and heat-treated in air under the following conditions: starting from room temperature (0~25 ◦C), after they had reached 600, 900, or 1200 ◦C, at a rate of 100 ◦C/h, each tube was maintained at the specified temperature for 1 h. The tubes were left to cool to 200 ◦C before being taken out of the furnace. Figure 3 shows a schematic of the heat treatment and the treatment conditions.

**Figure 3.** Schematic and conditions of heat treatment.

### *2.3. Analysis*

To compare the degrees of oxidation and the stability of the interfaces after exposure to high temperatures, 2 cm-long tubes were prepared as follows: each tube was immersed in epoxy and dried until it solidified. Subsequently, the tubes immersed in the solidified epoxy were cut, their cross-sections were polished, and the microstructures of the polished surfaces were analyzed by OM. An SEM analysis was performed to examine the interfacial structure obtained from the OM at high resolution. EDS was used to analyze the diffusion behavior of the metal components at the PST interface after high-temperature heat treatment. In addition, because the two tubes were adhered solely by physical force, the cross-sections were analyzed using transmission electron microscopy (TEM) to confirm the characteristics of the interface at a higher resolution. The cross-sections of the PSTs were prepared for the TEM analysis by focused ion beam (FIB) milling.

### **3. Results and Discussion**

Figure 4a shows the fabricated 90 cm-long PST and Figure 4b shows its cross-section before and after swaging. From Figure 4b it can be seen that some KNO3 filler remained inside the tube with the smaller internal diameter after swaging, which could be easily removed with water. Figure 5 shows the OM image of three points on the interface of the PST. This image confirms that there were no gaps or voids in all three points of the interface between the two tubes, which was brought about by applying stress radially inward and lengthwise. Furthermore, the bonding interface between the two tubes was determined to be the same along the entire circumference. The uniformity of the interface of the PST was further confirmed by SEM images of the metal interface, as shown in Figure 6. It is imperative to note that gaps in the interface can degrade the mechanical and chemical properties of the tube. An interface without three-dimensional defects, such as gaps or voids, does not undergo oxidation reactions; therefore, it does not experience mechanical and chemical decomposition, even in high temperatures or corrosive environments [19]. As shown in Figure 6, the PST fabricated by swaging exhibits an exceedingly dense interface without three-dimensional defects.

**Figure 4.** *Cont.*

(**b**)

**Figure 4.** (**a**) Photographs of the final pseudo-single tube (PST), (**b**) Comparison of size change before and after swaging.

**Figure 5.** Optical microscopy image of PST without heat treatment. Schematic of the PST manufacturing process by swaging.

A peculiarity was that the tube thickness increased by approximately 5% from 3.35 to 3.52 mm after swaging, which can be explained as follows: during swaging, powder fillers were subjected to compressive stress in the direction of the center of the tube. Owing to plastic deformation, the inner diameter of the tube did not return to its original state, even after removal of the compressive stress [20]. During the initial stages of swaging, the tube length remained unchanged, thus increasing the tube thickness.

**Figure 6.** Secondary electron microscope image of the PST interface without heat treatment.

Once the filling powder entered a high-density state in which it could not be compressed further during swaging, the metal tube was also no longer compressed under external compressive stress. After this step, the tube started becoming longer and thinner due to the tensile stress applied to the tube. When using the swaging process, the magnitudes of the compressive stress acting perpendicular to the axial longitudinal direction and the tensile stress acting in the longitudinal direction were the most important variables for determining the length and thickness of the final PST. It is crucial to design a PST such that the initial sum of the thickness of each tube is as close as possible to the thickness of the single Zircaloy-4 tube currently in use. In particular, after completion of the process, the outer tube should ideally be suited to possess a thickness that modifies the chemical properties of the surface with little effect on the mechanical properties of the Zircaloy-4 and the degree of neutron absorption. Therefore, it is preferable to start the process with an outer tube as thin as possible. However, as mentioned above, in this study, the swaging process was carried out using a commercially available thick SUS 316 tube for the purpose of proving that Zircaloy-4 and dissimilar metals can be produced in the PST form using the described tube process. We are currently investigating the relationship between compressive stress and tube length and thickness.

At room temperature, the physical pressure applied toward the center of the tube and in the longitudinal direction of the two tubes leads to the formation of an interface with a strong bonding force. This can be explained by the mechanism, as illustrated in Figure 2. During swaging, the tube was subjected to compressive and tensile stresses in the axial and longitudinal directions. However, the tensile stress acting in the axial direction of the tube was not generated by the tube being pulled from the outside in the swaging equipment. This stress is thought to have been generated by the KNO3 filler pressing toward the center. Tensile stress occurred in the longitudinal direction of the tube, as the inner diameter of the tube decreased, when a force of 4 t/cm<sup>2</sup> was applied to the center of the tube. This can be understood from the fact that the tube length increased after swaging, even if no tensile stress was applied. The most frequently mentioned cladding materials are Cr or Cr-based alloys and studies to optimize them are still incomplete. Because of this, it is not yet possible to measure or calculate the exact magnitude of the tensile stress generated in the longitudinal direction by the reduction of the inner diameter of the tube with the equipment used in this study. Clearly, however, the most important variable in the shaft process is the control of the compressive stress applied in the direction of the center of the shaft. We are designing and manufacturing new tube equipment that can accurately measure and control the magnitude of the tensile stresses created in the longitudinal direction induced by compressive stresses, along with the values of the compressive stresses.

Differences in physical properties, such as in the moduli of elasticity, shear moduli, and tensile strengths of the two metals after the tube process, result in different shrinkage behavior after stretching. Table 1 shows various physical properties of Zircaloy-4 and SUS 316. The total strain generated by swaging consists of a combination of elastic strain and plastic strain. An elastic deformation results in a return to the original state when the applied stress is removed; however, some permanent deformation, that is, plastic deformation, remained as the yield strength was exceeded. Finally, this plastic deformation resulted in a PST which is longer than the initial tube and has a smaller inner diameter than the initial tube [21]. As can be inferred from Table 1, Zircaloy-4 and SUS 316 tubes have different yield strengths, which can result in different final shrinkages. The first tube was 90 cm in length but post the tube process, the lengths of the Zircaloy-4 and SUS 316 tubes increased to 97.5 and 95.4 cm, respectively. Table 2 shows the sizes of the tubes before and after swaging. It can be seen that the outer and inner diameters decrease after the shaft.

**Table 1.** Physical properties of Zircaloy-4 and stainless use steel (SUS) 316.


**Table 2.** Outer diameters and inner diameters of the metal tubes before and after swaging.


Thus, SUS 316, which experienced a relatively small plastic deformation, shrank more than Zircaloy-4. It is likely that the residual stress generated at the interface between the metal surface with a low plastic strain (high shrinkage) and the metal surface with a high plastic strain (low shrinkage) caused the two metal surfaces to adhere strongly to each other. In other words, SUS 316, which experienced a high degree of shrinkage, generated compressive stress on the surface of Zircaloy-4, which experienced only a low degree of shrinkage. Meanwhile, tensile stress was generated on the surface of SUS 316. The mechanism of the strong physical adhesion of the tubes, owing to the different residual stresses generated by the shaft tube process, is shown in Figure 7. This mechanism leads to a stable interfacial structure between the metal surfaces, as shown in the cross-sectional SEM image of the tube after swaging (Figure 6). Upon close observation of the interface between the two metals in Figure 6, it is evident that there is a small corrugated interfacial structure. In other words, the compressive stress generated on the surface of Zircaloy-4 due to the different shrinkage amounts of the two metals caused the Zircaloy-4 surface to shrink again, thereby forming wrinkles on the surface. This corrugated structure may have resulted in the grasping of the surface of SUS 316, which in turn led to the bonding of the two different metals after room temperature condensation. Such strong adhesive interface behavior was well demonstrated with quenching heat treatment experiments after high-temperature exposure, as shown below.

**Figure 7.** Mechanism for strong physical adhesion of different tubes due to residual stresses generated by the shaft tube process.

To evaluate the thermal stability of the PST fabricated using the room-temperature condensation process, an experiment was performed where the PST underwent rapid cooling after exposure to high temperature under ambient conditions. Figure 8 shows an image of a Zircaloy-4 tube and a PST quenched to room temperature after being maintained at temperatures of 600, 900, or 1200 ◦C for 1 h. The color of the bare Zircaloy-4 cladding changed due to oxidation, even at 600 ◦C. Under 1200 ◦C, the PST disintegrated not only due to severe oxidation phenomenon but also experienced external surface peeling [22].

(**a**) **Figure 8.** *Cont.*

**Figure 8.** Appearance after maintaining (**a**) the Zircaloy-4-only tube and (**b**) the PST at 600, 900, and 1200 ◦C for 1 h.

(**b**)

Noticeably, the structural breakdown of the Zircaloy-4 tube after being exposed to high temperatures is exceedingly similar to the high-temperature disintegration of natural limestone [23]. Thermal shock tests of natural limestone concluded that its oxide layer is formed by an oxidation reaction at high temperatures. As the resulting oxide layer exhibits a large thermal expansion difference compared with the non-oxidized layer inside (owing to the differences in thermal shock received during cooling) and an immensely high amount of thermal stress accumulates at the interface between the oxide layer and the non-oxidized layer, the oxidized surface layer can be peeled off easily. In contrast, for the PST manufactured by the tube process, no surface peeling or structural failure was observed, even after heat treatment at 1200 ◦C. This is because the external material, SUS 316, which has a higher oxidation resistance than Zircaloy-4, inhibited the oxidation of Zircaloy-4 by preventing oxygen migration to the surface of Zircaloy-4. To observe the interface more clearly, the SEM analysis results of the PST cross-section before and after the heat treatment are shown in Figure 9. No collapse was observed at the interface of the quenched sample after exposure at 600 and 900 ◦C. In contrast, at 1200 ◦C, the interface of the quenched sample showed exfoliation and porosity due to volume expansion. It is assumed that this is caused by a large amount of oxygen diffusion from the center of Zircaloy-4 to the interface between the two metals. However, even when the PST was exposed to 1200 ◦C, the oxidation reaction of the Zircaloy-4 surface was suppressed, along with subsequent structural collapse.

**Figure 9.** *Cont.*

(**c**) 1200 °C

**Figure 9.** Scanning electron microscopy (SEM) images of cross-sections of the PST maintained at 600, 900, and 1200 ◦C for 1 h.

To confirm the diffusion behavior of the PST, an EDS analysis was conducted, as shown in Figure 10. After heat treatment at 600 and 900 ◦C, the diffusion of Zr did not occur in the SUS 316 direction and no metallic elements of SUS 316 were seen in the Zircaloy-4 direction. Compositional analysis by EDS line scanning showed that the interface between Zircaloy-4 and PST did not degrade even when exposed to high temperatures. However, after heat treatment at 1200 ◦C, a well-known eutectic reaction caused an interdiffusion reaction, in which an Fe-Zr complex with a melting point of 928 ◦C was formed [24]. Slight diffusion occurred at 1200 ◦C but it still exhibited high stability at operating and accident conditions. These results show that a PST consisting of a metal surface and Zircaloy-4 (e.g., 10Al-90Cr, etc.) [25], and fabricated using the tube process, has the potential to be used as ATF cladding. The interface of the sample, when exposed to 1200 ◦C, as well 600 ◦C and 900 ◦C, did not appear to peel off easily. From this, it can be concluded that an ATF cladding comprising the PST can maintain stability even at 1200 ◦C because the inside of the tube is isolated from external conditions during actual use. TEM analysis was performed to confirm the interfacial stability of the PST after exposure to high temperatures. Figure 11 is a cross-sectional TEM image that shows the interfacial structure of Zircaloy-4/SUS 316 of the PST when quenched to room temperature after being maintained at 900 ◦C for 1 h. Rapid cooling to room temperature after exposure to high temperatures is generally well-suited for inducing significant thermal stresses at an interface to observe the physical stability of the interface. Images were obtained at three different magnifications to observe the interfacial condition over a wide range of PST cross-sections and nanoscale. The width of the PST interface was 0.03 μm (30 nm) after heat treatment, as elucidated from the highest-magnification images. Moreover, no new interface was formed between the secondary phase and each tube due to the interdiffusion of the two tube components. This is in good agreement with the results of EDS, shown in Figure 10. Therefore, it is evident that the interface between two different tubes can maintain physicochemical stability, even when the PST formed by swaging is exposed to high temperatures. These results suggest that the development of several m-grade ATF claddings based on Zircaloy-4 is feasible, if conditions such as the type of tube used, its physical shape, and the optimum swaging pressure applied to the tube are established.

*Energies* **2020**, *13*, 3182

**Figure 10.** Energy-dispersive X-ray spectra of the PST before and after heat treatments at 600, 900, and 1200 ◦C for 1 h.

**Figure 11.** Cross-sectional transmission electron microscopy (TEM) image showing the interfacial structure of the Zircaloy-4/SUS 316 PST quenched to room temperature after being maintained at 900 ◦C for 1 h.

The results show that the tube process has a high potential for the development of an ATF cladding with a length of several meters, based on Zircaloy-4. The results suggest that controlling the compressive stress applied along the axial direction of the tube and the tensile stress occurring in the longitudinal direction will enable the facile and rapid development of ATF claddings, with their geometries calculated according to the design.

### **4. Conclusions**

A Zircaloy-4-based PST, over 90 cm in length, was successfully fabricated using a simple straight tube process that was carried out at room temperature. Zircaloy-4, a material exhibiting well-known physicochemical properties for nuclear power operation, was inserted into a SUS 316 tube. It was possible to fix the SUS 316 tube to the outer surface of Zircaloy-4 by applying compressive stress toward the center of the tube, such that it exhibited strong adhesion. PST cross-sectional observations using OM and SEM revealed no three-dimensional defects, such as gaps or voids, at the interfaces of the two metal tubes. To confirm the high-temperature stability and thermal shock behavior of the PST, it was quenched to room temperature after being kept for 1 h at 600, 900, and 1200 ◦C. It was inferred from the OM and SEM analyses that no macroscopic defects were formed at the interface between SUS 316 and the Zircaloy-4 cladding after heat treatment. Even if the prepared PST was maintained at 900 ◦C, an interdiffusion reaction between the elements constituting Zircaloy-4 and SUS 316 would not occur. The Zircaloy-4-based PST fabricated using a room-temperature condensation process has been found to have highly stable physicochemical properties, even on quenching from high temperatures. Although SUS 316 is not the optimal cladding material, PST using this showed better thermal stability than bare Zircaloy-4 tube and it was confirmed that the swaging process can be applied to ATF cladding production. In conclusion, if (1) a thin external dissimilar metal tube and (2) a conduit device capable of controlling the compressive stress applied during the condensation process and the tensile stress in the longitudinal direction are used, then ATF claddings, which are widely used in nuclear power plant operations, can easily be manufactured to a length of up to 4 m. In addition, the swaging method described herein is likely to find potential use for a variety of applications, especially in fields where PST structures are crucial, i.e., when specific surface properties and uniform physical properties throughout the tube are critical.

**Author Contributions:** Conceptualization, Y.S.Y.; methodology, Y.S.Y. and J.H.K.; software, J.H.K.; validation, Y.K. and J.H.K.; formal analysis, D.Y.K. and Y.S.Y.; investigation, D.Y.K. and Y.N.L.; resources, Y.K.; data curation, J.H.K. and Y.K.; writing—original draft preparation, D.Y.K.; writing—review and editing, Y.N.L. and Y.S.Y.; visualization, Y.S.Y.; supervision, Y.K. and Y.S.Y.; project administration, Y.S.Y.; funding acquisition, Y.K. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Engineering Research Center of Excellence Program of Korea Ministry of Science, ICT and Future Planning (MSIP)/National Research Foundation of Korea (NRF) grant number NRF-2008-0062609 and was funded by Korea Hydro and Nuclear Power Co. Ltd. (KHNP) grant number No. 2018-safety-10.

**Acknowledgments:** This work was supported by the Engineering Research Center of Excellence Program of Korea Ministry of Science, ICT and Future Planning (MSIP)/National Research Foundation of Korea (NRF) (Grant NRF-2008-0062609). This work was also supported by Korea Hydro and Nuclear Power Co. Ltd. (No. 2018-safety-10).

**Conflicts of Interest:** The authors declare no conflict of interest.

### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Application of the 2-D Tre**ff**tz Method for Identification of Flow Boiling Heat Transfer Coe**ffi**cient in a Rectangular MiniChannel**

**Mirosław Grabowski 1,\*, Sylwia Ho ˙zejowska 2, Beata Maciejewska 2, Krzysztof Płaczkowski <sup>1</sup> and Mieczysław E. Poniewski <sup>1</sup>**


Received: 25 May 2020; Accepted: 21 July 2020; Published: 2 August 2020

**Abstract:** The study presents the experimental and numeric heat transfer investigations in flow boiling of water through an asymmetrically heated, rectangular and horizontal minichannel, with transparent side walls. A dedicated system was designed to record images of two-phase flow structures using a high-speed video camera with a synchronous movement system. The images were analyzed with Matlab 2019a scripts for determination of the void fraction for each pattern of two-phase flow structures observed. The experimental data measured during the experimental runs included inlet and outlet temperature, temperature at three internal points of the heater body, volume flux of the flowing water, inlet pressure, pressure drop, current and the voltage drop in the heater power supply. The flows were investigated at Reynolds number characteristic of laminar flow. The mathematical model assumed the heat transfer process in the measurement module to be steady-state with temperature independent thermal properties of solids and flowing fluid. The defined two inverse heat transfer problems were solved with the Trefftz method with two sets of T- functions. Graphs were used to represent: the boiling curves, the local void fraction values, the boiling heat transfer coefficients and the errors of both of them for selected mass fluxes and heat fluxes.

**Keywords:** minichannel flow boiling; void fraction; inverse heat transfer problem; Trefftz method

### **1. Introduction**

A growing necessity of transferring very high heat fluxes from both miniaturized industrial equipment and home appliances generates demand for mini- or micro-scale cooling devices. Some cooling systems require low pumping power and consequently low Reynolds numbers for the flow [1–3]. In such cases, flow boiling heat transfer, which is characteristic of high heat transfer coefficients, appears to be the appropriate solution. The determination of the heat transfer coefficient requires knowledge of the parameters of the boiling fluid flowing in the minichannel. These parameters are in particular: temperature and pressure at characteristic points of the system, mass flux of the boiling fluid, heat flux and void fraction. The measurements of the void fraction were combined with photographic recording of boiling two-phase flow structures, characteristic for horizontal orientation of the minichannel [4,5].

The obtained experimental data for water flow boiling in a single rectangular minichannel at a low Reynolds number provided the basis for numeric analysis with the use of Trefftz method [6]. Direct heat conduction problems—DHCPs—Appear in mathematical modeling of steady-state heat transfer problems when the governing equation, domain, boundary conditions and physical properties

of the material are known. When any of these elements is unknown, we must deal with an inverse heat conduction problem—IHCP. High sensitivity to uncertainty of incoming data are a characteristic feature of many engineering inverse problems. This weakness can get enhanced when two or three succeeding inverse problems are taken into account. Therefore, solving inverse problems requires efficient and stable numeric methods. As the Trefftz method with a set of T-functions fulfills this requirement [7–11], it was used to solve the IHCP problem for flow boiling in the minichannel. Two sets of T-functions were applied to calculate: (a) two-dimensional temperature distributions of the heating copper block and the boiling water, (b) heating copper block temperature gradients, and (c) the heat transfer coefficient at the contact surface copper block–flowing fluid.

### **2. Experimental Facility**

*2.1. Design of the Flow Loop, Experimental Equipment and Data Collecting Procedure*

The flow loop and its components are presented in Figure 1.

**Figure 1.** Flow loop. 1—Measurement module with the minichannel (details in Figure 3), 2—Heating copper block, 3—Temperature and pressure sensors, 4 DC power supply, 5—Cooler with ventilator, 6—Rotameter, 7—Filter, 8—Pump, 9—Pressure control, 10—Compressed air valves, 11—Compressed air tank, 12—Preheater, 13—Facility control unit, 14—Computer for experiment control + LabView software, 15—High speed camera, Pca—Compressed air pressure sensor.

The minichannel, the basic element of the experimental stand (Figure 2), was constructed by glue bonding three transparent glass plates and a cuboid copper block, Figure 3. The cuboid copper block was a solid base for the minichannel and a mass heater as well. Four flat resistance heaters were placed on steel substrates and sintered to the external surface of the copper block. A TDK Lambda GEN 50–30 power supply provided direct current to the heaters, which generated heat required to initiate flow boiling inside the channel (Figures 1 and 2).

Optiwhite glass was selected for the minichannel walls. Optiwhite is colorless, super transparent float glass containing very limited amount of iron and having the highest light-transmission coefficients. Three transparent walls of the channel enabled proper lighting and recording the boiling two-phase structures with a high-speed camera (Phantom 711, Vision Research).

To prevent uncontrolled heating of the flowing fluid by the incandescent light, the LED system of our own design, based on LED diode Citizen CL-L233-HC13L1-C, was applied. LOCTITE® SI 5145 adhesive was used to glue the minichannel elements. The dimensions of the rectangular minichannel, length = 180 mm, width = 4 mm and depth = 1.5 mm, provided a 6 mm2 cross-section and a 2.18 mm

hydraulic diameter. Five thermocouples (Czaki TP-201) were mounted, one at the inlet, one at the outlet of the minichannel and three inside the body of the cuboid copper block, Figure 3. Two pressure sensors (Kobold 0–2.5 bar) were placed at the inlet and outlet of the channel. The flow of distilled water was generated by a precision gear pump (Tuthill Concord D) with a maximum capacity of 1.5 <sup>×</sup> <sup>10</sup>−<sup>7</sup> m3/s, featuring very stable laminar flow in the range 43 <sup>≤</sup> Re <sup>≤</sup> 229. Flows at low Reynolds number are quite frequently used in miniature cooling systems of electronic devices.

**Figure 2.** General view of the experimental stand; labels as in Figure 1.

**Figure 3.** Copper block with attached glass minichannel and location scheme of resistance heaters and thermocouples—Transverse section and axonometric view; dimensions in mm (pictorial view, not to scale).

In the experimental runs, it was necessary to record large numbers of data of various types. Those included video camera images of two-phase flow structures and numeric data of copper block and flowing fluid temperature, fluid pressure and volume flux. Data collection was done by the modular measurement-control system, Figure 4. The core of it was NI cDAQ-9178 chassis (National Instruments) that holds dependent modules, controlling both the measurements of boiling process parameters and the component assemblies, such as the pump, camera triggering, etc.

The following modules comprised the system, as shown in Figure 4: NI cDAQ-9178—the main module, NI 9214—Temperature measurement (CZAKI K-type TP-201 thermocouples), NI 9239—Voltage measurement (KOBOLD pressure sensors, 0–2.5 bar, NI 9203—Current measurement (KOBOLD pressure difference sensors, 0–2.5 bar), NI 9263—Voltage setting for pump control.

**Figure 4.** Block diagram of the data flow and control system.

In the supplementary measurements of heat loss to the environment, the entire minichannel module was treated as a cuboid with walls of different temperatures measured with the infrared camera. The calculations also accounted for thermal losses from the sight glass. The maximum value of the amount of heat released into the environment in the range of heat fluxes generated in the presented measurements did not exceed 2.1 W, which in relation to the total delivered thermal power was 0.86% to 1.9%. Since the heat losses were minor, they were disregarded in further calculations.

### *2.2. Procedure of Void Fraction Measurement and Computation in Rectangular Minichannel*

Photographs of the observed two-phase flow structures were taken with a high-speed video camera at the recording speed of 7000 fps. The recording speed was selected experimentally to maintain proper exposure and obtain sharp images of dynamic two-phase flow structures. Proper selection of the recording speed is essential for capturing geometric features of the vapor bubbles. The equations approximating spatial geometry of vapor bubbles were applied to convert flat images of bubbles into three dimensional ones, taking simultaneously into account the dimensions of both the bubbles and the minichannel. The shapes of sphere and ellipsoidal cylinder ended with one or two half-ellipsoids were used, depending on inter relations between the dimensions of the bubbles and those of the minichannel. Three cases I, II and III of the characteristic relation between the bubble and minichannel dimensions were specified for void fraction determination in the observed minichannel of length *L*, width *a* and depth *b,* Figures 5–7.

**Figure 5.** (**A**) Actual and (**B**) approximated shapes of vapor bubbles in longitudinal and transverse sections of the channel with flow direction, case I. Water, *Tin* = 63.1 ◦C, *Tout* = 123.8 ◦C, *pin* = 1.165 bar, Δ*p* = 0.0024 bar, *q* = 298.7 kW/m2, *G* = 7.2 kg/(m<sup>2</sup> s), *Re* = 55.0, *Lcam* = 20 mm.

**Figure 6.** (**A**) Actual and (**B**) approximated shapes of vapor bubbles in longitudinal and transverse sections of the channel with flow direction, case II. Water, *Tin* = 55.2 ◦C, *Tout* = 125.9 ◦C, *pin* = 1.0697 bar, Δ*p* = 0.0016 bar, *q* = 260.1 kW/m2, *G* = 8.6 kg/(m<sup>2</sup> s), *Re* = 66.9, *Lcam* = 20 mm

**Figure 7.** (**A**) Actual and (**B**) approximated shapes of vapor bubbles in longitudinal and transverse sections of the channel with flow direction, case III. Water, *Tin* = 72 ◦C, *Tout* = 120 ◦C, *pin* = 1.147 bar, Δ*p* = 0.0016 bar, *q* = 270.5 kW/m2, *G* = 5.9 kg/(m<sup>2</sup> s), *Re* = 43.5, *Lcam* = 20 mm.

Volume of the observed part of the minichannel was:

$$V = L \cdot a \cdot b,\tag{1}$$

### 1. **Small vapor bubbles: Rsb, i** ≤ **b**/**2.**

The void fraction for small bubbles in the minichannel was calculated from the equation [12]:

$$\varphi\_{sb} = \frac{4}{3Lab\sqrt{\pi}} \sum\_{i} A\_{sb,i'}^{\frac{3}{2}} \tag{2}$$

where *Asb*,*<sup>i</sup>* = π*R*<sup>2</sup> *sb*,*<sup>i</sup>* is the cross sectional area of single spherical bubble.

2. Large, elongated bubbles, fully visible: ellipse semi-axes P1 *lb, i*= *a*/2 and P2 *lb, i*= *b*/2.

The void fraction for large, elongated and fully visible vapor bubbles was [12]:

$$
\varphi\_{\text{lb}} = \sum\_{i} \varphi\_{\text{lb},i\prime} \tag{3}
$$

where <sup>ϕ</sup>*lb*,*<sup>i</sup>* <sup>=</sup> (*Llb*,*i*−*a*) <sup>π</sup>*ab* <sup>4</sup> <sup>+</sup> <sup>π</sup>*a*2*<sup>b</sup>* 6 *Lab* is the void fraction for a single elongated bubble, *Llb, i* is the bubble length, composed of the length of an ellipsoidal cylinder and two half-ellipsoids.

3. Large, elongated bubbles, partially visible: ellipse semi-axes P1 *lb, i*= *a*/2 and P2 *lb, i*= *b*/2.

The void fraction for large, elongated and partially visible vapor bubbles was [12]:

$$
\varphi\_{\text{lb}} = \sum\_{i} \varphi\_{\text{lb},i\prime} \tag{4}
$$

where <sup>ϕ</sup>*lb*,*<sup>i</sup>* <sup>=</sup> (*Llb*,*i*<sup>−</sup> *<sup>a</sup>* <sup>2</sup> ) <sup>π</sup>*ab* <sup>4</sup> <sup>+</sup> <sup>π</sup>*a*2*<sup>b</sup>* <sup>12</sup> *Lab* is the void fraction for a single elongated bubble, *Llb, i* is the bubble length, composed of the length of an ellipsoidal cylinder and one half-ellipsoid.

Three scripts were developed to analyze the void fraction for selected cases of two-phase boiling flow structures found in the recorded videos. MathWorks Matlab and two corresponding toolboxes, computer system vision and image processing, were used for that purpose. By applying the first plan method and the Gaussian model, cases I and II were detected and elaborated with the use of the Vision toolbox). The background subtraction method appeared suitable for case III. Sharpening and pixel multiplication applied for each video frame improved quality of the recorded images. These operations were finally followed by equalization of brightness and contrast of the images to obtain a larger number of details, as shown in Figure 8.

**Figure 8.** Picture of two-phase flow (**A**) before and (**B**) after the operations improving its quality.

Subsequently the perfected and binarized void frames underwent the following operations: morphologic opening and closing the picture, median filtration from artefacts and morphologic filling

of empty spaces in the observed objects—That is vapor bubbles. The application of the entire algorithm converted the image in Figure 8B into the black and white one, Figure 9.

**Figure 9.** Image of two-phase flow resulting from morphologic operations.

Indexing and measuring procedure was used to get geometric dimensions of each detected bubble. On the basis of these data and with the use of equations presented in this section, by converting the flat picture of the bubble into the expected spatial one, the sought void fraction was calculated in the observed part of the minichannel. The last action of the script was recording the obtained void fraction value for each video frame into a text file for further use.

The estimated values of the void fraction were compared to seven correlations to calculate the void fraction against vapor quality. The results were similar [12]. The maximum scatter did not exceed 16%.

### **3. Mathematical Model and Numeric Solution**

The mathematical model is a modification of the models described earlier in [2,13]. Analysis of the numeric calculation results presented in [13] showed that the change of the temperature along the width of the copper block body was negligible. In the two-dimensional model presented below, the copper block temperature depended only on two variables: *y* (along its height) and *x* (in relation to its length). Additionally, this assumption allowed considering heat conduction only in two elements of the measuring module: the copper block and the minichannel with flowing water.

By analogy to [13], we assumed that the heat transfer process in the experimental module was steady-state with constant copper and fluid properties. We also assumed that the temperature of the copper block in domain *<sup>D</sup>*<sup>1</sup> = <sup>+</sup> (*x*, *<sup>y</sup>*) <sup>∈</sup> *<sup>R</sup>*<sup>2</sup> : 0 <sup>&</sup>lt; *<sup>x</sup>* <sup>&</sup>lt; *<sup>L</sup>*, 0 <sup>&</sup>lt; *<sup>y</sup>* <sup>&</sup>lt; *<sup>H</sup>*<sup>2</sup> , satisfied the Laplace equation, i.e.,:

$$
\nabla^2 T\_{\mathbb{C}} = 0.\tag{5}
$$

For Equation (5), we assumed that the temperature of the copper block *TC* was known at three measuring points (*xi*, *H*1), Figure 10.

**Figure 10.** Scheme of the experimental stand with the adopted boundary conditions; (pictorial view, not to scale).

As in [13], in order to stabilize the numeric calculations two temperature values were added to boundary conditions at both ends of the fluid-copper heater contact surface:

$$\text{(a)}\ T\_{\text{C}}(0, H\_2) = 0.5 \left( T\_{\text{in}} + T\_{\text{approx}}(0) \right) \\ \text{b)}\ T\_{\text{C}}(L, H\_2) = 0.5 \left( T\_{\text{out}} + T\_{\text{approx}}(L) \right) \\ \tag{6}$$

where the broken line approximation of three measured temperatures inside the block, Figure 10, was used to obtain *Tapprox*. Both vertical walls of the copper block were insulated.

The heat flux *q* generated by each resistance heater was supplied to the copper block, i.e.,:

$$
\lambda\_{\mathbb{C}} \frac{\partial T\_{\mathbb{C}}(\mathbf{x}, 0)}{\partial y} = -q(\mathbf{x}). \tag{7}
$$

The heat flux *q*(*x)* assumed constant values *qi*, generated by a single heater in the segments, where four resistance heaters were located, Figure 10. Between these segments the heat flux equaled zero:

$$q(\mathbf{x}) = \begin{cases} \quad q\_i \text{ forx} \in D\_{i\prime} \, y = 0\\ \quad 0 \text{ forx} \in [0, L] - D\_{\prime} \, y = 0 \end{cases} \tag{8}$$

where *<sup>D</sup>* <sup>=</sup> <sup>4</sup> ∪ *i*=1 *Di* = [2.5; 41.5] ∪ [47.5; 86.5] ∪ [93.5; 132.5] ∪ [139.5; 178.5].

For fluid we assumed that, as shown in Figure 10:


$$\frac{1}{b} \int\_{H\_2}^{H\_2 + b} w\_x(y) dy = w\_{\text{arc}\_{\prime}} \tag{9}$$

where the *wave* was known


$$T\_f = \begin{cases} \begin{array}{c} T\_{\mathbb{C},s} \text{ if } T\_{\mathbb{C},s} < T\_{\text{sat}}\\ \begin{array}{c} T\_{\text{sat}} \text{ if } T\_{\mathbb{C},s} \ge T\_{\text{sat}} \end{array} . \end{cases} . \tag{10}$$

In domain *<sup>D</sup>*<sup>2</sup> = <sup>+</sup> (*x*, *<sup>y</sup>*) <sup>∈</sup> *<sup>R</sup>*<sup>2</sup> : 0 <sup>&</sup>lt; *<sup>x</sup>* <sup>&</sup>lt; *<sup>L</sup>*, *<sup>H</sup>*<sup>2</sup> <sup>&</sup>lt; *<sup>y</sup>* <sup>&</sup>lt; *<sup>H</sup>*<sup>3</sup> , , Figure 10, the fluid temperature satisfied the energy equation:

$$
\Delta\lambda\_f \nabla^2 T\_f = w\_x(y) c\_p \rho\_f \frac{\partial T\_f}{\partial \mathbf{x}} - \mu\_f \Phi - w\_x(y) \frac{dp}{d\mathbf{x}} + \Omega(\mathbf{x}),\tag{11}
$$

where function Φ = *d wx dy* 2 was the Rayleigh dissipation function, pressure gradient *dp dx* <sup>≈</sup> *pout*−*pin <sup>L</sup>* <sup>=</sup> <sup>Δ</sup>*<sup>p</sup> L* and Ω(*x*) was negative heat source. The heat flux absorbed by the bubbles, also called a negative heat source, was calculated from the formula, which used experimentally determined void fraction φ(*x*) [14–16]:

$$
\Omega(\mathbf{x}) = \frac{\theta \,\alpha\_{\rm com} \Delta T}{2 \mathcal{R}\_b} \varphi(\mathbf{x}). \tag{12}
$$

The bubble diameter 2 *Rb* in subcooled flow boiling was calculated using the correlation proposed in [17]:

$$
\Delta R\_b = \min\{1.4, 0.6 \cdot \exp\left(-\frac{\Delta T}{45}\right)\} \cdot 10^{-3}.\tag{13}
$$

*Energies* **2020**, *13*, 3973

The convective heat transfer coefficient α*con* in Equation (12) was given by Labuncov correlation [18] for laminar flow:

$$a\_{\rm con} = \frac{\lambda\_f}{2R\_b} 0.125 \text{Re}^{0.65} \text{Pr}^{\frac{1}{5}} \tag{14}$$

and Δ *T* was the difference between the fluid temperature *Tf* in the thermal sublayer δ *<sup>T</sup>* and the temperature inside the vapor bubble, which assumed to equal *Tsat*. Consequently, the fluid superheat Δ *T* was calculated from the dependence:

$$
\Delta T = T\_f - T\_{\text{sat}} = \frac{q \cdot \delta\_T}{\lambda\_f}.\tag{15}
$$

In (15) the thermal and hydraulic layer thicknesses were calculated from formulas [16], respectively:

$$
\Delta T = T\_f - T\_{\text{sat}} = \frac{q \cdot \delta\_T}{\lambda\_f}.
\
\delta\_T = \text{Pr}^{-\frac{1}{3}} \delta\_{\text{lt}},
\delta\_{\text{lt}} = \frac{2\mu}{f \, w\_b \, \rho\_f}.\tag{16}
$$

where *wb*(*y*) = *wx*(*y*) and the Fanning friction factor was calculated as in [19]:

$$f \text{Re}\_f = 24 \left( 1 - 1.3553 \frac{b}{L} + 1.9467 \left( \frac{b}{L} \right)^2 - 1.7012 \left( \frac{b}{L} \right)^3 + 0.9564 \left( \frac{b}{L} \right)^4 - 0.2537 \left( \frac{b}{L} \right)^5 \right). \tag{17}$$

The known copper block temperature distributions and the temperature gradient were applied to determine the heat transfer coefficient α(*x*) at the copper–fluid interface using the Robin boundary condition:

$$a(\mathbf{x}) = \frac{-\lambda\_{\mathbb{C}} \frac{\partial T\_{\mathbb{C}}}{\partial \mathbf{y}}(\mathbf{x}, H\_2)}{T\_{\mathbb{C}}(\mathbf{x}, H\_2) - T\_{f, \text{ave}}(\mathbf{x})}. \tag{18}$$

The reference temperature *Tf, ave* was the average fluid temperature along the depth of the fluid layer equal to the diameter of emerging vapor bubble:

$$T\_{f, \text{ave}}(\mathbf{x}) = \frac{1}{2\mathcal{R}\_b} \int\_{H\_2}^{H\_2 + 2\mathcal{R}\_b} T\_f(\mathbf{x}, y) dy. \tag{19}$$

The energy Equations (5) and (11) with a set of the adopted boundary conditions led to the solution of two succeeding inverse heat conduction problems (IHCPs) in two domains of different shapes and parameters. The Trefftz method with two sets of adequate T-functions was used to solve the problem.

The boundary conditions Equations (6)–(8) and the harmonic functions (T-functions), defined in [6], were employed to solve Equation (5). Linear combination of T-functions approximated the copper block temperature. These coefficients were found by minimization of the mean square error between the approximated temperature and the boundary conditions. T-functions for homogeneous energy equation, corresponding to energy Equation (11) were applied to find the flowing fluid temperature [20].

The fluid temperature was approximated by a sum of particular solution of Equation (11) and the linear combination of the T-functions shown in [9]. Similar to the case of copper block temperature, the boundary conditions for the fluid were adopted to compute the coefficients of that linear combination.

### **4. Results and Discussion**

The measurement procedure included recording the following quantities: temperature and pressure at the channel inlet and outlet, differences in inlet and outlet pressures, volume flux of the flowing medium, temperature at three selected points inside the copper block and also current and voltage of the heater power supply. The measurements were taken for the following ranges of the experimental parameters: total heat flux generated by four external flat heaters 188 <sup>≤</sup> *<sup>q</sup>* <sup>≤</sup> 340 kW/m2, inlet pressure 1.05 ≤ *p* ≤ 1.17 bar (practically constant), inlet fluid subcooling 3.6 ≤ Δ *T* ≤ 70.7 K and mass flux 2.2 <sup>≤</sup> *<sup>G</sup>* <sup>≤</sup> 8.6 kg/(m2 s), 43 <sup>≤</sup> *Re* <sup>≤</sup> 229.

At the same time, flow boiling in the minichannel was filmed using high-speed video camera and local void fractions ϕ*(x)* were measured on the basis of recorded images. The scripts developed in Section 2.2 allowed the calculation of the void fraction for each pattern of two-phase flow structures observed in this study. The experimental local values of the void fraction, measured along the minichannel length, were approximated with broken lines, as shown in Figure 11.

**Figure 11.** Examples of recorded flow structures presented in (**A1**): *q* = 298.8 kW/m2, *G* = 8.6 kg/(m<sup>2</sup> s); (**B1**): *q* = 223.2 kW/m2, *G* = 2.2 kg/(m2 s). Corresponding void fractions (blue points), approximated with broken lines (red lines), are shown in (**A2**) and (**B2**).

The results of numeric calculations presented in Figure 12 were obtained for experimental parameters given in Section 4. In the first step, nine T-complete functions were used to calculate the approximate two-dimensional temperature of the copper block *TC* for Equation (5). Next, by combining four T-complete functions for Equation (20) and particular solution of Equation (11), the field of the fluid temperature *Tf* was found.

**Figure 12.** Two-dimensional temperature fields of (**a**) the copper block and (**b**) flowing fluid obtained by the Trefftz method for (**A**) heat flux *q* = 223.2 kW/m<sup>2</sup> and mass flux *G* = 5.9 kg/(m<sup>2</sup> s) and (**B**) heat flux *q* = 340 kW/m2,mass flux *G* = 5.9 kg/(m<sup>2</sup> s), (pictorial view, not to scale).

Fluid flowing into the minichannel substantially reduces the temperature of the copper block across the contact area, Figure 12. The fluid motion in the minichannel causes heat transfer in the direction of its axis. As a result, the temperature in the entire channel increases as does the void fraction, which impairs heat transfer between the fluid and the heater. For this reason, in the final section of the module, both the liquid temperature and the heater temperature increase significantly.

The temperature distribution in both the fluid and the heating block, Figure 12, is as expected, the isotherms run in a manner likely for the direction and intensity of flow and location of the heaters. The charts are not fully comparable, because both parameters, *G* and *q*, are different for each case, but the impact of the parameter change on the temperature field is clear.

The results of numeric calculations presented in Figure 12 were obtained for experimental parameters given in Section 4. In the first step, nine T-complete functions were used to calculate the approximate two-dimensional temperature of the copper block TC for Equation (5). Next, by combining four T-complete functions for Equation (20) and particular solution of Equation (11), the field of the fluid temperature *Tf* was found.

The boundary conditions adopted for the copper block and liquid temperatures were satisfied with high accuracy. The maximum differences between the calculated and measured temperatures of the copper block did not exceed 24 K (5.93%) for the first measuring point (from the inlet to the minichannel), 19 K (4.45%) for the second and 16 K (3.72%) for the third measurement point. Throughout the study, the largest errors were always at the first measuring point and the smallest at the last. When identifying the liquid temperature, the maximum differences between the calculated and measured temperature did not exceed 13 K (3.71%) at the inlet to the minichannel and 15 K (4.65%) at the outlet from the minichannel. As in the case of thermocouples placed in the copper block, the differences between the measured and calculated temperatures for liquids at the inlet to the minichannel were in most cases higher than the differences at the outlet from the minichannel.

Variation of the heat transfer coefficient for the flowing two-phase mixture is shown in Figure 13 as the function of the distance from the minichannel inlet. It is presented for three selected mass fluxes and varying heat flux. The heat transfer coefficient decreases as the distance from the inlet grows, which is related to increasing void fraction. Some variation in this trend can be observed both at the channel inlet and outlet. The growing heat flux also increases the heat transfer coefficient. Analysis of graphs in Figure 13 shows a known pattern—an increase in flow intensity causes an increase in heat transfer coefficient.

**Figure 13.** Heat transfer coefficient as a function of minichannel length (**a**) for *G* = 8.6 kg/(m2s), (**b**) *G* = 5.9 kg/(m<sup>2</sup> s) and (**c**) *G* = 2.2 kg/(m2s).

In the initial segment of the minichannel, changes in the calculated values of heat transfer coefficients as a function of distance from the inlet, Figure 13, are non-physical in nature. This is the result of boundary condition (6) adopted to stabilize the numeric procedure.

The shape of boiling curves presented in Figure 14 for various distances from the minichannel inlet is similar to the initial section of standard boiling curve. The curves originate with segments of steep dependence of the heat flux *q* versus the temperature difference *T*C, s − *T*sat and flatten with further increase of that difference. This results from the growing vapor phase (void fraction) and approaching first boiling crisis. In the discussed experiment, the vapor was not superheated to prevent possible thermal damage of the research facility. Therefore, the boiling curves are limited to their initial segments.

**Figure 14.** Boiling curves for selected points *x* along the minichannel length for (**a**) *G* = 2.2 kg /(m2 s), (**b**) *G* = 5.9 kg /(m2 s) and (**c**) *G* = 8.6 kg /(m2 s).

The smaller the distance from the inlet, the greater is the heat flux q for the same temperature difference between the heating surface and the boiling liquid. This result agrees qualitatively with the results shown in Figure 13, because the heat transfer coefficient increases with the decreasing distance to the inlet.

With increasing distance from the channel inlet, the temperature difference *T*C, s − *T*sat must increase to induce heat transfer of expected value. This is due to the increase in the amount of gaseous phase in the two-phase mixture and thus in its insulating properties.

The uncertainties of the experimental parameters determined in study [2] were applied to calculate the mean relative error (*MRE)* of the heat transfer coefficient α (*x*). Figure 15 shows the *MRE* as a function of both heat flux and mass flux. When the heat flux supplied to the copper block increases, the *MRE* decreases from the value of 4.2% for heat flux *q* = 188 kW/m2 down to 1.0% for *q* = 340 kW/m2. For smaller mass fluxes, the *MRE* values are lower.

**Figure 15.** Mean relative error vs. (**a**) heat flux *q*, and (**b**) mass flux.

### **5. Conclusions**

The experimental data gained at the Warsaw University of Technology, Plock Campus created a ground for the completed numeric analysis. The design of the experimental setup allowed setting and automatically acquiring temperature of the copper block and flowing fluid, fluid pressure, volumetric flow rate, heat flux and video recording for various positions of the video camera along the minichannel length. During the experimental run the research schedule of the setup was under control of the program created in the LabView environment. Data acquisition was controlled by the same script as well. The configuration of the setup is flexible enough to provide the potential for further expansion of the research program for new settings of thermal and flow parameters, channel space orientations and channel geometry.

The unique feature of the photographic method developed in this research is the video registration of the observed two-phase flow structure coupled with the measurement of the corresponding local value of the void fraction. This idea is based on the conversion of the flat image of the vapor bubble into the expected spatial shape. The obtained values of the void fraction were successfully applied in boiling two-phase flow modeling and numeric calculations. Disadvantages of the photographic method include 1) the requirement of quick collection of a large amount of data and 2) laborious conversion and refinement of the data into the required void fraction values with dedicated Matlab scripts.

The mathematical approach to solving two sequential IHCPs in flow boiling in the asymmetrically heated rectangular and horizontal minichannel was validated against experimental data. The employed meshless Trefftz method, making use of two types of T-complete functions, allowed obtaining approximate solutions of two coupled energy equations with their boundary conditions. The method provided stable solution of IHCP and consequently granted determination of: 1) the two-dimensional temperature distributions in two contacting domains, having different shapes and thermal features (copper block and flowing fluid), 2) the heat flux transferred from the block to the flowing two-phase mixture, and 3) the heat transfer coefficient on the contact surface between the copper block and flowing fluid. The Trefftz method is not limited by the number or character of boundary conditions. They can be established directly from temperatures or temperature gradients and they may have continuous or discrete form as well. It is worth noting that the approximate solutions satisfy the governing partial differential equations exactly and numeric calculations are not very much complicated because T-complete functions are polynomials in the investigated case.

In the near future the research program will be extended to include water, ethanol and FC72 and various minichannel space orientations [21].

**Author Contributions:** Conceptualization, S.H., M.E.P.; methodology, M.G, K.P., S.H. and B.M.; validation, S.H., B.M.; formal analysis, S.H.; investigation, M.G. and K.P.; writing—Original draft preparation, M.G., S.H. and M.E.P.; writing—Review & editing, S.H., B.M. and M.E.P.; funding acquisition, M.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** The research was funded by the Warsaw University of Technology, Plock Campus, grant no. 504/04480/7193/44.000000 and the Poland National Science Center, grant no. UMO-2018/31/B/ST8/01199.

**Conflicts of Interest:** The authors declare no conflicts of interest.

### **Nomenclature**

#### **Nomenclature**



### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
