Air-source ammonia heat pump for district heating: a field modeling approach with focus on frosting-defrosting cycles

News

Air-source ammonia heat pump for district heating: a field modeling approach with focus on frosting-defrosting cycles

Products | R&D

LUVE Leading the low charge NH3 revolution – ammonia evaporators

Get rights and content

Under a Creative Commons license

Publication: Applied Thermal Engineering

Publisher: Elsevier

Date: 1 January 2026

Original publication link

Highlights

Evaporator field optimization for district heating heat pumps.
Dedicated model for frost formation in fin-and-tube evaporators.
Experimental validation with transient tests in frosting and defrosting.
Numerous field configurations and defrosting strategies considered.
Best design options largely different when including or excluding frost formation.

Abstract

Frost formation on air-source evaporators in cold and humid climates significantly increases energy consumption due to airflow obstruction and thermal insulation. While this issue is well recognized, effective control strategies for multi-evaporator fields remain underexplored. This paper presents a novel evaporator model that combines experimental validation under transient frosting conditions with a new algorithm for solving frost accumulation and defrosting process across an entire evaporator field coupled to the heat pump cycle.

The model is applied to simulate a large-scale ammonia heat pump serving a Danish district heating network over an annual period. Results show that neglecting frosting in design analyses can lead to strongly overestimated performance, with total costs increasing by up to 70 % once frost and defrosting cycles are considered. The study also identifies an optimal field configuration of 80 evaporators with a fin pitch of 9 mm and an airflow velocity of 3.5 m/s, and shows that defrosting 20 % of the field simultaneously is more effective than defrosting 10 %. Overall, the outcomes emphasize the critical importance of accounting for frost and defrosting cycles in thermo-economic optimization, showing that the optimal design depends on local climatic conditions and should be supported by dedicated modelling tools.

Nomenclature

Symbols

Specific heat, J/(kgK)

Diameter, m

Heat transfer coefficient, W/(m²K)

Specific enthalpy, kJ/kg

i_sv

Latent heat of desublimation, J/kg (2830 kJ/kg)

K_m

Mass transfer coefficient, kg/(m²s)

Mass flow rate, kg/s

N_r

Number of tube rows

Pressure, Pa

Pitch, m

Thermal power, kW

Radius, m

Thickness, m

Internal energy, J

Uncertainty, [%]

Overall heat exchanger conductance (W/K)

Velocity, m/s

Greek lettersα

Thermal diffusivity, (m²/s)

δ_f

Frost layer thickness, m

Efficiency, [-]

Thermal conductivity, W/(mK)

Density, kg/m³

Time, s

Specific humidity, kg_H2O/kg_air

Subscripts $co$

Condensing

Evaporation

Frost

int

Intermediate

Iso-entropic

Liquid-to-vapor

ref

Refrigerant

sur

Surface

Tube

Water vapor

1. Introduction

District heating networks (DHN), first introduced in the late 19th century, have consistently played a crucial role in meeting the heating demand of the building sector across numerous countries [1]. Current international policy anticipates that a significant portion of residential heating demand (approximately 20 %) will be supplied by sustainable and carbon-free district heating systems [2]. In this evolving energy landscape, large-scale vapor compression heat pumps emerge as an efficient technology with the potential to significantly contribute to the decarbonization of heating systems and foster greater integration between the heating and power sectors by increasing the utilization of renewable energy sources [3]. Air-source heat pumps, in particular, offer the notable advantage of flexible installation, typically employing extensive fin-and-tube air-refrigerant heat exchangers as evaporators [4].

However, a key operational challenge arises as district heating networks are frequently situated in geographical areas characterized by cold and wet weather conditions. Such environments can negatively impact the proper functioning of the evaporator field due to the occurrence of frost formation [5]. Frost formation on air-refrigerant fin-and-tube heat exchangers has been a widely studied issue for decades in both refrigeration and air conditioning systems [6]. Experimental studies have consistently shown that the reduction of the air-side flow area and the insulating properties of frost are primary contributors to evaporator performance degradation [7]. The properties of frost and its formation dynamics are highly dependent on environmental and surface boundary conditions [8]. Consequently, detailed experimental investigations have been conducted on various surfaces, including plates [9], cylinders [10], and surfaces with specific characteristics like interaction with water drops or finning [11,12]. Beyond experimental efforts, diverse methods have been developed to predict frost behavior, including its growth and properties during the formation transient [6]. These methods range from empirical correlations based on reduced parameters [13] to more localized approaches [14]. The latter are typically founded on transient convective and diffusive heat and mass transfer equations, enabling the prediction of property distribution within the frost layer. Notably, local approach analyses have served as the foundation for developing detailed models of fin-and-tube heat exchangers that account for frost formation [[15], [16], [17]].

Once frost accumulates on the evaporator surface, an appropriate defrosting technique and strategy become essential to ensure the continuous operation of the system. The primary objectives of defrosting are twofold: achieving optimal timing to prevent both excessive frosting and unnecessary defrosting, and implementing an effective technique to clean the surface as efficiently as possible. For large-scale heat pumps, hot gas bypass defrosting is a benchmark process [18]. This method utilizes a portion of the hot refrigerant vapor from the compressor outlet to alternately defrost one or a group of evaporators in the field. Extensive experimental and modeling activities have been dedicated to predicting the time required to melt a given amount of frost [19,20], particularly with the aim of analyzing subsequent frosting/defrosting cycles.

The combined effect of frost formation and the defrosting process invariably leads to a deviation from the nominal performance of the heat pump. This deviation is more pronounced under critical weather conditions. While this performance degradation has been thoroughly investigated for small and medium-sized heat pumps [21], there remains a notable gap in knowledge for large-scale systems. The complexity of multiple evaporators operating in parallel within a field, mutually influencing each other’s airside flow, makes this problem particularly challenging from both experimental and modeling perspectives. To date, this specific topic for large systems has been addressed only by experimental studies. For instance, the work by Li X. [22] et al. or Liang et al. [23] evaluates the effect of alternate defrosting of evaporators for both direct and indirect expansion systems. Few evaporators are involved in these studies compared to dozens of modules making up a field. At this purpose, on field experimental studies have been conducted by Hansen et al [4] or Liang S. et al. [24] but the latter regards a field of air-source heat pumps instead of a single system. To the best of our knowledge, no modelling studies focusing on evaporator fields are available in literature. Although detailed modeling approaches for frost formation and melting processes are available for single fin-and-tube heat exchangers, they are often not directly applicable to large, complex multi-evaporator systems. Current methodologies for such systems frequently rely on simplified SCOP or COP penalization terms to account for frost formation. This is the case of the work by Wei et al. [25], in which the penalty associated with frosting and defrosting cycles is summarized by a correlation based on ambient temperature. As indicated by Cabrol and Rowley [26], the use of this method requires correction coefficients based on air humidity calibrated for the specific machine. Such methods may not accurately capture the intricacies of this complex phenomenon. For instance, evaporator surface conditions are fundamental for estimating frost formation. Furthermore, other methods of COP evaluation, as investigated by Pieper et al. [27], do not consider frost formation.

Against this background, the present work aims to develop a comprehensive model for a large-size air-source heat pump. The novelty of the study is related to several aspects. Firstly, the evaporator model is validated with dedicated experiments under transient frosting conditions, which are quite rare in literature and provide a solid basis to assess the model reliability. Secondly, a new algorithm is introduced to solve the frost accumulation and defrosting process across the entire evaporator field, explicitly coupled with the heat pump cycle, thus having realistic simulations at system level. Finally, the model is applied to perform seasonal analyses that quantify the penalization on performance caused by recurrent frosting and defrosting cycles, allowing the identification of optimal field configurations and defrosting strategies. The impact of fin and tube geometrical characteristics is therefore assessed to provide practical guidelines for the thermos-economic design of large evaporator fields.

The paper is organized as follows: first, the system model is described, with particular attention to local heat transfer within the evaporator under frosting and melting conditions. Subsequently, the model is calibrated and validated using dedicated experimental tests. Finally, a broad domain of evaporator field configurations (in terms of geometrical and field size values) is investigated under frost-free conditions and evaluated based on total costs, revealing a more constrained area of thermo-economic convenience in design conditions. These identified optimal configurations are then simulated, incorporating frosting/defrosting cycles, to determine the true minimum total costs.

2. System description and general equations

This section presents the description of the system model, starting from the heat pump model and focusing on the detailed modelling of the evaporator field with frosting and defrosting cycles.

2.1. District heating network

The heat load of a district heating network is related to two main factors: the utilization factor

f_{ut}

(i.e. the number of buildings out of the total requiring heat at a certain time of the day) and the ambient temperature that is the main variable involved in the heat thermal load evaluation for a building in the heating season [28]. Therefore, the heat rate required by the network is:(1)

{\dot{Q}}_{DHN} = f_{ut} \dot{Q} (T_{amb}) = f_{ut} [(\frac{{\dot{Q}}_{des}}{T_{amb, d e s} - 16}) (T_{amb} + 16)]

The heat thermal load

\dot{Q} (T_{amb})

is a linear function of the ambient temperature, defined by the design heat thermal load

