Dynamic simulation of indirect air conditioning systems with optimized computational time

The simulation of HVAC systems is a powerful tool to improve the energy efficiency in buildings. The modelling of such systems faces several obstacles due to both the physical phenomenology present and the numerical resolution difficulties. The present work is an attempt to develop a robust, fast, and accurate model for HVAC systems that can interact with the other relevant systems involved in buildings thermal management. The whole system model has been developed in the form of libraries under the Modelica language to exploit its advantageous characteristics: object-oriented programming, equationbased modelling, and handling of multi-physics. The global resolution is carried out dynamically so that not only steady-state predictions can be conducted but also control strategies can be studied over meaningful periods of time. This latter aspect is crucial for optimizing energy savings. The libraries include models for all the system individual components such as pumps, compressors or heat exchangers (operating with twophase flows and/or moist air) and also models assemblies to account for vapour compression units and liquid circuits. An illustrative example of an indirect air conditioning system is detailed in the present work in order to highlight the model potential.


Introduction
According to the International Energy Agency, the buildings sector is one of the world main energy consumers together with the other two key sectors, namely, industry and transportation [1]. In general, Heating Ventilation and Air Conditioning (HVAC) accounts for a significant proportion of the buildings total energy usage. Both the running cost and the indirect negative impact on the environment produced by HVAC units can be reduced by enhancing the systems efficiencies. On the contrary, inefficient operation or poor design could lead to undesirable waste of energy [2]. Furthermore, the need for more sustainable systems becomes evident due to the rising demand for maintaining comfortable indoor air conditions (e.g. space cooling [3]).
The global efficiency of HVAC systems can be maximized through the optimization of its operating parameters and/or the implementation of appropriate control procedures. In this regard, an accurate simulation model for HVAC systems is the perfect tool to conduct studies on both optimization and control strategies to achieve better cost-effective solutions. However, the simulation of such systems is a challenging assignment due to their structure complexity (interaction between different fluid circuits), to the amount of components included (different types of heat exchangers, turbomachinery, ducting elements, sensors, controllers), and to the physical phenomena involved (heat and mass transfer, fluid phase changes).
Many models to predict the operating behaviour of HVAC systems have been proposed in the technical literature. The modelling techniques can be classified into three different methods, namely, data driven, physics-based and grey-box [4]. The goal of this work is to show the progress done so far of an on-going project to develop a reliable tool to simulate HVAC systems. The model is being developed using the Modelica language and consists of several libraries that account for all components and fluids present in such systems.
Modelica is an equation-based language (models are based on equations where assignments are not needed unlike imperative programming languages). The whole system of equations is systematically rearranged and algebraically manipulated to generate an optimized code prior to simulation so that Modelica users or programmers do not have to deal with this task. Modelica is also an object-oriented language where models can be associated following real physical connections and where multi-physics simulations can be handled with ease. The use of such approach for building energy modelling, simulation and optimization has been recommended due to its important advantages over approaches based on imperative languages [5]: • Better models extensibility to new use cases, • More efficient simulations in different use cases, • Easier implementation of control sequences.
The development of the aforementioned libraries for HVAC systems is currently focused on modelling the vapour compression system (VCS) and the auxiliary liquid circuits. Three main attributes are required for individual components and components assemblies: robustness (numerical reliability within a wide range of operating conditions), accuracy (numerical predictions close to measured values), and minimized computational time. This latter attribute is crucial to guarantee reasonable computational time in HVAC systems with large and complex configurations but also to conduct control procedures over large periods of time.
Several difficulties were faced during the library development, namely, the initialization of closed circuits, the management of the refrigerant charge, the handling of near-zero mass flow rates, and the implementation of low time consuming heat exchangers with two-phase flows and moist air. All these aspects were carefully tackled to meet the numerical sturdiness and computational time requirements.
This document presents the simulation of an air conditioning system made up of a vapour compression cycle linked to a liquid loop that cools down a room through an air-to-liquid heat exchanger. In the following sections the mathematical formulation of components and systems is briefly described and the most relevant numerical results of the studied system are shown.

