Switching criteria analysis for a condenser moving boundary model

The switching moving boundary (SMB) scheme is widely used for numerical simulation of heat exchangers, especially in vapour compression refrigeration system. The method divides the heat exchanger into up to three zones (superheated vapour, two-phase fluid and subcooled liquid), depending on the operating conditions and treats each as a lumped-parameter system. During transients the number of zones may vary, and suitable switching criteria and related threshold values must be given, for this change to occur. Also, variables for the zones momentarily phased out need to be tracked, to allow smooth operation of the model. In this work, a switching moving boundary model of a brazed plate condenser in counter-flow arrangement is studied to obtain some useful guidelines for the selection of the optimal switching thresholds. It is shown that these values influence both the accuracy and the numerical robustness of the model; moreover an appropriate selection of these parameters may allow the usage of a fixed-step solver, which is more suitable for real-time simulations.


Introduction
Air conditioning and refrigeration are responsible for about 20% of the global energy demand, and this figure is expected to more than double by 2050 [1].Vapour compression systems (VCSs) are the most widespread technology in Heating, Ventilation, Air-Conditioning and Refrigeration (HVAC&R).In these systems a working fluid (refrigerant) is recirculated between two sections in which it exchanges thermal energy with two environments at different temperatures; these sections are kept at two different pressure levels through the action of a compressor and an expansion valve.In order to contain energy consumption and greenhouse gases (GHG) emissions, HVAC&R stakeholders strive, among others, to improve the energy efficiency of VCSs, through better system design and to implement suitable control strategies.Numerical modelling of VCSs can be the best way to achieve these goals: models can be used in design, system analysis, control tuning and commissioning phases, but they can also be applied as software-in-the-loop controls for performance optimization or for automatic fault detection and diagnosis (AFDD) [2].
In VCS dynamic modelling, the largest effort must be devoted to the heat exchangers for two reasons: -the dynamic behaviour of the VCSs strongly depends on the thermal transients, indeed the dominant time scale of the process is related to the heat transfer [3]; -the working fluid undergoes a non-instantaneous phase transition, in both heat exchangers for subcritical systems, or in the evaporator only for supercritical ones; this implies that the heat transfer rates and the pressure levels are coupled.According to [4], heat exchanger models can categorised as lumped parameter, finite volume and moving boundary.In the first approach, the process is described employing a single control volume, for which the balance equations, for mass, energy and momentum, are written; heat transfer is therefore based on the physical properties averaged over the control volume.This assumption is liable to introduce non-negligible inaccuracies since physical properties of the working fluid vary significantly during the transport process, especially in the two-phase region, and so does the heat transfer coefficient.At the other end, the finite volume technique divides the heat exchanger volume into a given number of cells, for each of which the same balance equations are written; in this way, the properties of the fluids can be evaluated more precisely, reaching an almost local estimate when the grid becomes fine enough, and capturing the spatial gradient of the process variable; the drawback of this approach is that it requires an increasingly demanding computational effort if discretisation needs high degrees of grid refinement.
In the last decades, the moving boundary (MB) approach has been developed, with the aim of combining the accuracy and flexibility of the FV scheme and the low computational cost of a lumped parameter approach.In the early realizations of the MB scheme [5], the heat exchanger volume is divided into a maximum of three lumped systems, depending on the phase of the working fluid (superheated vapour, two-phase mixture or subcooled liquid).The same balance equations can be used for each subsystem, but in this approach, it is possible to deal with the discontinuities associated with phase transition; indeed, the physical properties are averaged over a volume of fluid under uniform conditions and a suitable integral correlation to evaluate heat transfer is employed for each zone.The drawback is that, during large transients, the conditions of the working fluid may vary sharply, causing zones to disappear or reappear, thus leading to numerical issues during computation, [2].To avoid such complications, the MB scheme has been improved, allowing the interchange between different representations, [6], or modes, [7], defined on the basis of the number of active zones [8].Transition is controlled by the definition of some switching criteria, which require the choice of a monitored variable and of its corresponding threshold; moreover, when one zone disappears, the corresponding equations are no longer available, so they must be substituted with some pseudo state equations [6], in order to avoid discontinuities at the zone's reappearance.The resulting scheme, called switching moving boundary, (SMB), has performed well when employed in the modelling of finned-tube condensers and evaporators for the simulation of on-off transient of a complete VCSs [6,7].
The selection of the switching criterion has considerable influence on the accuracy and the stability of the model, as demonstrated in [9] in the case of an evaporator, yet no such analysis exists for condensers, nor for types of heat exchangers other than the fin-and-tube.In the present work, the effects of the activation and deactivation thresholds are numerically investigated for the case of a brazed-plate (BP) condenser; among the various formulations of the SMB scheme, the one based on the mean void fraction is chosen, in accordance with [6,7].The aim is to investigate the effects of the switching criteria on the accuracy of the model, analysing the behaviour of the predicted values of some variables of interest for the process as the thresholds vary, and on the stability and robustness of the model, exploring the numerical performance in terms of simulation time and minimum integration step size.Some useful guidelines are provided to select the optimal settings for the switching thresholds, so as to restrain the reduction in the integration step caused by the increase of both the stiffness and the ill conditioning of the system of differential equations, which occurs especially close to mode transition.The ultimate goal is to propose an SMB scheme suitable for fixed-step integration with optimal settings for the switching thresholds, to allow heat exchanges modelling to be included into the larger framework of the simulation of energy systems for control design purposes.