{\dot{Q}}_{des}

at the design temperature

T_{amb, d e s}

(i.e. the ambient temperature covering the 97 % of the possible winter season temperature [28] and the off temperature (16 °C). The utilization factor changes along the day and its magnitude is different depending on the month of the winter season as shown in Fig. 1 that is an example pattern taken from the work of Noussan et al. [29].

2.2. Heat pump scheme and model

The architecture of the large-size air-source ammonia heat pump considered in this work is one of the common configurations for this kind of systems [27]: a two-stage compression thermodynamic cycle with an intermediate pressure liquid receiver and a second liquid receiver on the low-pressure side. The low-pressure receiver connects the low-stage compressor with the air-refrigerant evaporator field, allowing a wet direct expansion working mode. The system scheme is shown in Fig. 2 where all the main components are shown including also the branch related to the hot-gas bypass defrosting system (highlighted in red) whose details are described in the following section.

The main assumptions of the thermodynamic cycle, whose specific enthalpy-pressure diagram is shown in Fig. 3, are the following:

•
The inlet vapor quality for both compression stages (1 and 3 in Fig. 2) is equal to 1.
•
Refrigerant in condition (5) is subcooled or saturated liquid, while it is saturated liquid for condition (7).
•
The intermediate pressure $p_{int}$ is evaluated as known in literature to get the same pressure ratio for the two stages $p_{int} = p_{3} = \sqrt{p_{1} \cdot p_{co}}$ [3].

Apart from the evaporator field, the main components of the system are modelled in a simplified way. Expansion valves are modelled with iso-enthalpic transformations. The compressors’ outlet conditions (2), (4) in Fig. 2) are evaluated through the global efficiency

η_{c}

and the iso-entropic outlet conditions as stated in Eq. (2).(2)

i_{2 (4)} = i_{1 (3)} + \frac{(i_{2 (4) i s} - i_{1 (3)})}{η_{c_{1 (2)}}}

Global efficiency curves for screw compressors are taken from literature [27]. The condenser is not modelled in this work, and the condensation temperature is assumed as a constant, having a regulation on the water flow rate to the district heating network. Since the condenser must balance the heat network load (

{\dot{Q}}_{DHN} = {\dot{Q}}_{co})

, the condenser mass flow rate

{\dot{m}}_{co}

is evaluated as depicted in Eq. (3).(3)

{\dot{m}}_{co} = \frac{{\dot{Q}}_{co}}{(i_{4} - i_{5})}

The energy balance on the intermediate liquid receiver leads to the intermediate mass flow rate

{\dot{m}}_{int}

of Eq. (4).(4)

{\dot{m}}_{int} = \frac{{\dot{m}}_{co} (i_{3} - i_{6}) + {\dot{m}}_{hgbp} (i_{3} - i_{2})}{(i_{2} - i_{7})}

Furthermore, the thermodynamic condition on the return line from the evaporator field is evaluated through the energy balance on the low-pressure liquid receiver through Eq. (5):(5)

i_{11} = \frac{{\dot{m}}_{ev} i_{9} + {\dot{m}}_{int} (i_{1} - i_{8}) + {\dot{m}}_{hgbp} (i_{1} - i_{15})}{{\dot{m}}_{ev}}

The refrigerant mass flow rate

{\dot{m}}_{ev}

circulating from the low-pressure liquid receiver to the evaporator field is defined through the coupling of a fixed quadratic characteristic curve for the parallel set of recirculating pumps considered and the pressure drop calculation along the low-pressure circuit. Given the mass flow rates, the part load regulation factor

(\frac{n}{n_{0}})

is evaluated for screw compressors from the following equation:(6)

{\dot{m}}_{co (i n t)} + {\dot{m}}_{hgbp} = ρ_{1 (3)} {\dot{V}}_{\max}_{1 (3)} η_{V_{C_{1 (2)}}} {(\frac{n}{n_{0}})}_{C_{1 (2)}}

where volumetric efficiency

η_{V}

presents a decreasing linear relationship with respect of pressure ratio and

{\dot{V}}_{\max}_{1 (3)}

is the ideal volumetric flow rate according to the compressor displacement. Therefore, compressors consumption, the evaporator field cooling capacity and the COP are defined as follows.(7)

\begin{matrix} {\dot{L}}_{C_{1}} + {\dot{L}}_{C_{2}} = ({\dot{m}}_{int} + {\dot{m}}_{hgbp}) (i_{2} - i_{1}) + ({\dot{m}}_{co} + {\dot{m}}_{hgbp}) (i_{4} - i_{3}) \\ {\dot{Q}}_{ev} = {\dot{m}}_{ev} Δ i_{lv} x_{11} \\ C O P (θ) = \frac{{\dot{Q}}_{DHN}}{{\dot{L}}_{C_{1}} + {\dot{L}}_{C_{2}} + {\dot{L}}_{aux}} \end{matrix}

The consumption from auxiliary components (

{\dot{L}}_{aux}

) refers to recirculating pumps and axial fans of the evaporator field.

2.3. Hot gas bypass system

As stated in the introduction, the heat pump modelled uses a hot gas bypass system to realize defrosting. The common technique to define the defrosting process start time is to use a time-based strategy that is active only when climatic conditions are favourable for frost formation through air temperature and humidity measurement combined with evaporation temperature one [30]. Furthermore, a practical method used to define the end of the defrosting cycle is through a set-point temperature of the coil to be reached after the melting of frost. According to the ASHRAE technical specifications regarding ammonia systems [18], the basic schematization of the hot gas bypass system is showed in Fig. 4.

The defrosting process happens at an intermediate pressure

p_{relief}

fixed by the DRR (“defrost relief regulator”) while the HGSV (“hot-gas solenoid valve”) adjusts the mass flow rate

{\dot{m}}_{hgbp}

according to the subcooling control at the evaporator outlet. The vapor quality related to the hot-gas bypass mass flow rate

x_{hgbp}

is evaluated through an iso-enthalpic transformation from the relief pressure to the evaporating pressure. This mass flow rate is provided through an increment of the compressor speed in order not to reduce the mass flow rate needed by the condenser to meet the district heating network load. For this reason, there is an additional energy consumption

E_{hgbp}

to be considered in the whole system energy balance as described in Eq. (8).(8)

E_{hgbp} = {\dot{m}}_{hgbp} Δ i_{1 - 4} θ_{def}

The total time of the defrost process

θ_{def}

takes into account the time needed to melt the frost and heat the coil resulting from the model of the defrost transient described in the following sections. Moreover, an additional fixed time is considered for the starting transient related to the pressurization of the evaporator at the relief pressure and the final dripping phase realized through the activation of fans before working conditions are restored [18].

2.4. Evaporator field

The evaporator field is made by several fin and tube heat exchangers working in parallel as shown in Fig. 5(a) having a wet direct expansion on the refrigerant side (outlet vapor quality less than 1). Each evaporator can be in one of the three statuses: free surface, frost formation or defrosting. For this reason, a global conductance of the field

U A_{field}

could be defined as the sum of the conductance of the active evaporators

N_{ev, a c t}

(i.e. those that are not defrosting):(9)

\begin{matrix} U A_{field} = \sum_{n = 1}^{N_{ev, a c t}} (U {A)}_{n} \\ {\dot{Q}}_{ev} = U A_{field} Δ T_{air - r e f .} \end{matrix}

This global conductance is the key parameter to define the coupling evaporation temperature (i.e. saturation temperature within low-pressure receiver) between the field and the rest of the system. Each evaporator of the field is equipped with axial fans. In case of partial load conditions, a 3-step regulation for axial fan rotating speed is considered to keep at least a minimum temperature difference between the inlet ambient air temperature and the evaporating refrigerant to 3 °C. Also, the refrigerant mass flow rate

{\dot{m}}_{ev}

regulation is made de-activating some of the pumps in parallel after the low-pressure liquid receiver.

The use of horizontal evaporators can lead to air recirculation, especially when the air is blown upwards. Also, their proximity in a field can cause mutual interference in the air flows. This results in the so-called cold island effect [24] (i.e. lower average inlet temperature and higher relative humidity than undisturbed air) with air flows that differ from the nominal values. Furthermore, this effect is strongly dependent on external factor such as wind [31] or the occurrence of frost formation alters the air flow opening angle at the fan outlet, modifying the potential recirculation. Given the considerable uncertainty regarding the recirculation factors to be implemented and the lack of experimental data, this effect is not considered.

2.5. Matching of components defining the current balancing point

Evaporating conditions within low-pressure receiver must be such as to satisfy the cooling capacity defined by both the thermodynamic cycle and the heat transfer modelling of the evaporator field.

Assigned a first guess saturation temperature within low-pressure receiver (

{T_{8} = T}_{9} = T_{11} = T_{1}

), the thermodynamic cycle is established by the set of equations described within Section 2.2. Modelling the evaporator field leads to an evaporator field cooling capacity depending on the global conductance expressed in Eq. (9). The balancing point is the one that allows to have the same cooling capacity from Eq. (9) and Eq. (7) for the guess saturation temperature.