System configuration
The system studied herein consists of two different flow circuits interacting between each other (see Figure 1). At the system core we find the vapour compression cycle which is made up of several major elements, namely, a reciprocating compressor (to propel the fluid flow), a refrigerant-to-air condenser (to eject heat towards the environment), a refrigerant-to-liquid evaporator (to extract heat from the liquid loop), and some additional minor elements.
The liquid circuit interacts with the vapour compression loop through the refrigerant-to-liquid heat exchanger which is a shared component. It also contains a pump (to drive the liquid flow), a liquid-to-air heat exchanger (to cool down the indoor environment), and other additional components like valves, tubes and tanks.

Components
This section is devoted to introduce the models of components that have been developed and included in the Modelica libraries created. The primary goal was to provide low time consuming components to prevent from having bottlenecks during the numerical resolution. Therefore, the modelling of most components has been based on simplified mathematical formulations and fed with performance maps. The mathematical formulation of the most determinant components is presented in this section.

Turbo machines
The vapour compression system and the liquid loop are driven by a compressor and a pump, respectively.

Compressor
The compressor used in the vapour compression system is of the hermetical reciprocating type. The compression process is considered instantaneous as its time scale is significantly smaller compared to the whole vapour cycle thermal and fluid-dynamic time scale. Both the mass flow rate produced and the compressor work done are expressed via the volumetric and the isentropic efficiencies: These efficiencies are also defined as a function of the compression ratio from experimental data or from high detailed simulations. In addition, a dynamic calculation for the fluid contained inside the compressor shell is carried out:

Pump
The pump model calculates the pressure gain together with the consumed power by means of performance maps. These two variables are expressed as a function of the volumetric flow and the working frequency: The maps are previously generated from empirical data, or alternatively, from more advanced models.

Heat exchangers
Heat exchangers are probably the most challenging components of the whole system. Three types of heat exchangers were implemented, namely, a condenser (refrigerant vs. dry air), an evaporator (refrigerant vs. liquid), and the air-to-liquid heat exchanger (liquid vs. moist air). These components share a common internal structure where two fluid flow models are thermally linked through a specific component that stands for the solid interface between them (see Figure 2). This simple approach allows to easily build a large variety of heat exchangers by combining different fluid flow models. The most relevant aspects of these crucial components are detailed.

Fig. 2. Heat exchanger internal scheme
Fluid flows are calculated from lumped heat transfer analyses in order to optimize the computational time, which is a critical requirement for the simulation of large and complex systems. The model accuracy relies on the empirical information provided (e.g. heat transfer coefficients) due to its lumped nature. The acquisition of appropriate heat transfer coefficients is essential for the model accuracy.

Liquid flow model
The prediction of the heat transfer between the liquid flow and the solid interface is carried out considering a single control volume. The method implemented is based on an ε-NTU approach in order to optimize the calculation speed and also to prevent the calculation of unreal temperature values. The sensible heat is calculated by:

Refrigerant flow model
The refrigerant model implemented herein is based on the switching moving boundary approach SMB [6] as it has low computational time (only three control volumes are used) unlike the traditional fixed volumes approach (where a large number of control volumes are needed).
However, SMB models possess added numerical complexity not only because a mean void fraction has to be calculated but also because volumes have variable lengths and may be activated or deactivated. The successful implementation of the SMB model was the outcome of an extensive programming and verification effort. The flow domain scheme for an evaporator is depicted in Figure 3 where three different zones are distinguished: sub-cooled, two-phase, and super heating. The conservation equations of mass and energy are calculated at each zone. Depending on the current flow conditions, the model switches between different operating modes so that zones are activated or deactivated as shown in Figure 4. The switching between modes (activation or deactivation of zones) leads to major modifications of conservative equations at each individual control volume. During transitions some local boundary values are rearranged and some equations are fully replaced. This procedure is necessary to retain the equations structure solved by Modelica. The switching criteria between modes has been carefully addressed. A zone is deactivated whenever its length becomes smaller than a reference minimum length value, and on the contrary, a zone is activated whenever its local outlet enthalpy crosses any boundary saturation enthalpy value.
The SMB numerical stability has been extensively tested at all possible transitions and operating modes. For instance, Figure 5 depicts the evolution of the zone lengths for a particular robustness test conducted on an evaporator. In this case the inlet specific enthalpies of both fluids were modified independently over time in order to force the heat exchanger to switch between modes 1, 2, 3 and 4 as the sub-cooled and the superheated zones were activated or deactivated (the CPU time to run this particular case was about 1 second).
The pressure drop throughout the domain is calculated following an approach based on the assumption that the variation is linear.