Model definition
In the MB description, a maximum of three control volumes can be defined in the heat exchanger, depending on the local thermodynamic state of the working fluid.Each control volume includes three elements, the working fluid, the secondary fluid and a portion of the wall of the heat exchanger, for each of which the balance equations may be written.
Consistently with previous works on SMB [5,6,8], some simplifying assumptions are made: -one-dimensional, compressible and unsteady flow for the working fluid; -one-dimensional, incompressible steady flow for the secondary fluid; -cross-sectional area constant for both fluids; -negligible pressure drops for both fluids; -negligible axial conduction; -negligible wall thermal resistance; -thermal storage in the secondary fluid is neglected.
As a consequence of the assumptions above, the momentum equation is unnecessary for both fluids; moreover, since flow and heat transfer to the secondary fluid are stationary, no differential equation is needed to describe its behaviour, the algebraic equation for heat transfer in each control volume suffices.
In the following, the integral form of the balance equations applied to each control volume for the derivation of the MB formulation are presented.
The refrigerant mass balance equation: where ρ is the density of the working fluid, V is volume occupied by the fluid, and ṁin e ṁout are the inlet and outlet mass flow rates respectively.The working fluid internal energy balance: where u is the specific internal energy of the fluid, h in and h out are the specific enthalpies of the fluid at the inlet and at the outlet and Qr is the heat transfer rate between the heat exchanger wall and the refrigerant.The wall energy balance where C w is the thermal capacity of the wall in the control volume, T w is the wall temperature and Qs is the heat transfer rate between the secondary fluid and the wall.Selection of the state variable has a key role in the realization of the model; although this choice is arbitrary, given the redundant nature of thermodynamic properties [10], it may have consequences on the switching criteria employed and on mass conservation during the transients, [11].

SMB equations
In this work, the SMB model of a BP condenser is obtained using the mean void fraction scheme; the detailed procedure for the derivation of the governing equation, from the balance equations, eqs.( 1) to (3), can be found in [12].Expanding the balance equations for each zone, a system of nine differential equations is obtained.The working fluid is characterized by two state variables for each zone; the pressure, assumed as constant throughout the heat exchanger, must be added to these six state variables.Finally the condenser wall is described by its lumped temperature in each zone.Given that the ten time derivatives of the state variables are the unknown, another equation must be integrated to close the system.According to [6], it is convenient to add a tracking equation for the mean void fraction: where γtot represents the mean void fraction corresponding to a complete phase change, complete condensation in this case, and K γ is a tunable relaxation gain.The mean void fraction is calculated employing the Zivi correlation [13].
The general form of the system of differential equations is Z(x, u) ẋ = f (x, u), where the state vector is

Refrigerant heat transfer rate
The heat transfer rate Qr,j for the working fluid in each zone is calculated through: where α r,j is the heat transfer coefficient for the j th zone, and A r is the area of the heat transfer surface in the refrigerant side.
The physical properties of the working fluid averaged for each zone are evaluated employing a suitable table compiled using the CoolProp library [14].The heat transfer coefficient for the refrigerant side is calculated for each zone employing a set of two correlations depending on its phase: for single phase in BP heat exchangers the correlation suggested by Martin, [15], is employed, while for condensation in BP heat exchangers the correlation proposed by Longo et al., [16] is adopted.

Secondary fluid heat transfer rate
Considering that the two fluids involved in the process are counter-flowing, the heat transfer rate Qs of the secondary fluid in each zone can be calculated with: where ṁs is the mass flow rate of the secondary fluid, c ps is its specific thermal capacity at constant pressure, T si,j is the temperature of the fluid at the inlet of the control volume of the j th zone , T w,j is the mean temperature of the wall, and N T U s,j is the number of transfer units.
The physical properties of the secondary fluid are evaluated at the inlet temperature.The heat transfer coefficient is calculated with the correlation recommended by Martin, [15].