For the guess value there will be a mismatching between the two mentioned values of the cooling capacity. A new guess value can be iteratively assigned in order to reduce this error below a threshold value established for the present work at 1 %.

3. Evaporator modelling

The modelling approach used in this work is based on dividing the heat exchanger coil in several elements

i

along each tube

j

of the refrigerant circuit for all the tube rows

k

as depicted in Fig. 5(b). The total evaporator cooling capacity

{\dot{Q}}_{ev}

is equal to the sum of the elementary heat transfer contributions

{\dot{Q}}_{i, j, k}

. Each element is modelled according to a finite volume approach, defining a suitable calculation mesh.(10)

{\dot{Q}}_{ev} = \sum_{k} \sum_{j} \sum_{i} {\dot{Q}}_{i, j, k}

The heat transfer domain of a single evaporator element is depicted in Fig. 6 using a radial symmetry with respect of the tube axis and an axial symmetry with respect of the fin axis. Apart from the fin and the tube, a growing frost thickness is considered whenever the surface temperature is below 0 °C and the dew point of moist air. According to the finite volume approach, it is defined a grid with a single row of M volumes for the tube and a single column of N volumes for the fin, since fin and tube pitches are at least an order of magnitude more than their thicknesses. Instead, the frost thickness is divided in multiple volumes P that adapt their dimension according to the rate of growth defined by frost formation model. The heat transfer for the frost is characterized by convection with air and conduction with the coil (radiation is neglected), except for the overlapping zone at the elbow between the tube and the fin base where only conduction is present (darker blue square in Fig. 6). For each volume V of the grid, the energy balance given by the sum of the heat fluxes on the s edges is defined as follows:(11)

\sum_{s} {\dot{Q}}_{s} = \frac{d U_{V}}{d θ}

The accumulation term in case of frost volumes takes into account only the mass of water, neglecting the air contribution. The details about the characterization of the heat rates

{\dot{Q}}_{s}

are provided in Appendix A1. According to the symmetry, the fin axis, the tip of the fin and the tube end are considered as adiabatic surfaces, while the convective heat transfer with air and refrigerant is evaluated through heat transfer coefficient calculation depicted in Appendix A2.

The sum of the heat rates on the volume edges in contract with the airflow define the heat rate

{\dot{Q}}_{ai r_{i, j, k}}

for the element of the evaporator shown in Fig. 5(b). The same is valid for the volume edges in contact with the refrigerant flow, resulting in the heat rate

{\dot{Q}}_{{ref}_{i, j, k}}

. These contributions are used to estimate the enthalpy variations across the element as depicted in Fig. 5(b).

The refrigerant within each tube-element experiences a pressure drop

{Δ p}_{i, j, k}

that is equal to the sum of a frictional

({Δ p}_{f r_{i, j, k}})

, momentum

({Δ p}_{m_{i, j, k}})

and the gravitational

({Δ p}_{g_{i, j, k}})

contributions.(12)

{Δ p}_{i, j, k} = {Δ p}_{f r_{i, j, k}} + {Δ p}_{m_{i, j, k}} + {Δ p}_{g_{i, j, k}}

Instead, air side pressure drops are evaluated with proprietary correlations provided by the manufacturer. The actual air flow is obtained by matching the air side pressure drop to the axial fan characteristic curve. For this operation, a blockage ratio of the free-flow area due to the frost presence is considered. Details about pressure drop correlations for both air and refrigerant side are provided in Appendix A2.

3.1. Frost formation model

When the surface temperature of the evaporator is both below 0 °C and the dew point of moist air, the model takes into account frost formation on the evaporator surface. According to the state of the art about frost formation models on fin and tube heat exchangers, the sub-model here presented counts the following assumptions [32]:

•
The frost layer growth is one-directional and orthogonal to the evaporator surface, but due to the finite volume approach, different layer thickness will form on the fin and the tube for the single evaporator element.
•
Frost thermophysical properties are considered as constants within the frost layer.
•
Heat transfer within frost layer is only conductive and its thermal conductivity is only dependent on frost density.

The first assumption neglects the portion of frost that is growing on the leading edge of the fins on the airflow frontal area of the evaporator. According to Padhmanabhan et al. [15], the growth of this portion of frost is uninfluential on the airflow chocking until the frontal area is quite completely occupied by the frost having its propensity to grow out of the coil envelope. According to the second hypothesis, frost layer is regarded as a homogenous porous layer. This assumption has been made in line with the numerous models developed in literature, as resumed in the review by Li and Wang [32]. The third assumption comes from the mathematical model proposed by Tao et al. [33] where the thermal conductivity is expressed as a series expansion of frost density.

Frost formation could be assumed as a phenomenon consisting of convective and diffusive heat and mass transfer involving a total mass flux of water vapor

G_{f}

toward a cold surface that leads to frost accumulation [9] along the z direction orthogonal to the surface.(13)

G_{f} = \frac{d}{d θ} \int_{0}^{δ_{f}} ρ_{f} (z) d z

Starting from an initially clean surface, frost growth could be divided in an early and short crystal growth period and a fully developed frost layer growth period that starts when frost crystals branch out and join to create a layer that could be assumed as a macroscopic homogeneous porous medium [33]. In this stage, the total mass flux of water vapor consists of a “growth” contribution

G_{g}

that is related to the convective heat and mass transfer and of a “densification” contribution

G_{d}

that is related to the diffusive mass transfer within the frost layer where a water vapor concentration gradient occurs.(14)

G_{f} = G_{g} + G_{d}

Fig. 7 shows the schematization of the main parameters involved in frost formation for a generic volume of frost of the evaporator element (the z direction becomes the x axis for the frost growing on the fin or the r axis for the frost growing on the tube). The heat continuity boundary condition applies at the frost surface in contact with the airflow as follows:(15)

k_{f} {\frac{d T_{f}}{dz}|}_{z = δ_{s}} = \dot{q} + G_{g} ∙ i_{sv}

Where the sensible contribution is given by the convective heat flux

\dot{q}

and the latent contribution is related only to the growth contribution. Through the analogy between convective heat and mass transfer coefficient using the Lewis number, the total water vapor mass flux could be evaluated as follows:(16)

G_{f} = K_{m} (ω_{air} - ω_{sur})

Additional details about the mass transfer coefficient calculation are provided in Appendix A3. According to Hermes et al. [9], using the boundary condition of Eq. (15), the densification contribution could be expressed as follows:(17)

G_{d} = b ∙ G_{g} ∙ δ_{f} {\frac{d T_{f}}{dz}|}_{z = δ_{s}}

This equation is obtained assuming that, according to the empirical evidence from the first studies of Hayashi et al [34] to the latest ones [6], the density of the frost layer can be expressed with an exponential function of the frost temperature and other boundary condition temperatures in the form

ρ_{f} = a e^{(b T_{f} + T)}

. Through Eq. (17), the overall mass balance in Eq. (14) becomes a second order equation that can be solved with respect of the growth contribution.(18)

\frac{b ∙ δ_{f} ∙ i_{sv}}{k_{f}} G_{g}^{2} + (1 + \frac{b ∙ δ_{f} ∙ \dot{q}}{k_{f}}) G_{g} - G_{f} = 0

Therefore, Eq. (14) can be integrated using an integral approach with respect of time for both frost density and layer thickness, following a common approach in literature, of which the work by Lee et al. [35] is an example.(19)