Moist air flow model
The dry air flow calculation is conducted in a similar way to the liquid case presented in section 3.2.1. However, when moist air is considered the heat transfer due to water condensation should be also considered. The total heat transferred by the fluid is now calculated from both its sensible and latent terms: The sensible heat term is calculated as in section 3.2.1, while the humidity at the exit is expressed as follows:

Empirical information
The heat exchangers level of accuracy depends on the calibration process done prior to simulations. Both the thermal performance and the pressure drop calculation are based on performance maps. These maps are previously obtained from a more advanced model, or alternatively, from experimental tests.

Minor components
In addition to the most important components other devices such as valves, ducts, tanks, sensors, and controllers are also included in the libraries.

Expansion devices
The valve models used to regulate the mass flow rate through fluid circuits are based on lumped parameters. The mass flow rate value is calculated as a function of the pressure difference. The main equation includes a flow coefficient but also an effective cross-sectional area which is adjusted from an external control signal.

Tanks
Reservoirs are typically found in vapour compression cycles (to store or release refrigerant as needed) and inside liquid loops (to allow fluid expansion). The volume pressure value is unique for the whole tank. The mass and energy conservation equations are expressed as follows: The particular refrigerant reservoir used in the vapour compression system handles two-phase flows and assumes complete separation of phases. In that case, the mean density and total enthalpy are calculated from the liquid volume ratio: The outlet enthalpy calculation will depend on the tank geometry and characteristics.

Ducting
Ducts are minor elements that account for the pressure losses occurring through the linking tubes located between the VCS or liquid loop components. The pressure drop calculation is based on empirical parameters. Heat transfer with the environment can also be considered in this component.

Systems
Several difficulties arise when complex systems are numerically solved. The complications increase when some specific characteristics are present such as closed loops or fluids experiencing phase changes. For instance, many aspects of the vapour compression cycle demand special attention, namely, the self-regulation thermostatic valve, the system initialization stability, the strategy to reach a steady-state operating point, the handling of start-up/shut-down events, the flow at low/null mass flow rate conditions, and the management of the system refrigerant charge. These critical aspects were tackled in the present work focusing on achieving robustness and minimum computational time.

Initialization
The system resolution at time equal to zero is critical for models with dynamics as no information on the previous time step is available. The initialization goal is to assign consistent initial values for all variables and derivatives but also to provide suitable initial constraints to the equations. In this sense, a global component dedicated to this particular purpose has been developed in order to automatically set the initial values of all variables. In addition, both homotopy transformations and dynamic relaxation of certain variables were implemented to facilitate the resolution of components during the initialization [7].

Steady-state resolution
A direct calculation of the initial steady-state cannot be achieved by setting all derivatives equal to zero because this leads to a singularity. This issue could be addressed by reaching the steady-state with an initial meaningless transient relaxation starting from a set of guessed initial values. The criteria for the steady-state convergence is defined by the user and can involve one or more variables from any component from any part of the system. A comprehensive set of tests to assess both the model initialization and the steady-state convergence robustness has been carried out.

Start-up/shut-down operations
The system simulation at (near) zero mass flow rates is a considerable challenge. At such conditions oscillations appear due to the fast dynamics in the model that force the solver time step to be significantly reduced. Several numerical aspects were considered in order to prevent the simulation from failing at such conditions. On the one side, equations were rearranged to avoid any possible division by zero, and on the other side, regularization was added to all pressure drop equations. Additional improvements were obtained by merging the heat transfer coefficient values of single-and two-phase flows when the mass flow rate approaches zero. The resolution of consecutive shut-downs and start-ups is essential for control purposes in some applications.