Model implementation
The SMB model of a BP model presented in this work is realized in the Simulink ® environment, employing the S-function block, that is a powerful mechanism for extending the capabilities of this system implementing code blocks written in Matlab ® .
The core of the block is constituted by the system of equations which includes the ten differential equations and the algebraic heat transfer equations, eqs.( 5) and (6).This system takes different forms depending on the number of active zones, i.e. on the phase of the working fluid.
The boundary conditions for the system of equations constitute the inputs to be provided to the S-function block; these are the refrigerant mass flowrate at the inlet and outlet of the condenser, the refrigerant temperature at the inlet and the stationary mass flow rate and the inlet temperature of the secondary fluid.Further details can be found in [12].

Analysis of the switching criteria
The system of ten differential equations, is employed only when all zones are active, that is, in the heat exchanger the working fluid is present as superheated vapour, two phase mixture and subcooled liquid.This representation or mode is called V-TP-L, consistently with the convection adopted in [11,17].
In this mode the matrix Z(x, u) is full rank only if ζ j > 0 for each zone, otherwise it is no longer invertible, [2].Moreover, when a the normalized length tends to zero ζ → 0, the condition number of the matrix Z(x, u) skyrockets, causing numerical issues in the solution.The SMB scheme overcomes these problems, allowing the disappearance and the reappearance of a zone by changing the set of state equations employed according to same criteria.
In order to avoid singularities, disappearance of a zone is triggered when the corresponding normalized length reaches a minimum value, ζ min , rather than allowing the length to reach zero.The first switching criterion can be thus expressed as ζ j < ζ min .
From mode V-TP-L, the model can switch to V-TP mode if the subcooled liquid zone is extinguished; during normal operation of the condenser, disappearance of the super-heated vapour zone is unphysical, since the heat exchanger is directly fed by the compressor.When the subcooled liquid zone disappears, the corresponding differential equations are no longer available, and eq.( 4) is no longer needed.In this case, the three state variables (ζ 3 , h 2 , T w3 ) corresponding to the subcooled liquid zone must be tracked in some way in order to avoid discontinuities when the region reappears, as pointed out in [6,7,8].The normalized length ζ 3 tracks the threshold value for the minimum length ζ min , the enthalpy h 3 tracks the reference value corresponding to the saturated liquid state h l and wall temperature T w3 tracks the wall temperature in the two-phase zone, as shown in eq.(7).
When operating conditions lead to the reappearance of the sub-cooled liquid zone, the system of state equations must be changed.This is achieved again using another switching criterion, based on the mean void fraction in the two-phase zone.Indeed, if complete condensation occurs in the TP zone, and sub-cooled liquid starts to accumulate towards the outlet of the heat exchanger, the mean void fraction in the TP region is less than the value corresponding the complete phase transition γtot .So the sub-cooled zone formation can be identified by γ < γtot .
To make the two criteria more consistent, it is possible to convert comparison of the mean void fraction into comparison of the normalized lengths, since the same condition may be formulated defining a maximum value of normalized length ζ max corresponding to the volume occupied by the subcooled liquid inside the two-phase region.In this way the second criterion can be written as ζ 2 (γ − γtot ) > ζ max .

Analysis
The effects of the thresholds have been studied considering the accuracy in the evaluation of the process variables and the numerical performance of the model.The numerical test consists of a mass flow rate transient: all boundary conditions are kept constant, except the inlet mass flowrate of the refrigerant, which is subject to an increase of 10% of its initial value and then to an equal decrease before being brought back to its initial value.The input values used for the variables are ṁr,in,ref = 0.05 kg • s −1 , ṁr,out = 0.05 kg • s −1 , ṁs = 0.28 kg • s −1 , ,T f,in = 36 • C while the refrigerant mass flow rate at the condenser's inlet varies by 10% around the equilibrium value ( ṁr,in,ref ).To the Authors' best knowledge, no experimental data are available for direct comparison, since imposing the mass flowrate would imply to be able to test the heat exchanger alone, on one hand and the evaporator should also allow optical access to determine the extent of the regions (although this could be bypassed by e.g.thermographic mapping of outer surfaces), for similar reasons no experimental comaprison was carried out in [9].What can be compared is the influence of correctly assessing the switching threshold has on the overall simulation of a compression vapour system, which is beyond the scope of the paper.The total heat transfer rate, the pressure of the working fluid, and the outlet temperatures of the refrigerant have been chosen as monitored variables in order to evaluate the accuracy.
The numerical integration is carried out via a variable-step solver available in Matlab ® , the commonly used ode45 with a relative tolerance of 10 −3 .
On the other hand, a variable-step solver should not be employed in a control-oriented application involving several subsystems, which may require different integration step-size at any given time.Also, these methods are often unsuitable for use in a real-time environment due to requirements on sample inputs and update outputs at precise time intervals.This issue may be overcome using a fixed-step solver: in this case an appropriate integration step-size must be chosen to guarantee the accuracy required, stability of the solution and compatibility with larger frameworks [18].
For these reasons, the numerical performance of the model when the thresholds vary is evaluated in terms of simulation time and minimum step size; the latter is automatically chosen by the solver in order to achieve the prescribed accuracy, and it depends mostly on the condition number and the stiffness of the matrix Z(x, u) at certain point of the simulation [19].