\{\begin{matrix} \frac{d ρ_{f}}{d θ} = \frac{G_{d}}{δ_{f}} \\ \frac{d δ_{f}}{d θ} = \frac{G_{g}}{ρ_{f}} \end{matrix}

According to Eq. (19), a couple of starting conditions for both initial frost thickness and density are needed. Despite initial frost thickness could be set to very small values (1

μ

m), the initial value for density is a challenging problem. Indeed, this value has a close relationship with the characteristic of ice crystal growth which is strongly dependent from all the boundary conditions of the problem [36]. For these reasons, the approach used in some model with a fixed starting value for density [37,38] cannot be considered the right solution for the prediction of frost formation in a wide range of boundary conditions, while a general correlation should be used. To the best of our knowledge, among all the correlations developed in literature, the one that better, but not completely, covers the range of application of this model is the one by Hermes et al. [39], valid for air temperature between 2 and 7 °C, relative humidity between 75 % and 86 % and a surface temperature between −5 and −10 °C.(20)

ρ_{f} = a e^{(b T_{f} - c T_{air, d e w})} = 494 e^{(0.11 T_{f} - 0.06 T_{air, d e w})}

3.2. Defrost process

The dual process to frost formation involves a source of energy that is able to melt the frost on the surface to restore the free working condition of the evaporator. This process is a complex combination of heat and mass transfer mechanisms [20], involving sensible and latent heat exchange with permeation of melted water within the remaining frost with dragging and dropping of frost pieces along the coil surface. However, according to literature [19,20], some models have been developed to estimate the time needed to melt a certain amount of frost under some simplifying assumptions. In this model, the defrost process is solved with the finite volume approach previously described. In this case, frost volumes do not present a moving boundary, having the melting process leading to a progressive increase to a porosity equal to 1.

The defrosting transient consists in a first stage of sensible heating of the coil and of the frost thickness. When the melting temperature is reached by the frost, the melting process starts. The frost mass variation is the following:(21)

\frac{d m_{f, V}}{d θ} = \frac{d U_{V}}{d θ} \cdot \frac{1}{i_{sv}}

The permeation of melted water and its sensible heat transfer or vaporization of retained water are not considered. However, these contributions account for about 5 % of the heat consumption during defrosting according to the experimental evidence of Song et al. [40]. Instead, the thickness of liquid that forms on the wall can cause entire portions of frost to fall or detach in cases of very intense frost formation. In these cases, the total energy required to clean the evaporator would be less than that needed to melt the deposited mass. Therefore, a potential conservative estimation of defrosting time may result from the use of this model. When all the frost is melted, a dry heating stage of the coil occurs until the set-point temperature that ends the defrost process is reached. The defrosting process transient ends when the total amount of frost has been melted and the set point temperature is reached for all the heat exchanger elements.

The refrigerant side inlet conditions are defined by the hot gas bypass system model. The heat transfer toward the frosted coil leads to the decrease of the refrigerant specific enthalpy having a de-superheating and a subsequent condensation with the final subcooling along the circuit.(22)

{\dot{m}}_{hgbp} Δ i_{ref} + {\dot{Q}}_{re f_{i, j, k}} = 0

For the evaluation of the refrigerant heat rate, the refrigerant side heat transfer coefficient is calculated using Shah correlation while condensing [41] and Dittus-Boelter one [42] while de-superheating or sub-cooling. The heat transfer with ambient air is natural convection since the evaporator axial fans are idle during the defrosting process.

4. Finite volume grid and time-step analysis

The finite volume approach is solved using an implicit method. A proper grid and time-step analysis have been conducted to reduce the computational effort to run seasonal simulations. Some critical combinations of boundary conditions (high temperature difference between air and refrigerant, 98 % relative humidity and maximum fan speed), have been evaluated in terms of time derivative of the frost density and thickness to choose the one to be used for the analysis as shown in Fig. 8.

The simulation with a temperature difference of 20 °C presents the highest magnitude for time derivatives. The minimization of numerical error in this condition ensures that the numerical error is minimized for all other, less critical, conditions. Using the set of conditions in green in Fig. 8, several frost formation transient simulations have been evaluated changing the mesh size and the time-step control logic. The mesh size has been varied using 3 to 30 volumes for the tube, 3 to 110 volumes for the fin and 3 to 9 volumes to divide the frost layer. Fig. 9 presents the numerical error of the predicted mass of frost with respect to the one obtained with the denser mesh for a single evaporator element. The colorbar is related to the number of volumes for the fin since it resulted the most significant dimension of the grid. For this reason, the chosen mesh presents around 50 volumes for the fin and a combination of tube and frost volumes that allows to maintain the elapsed time for the single timestep of the evaporator element below 1e-4 s to have a seasonal simulation completed in several hours.

Furthermore, Table 1 resumes the results for different auto-adaptive time-steps. In this case, the numerical error

E R_{sim}

is defined as the percentage error between the time (

θ_{pen})

needed to accumulate a quantity of frost that leads to a −50 % in the evaporator cooling capacity with respect of the reference timestep of 0.1 s.(23)

E R_{sim} = \frac{{θ_{pen}|}_{auto - a d a p t i v e} - {θ_{pen}|}_{0.1 s}}{{θ_{pen}|}_{0.1 s}} \cdot 100

Apart from a fixed timestep of 0.5 s, auto-adaptive timesteps based on the variation of the evaporator cooling capacity or the mean frost thickness on the evaporator (control variable c) have been considered with different thresholds

Z

needed to double the previous timestep as following:(24)

Δ θ_{θ + 1} = 2 \cdot {Δ θ}_{θ} \leftrightarrow \frac{c_{θ} - c_{θ - 1}}{c_{θ - 1}} \cdot 100 < Z

Therefore, according to Table 1 the auto-adaptive timestep based on the variation of the evaporator heat rate with a 1 % threshold has been chosen for all the simulations executed to get the presented results.

Table 1. Timestep sensitivity analysis for one frosting cycle on a single evaporator.

Timestep type	0.1 s (fixed)	0.5 s (fixed)
N. of iteration	80,387	16,120	58	23	466	117
Running time [min]	1742	349	1.25	0.49	10.1	2.53
$θ_{pen}$ [s]	8040	8061	7911	7783	8083	7956
$E R_{sim}$ [%]	0	0.3	−1.6	−3.2	0.5	−1.0

Numerical errors related to the settings chosen for simulations are much lower than uncertainties related to experimental results. Therefore, their effect on results obtained can be considered contained within the expanded uncertainties mentioned in Appendix B.

5. Calibration and experimental validation

Given the size of these plants, it is very complex to run an experimental campaign on the whole system. However, evaporators are the most critical part of the system, and they are the crucial component of this work. For these reasons, an experimental campaign involving a wavy fin and tube evaporator arranged horizontally as done in DHN evaporator fields has been carried out. Typical environmental conditions of cold climatic locations have been replicated in a climatic chamber to analyse frost formation. Tab. 2 shows the summary of the experimental conditions tested: room temperature has been varied from 15 to −15 °C with high relative humidity values to test both dehumidifying conditions and frost formation ones. These conditions have been repeated changing air velocity values according to the rotating speed of the axial fan that is coupled with the evaporator. The test facility description and all the details about the experimental procedure are resumed in Appendix B1.

Table 2. Summary of experimental campaign test conditions.

Test number	Air velocity [m/s]	Room temperature [°C]	Room RH [%]
1	1.4	15	80
2	1.4	10	80
3	1.4	5	80
4	1.4	−5	75
5	1.4	−15	75
6	2.1	15	80
7	2.1	10	80
8	2.1	5	75
9	2.1	−5	75
10	2.1	−5	80
11	2.1	−15	75
12	3	−5	80

The results coming from these tests have been used to run a comparison with the evaporator model with frost formation. Each test lasted a few hours having subsequent frosting/defrosting cycles as depicted in Fig. 10. In particular, Fig. 10(c) shows the evaporator surface right after the defrosting process is completed. Before restarting the evaporator, the axial fan is switched on to drip off water droplets. This dripping phase is usually effective with a negligible retained liquid mass. Significative data to be used for the calibration and the experimental validation were obtained from the data reduction and uncertainty analysis reported in Appendix B2. The experimental boundary conditions for both air and refrigerant side for each test were given as input to the model to make an experimental comparison regarding the thermal power, air-side pressure drop and the total mass of frost. The overall statistical comparison has been conducted using the mean absolute percentage error (MAPE) and the mean relative percentage error (MRPE) of the predicted values with respect of the experimental one. Fig. 11 shows the overall results for the thermal power and pressure drop. These results have been obtained through a calibration process involving the air-side heat transfer coefficient and pressure drop correlations that is described in Appendix B3.

Furthermore, Fig. 12 presents the comparison in terms of the total mass of frost accumulated on the evaporator surface, showing a good degree of accordance between the model and the experimental results. The points shown refer only to the end of the frosting cycles, which is when the mass was measured.

The frames obtained by the camera have been analysed to obtain an estimate of frost thickness on the first row of the evaporator as shown in Fig. 13. In this figure, frost formation on a single element of the evaporator is shown at the end of the transient with the relative frame from the camera. Instead, the smaller box shows the result for the whole transient.

6. Case study and methodology

The evaporator field investigated in this work is associated to a 7.35 MW heat pump serving a district heating network in Aarhus (Denmark) providing hot water at 62 °C with a return temperature of 30 °C. This field can be realized in a huge number of ways related to the several geometric features of a single fin and tube heat exchanger, the number of evaporators and the axial fan choice, having a potential domain of million design solutions. However, the mean computational time needed to run detailed seasonal simulation does not allow to explore all this domain. At this purpose, the first part of results involves the search for a restricted domain of optimal thermo-economic solutions in design conditions to be investigated on a seasonal base. The idea behind this part is the assumption that the poorest design solutions cannot overcome the best ones in a seasonal performance assessment that considers frost formation which is a cause of performance decay for all the solutions. This restricted domain is chosen defining the best solution in terms of total costs and using a threshold (5 %) in terms of total cost percentage difference with respect of this optimal solution. After that, a group of field configuration is chosen from this subdomain to run seasonal simulations with a sensitivity analysis on field design parameters. Also, different time-based defrosting strategy options will be investigated. According to that, the total number of seasonal simulations is the combination of the field configuration extracted from the optimal subdomain in design conditions tested with each defrosting options. The number of simulations is related to the results obtained.

Total costs

C_{tot}

are evaluated as the sum of the total energy consumption

E_{cons}

and the field investment cost

C_{field}

(provided by the manufacturer) in a 10-year scenario considering a fixed energy price

C_{en}

(0.3 €/kWh).(25)

C_{tot} = E_{cons} \cdot C_{en} \cdot 10 + C_{field}

7. Simulation results

7.1. Preliminary results in design conditions

According to methodology, Fig. 14 shows the starting domain of design solutions investigated, spanning the main variables affecting both the global conductance of the field, such as the air flow velocity and the fin pitch, and its total heat transfer surface such as the field size in terms of the number of evaporators. These variables have been also considered due to their significant influence on frost formation. Other geometrical features leading both to similar effects on heat transfer surface such as the air-side frontal area or the number of tube rows, and on global conductance such as tube and row pitches are assumed as constants. Fig. 14 shows the results in terms of a colormap with respect of the three chosen design variables. The colormap refers to the relative percentage difference in total costs

Δ C_{tot}

, taking the one with the lowest total costs as a reference. The best solution has a field with 80 evaporators, a 9 mm fin pitch with 2.5 m/s for air velocity. The energy consumption for the evaluation of total costs is evaluated considering the COP at the design conditions for these preliminary results.

From the analysis of Fig. 14, it is possible to identify a zone of optimal thermo-economic solutions (within 5 % range) around the best one, having a subset of 48 out of 150 total solutions. The COP for this subset spans from 2.81 to 3.18. This range is comparable to design COP available in literature for installed systems. For instance, a COP of 2.1, 3.4 and 3.5 is registered for Danish installations of 1.2 MW, 2.2 MW and 5 MW respectively [43]. Starting from the relative optimal solution of this preliminary analysis, a restricted simulation domain, depicted in Fig. 14, is defined to investigate how the change of design variables affects frost formation and consequently the seasonal performance. Indeed, 13 solutions are selected considering the variation of one variable at time and moving on the main diagonals of the subset (bubbles circled in purple in Fig. 14).

7.2. Seasonal simulation plan

According to the previous preliminary results, Table 3 resumes the seasonal simulation domain investigated in this work in terms of selected field geometrical configurations and defrosting strategy. The temperature–time based defrosting strategy needs two main input setting that have been included in the variable parameters of the simulation domain: the number of evaporators to be defrosted simultaneously and the defrosting starting time. In particular, a fraction of 10 % and 20 % of the field have been selected for the first input. The second input is based on the mean time needed to reach a −30 % or a −60 % decay of the evaporator cooling capacity for the selected configurations in a reference simulation where the worst combination of boundary conditions in terms of DHN load and weather are given as fixed input during all the frost formation transient. The total permutations of Table 3 (4 defrosting strategies for each of the field configurations) led to a group of 52 seasonal simulations. The model has been developed and executed in MATLAB, using parallel computing techniques for the evaporator field.

Table 3. Seasonal simulation domain as a combination of field configurations extracted from the preliminary analysis and defrosting strategy parameters.

Field configurations				Defrosting strategy parameters
Evaporator field size	Airflow Velocity [m/s]	Fin pitch [mm]	Fin type	Fraction of the field to be simultaneously defrosted	Defrosting start time definition
80	1.5	9	Wavy	10 % and 20 %	20 h and 40 h
80	2.5	9
80	3.5	9
80	2.5	6
80	2.5	12
40	2.5	9
120	2.5	9
40	2.5	3
40	4.5	9
120	1.5	6
60	4.5	15
120	2.5	15
120	1.5	12

An example of the operation of the evaporator field during a seasonal simulation for one of the permutations of Table 3 is presented in Fig. 15 where a 100-hour window is shown. When the cooling capacity of an evaporator is equal to zero the defrosting process is going on. The red dashed box highlights a complete cycle of defrosting for the whole field, that means all evaporator groups have been defrosted, from the first one (blue line) to the last one (grey line). According to defrosting strategy, an interval of 40 h between two consecutive defrosting for the same evaporator is visible. During this time interval, all evaporators are defrosted using a proper temporal shift between each other to avoid simultaneous defrosting.

Seasonal simulation results for all the permutations are reported in the table provided in Appendix C. Instead, the main outcomes are here commented. Fig. 16 shows the same plot of Fig. 14 but related to the subgroup of seasonal simulations characterized by the same defrosting strategy that contains the overall optimal solution. Compared to previous results, the optimal solution (whose edge is highlighted in red) has moved from 2.5 m/s to 3.5 m/s with the same fin pitch and evaporator field size. Apart from that, the main outcome is that the variation of total costs has changed from the preliminary 5 % to a maximum of 70 %, proving the importance of taking into account frosting/defrosting cycles to analyse design conditions that seems to be equivalent to each other in a preliminary analysis. The effect of the design parameters is evaluated changing one parameter at time following the red arrows in Fig. 16.

The effect of the airflow velocity with the variation of the defrosting strategy is shown in Fig. 17(a). Total costs decrease with the increase of airflow velocity that leads to a higher field mean evaporation temperature as shown in Fig. 17(b) for a 100-hour window of the seasonal simulation. The 20 % strategy is more convenient for all the velocities while a lower defrosting frequency is better when the velocity increases.

The higher evaporation temperature leads to a reduction in the compression consumption with an increase for the axial fan power. This trend is shown in Fig. 18 where the total seasonal energy consumption varying the airflow velocity is represented. The energy fraction accounted to HGBP system is related to the additional compressor consumption described in Eq. (8). The solution at 3.5 m/s is slightly better than the one at 2.5 m/s but going beyond this velocity is not convenient anymore since the solution at 4.5 m/s has been discarded in the preliminary analysis.

The effect of the evaporator field size with the variation of the defrosting strategy is presented in Fig. 19(a). In this case the higher defrosting frequency is more convenient when the field size is the lowest because of a more intense frost formation as shown in Fig. 19(b) that represents the field cumulative mass of frost during the heating season. Despite the highest field size presents the lowest accumulation of frost, total costs are lower for the intermediate field size.

At this purpose, Fig. 20(a) shows that increasing the field size brings to a higher evaporation temperature, but the energy consumption related to the defrosting system increases for 80 to 120 evaporators as shown in Fig. 20(b). Indeed, the energy consumption related to hot-gas bypass system is directly proportional to the field size and inversely proportional to the evaporation temperature. Therefore, the small gain in the evaporation temperature passing from 80 to 120 evaporators do not overcome the higher consumption due to a wider field.

Lastly, the analysis considering the fin pitch variation is shown in Fig. 21. According to Fig. 21(a), the effect of the defrosting strategy seems more impactful at the narrower fin pitch. Furthermore, Fig. 21(b) presents the field airflow blockage ratio in terms of mean and maximum values for the whole heating season. The narrower fin pitch presents a higher maximum blockage ratio. Instead, the higher fin pitch presents a lower maximum value but a higher mean value. This is due to a stronger accumulation of frost that is related to a lower mean evaporation temperature. Therefore, these two opposing effects justify the position of the minimum values in total cost difference for the 9 mm fin pitch.

Analysing overall results, defrosting 20 % of the field is more cost-effective than defrosting 10 % in almost all situations. Defrosting 20 % of the field rather than 10 % results in a lower total number of defrosts at the same frequency (i.e. in a 40-hour frosting cycle, defrosting happens 10 times for 10 % fraction while 5 times for the 20 % fraction). However, defrosting a larger fraction of the field means a lower number of active evaporators serving the system, having a stronger, but less frequent, reduction of evaporation temperature. A higher number of defrosting means a lower mean amount of frost on the evaporator, but time-based strategies may run into unnecessary defrosting. The analysis of defrosting strategies is based on these aspects, and their combination leads to the result obtained that is strictly related to this case study.

The sensitivity analysis on electricity price variability is shown in Fig. 22 in terms of total cost percentage difference. Despite the optimal solution is the same for each price, increasing the cost of energy from 0.05 €/kWh to 1 €/kWh does not always lead to an increase in cost difference. This is related to the relative weight of investment costs, although they represent on average 3 % of total costs. In fact, in cases where the investment cost of the field is lower than that of the optimal solution, there is a maximum increase in the difference in total costs equal to 10 %. In the opposite cases, the difference in total costs decreases, up to −11 %.

8. Conclusions and further comments

This work has presented a comprehensive model to evaluate the seasonal performance of a large-size air-source heat pump served by an evaporator field. A detailed modeling approach for the parallel evaporators, accounting for frost formation and defrosting transients, has been described. This model was subsequently tested and calibrated through a dedicated experimental campaign conducted on a horizontal wavy fin-and-tube evaporator. According to that, these results refer only to the case study analysed and not generalizable to other climatic conditions or system characteristics such as size or refrigerant used. The validated model was then employed to evaluate various evaporator field configurations, including variations in the number of evaporators, airflow velocity, and fin pitch, for a specific Danish district heating case study, aiming to identify an optimal solution in terms of total costs. A preliminary analysis, investigating a broad design domain, successfully identified a region of thermo-economic convenience for further seasonal simulations.

The primary outcome of this study highlights a critical discrepancy: a preliminary design analysis that neglects frost formation can suggest a domain of nearly equivalent thermo-economic solutions (showing only up to 5 % deviation in total costs compared to the best solution). However, when frost formation and the subsequent necessity of defrosting are considered, there is a substantial difference among solutions within this equivalent domain (up to 70 % variation), leading to a significant reduction of the convenient solution space.

Among the investigated design variables, the optimal evaporator field size was found to be 80 evaporators. Increasing the airflow velocity up to 3.5 m/s leads to a reduced compressor consumption without exceeding the additional fan energy penalty, while an intermediate fin pitch of 9 mm was most effective in balancing airflow obstruction and heat transfer. Regarding the defrosting strategy, the simultaneous defrosting of 20 % of the evaporator field consistently was more convenient than defrosting 10 %.

A key contribution of this work is the assessment of quantitative seasonal results that directly support design and optimization. The outcomes demonstrate how apparently minor design choices, such as airflow velocity, fin pitch, or defrosting fraction, lead to significant differences in seasonal efficiency and total costs, up to 70 %. These insights provide practical guidance for engineers and decision-makers, enabling the development of more robust and cost-effective heat pump systems under real operating conditions.

This work has therefore demonstrated the importance of including frosting and defrosting cycles into the design and optimization process for evaporator fields serving large-scale heat pumps. The model has proven effective in reproducing the performance of a single evaporator under controlled laboratory conditions, providing a reliable foundation for field-scale simulations. We acknowledge that the accuracy of the predictions reflects the assumptions adopted, and further comparisons with experimental data could help to quantify their influence on deviations between simulations and real operation. Given the complexity and duration of such experimental investigations, the results presented here already offer a solid basis to guide design choices and highlight the thermo-economic relevance of accounting for frosting and defrosting phenomena in large-scale heat pump applications. Clearly, the comparison of numerical predictions with experimental results on a real system could further assess the effect of the assumptions on simulation deviations and allow the improvement of the model predictive capability.

Credit authorship contribution statement

A.F. Passarelli: Writing – original draft, Investigation, Data curation. U. Merlo: Validation. F. Pelella: Formal analysis, Data curation. L. Viscito: Writing – original draft, Formal analysis, Data curation. S. Filippini: Supervision. A.W. Mauro: Writing – review & editing, Supervision, Conceptualization.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements

Acknowledgment is made to the Ministry of University and Lu-Ve S.p.A. for co-funding the DM-352 PhD scholarship in Industrial Engineering at the University of Naples Federico II and for providing the experimental facilities.

Appendix A. Evaporator model closure equations

A1. Finite volume equations characterization

The heat transfer problem shown in Fig. 5 can be solved with a system of equations involving conductive heat fluxes among the tube, fin and frost contact surfaces and convective heat fluxes with the refrigerant inside tubes and the external airflow. These heat rates

{\dot{Q}}_{s}

are described with the following general equation:(A1.1)

{\dot{Q}}_{s} = C_{s} A_{s} Δ T_{s}

Where the conductance

C_{s}

are expressed using the expression of the series of conductance since the domain is composed by different materials:(A1.2)

C_{s} = \frac{k_{1} k_{2}}{k_{1} s_{2} + k_{2} s_{1}}

In case of conductance along the x axis, the thicknesses involved in Eq. (A1.2) are equal to the half of the distance between the centres of the adjacent volumes considered, while for the r axis the thicknesses refer to the logarithm of the radii ratio. Consistently with the radial symmetry of the domain, the edge surfaces

A_{s}

are circular crowns when parallel to the r axis or cylindrical when parallel to the x axis.

A2. Convective heat transfer coefficients and pressure drop correlations

The refrigerant convective heat transfer coefficient is evaluated using the correlation of Gungor-Winterton [44], while proprietary correlations provided by the manufacturer are involved in the air-side heat transfer coefficient calculation. The general expression of the correlation is the following:(A2.1)

h_{air} = a \cdot w^{b}

Where appropriate coefficients

a

and b are used depending on the geometrical features of the fin and tube coil. Alternatively, the model can be used with Wang et al. correlations [45,46].

Regarding pressure drop on the refrigerant side the frictional contribution is evaluated with Muller-Steinhagen & Heck correlation [47], the momentum contribution is evaluated according to [48] and the gravitational one is equal to:(A2.2)

{Δ p}_{g_{i, j, k}} = ρ_{ref} g P_{fin}

To take into account tube bends, a correction factor is applied to the frictional contribution according to [49]. Instead, total air side pressure drops are evaluated with proprietary correlation provided by the manufacturer with the general form expressed in Eq. (A2.2) whose coefficients c, d and e depends on the evaporator geometrical features.(A2.3)

Δ p_{air} = c \cdot w^{d} \cdot N_{r}^{e}

Alternatively, Kays and London correlation [50] can be used with the friction factor evaluated with Wang et al. correlations [51], calculating the geometrical factors of the evaporator element considering frost thickness.

A3. Mass transfer coefficient

The Lewis analogy allows to evaluate the mass transfer coefficient as follows:(A3.1)

K_{m} = \frac{h_{air}}{(c_{p, a i r} + ω_{air} c_{p . v}) L e^{\frac{2}{3}}}

Since the molecular diffusion is weaker within a porous solid than a gas due to the impediment of pores, the Lewis number

Le

taken into account in frost formation models is expressed as the ratio of air thermal diffusivity and effective vapor diffusivity

D_{eff}

:(A3.2)

Le = \frac{α}{D_{eff}}

Starting from the correlation for vapor diffusivity created from the data presented in the work by Lee et al. [35], the effective diffusivity is a function of frost porosity

γ

and tortuosity

τ

according to Brian et al. [52]:(A3.3)

D_{eff} = \frac{γ}{τ} D_{v} = \frac{γ}{τ} (0.0013 T_{air} + 0.2219) \cdot 10^{- 4}

With:(A3.4)

γ = \frac{(ρ_{f} - ρ_{ice})}{ρ_{v} - ρ_{ice}}

While tortuosity is evaluated using the equation proposed by Zehnder et al. [53]:(A3.5)

τ = \frac{γ}{1 - 1 - \sqrt{1 - γ}}

Finally, the air thermal diffusivity

α

is equal to:(A3.6)

α = \frac{λ_{air}}{ρ_{air} c_{p, a i r}}

Hermes et al. correlation [9] is used for frost thermal conductivity that is expressed as a function of only frost density:(A3.7)

k_{f} = 0.132 + 3.13 \cdot 10^{- 4} ρ_{f} + 1.6 \cdot 10^{- 7} ρ_{f}^{2}

Appendix B. Experimental campaign

B1. Test facility description and experimental procedure

The experimental campaign has been conducted in the climatic chamber of Lu-Ve S.p.A. whose schematization is showed in Fig. B1 with all the measurement along the three main parts: liquid and suction line and the branch for the hot gas defrosting process. The testing evaporator is connected to a direct expansion system working with R507A, so a superheating control is set at the evaporator outlet. The experimental conditions of temperature and relative humidity in the chamber are controlled by a PID system connected with the chamber measurement and the set point values given at the beginning of each test. The sequence of frosting and defrosting phases within each test has been controlled with a timer, fixing the starting time of defrosting after frosting phases lasting from 3 to 6 h. The defrosting was carried out activating the hot gas line in Fig. B1.1 when the ambient temperature was below 0 °C, instead a forced convection defrosting was activated in the opposite case. In both cases, defrosting lasted about 25 to 30 min. After a single frosting/defrosting cycle, a mass measurement of the drainage water collection was carried out to evaluate the total mass of frost accumulated on the evaporator surface.

The tested evaporator is characterized by a staggered configuration with a 55 mm tube pitch and a 27.5 mm row pitch. Each row has 18 internal micro-finned copper tubes with an internal diameter of 12.61 mm and a length of 1215 mm divided in 9 circuits. The fin typology is wavy with a 5 mm pitch and a 0.3 mm thickness. The evaporator is positioned horizontally with a suction fan on its upper side as shown in Fig. B1.2(a), having a counter-current flow displacement as schematized in Fig. B1.2(b) showing also the ten thermocouples installed (six on the suction side and four on the fan outlet), the differential pressure transducer and an illuminated plexiglass window realized on the evaporator to install a camera that frames part of the first row on the refrigerant side. Camera frames have been used for the estimation of frost thickness. Considering the fin with the highest degree of orthogonality in the frame, frost thickness refers to the mean value along the fin. This value has been collected using a RGB pixel analysis to define the position of frost surface.

B2. Data reduction and uncertainty analysis

Tab. B2.1 shows the summary of all the measurements with their respective uncertainties. The thermal power of the testing evaporator was evaluated on the refrigerant side as follows:(B2.1)

{\dot{Q}}_{ref} = {\dot{m}}_{ref} (i_{s} - i_{l})

The superheated vapor enthalpy

i_{s}

is evaluated using the measured temperature

T_{rs}

at the measured pressure

p_{ev}

at the evaporator outlet. The inlet enthalpy

i_{l}

is evaluated, assuming an iso-enthalpic lamination, through the subcooled temperature at condenser outlet

T_{liq}

and using the compressor outlet pressure measurement

p_{co}

. Therefore, according to the expanded uncertainty calculation [54], the mean relative uncertainty for the thermal power is evaluated using the following equation:(B2.2)

u_{{\dot{Q}}_{ev}} = K \sqrt{{\dot{m}}_{ref}^{2} u_{(i_{s} - i_{l})}^{2} + Δ i^{2} u_{{\dot{m}}_{ref}}^{2}}

For the whole experimental campaign, assuming a coverage factor K equal to 2, this expanded uncertainty is equal to

\pm

10 %. The experimental evaporator cooling capacity has been used for the calibration of the air side heat transfer coefficient. Therefore, its expanded uncertainties can be applied to this calibration process, resulting in a calibration uncertainty that has been used in seasonal simulations. A mean

\pm

4.5 % variation for the seasonal energy consumption has been obtained as quantification impact of experimental uncertainties on seasonal simulation results.

Table B1. Measurement uncertainties.

Measurement type	Uncertainty
Inlet air-side thermocouples	$\pm$ 0.2 °C
Outlet air-side thermocouples	$\pm$ 0.3 °C
Evaporator outlet superheating, $T_{rs}$	$\pm$ 0.2 °C
Condenser outlet thermocouple, $T_{liq}$	$\pm$ 0.2 °C
Evaporator outlet pressure, $p_{ev}$	$\pm$ 0.01 bar
Compressor outlet pressure, $p_{co}$	$\pm$ 0.01 bar
Air side pressure drop
Refrigerant mass flow rate, ${\dot{m}}_{ref}$	$\pm 0.1$ % of read value

Fig. B2.1 shows the subsequent frosting cycles for the experimental room conditions of −5°C and 75 % for relative humidity. Taking from this set of cycles the one with the lower uncertainty, the values of thermal power and pressure drop related to the start, the half and the end of the selected cycle, calculated with a moving average on 10 values, are extracted to run a comparison with the model.

B3. Calibration of the air side heat transfer coefficient and pressure drop with frost formation

The mentioned correlations for the air-side heat transfer coefficient (Eq. A2.1) and pressure drop (Eq. A2.2) have been developed for evaporators working in free surface conditions. Therefore, a calibration of these correlations in case of frost formation has been conducted using experimental data. The heat transfer coefficient correlation has been modified as follows:(B3.1)

h_{air} = 0.86 \cdot a \cdot w^{b}

While a factor that takes into account the experimental blockage ratio has been added to the pressure drop correlation to consider the effect of frost formation:(B3.2)

Δ p_{air} = a \cdot v^{b} \cdot N_{r}^{c} \cdot {BR}^{2.5}

Appendix C. Seasonal simulations results

Table C1 presents the results for all the set of seasonal simulations ordered from the best solution to the worst one in terms of percentage difference in total costs with respect to the optimal solution.

Table C1. Seasonal simulation results.

N° simulation	Field size	Fin pitch [mm]	Airflow velocity [m/s]	Defrosting fraction of the field	Defrosting start time [h]	ΔC_tot [%]
12	80	9	3.5	20 %	40	0.00 %
10	80	9	3.5	20 %	20	1.51 %
8	80	9	2.5	20 %	40	1.55 %
6	80	9	2.5	20 %	20	1.68 %
40	120	6	1.5	20 %	40	2.03 %
18	80	12	2.5	20 %	20	3.45 %
28	120	9	2.5	20 %	40	3.48 %
20	80	12	2.5	20 %	40	3.51 %
46	120	15	2.5	20 %	20	3.66 %
16	80	6	2.5	20 %	40	4.17 %
52	120	12	1.5	20 %	40	4.76 %
48	120	15	2.5	20 %	40	4.90 %
50	120	12	1.5	20 %	20	5.47 %
11	80	9	3.5	10 %	40	5.49 %
2	80	9	1.5	20 %	20	6.23 %
7	80	9	2.5	10 %	40	7.88 %
9	80	9	3.5	10 %	20	8.05 %
5	80	9	2.5	10 %	20	8.79 %
19	80	9	2.5	10 %	40	9.71 %
26	120	9	2.5	20 %	20	10.07 %
17	80	12	2.5	10 %	20	10.11 %
34	40	9	4.5	20 %	20	10.18 %
14	80	6	2.5	20 %	20	10.20 %
4	80	9	1.5	20 %	40	10.47 %
47	120	15	2.5	10 %	40	10.56 %
15	80	6	2.5	10 %	40	10.58 %
39	120	6	1.5	10 %	40	10.66 %
27	120	9	2.5	10 %	40	11.35 %
38	120	6	1.5	20 %	20	11.67 %
45	120	15	2.5	10 %	20	11.81 %
51	120	12	1.5	10 %	40	12.05 %
42	60	15	4.5	20 %	20	13.15 %
1	80	9	1.5	10 %	20	14.68 %
33	40	9	4.5	10 %	20	16.05 %
49	120	12	1.5	10 %	20	16.07 %
44	60	15	4.5	20 %	40	16.25 %
3	80	9	1.5	10 %	40	17.98 %
22	40	9	2.5	20 %	20	18.98 %
23	40	9	2.5	10 %	40	18.98 %
41	60	15	4.5	10 %	20	19.68 %
43	60	15	4.5	10 %	40	21.30 %
13	80	6	2.5	10 %	20	24.88 %
36	40	9	4.5	20 %	40	25.05 %
21	40	9	2.5	10 %	20	26.31 %
25	120	9	2.5	10 %	20	27.17 %
35	40	9	4.5	10 %	40	30.42 %
37	120	6	1.5	10 %	20	30.63 %
30	40	3	2.5	20 %	20	31.31 %
31	40	3	2.5	10 %	40	31.31 %
24	40	9	2.5	20 %	40	44.04 %
29	40	3	2.5	10 %	20	52.53 %
32	40	3	2.5	20 %	40	70.83 %

Data availability

Data will be made available on request.

References

[1]
Y. Gong, G. Ma, Y. Jiang, L. Wang
Research progress on the fifth-generation district heating system based on heat pump technology
J. Build. Eng., 71 (2023)
Google Scholar
[2]
IEA 2021, Net Zero by 2050, A roadmap for the Global Energy Sector.
Google Scholar
[3]
P.H. Jorgensen, T. Ommen, B. Elmegaard
Quantification and comparison of COP improvement approaches for large-scale ammonia heat pump systems
Int. J. Refrig, 129 (2021), pp. 301-316
Google Scholar
[4]
S.O.K. Hansen, J.K. Jensen, M. Dominato, L.O. Reinholdt, W.B. Markussen, Results from measurement of heat pump for district heating using ambient air as heat source, 2022 15th IIR-Gustav Lorentzen conference on Natural Refrigerants, Trondheim Norway.
Google Scholar
[5]
J.J. Aguilera, W. Meesenburg, T. Ommen, W.B. Markussen, J.L. Poulsen, B. Zuhlsdorf, B. Elmegaard
A review of common faults in large-scale heat pumps
Renew. Sustain. Energy Rev. (2022)
Google Scholar
[6]
M. Song, C. Dang
Review on the measurement and calculation of frost characteristics
Int. J. Heat Mass Transf., 124 (2018), pp. 586-614
View PDF View article Google Scholar
[7]
F. Zhaosong, L. Lafang, Y. Cheng, H. Wang, M.O. Oladokun
Investigation on thermal-hydraulic performance of outdoor heat exchanger in air source heat pump
Exp. Therm Fluid Sci., 84 (2017), pp. 28-38
Google Scholar
[8]
L. Shangwen, S. Mengjie, P. Libor, J. Shen
A numerical study on frosting and its early stage under forced convection conditions with surface and environmental factors considered
Sustain. Energy Technol. Assess., 45 (2021), Article 101202
Google Scholar
[9]
C.J.L. Hermes, R.O. Piucco, J.R. Barbosa Jr, C. Melo
A study of frost growth and densification on flat surfaces
Exp. Therm Fluid Sci., 33 (2) (2009), pp. 371-379
View PDF View article Google Scholar
[10]
A.R. Tahavvor, M. Yaghoubi
Experimental and numerical study of frost formation by natural convection over a cold horizontal circular cylinder
Int. J. Refrig, 33 (7) (2010), pp. 1444-1458
View PDF View article Google Scholar
[11]
C.J.L. Hermes, V.S. Nascimento Jr, F.R. Loyola, R.P. Cardoso, A.D. Sommers
A study of frost build-up on hydrophilic and hydrophobic surfaces under forced convection conditions
Exp. Therm Fluid Sci., 100 (2019), pp. 76-88
View PDF View article Google Scholar
[12]
X. Wu, S. Hu, F. Chu
Experimental study of frost formation on cold surfaces with various fin layouts
Appl. Therm. Eng., 95 (2016), pp. 95-105
View PDF View article Google Scholar
[13]
H.W. Schneider
Equation of the growth rate of frost forming on cooled surfaces
Int. J. Heat Mass Transf., 21 (1977), pp. 1019-1024
Google Scholar
[14]
A. Saygin, A.M. Basol, M. Arik
An Euleriam multiphase frost model based on heat transfer measurements
Int. J. Heat Mass Transf., 216 (2023)
Google Scholar
[15]
S.K. Padhmanabhan, D.E. Fisher, L. Cremaschi, E. Moallem
Modeling non-uniform frost growth on a fin-and-tube heat exchanger
Int. J. Refrig, 34 (8) (2011), pp. 2018-2030
View PDF View article Google Scholar
[16]
R. Mastrullo, A.W. Mauro, L. Menna, A. Palma, G.P. Vanoli
Transient model of a vertical freezer with door openings and defrost effects
Appl. Energy, 121 (2014), pp. 38-50
View PDF View article Google Scholar
[17]
H.-Y. Ye, K.-S. Lee
Performance prediction of a fin-and-tube heat exchanger considering air-flow reduction due to the frost accumulation
Int. J. Heat Mass Transf., 67 (2013), pp. 225-233
View PDF View article Google Scholar
[18]
ASHRAE Handbook, Refrigeration, 2018.
Google Scholar
[19]
J.A. Dopazo, J. Fernandes-Seara, F.J. Uhia, R. Diz
Modelling and experimental validation of the hot-gas defrost process of an air-cooled evaporator
Int. J. Refrig, 33 (2010), pp. 829-839
View PDF View article Google Scholar
[20]
N. Hoffenbecker, S.A. Klein, D.T. Reindl
Hot gas defrost model development and validation
Int. J. Refrig, 28 (2005), pp. 605-615
View PDF View article Google Scholar
[21]
S. Mengjie, D. Shiming, D. Chaobin, M. Ning, W. Zhihua
Review on improvement for air source heat pump units during frosting and defrosting
Appl. Energy, 211 (2018), pp. 1150-1170
Google Scholar
[22]
X. Li, G. Ma, T. Lu, L. Gao, W. Rong, Y. Gong, S. Xu
Experimental study on defrosting multi outdoor units in turn for air source heat pump using hot gas
Int. J. Refrig, 174 (2025), pp. 76-85
View PDF View article Google Scholar
[23]
C. Liang, X. Li, X. Meng, W. Shi, J. Gu, B. Wang, Y. Lv
A year-round efficient air source heat pump with separate and distant heat extraction for non-decaying heating capacity during defrosting
Energy, 315 (2025)
Google Scholar
[24]
S. Liang, H. Wang, X. Bai, G. Wang, X. Gao, X. Tian
A field study on the impacts of cold island effect on operating performances of air source heat pump array for space heating
J. Build. Eng., 82 (2024), Article 108266
View PDF View article Google Scholar
[25]
W. Wei, B. Wang, H. Gu, L. Ni, Y. Yao
Investigation on the regulating methods of air source heat pump system
Energy Build., 235 (2021)
Google Scholar
[26]
L. Cabrol, P. Rowley
Towards low carbon homes – a simulation analysis of building-integrated air-source heat pump systems
Energ. Buildings, 48 (2012), pp. 127-136
View PDF View article Google Scholar
[27]
H. Pieper, T. Ommen, J.K. Jensen, B. Elmegaard, W.B. Markussen
Comparison of COP estimation methods for large-scale heat pumps used in energy planning
Energy, 205 (2020)
Google Scholar
[28]
J.L. Threlkeld
Thermal Environmental Engineering, Prentice-Hall Inc, New York NY (1970)
Google Scholar
[29]
M. Noussan, M. Jarre, A. Poggio
Real operation data analysis on district heating load pattern
Energy, 129 (2017), pp. 70-78
View PDF View article Google Scholar
[30]
M.H. Kim, K.S. Lee
Determination method of defrosting start-time based on temperature measurements
Appl. Energy, 146 (2015), pp. 263-269
View PDF View article Google Scholar
[31]
B. Rogie, J.K. Jensen, S.O.K. Hansen, W.B. Markussen
Analysis of cold air recirculation in the evaporators of large-scale air-source heat pumps using CFD simulations
Fluids, 5 (4) (2020), p. 186
Google Scholar
[32]
S.M. Li, C.C. Wang
Predictive models on the frost formation for plain surface – a review and comparative study
Int. Commun. Heat Mass Transfer, 129 (2021)
Google Scholar
[33]
Y.X. Tao, R.W. Besant, K.S. Rezkallah
A mathematical model for predicting the densification and growth of frost on a flat plate
International Journal of Heat and Mass Transfer, 36 (2) (1993), pp. 353-363
View PDF View article Google Scholar
[34]
Y. Hayashi, A. Aoki, S. Adachi, K. Hori
Study of frost properties correlating with frost formation types
J. Heat Transfer, 99 (1977), pp. 239-244
Google Scholar
[35]
K.-S. Lee, W.-S. Kim, T.-H. Lee
A one-dimensional model for frost formation on a cold flat surface
Int. J. Heat Mass Transf., 40 (18) (1997), pp. 4359-4365
View PDF View article Google Scholar
[36]
Y.B. Lee, S.T. Ro
Analysis of the frost growth on a flat plate by simple models of saturation and supersaturation
Exp. Therm Fluid Sci., 29 (2005), pp. 685-696
View PDF View article Google Scholar
[37]
L. Zhang, M. Song, N. Mao, J. Dong
Temporal and spatial frost growth prediction of a tube-finned heat exchanger considering frost distribution characteristics
Int. J. Heat Mass Transf., 183 (2022)
Google Scholar
[38]
F. Breque, M. Nemer
Modeling of a fan-supplied flat-tube heat exchanger exposed to non-uniform frost growth
Int. J. Heat Mass Transf., 75 (2017), pp. 129-140
View PDF View article Google Scholar
[39]
D.L. da Silva, C.J.L. Hermes, C. Melo
First-principles modeling of frost accumulation on fan-supplied tube-fin evaporators
Appl. Therm. Eng., 31 (2011), pp. 2616-2621
View PDF View article Google Scholar
[40]
M. Song, X. Xu, N. Mao, S. Deng, Y. Xu
Energy transfer procession in an air source heat pump unit during defrosting
Appl. Energy, 204 (2017), pp. 679-689
View PDF View article Google Scholar
[41]
M.M. Shah
A general correlation for heat transfer during film condensation inside pipes
Int. J. Heat Mass Transf., 22 (4) (1979), pp. 547-556
View PDF View article Google Scholar
[42]
F.W. Dittus, L.M.K. Boelter
Univ California Publ Eng., 2 (1930), pp. 443-461
Google Scholar
[43]
Store Varmepumper i fjernvarmen driftserfaringer, 2019 Dansk Fjernvarme. https://danskfjernvarme.dk/media/iuujntqk/rapport-om-driftserfaringer-med-store-varmepumper-i-fjernvarmen-master_final.pdf.
Google Scholar
[44]
K.E. Gungor, R.H.S. Winterton
A general correlation for flow boiling in tubes and annuli
Int. J. Heat Mass Transf., 29 (3) (1986), pp. 351-358
View PDF View article Google Scholar
[45]
C.C. Wang, W.L. Fu, C.T. Chang
Heat transfer and friction characteristics of typical wavy fin-and-tube heat exchangers
Exp. Therm Fluid Sci., 14 (2) (1997), pp. 174-186
View PDF View article Google Scholar
[46]
C.C. Wang, C.J. Lee, C.T. Chang, S.P. Lin
Heat transfer and friction correlation for compact louvered fin-and-tube heat exchangers
Int. J. Heat Mass Transf., 42 (1992), pp. 1945-1956
Google Scholar
[47]
H. Muller-Steinhagen, K. Heck
A simple friction pressure drop correlation for two-phase flow in pipes
Chem. Eng. Process., 20 (6) (1986), pp. 297-308
View PDF View article Google Scholar
[48]
A.W. Mauro, J.M. Quiben, R. Mastrullo, J.R. Thome
Comparison of experimental pressure drop data for two phase flows to prediction methods using a general model
Int. J. Refrig (2007), pp. 1358-1367
View PDF View article Google Scholar
[49]
M. Tammaro, A.W. Mauro, J. Bonjour, R. Mastrullo, R. Revellin
Curvature ratio effect on two-phase pressure drops in horizontal return bends: Experimental data for R134a
Heat Transfer Eng., 34 (13) (2013), pp. 1124-1132
Google Scholar
[50]
W.M. Kays, A.L. London
Compact Heat Exchangers (3rd ed), McGraw-Hill, New York (1984)
Google Scholar
[51]
Y.M. Ko, J.Y. Song, J.W. Lee, S. Sohn, C.H. Song, M. Khoshvaght-Aliabadi, Y. Kim, Y.T. Kang
A critical review on Colburn j-factor and f-factor and energy performance analysis for finned tube heat exchangers
Energy, 287 (2024), Article 129609
View PDF View article View in Scopus Google Scholar
[52]
P.L.T. Brian, R.C. Reid, Y.T. Shah
Frost Deposition on Cold Surfaces
Industrial & Engineering Chemistry, 9 (1970), pp. 375-380
View at publisherCrossref View in Scopus Google Scholar
[53]
B. Na, R.L. Webb
Mass transfer on and within a frost layer
Int. J. Heat Mass Transf., 47 (2004), pp. 899-911
View PDF View article View in Scopus Google Scholar
[54]
R.J. Moffat
Describing the uncertainties in experimental results
Exp. Therm Fluid Sci., 1 (1988), pp. 3-17
View PDF View article Google Scholar

In this section

Design optimization of NH3 evaporators for heat pumps in district heating

Air-source ammonia heat pump for district heating: a field modeling approach with focus on frosting-defrosting cycles

URL has been copied successfully!