Refrigerant charge management
The refrigerant mass amount conservation is a critical aspect to consider for simulations. The total system refrigerant amount must be provided during the initialization in order to have a unique solution, otherwise the system will converge to different solutions depending on the initial guessed values. The refrigerant charge amount value can experience non-negligible variations during the simulation execution as long as the mass conservation undergoes small deviations in time due to the lumped nature of the two-phase flow heat exchanger models used here.
The approach used in the present model consists in correcting the refrigerant system charge whenever its value deviates too much from the reference value. A global component dedicated to manage the refrigerant charge has been developed. The idea was to calculate the total refrigerant mass during time steps and compare it to the reference value. When the defined deviation criterion is not fulfilled mass is artificially injected or extracted into or from the system, respectively.  Figure 6 shows the evolution of the refrigerant amount when the charge correction is applied and when it is not.
These results are obtained from the same simulation detailed in section 5.2 where the system has been subjected to several transient events.

Transient resolution
The system transient response throughout time can produce numerical issues when extreme modifications of the boundary conditions occur. The system model has been tested at these conditions to assure robustness. For instance, it has been subjected to multiple external inputs changes during test simulations: the compressor frequency, the by-pass valve opening degree, the heat exchangers secondary fluids conditions, and so on.

Thermostatic valve operation
The vapour compression cycle is controlled by means of a thermo-static valve located between the condenser and the evaporator. This type of self-controlling systems can generate several numerical instabilities. The proposed approach consisted in using the valve formulation described in section 3.3.1 combined with a PID controller with limited output and set point weighting (it compares the reference superheating value against the measured one and provides a suitable opening degree parameter for the valve). The PID parameters should be adjusted from experimental data.

Results
The air conditioning system detailed in sections 2, 3 and 4 has been assembled in the Dymola framework using the recently developed libraries. A comprehensive study was carried out to test the system numerical robustness and computational time at different operating conditions. The system graphical scheme is shown in Figure 5 where the two main fluid closed circuits are observed. Firstly, the VCS which is made up of two heat exchangers, an hermetical reciprocating compressor, two valves, connecting tubes and a reservoir. The VCS role is to extract heat from the evaporator and eject it to the environment -whose conditions are defined at the condenser air flow inlet-. The VCS performance can be regulated either by changing the by-pass valve opening degree or by modifying the compressor working speed. The evaporator, which is shared by both circuits, sucks heat from the liquid loop. The VCS is self-regulated as its operating point is controlled by the expansion valve whose level of aperture is determined by the super heating degree at the evaporator outlet. And secondly, the liquid circuit which is composed of two heat exchangers, one centrifugal pump, a tank, and a regulation valve. The liquid mass flow rate is regulated through both the valve level of contraction and the pump working speed. The role of the liquid-to-air heat exchanger is to cool down the room, which is in turn characterized with the inlet air flow boundary conditions (water condensation is considered).
The inputs needed for the model are the external conditions (environment and room return air temperatures) and the system operating conditions (compressor and pump speeds, valves opening level, superheating degree). If a particular study is required, the environment data should be obtained from a weather file while the room return air should be calculated from an air circuit.

Steady-state simulations
A full set of steady-state simulations has been conducted to analyse the influence of the AC system main external variables, namely, the environment air temperature, the room return air temperature, and the VCS compressor working frequency. These variables were modified over a wide range starting from a reference case defined in the following way: room return air temperature of 24 ºC, environment air temperature of 35 ºC, and a defined compressor working frequency.  Figure 6 shows the system response depending on the room return air temperature value. The variables shown include the room supply air temperature and three dimensionless values: the heat extracted, the VCS coefficient of performance, and the VCS mass flow rate.
It is seen that higher room return air temperatures entail higher supply air temperatures. The increase of the return air temperature forces the VCS evaporation temperature to rise and the VCS mass flow rate to increase to cope with the heat demand. Both the air-toliquid heat exchanger and the compressor power increase but in a different rate so that the COP is enhanced. The influence of increasing the environment temperature is shown in Figure 7. In this case the VCS high pressure must increase to enable heat rejection to the ambient. Therefore the compressor power consumption grows and the VCS coefficient of performance deteriorates. The overall cooling power diminishes and the room supply temperature rises. Figure 8 presents the system behaviour as the compressor frequency is increased. In this case, the power consumption and the VCS mass flow rate increase as well as the heat extracted in the evaporator. However, the latter parameter increases at a lower rate so that the coefficient of performance diminishes. To produce these studies more than 80 steady-state simulations were carried out. The computational time for each particular test ranged between 4 and 5 seconds. The time values mentioned in this document stand for the CPU time for integration provided by Dymola (version 17). The specifications of the PC used to run all cases presented herein are: Intel(R) Core (TM) i5-2540M CPU @ 2.60GHz with 8192MB RAM.