Results ans discussion
4.1.Effect of thresholds on the accuracy of the model Reduction in thresholds' value produces a convergence of the predicted values; in fig. 1 the convergence of the dynamic behaviour of the monitored variables, pressure, heat transfer rate, outlet temperature of refrigerant when ζ min drops from 10 −1 to 10 −5 can be appreciated qualitatively.
The normalized root mean squared error (NRMSE) is chosen as a metric for convergence.The NRMSE for each monitored variable (x) is calculated with eq. ( 8) using as reference signal (x ref ), corresponding to the minimum value of both thresholds (10 −5 ) with the relative tolerance set to 10 −4 .
The threshold ζ min is demonstrated to be more influential on convergence for pressure than ζ max .A good level of of accuracy is reached for ζ min and ζ max lower than 0.5 × 10 −2 .

Effect of thresholds on the numerical performance of the model
The numerical performance of the model was evaluated through the time needed to complete the simulation and the minimum integration step-size required by the solver.The variable-step solvers adapt the integration step-size to achieve the prescribed accuracy, indeed absolute or relative tolerance must be included in their settings.If the matrix of the system Z(x, u) is ill-conditioned, accuracy of the solution deteriorates, due to numerical reasons; the solver tries to compensate for this effect reducing the integration step-size.Because of this, the minimum step-size used by a variable-step solver is a reliable indicator of the ill-conditioning and of the stiffness of the numerical problem at any time during the simulation.The stiffness of the system of differential equations Z(x, u) is evaluated through a stiffness ratio (SR), i.e. the ratio between the absolute values of the largest and the smallest eigenvalues, in modulus, λ.This is done at each time step: the SR is quite large for the whole duration of the simulation because the SMB system involves dynamics with very different time scales, which correspond to eigenvalues liable to vary of orders of magnitude.Moreover, when the normalized zone length tends to zero, the SR skyrockets, as does the condition number.The presence of such peaks in the SR time signal is evaluated using as smoothness factor (SF) the non-dimensional jerk metric defined as follows [20]: Higher values of SF correspond to smoother signals.Using this metric, the presence of peaks close to the switching points can be conveniently detected.
When the thresholds become higher, the step size increases; this behavior is shown in fig. 2  (b).For larger values of ζ max , the model is unable to capture the transition in mode; in these cases, the larger step-size corresponds to an non-physical description of the state of the working fluid in the condenser.For ζ max < 0.5 × 10 −3 , the minimum step-size shows a jump for larger

Optimal thresholds
The analysis above evidences two opposite trends.The accuracy is enhanced by the reduction of the thresholds ζ min and ζ max , yet, the stiffness of the system increases, and a smaller integration step-size is needed to satisfy accuracy requirements.A variable-step solver can do the task, while keeping simulation time within reasonable bonds.Yet, to include the SMB model into a more complex energy system framework, a fixed-step solver is preferable.The choice of the switching thresholds determines the accuracy in the computed values and the stiffness of the system of differential equations and, thus, the optimal integration step-size.This work proves that a value of the switching thresholds ζ min and ζ max between 0.5 × 10 −2 and 10 −3 , is a satisfactory compromise; the model yields a solution with a relative tolerance of 10 −3 implementing a fixed step-size in the order of 0.05s, which appears suitable for control design purposes.

Conclusions
In this work a numerical analysis of a SMB model based on the mean void fraction scheme for a BP condenser in counter-flow arrangement is carried out.The switching thresholds are checked with regard to the accuracy and the numerical performance of the simulation, studying a transient in the mass flow rate of the refrigerant.The general trend is that the smaller the thresholds, the higher the accuracy in the estimation of the process variables.On the other hand, if the threshold becomes too small, the system of differential equations becomes stiffer and may show ill conditioning.Some guidelines for the selection of the optimal thresholds are provided, with the ultimate goal to support the usage of a fixed-step solver, which is more suitable to include this condenser model in a real-time simulation framework for control design purposes.
40th UIT International Heat Transfer Conference (UIT 2023) Journal of Physics: Conference Series 2685 (2024) 012040 IOP Publishing doi:10.1088/1742-6596/2685/1/0120408 values of ζ min .This behaviour corresponds to a reduction of the peaks in the stiffness factor as shown in fig. 2 (c); indeed the smoothness factor of the SR increases for large values of ζ min .