Transient simulations
In this section an illustrative test of the AC system response to several changes of the internal operating conditions is detailed with the aim of showing the flexibility of the developed model (the external boundary conditions are maintained constant for this test: environment temperature 35 ºC and room return air temperature 24 ºC).
The system transient behaviour has been studied for a period of 2500 seconds and has been subjected to multiple adjustments/events at different moments of the simulation. The system events sequence is listed as follows: • System is initialized (see section 4.1), • 400 seconds: system reaches steady-state, • 500 seconds: compressor frequency is reduced (75% of the initial working frequency), • 1000 seconds: the VCS by-pass valve is opened up to a certain degree, • 1500 seconds: the pump rotating speed is increased (50% more than the initial rotating speed), • 2000 seconds: VCS superheating value increased 5 ºC.
The complete simulation took about 16 seconds and the evolution of the VCS dimensionless mass flow rate, the heat of the air-to-liquid heat exchanger, and the room supply air temperature are plotted in Figures 9, 10 and 11, respectively.   From Figures 9, 10 and 11 it is observed how the system initializes from unreal values and reaches the steady-state condition after about 400 seconds. At 500 seconds, when the compressor working frequency is decreased, the system mass flow rate is reduced together with the heat extracted at the air-to-liquid heat exchanger while the supply air temperature grows accordingly. At 1000 seconds the opening of the by-pass valve leads to a recover on the mass flow rate. However, the cooling capacity of the VCS is degraded and the room supply air temperature is increased further. At 1500 seconds the pump operates at a higher speed and enhances the heat exchange in both the air-to-liquid heat exchanger and the evaporator resulting in a lower room supply air temperature. Finally, at 2000 seconds, when the VCS super-heating value is incremented, the heat transferred through the evaporator is reduced, leading to an even higher room supply air temperature.

Start-up/shut-down simulations
This section is devoted to test more challenging transient situations faced by the VCS: starts-ups and shut-downs. The simulation of such conditions pose significant numerical difficulties that must be carefully handled. The illustrative case presented herein serves to highlight the robustness attained by the model.
The simulation consists of a transient case lasting 4000 seconds and subjected to consecutive start-ups and shut-downs. The computational time used for this simulation was about 49 seconds which is relatively low considering the stress faced by the solver at the start-up and shut-down events.    After about 500 seconds the system reaches the steadystate and the compressor is immediately turned off. At that moment, the mass flow rate falls to zero while the heat flowing through the air-to-liquid heat exchanger decreases progressively towards zero. The room supply air temperature increases progressively towards the upper limit of 24 ºC which corresponds to the return air temperature. This sequence is repeated two more times until the end of the simulation.

Conclusions
The focus of the present work was to develop a robust and accurate model to simulate HVAC systems based on Modelica libraries with optimized time consumption. The model has been comprehensibly tested to ensure the aforementioned capacities.
The tool has shown to consume low computational time when calculating steady-state simulations. For instance, the time needed to predict the steady-state of the system configuration studied in the present work was about 5 seconds.
The tool has proven to be numerically robust at all levels. No issues were reported during the initialization of any case conducted. The tool has been satisfactorily tested at demanding transient conditions including startups and shut-downs. In addition, simulations at (near) zero mass flow rate were also carried out with success.
The libraries are now ready to be extended by adding more components (e.g. chambers, fans) and by building up more complex system configurations (e.g. additional and/or more complex fluid loops). The tool is also ready to study control procedures during short or long periods of time due to the robustness and quickness achieved. This work has been developed within the research project "MALET -Development of Modelica Libraries for ECS and Thermal Management Architectures" within the CleanSky2 framework.