Optimization and multivariable control of refrigeration systems

Proyecto Fin de Carrera Ingeniería de Telecomunicación Formato de Publicación de la Escuela Técnica Superior de Ingeniería Autor: F. Javier Payán Somet Tutor: Juan José Murillo Fuentes Dep. Teoría de la Señal y Comunicaciones Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2013 Tesis Doctoral Ingeniería Industrial Optimization and Multivariable Control of Refrigeration Systems Autor: Guillermo Bejarano Pellicer Directores: Manuel G. Ortega Linares Francisco Rodríguez Rubio Ingeniería de Sistemas y Automática Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2017 Tesis Doctoral Ingeniería Industrial Optimization and Multivariable Control of Refrigeration Systems Autor: Guillermo Bejarano Pellicer Directores: Manuel G. Ortega Linares Profesor Titular Francisco Rodríguez Rubio Catedrático de Universidad Ingeniería de Sistemas y Automática Escuela Técnica Superior de Ingeniería Universidad de Sevilla 2017 Tesis Doctoral: Optimization and Multivariable Control of Refrigeration Systems Autor: Guillermo Bejarano Pellicer Directores: Manuel G. Ortega Linares Francisco Rodríguez Rubio El tribunal nombrado para juzgar la Tesis arriba indicada, compuesto por los siguientes doctores: Presidente: Vocales: Secretario: acuerdan otorgarle la calificación de: El Secretario del Tribunal Fecha: To my family and Patri, the real support of my life Acknowledgements The successful completion of this Thesis is due to the effort, cooperation, support, help, understanding, and dedication of many people. May the reader forgive me for writing the acknowledgements in my mother tongue. En primer lugar, quisiera mostrar mi agradecimiento a mis directores de Tesis, Manolo Gil y Paco Rodríguez Rubio. Desde mis comienzos en tareas de investigación creyeron en mí mucho más que yo mismo, dándome múltiples oportunidades para desarrollar mis capacidades y mis ideas. Sin su ayuda y apoyo constante ninguno de mis logros habría sido posible, y esta Tesis no es más que el resultado de varios años de formación, seguimiento y guía acertada por su parte. En especial quisiera agradecer a Manolo Gil su predisposición a ayudar en cualquiera tarea, relacionada o no con la Tesis, y la confianza que me dio desde que finalicé mi Carrera, contando conmigo en múltiples proyectos de muy diferente índole, en los cuales nos embarcamos juntos sin saber muy bien cómo llegaríamos al final, pero que han resultado sin duda exitosos. Gracias por todo. También debo agradecer su intenso trabajo a todos aquellos que han participado conmigo en el desarrollo de los proyectos sobre refrigeración en los que se enmarca la Tesis. Profesores como Fernando Castaño, Manolo Vargas, y Carlos Vivas, que ya fue tutor de mi Proyecto Fin de Carrera, con diversas ocupaciones docentes e investigadoras y poco tiempo disponible, han hecho lo posible por ayudar en el desarrollo del proyecto, permitiendo con sus consejos y paciencia que hoy presente este trabajo. Asimismo, debo agradecer a Julio Normey-Rico y Alexandre Trofino-Neto su magnífica atención durante mi estancia de investigación en Brasil, así como su valiosa ayuda en todos los aspectos de una estancia que ha resultado a nivel personal extremadamente enriquecedora, además de productiva y provechosa a nivel profesional. Cómo no, tengo mucho que agradecer también a mis compañeros de trabajo. A David por su gusto artístico para las gráficas y su intensa dedicación a la planta experimental, a José Joaquín por sus conceptos de transferencia de calor que van a resultar extremadamente útiles en el proyecto sobre almacenamiento de frío, y sobre todo a José Enrique, mi II Contents 3.1.2. Sensors 58 3.1.3. Programmable Logic Controller (PLC) 59 3.1.4. Supervision software and OPC communication 61 3.2. Identification 63 3.2.1. Expansion valves 63 3.2.2. Compressors 65 3.2.3. Evaporators 68 3.2.4. Condenser 74 3.3. Validation 85 3.3.1. One-compression-stage, two-load-demand configuration 85 3.3.2. Two-compression-stage, two-load-demand configuration 91 3.4. Final remarks 95 4. Optimization and controllability analysis 97 4.1. Global optimization 98 4.1.1. Problem statement 98 4.1.2. Nonlinear steady-state model 99 4.1.3. Optimization results 100 4.2. Controllability analysis 104 4.2.1. Analysis based on linear theory 106 4.2.2. Nonlinear analysis based on the phase portrait method 108 4.3. Final remarks 115 5. Control 117 5.1. Robust H∞tracking controller 119 5.1.1. Modelling 119 5.1.2. Controllability analysis 122 5.1.3. H∞control synthesis 128 5.1.4. Tracking controller comparison 133 5.2. Suboptimal hierarchical control strategy 146 5.2.1. Overview 146 5.2.2. Optimizer 148 5.2.3. Practical Nonlinear Model Predictive Control (PNMPC) 152 5.2.4. Energy-efficiency-aware controller comparison 155 5.3. Final remarks 164 6. Conclusions 165 6.1. Thesis contributions and conclusions 165 6.2. Future work 169 Appendix A.Signal distribution and PLC bus cycle 173 A.1. Signal wiring 173 A.2. PLC bus cycle 174 Contents III List of Figures 179 List of Tables 183 Bibliography 185 Acronyms CV configuration valve DMC Dynamical Matrix Control EEV electronic expansion valve ESC extremum-seeking control FB+FF feedback-plus-feedforward FV finite-volume GPC Generalised Predictive Control IP initial point ISE Integral of Square Error ITAE Integral of Time multiplied by Absolute Error LQG linear-quadratic-Gaussian LQR linear-quadratic regulator MB moving boundary MIMO multiple-input-multiple-output MPC Model Predictive Control OLE Object Linking and Embedding OP operating point OPC OLE for Process Control P-h pressure-specific enthalpy PCM Phase change material PFC Predictive Functional Control PI proportional-integral PID proportional-integral-derivative PLC Programmable Logic Controller PNMPC Practical Nonlinear Model Predictive Control PO Percentage Overshoot PWM pulse-width-modulation QP quadratic programming V VI Acronyms RECS Residential Energy Consumption Survey RHP right-half plane RMS Root Mean Square RNGA Relative Normalized Gain Array SCADA Supervisory Control And Data Acquisition SHC suboptimal hierarchical control SISO single-input-single-output SMB switched moving boundary SVC superheat valve controller TES Thermal Energy Storage TV-ESC time-varying ESC VFD variable frequency drive VSC variable speed compressor Notation Latin symbols ALinear model state matrix AHeat transfer perimeter [m] AvExpansion valve opening [%] BLinear model input matrix bCompressor parameter bi j Element {i,j} in matrix B CControllability matrix CThermal capacity ratio CH Control horizon COP Coefficient of performance cCompressor parameter [m3] ceev Electronic expansion valve coefficient [m2] cpSpecific heat at constant pressure [J kg-1 K-1] cpuv Pressure upholding valve coefficient [m2] cvSpecific heat at constant volume [J kg-1 K-1] DLinear reduced model input vector DF(z)Low-pass-filter denominator (discrete) diElement {i} in vector D ˆ EOPi(s) Multiplicative output uncertainty transfer matrix at the operating point OPi(continuous) EMechanical energy [W h] eError vector F(z)Discrete low-pass filter fForce function in the heat exchanger model fiElement {i} in force function f GJacobian matrix G(s)Continuous transfer matrix VII VIII Notation ˆ G(s)Scaled continuous transfer matrix ˆ GOPi(s) Scaled continuous transfer matrix at the operating point OPi GP (s)Generalised plant for the S/KS/T Problem (continuous) gHeat exchanger characteristic function hSpecific enthalpy [J kg-1] JCost function K(s)Controller continuous transfer matrix ˆ K(s)Scaled controller continuous transfer matrix KStatic gain kDiscrete time step kdrop Pressure drop factor [kg-1 m-1] LHeat exchanger length [m] MI Number of manipulated inputs ˙mMass flow [kg s-1] NLF Generic nonlinear function NCompressor speed [Hz] NTU Number of Transfer Units niParameters of ceev correlation (∀i=0,1,...,4) PPressure [Pa] PH Prediction horizon QTracking error weighting matrix ˙ QCooling power [W] qVapour quality RControl effort weighting matrix RThermal resistance [m2K W-1] RE Relative error [%] SControllable subspace S0(s)Output sensitivity continuous transfer matrix SM Scaling matrix StCompressor parameter [m3] sSpecific entropy [J kg-1 K-1] T0(s) Output complementary sensitivity continuous transfer matrix Tzω(s)Generalised plant closed-loop continuous transfer matrix TTemperature [K] TSC Degree of subcooling [K] TSH Degree of superheating [K] tTime trRise time [s] tsSettling time [s] UA Global heat transfer coefficient [W K-1] uManipulable input vector VRHeat exchanger volume [m3] Notation IX vGeneralised plant measured vector vSpecific volume [m3kg-1] WKS(s)Control sensitivity weighting matrix (continuous) WS(s)Output sensitivity weighting matrix (continuous) WT(s) Output complementary sensitivity weighting matrix (continuous) ˙ WMechanical power [W] wInput vector of the heat exchanger model xState vector yOutput vector ZCoefficient matrix ZTotal number of experimental data zGeneralised plant error vector zTransfer function zero zi j Element {i,j} in coefficient matrix Z Greek symbols α Overall heat transfer coefficient between the refrigerant and the secondary flux [W m-2 K-1] βOutput sensitivity weighting matrix parameter within the H∞controller design ΓExperimental output vector in the identification process γRatio between the energy of zand ω ¯ γMean void fraction at the heat exchanger δDistance to the controllable subspace εHeat exchange effectiveness ζNormalised heat exchanger zone length ηCorrection to predicted output θParameter vector in the identification process κBandwidth parameter within the H∞controller design ΛProjection matrix λConvective heat transfer coefficient [W m-2 K-1] µOutput sensitivity weighting matrix parameter within the H∞controller design νPrediction error ξResidue vector in the identification process ρDensity [kg m-3] σSingular value τTime constant [s] υDifferential increment φMeasurable variable set ϕExperimental input matrix in the identification process χDeviation with respect to the desired state XNotation χiElement {i} in vector χ ψOptimization variable set ωExogenous signal vector ωFrequency [rad s-1] ωB Crossover frequency with -3 dB of output sensitivity function [rad s-1] ωCN Frequency where the condition number reaches its minimum [rad s-1] ωcr Achievable control bandwidth [rad s-1] ωmin,SV Crossover frequency of minimum singular value [rad s -1 ] ωT Crossover frequency of output complementary sensitivity weight [rad s-1] Subscripts and superscripts boost booster compressor ccondenser cntrl control comp compressor cond conduction cycle vapour-compression cycle diag diagonal dis discharge eevaporator e5evaporator at 5°C e20 evaporator at -20°C eev electronic expansion valve err tracking error FB feedback FF feedforward fsaturated liquid f ast fast f orced forced response f ree free response gsaturated vapour he heat exchanger in inlet int integral effect inv inverse is isentropic it iterative m1condenser mode 1 m2condenser mode 2 main main compressor Notation XI max maximum meas measured mid middle min minimum nom nominal operating point out outlet PNMPC Practical Nonlinear MPC past past and current predict predicted RHP right-half plane re f reference sat saturation sc subcooled liquid zone sec secondary flux settle settling sh superheated vapour zone slow slow start control start suc suction surr compressor surroundings tot complete phase change t p two-phase zone track tracking trns f heat transfer uinput valve expansion valve wwall wc without correction zzero 6Chapter 1. Introduction (specifically from their secondary fluxes) by evaporating the refrigerant at low pressure and temperature. The temperature and pressure of the refrigerant are increased at the compression stage. Then, heat is transferred to the secondary flux at the condenser by condensing and subsequently subcooling the initially superheated refrigerant. Pressures are usually quite different at the cooling and freezing chamber, which involves different evaporation temperatures. Thus, the pressure upholding valve just maintains the pressure difference between Evaporator 1 and Evaporator 2. The less restrictive it is, the less the pressure difference is between both evaporators, which implies that the temperature difference between the refrigerant and the secondary flux at Evaporator 1 is greater. At the output of the pressure upholding valve, two refrigerant flows are mixed together: that circulating through the Evaporator 2 at low temperature, and that expanded through the pressure upholding valve which comes from the Evaporator 1. The expansion valves allow the cycle to be closed by holding up the pressure difference between the condenser and each evaporator. Expansion valve 1 Secondary flux Expansion valve 2 Secondary flux Figure 1.4 One-compression-stage, two-load-demand refrigeration cycle. In these applications the availability of multi-compression-stage refrigeration cycles is useful to achieve greater energy efficiency. For instance, a two-compression-stage, twoload-demand refrigeration system is represented in Figure 1.5. In this case the refrigerant flow is compressed in two stages, using the main compressor and the booster one. Since there are different temperature requirements at the freezing and cooling chambers, the 1.1 Context 7 mentioned extra valve is required in order to hold up the pressure difference between the Evaporator 1 output and the discharge of the booster compressor. Expansion valve 1 Secondary flux Expansion valve 2 Secondary flux Figure 1.5 Two-compression-stage, two-load-demand refrigeration cycle. In addition to multi-stage configurations, automatic control and optimization have potential to improve the performance of current refrigeration systems and they are expected to contribute to energy efficiency enhancement. However, it requires a deep knowledge of the process and modelling of such systems with required accuracy plays a central role, since many of the developed optimization techniques for vapour-compression refrigeration cycles involve components and system modelling, at least static models. Moreover, availability of a purpose-built breadboard refrigeration plant which allows to test such optimization techniques, along with novel control strategies, is of the utmost importance to obtain real enhancements regarding energy efficiency and control performance. Some examples of experimental plants found in the literature are those described by Rasmussen et al. [ 1 ], Larsen [ 10 ], and Schurt et al. [ 11 ], but they are usually one-compression-stage refrigeration cycles, with one or more cooling demands. That is why an experimental refrigeration plant, which can be configured to work with up to two compression stages and up to two load demands, has been set up at the Department of System Engineering and Automatic Control of University of Seville. Optimization and control issues which arise when regulating actual multi-stage, multi-load-demand refrigeration systems are intended to be simulated using this experimental facility. 8Chapter 1. Introduction 1.2 State of Art 1.2.1 Problem statement In this Thesis a particular application of refrigeration systems is considered. The cycle is expected to provide a certain cooling power ˙ Qe to a continuous flow entering the evaporator as secondary flux. Neither the mass flow ˙me,sec nor the inlet temperature Te,sec,in of such secondary flux are to be controlled, since they may be managed by another high-level controller. Therefore, the cooling demand can be expressed as a reference on the outlet temperature of the evaporator secondary flux Te,sec,out , where the mass flow and inlet temperature behave as measurable disturbances to the refrigeration system. The difference with respect to the conventional case analysed in the literature lies in not considering the secondary mass flow as manipulable, since conventionally the cooling demand is merely a certain thermal power to be provided to the secondary flux for any mass flow, which is used as an additional manipulated input. Regarding the condenser, the inlet temperature Tc,sec,in and mass flow ˙mc,sec of the secondary flux are also considered as disturbances. It is shown in the literature that adding controllable inputs to the problem does not remove its underactuated features. Some works study the degrees of freedom of a vapourcompression cycle [ 12 , 13 , 14 ]. Particularly, Jensen and Skogestad consider more manipulated variables (four, including variable heat transfer coefficients at the evaporator and condenser, or equivalently considering adjustable secondary fluxes) than those considered in this Thesis [ 13 ]. Jensen and Skogestad state that five design specifications must be proposed to define a cycle given equipment, while only four controllable inputs are available. The fifth variable to be manipulated in operation is the so-called active charge, related to the refrigerant mass flow. One of the main conclusions of the work by Jensen and Skogestad is that this additional degree of freedom related to the active charge cannot be manipulated without introducing a variable-level liquid tank in the cycle, which implies modifying the topology of the system. Therefore, adding the secondary mass flows as manipulated inputs only increases the complexity of the control system and it does not help reduce the degree of underactuation, neither solves it the controllability issues, as described later in this Thesis. Then, only the compressor speed N and the valve opening Av are considered as manipulated variables, whereas both secondary flux mass flows and inlet temperatures act as measurable disturbances, as depicted in the schematic drawing shown in Figure 1.6. 1.2.2 Modelling Refrigeration systems are, as generally known, closed cycles, whose components are connected through various pipes and valves, which causes strong nonlinearities and high coupling between variables. This is why dynamic modelling of vapour-compression refrigeration systems is not definitely trivial matter. The most important elements regarding the dynamic modelling are the heat exchangers, since the expansion valve, the compressor, and the thermal behaviour of the secondary fluxes are statically modelled, because their dynamics are usually at least one order of magnitude faster than those of the evaporator and condenser. 1.2 State of Art 9 COMPRESSOR EXPANSION VALVE EVAPORATOR CONDENSER Tsurr N Av ˙mc,sec Tc,sec,in ˙me,sec Te,sec,in hc,in ˙mc,in Pc ˙me,in he,in Pe he,out hc,out Figure 1.6 Diagram of variables involved in the cycle (blue colour refers to manipulated variables, red to disturbances and green to internal variables). A very detailed model of a heat exchanger is based on mass, energy, and momentum balances of the refrigerant, the secondary flux, and the material separating them which comprises the heat exchanger itself. This approach, dating to MacArthur [ 15 ], involves spatially discretizing the heat exchanger into an arbitrary number of control volumes [ 16 , 17], and thus leading to a numerical solution of a set of differential equations discretized into a finite difference form [ 18 , 19 ]. The finite-volume (FV) approach provides very detailed knowledge about the system statics and dynamics, but due to its computational cost and complexity it is inappropriate for identification purposes and model-based control strategies. A simpler model with better balance between accuracy and computational cost may be obtained using the moving boundary (MB) approach. This methodology divides the heat exchanger into a number of zones corresponding to different refrigerant states: superheated vapour, two-phase fluid, and/or subcooled liquid [ 1 , 20 ]. Mass and energy balances are applied to each zone and, taking into account theses balances as well as other system constraints, the refrigerant variables at the heat exchanger outlet are obtained. The zone lengths are state variables, since they can vary with time depending on inputs and disturbances. A step further is developed by McKinley and Alleyne: some control volumes are allowed to completely disappear and reappear without simulation issues, giving rise to the switched moving boundary (SMB) model [ 21 ]. Different representations of the heat exchanger model, also known as modes, are defined depending on the existence or absence of each zone. For example, in the case of the evaporator, two modes are defined depending on the amount of superheated vapour, as seen in Figure 1.7. In the case of the condenser, 10 Chapter 1. Introduction up to five different modes are considered, as shown in Figure 1.8. This model is validated and extended to startup and shutdown processes by Li and Alleyne [22, 23]. EVAPORATOR Superheated Two-phase Two-phase Mode 1 Mode 2 Figure 1.7 Evaporator modes according to the SMB approach [23]. Mode 1 Mode 3 Mode 4 Mode 5 Superheated Superheated Superheated CONDENSER Subcooled Two-phase Two-phase Mode 2 Subcooled Two-phase Two-phase Figure 1.8 Condenser modes according to the SMB approach [23]. The complexity, computational load, and accuracy of the SMB model have been recently compared to those of the FV formulation [ 24 ]. It has been concluded that while the SMB approach can execute much more quickly in simulation than the FV approach, there is little difference in the achievable accuracy with respect to experimental data. 1.2 State of Art 11 In addition to the FV and MB approaches, black-box models can be found in the literature, where any variable defining the cooling power is to be predicted from the manipulated variables. For instance, Bittanti and Piroddi propose a neural network approach to identification of a heat exchanger [ 25 ], while Romero et al. propose a simplified black-box model to predict accurately the chilled water temperature dynamic response of a vapourcompression chiller, where the Box-Jenkins structure gets the best fit to experimental results [ 26 ]. This modelling approach is suitable for control purposes, but only the identified output variables can be controlled. 1.2.3 Optimization and Control To achieve high energy efficiency while satisfying the cooling demand, it must be taken into account that heat transfer at the evaporator is determining for the overall efficiency. Heat transfer is widely recognised to be much higher when the refrigerant flow is twophase. Thus, the highest evaporator efficiency would be achieved if the refrigerant at the evaporator outlet was saturated vapour. This ideal behaviour is not advisable nor applicable in practice, since the risk of liquid droplets appearing at the evaporator outlet is very high in transient, which must be definitely avoided because the evaporator outlet matches the compressor intake. Therefore, the approach conventionally applied in industry consists in operating the cycle with a certain degree of superheating of the refrigerant at the evaporator outlet ( TSH ), which is held low to approximate to the ideal behaviour previously described. Therefore, the conventional control scheme is very simple: in addition to the reference on Te,sec,out imposed by the cooling demand, a low but constant set point on TSH is applied and the controller is designed to get these two variables to track their references as efficiently as possible in presence of disturbances by manipulating Nand Av. To design the tracking controller used within the conventional control scheme, it is important to take into account that the difficulty in controlling this process lies in high thermal inertia, dead times, high coupling between variables, and strong nonlinearities. The most used linear techniques which can be found in the literature regarding tracking control are decentralised control [ 27 , 28 , 29 , 30 ]; decoupling multivariable control [ 31 ], LQG control [ 11 , 32 , 33 ]; model predictive control (MPC) [ 34 , 35 , 36 , 37 ], and robust H∞control [38]. Regarding decentralised control, Marcinichen et al. have proposed a SISO dual strategy for N and Av simultaneous control using PI controllers and the most used matching in industry: regulating the cooling capacity (namely Te,sec,out ) through N and the degree of superheating with Av . Wang et al. implement an hybrid PID-Neural Network controller, which consists in a neural network which tunes the PID parameters online [ 28 ]. Underwood proposes a control strategy where two PID controllers (one for each manipulated variable) are jointly tuned through optimisation techniques [ 29 ]. Meanwhile, Salazar and Méndez apply a decentralised PID control over a single-stage transcritical CO 2 refrigeration cycle, where a lumped energy balance model to derive a set of non-linear first order differential equations for the gas cooler and the evaporator is proposed [ 30 ]. Two insights are considered: linear PID control (by means of linearisation of the obtained models) and PID feedforward compensation techniques. 12 Chapter 1. Introduction Furthermore, Shen et al. propose a decoupling control system design based on the Relative Normalized Gain Array (RNGA) [ 31 ]. Using the gain and phase information provided by the RNGA, an equivalent transfer function matrix for a closed-loop control system is obtained, and its relationship with the inverse of the process transfer matrix is established. Moreover, the criterion to determine a stable, proper, and causal ideal-diagonal decoupler is established. The main advantage of this method is its simplicity, since it does not require extensive calculation effort. As fas as LQG control is concerned, similar approaches have been proposed by He and Schurt et al. [ 11 , 32 , 33 ]. The three works design a MIMO controller based on the LQG method using a Kalman filter for the estimator design. An integrator is incorporated to the design to reject disturbances and also to keep the time response of the closed-loop system approximately the same as that of the open-loop system. Schurt et al. use only a linearised model to derive the controller [ 11 , 33 ], while He linearises the system in a few operating points and a gain-scheduling scheme is proposed [32]. Concerning MPC, Razi et al. propose a neuro-predictive control algorithm, where a neural network is used to estimate the state and the MPC calculates the control actions, compensating the delays of the model [ 34 ]. The performance of the neuro-predictive control algorithm is compared with that of a fuzzy PID controller, showing better performance although more computational load. Furthermore, Sarabia et al. study a hybrid MPC controller applied to a supermarket refrigeration system, where the compressor speed and valve opening are only on-off variables [ 35 ]. This study is extended by Ricker and the approach suggested turns out to be scalable to large systems, since the computational effort is modest [ 36 ]. It increases linearly with the number of display cases, but it is independent of the number of compressors, since online optimization is avoided. Fallahsohi et al., in turn, propose a Predictive Functional Control (PFC) to regulate the degree of superheating using the expansion valve opening as manipulated variable [ 37 ]. This method is shown to be more accurate than the conventional PID control. As the PFC improves disturbance rejection compared to a PID control, it is possible to reduce the reference on TSH and to prevent any unevaporated refrigerant liquid from reaching the compressor. As a consequence, the use of PFC leads to an increase of COP which depends on operating conditions. Regarding the robust control strategy, Larsen and Holm have designed a MIMO H∞ controller which solves the S/KS problem [ 38 ]. For the considered system the coupling between Te,sec,out and Av is weak, which is taken into account when designing the reducedorder controller. A comparison between the designed controller and a SISO one (calculated by easing the closed-loop system bandwidth but following the same procedure) is also addressed. In this work several operating points are not taken into account when solving the S/KS problem, but only uncertainties along the frequency are considered. However, optimal steady-state operation is not explored through the references in the conventional control scheme, but the low degree of superheating set as reference is believed to lead the cycle to high energy efficiency. Nevertheless, some studies have shown that there are plenty more scope for efficiency improvement if optimal references are imposed on the evaporator and condenser pressures [ 39 , 40 ]. Alternatively to conventional set point selection, it is possible to perform a global optimization in order to obtain the optimal operating point, given a certain cooling demand, a specific facility and some 1.2 State of Art 13 constraints regarding the secondary fluxes. A steady-state model of the cycle is required, being the solution of the global optimization a cycle which generates the desired cooling power with as high as possible energy efficiency. This efficiency is measured by any metrics. The solution should also observe all physical and technological constraints of the components, in addition to some operating constraints, for example that of minimum degree of superheating for safe operation of the compressor. Within this line of research Jain and Alleyne implement a exergy-based global optimization [ 41 ]. The destroyed exergy used as metrics presents some advantages over other first-law metrics, such as energy consumption, since it is a measure of how effectively the exergy supplied to a system, in this case work done on the system by the compressor, is used. It also enables the objective function to be scalable with respect to subsystem configuration and subsystem capacity. The optimal cycle and simultaneously the optimal values of the manipulated variables are calculated. Regarding the optimization procedure, Zhao et al. develop a modified genetic algorithm combined with a solution strategy for a group of nonlinear equations to minimize energy consumption [ 42 ]. Some experimental results show that the calculated set points may reduce energy consumption compared to traditional on-off control. Zhao et al. also present a steady-state decentralised optimization, where total energy consumption is minimized [ 43 ]. The novelty consists in the decentralised algorithm, which divides the global problem into three optimization subproblems. They are constrained by some shared variables which are updated online. The decentralised procedure achieves good accuracy compared to a global optimization while dramatically reducing the computation time. Regarding the energy-efficiency-aware control strategies, Jain presents a multivariable feedback-plus-feedforward (FB+FF) control strategy to achieve the set points generated by the offline global optimizer previously commented [ 12 ]. The latter also generates the values of the steady-state control actions, which are used as feedforward contribution to the control law. The feedback controller is a linear-quadratic regulator (LQR) based on a linear model identified around the desired point. The optimal cycle generated by the optimizer is defined by three variables, but only two control actions are available, thus the optimal solution is projected from the three-degree-of-freedom optimization space onto the two-degree-of-freedom control space, in such a way that the achieved steady-state cycle is no longer optimal, but suboptimal. Given this limitation, Jain and Alleyne propose a MPC which uses a dynamic exergybased cost function to determine the optimal control actions. The objective is to maximize the exergetic efficiency while achieving the desired cooling load [ 44 ]. That is, online optimization and control are achieved at once. Simulation results comparing the exergybased MPC and a first-law MPC are presented, showing that optimal cycles for a given facility and a desired cooling demand are slightly different when optimizing the COP and the exergetic efficiency. Jensen and Skogestad deal with the selection of controlled variables which allow to achieve in practice close-to-optimal operation with a constant set point policy, considering some modifications on the basic cycle such as a liquid receiver on the low pressure side which ensures that the vapour entering the compressor is saturated [45]. Larsen et al. propose an online optimization where no dynamic model is needed [ 46 ]. The authors consider four control actions and two major control objectives: to 14 Chapter 1. Introduction satisfy the cooling demand and to hold the degree of superheating low and constant. The two remaining control actions are devoted to controlling the condenser and evaporator pressures independently, in such a way that the total power consumption is minimized. The optimization is based on the power gradient estimation, which must be convex, and the derivation of a simple static model. Other line of research concerning energy-efficiency-aware control is extremum-seeking control (ESC). Burns and Laughman propose an online optimization of energy consumption using ESC, that means, without relying on a dynamic model of the system [ 47 ]. It requires that the relationship between a system input and the output to be optimized is convex. In the context of this application, convexity is required in the input-output map between one of the vapour compression system inputs and the electric power consumption. The evaporator fan speed is selected as the variable to be disturbed, whereas the remaining inputs are devoted to regulating the cooling demand and the system states by means of a feedback controller. The ESC is designed so that the adaptation is slow compared to dominant dynamics of the cycle and the regulating controller. Guay has recently proposed a novel technique called time-varying ESC (TV-ESC) with greater convergence rate than traditional perturbation-based ESC [ 48 ], which has been successfully applied to refrigeration systems [ 49 ]. The control architecture is also modified, since the ESC deals with the adjustment of the set point of the compressor discharge temperature, whose relationship with power consumption is shown to be convex, while the compressor speed regulates the cooling power. The temperature at the compressor discharge is chosen instead of the degree of superheating, since the first one can be easily measured because the refrigerant at this location in the cycle is always superheated vapour, whereas the last one is often difficult to measure, it is not defined for values less than zero and produces no change in sensible temperature when two-phase refrigerant exits the evaporator. 1.3 Objectives In recent years energy efficient control of existing vapour-compression systems has attracted great attention. Several elements are required to achieve optimal operation of a refrigeration system from the point of view of energy efficiency. If a classical hierarchical control structure is proposed, as shown in Figure 1.9, where a reference generator calculates the set points for the controlled variables and the controller merely gets the outputs to track their references, three key factors are involved: •Modelling : the static and dynamic behaviour of the whole cycle must be modelled for optimization purposes and controller design, thus modelling plays a key role. •Optimization : a model-based optimization stage is intended to generate optimal references to the control system, considering energy efficiency and power consumption, and using a steady-state model of the cycle. •Control : the control stage must be designed based on the dynamic model in such a way that the calculated references are achieved despite uncertainties and disturbances. The main objective of the Thesis is to analyse in depth these three factors. 1.3 Objectives 15 Cooling demand Optimal set points Control actions Measured variables Disturbances SYSTEM OPTIMIZER Efficiency metrics Static model CONTROLLER Dynamic model Figure 1.9 Hierarchical optimization and control structure. Regarding modelling, different approaches have been developed in the literature, focusing on the heat exchangers, because their dynamics are usually at least one order of magnitude slower than those of the remaining elements of the cycle. The modelling technique achieving best balance between accuracy and computational cost is the MB approach. However, although it has been verified that the complexity and computational load of the MB model is minor compared to the FV formulation while achieving comparable accuracy, the model order remains high enough to limit its use to simulation and hinders its integration in model-based control strategies, since the simplest refrigeration cycle is made up of two heat exchangers and each one generates a high-order state vector. The first objective of the Thesis is to develop a simplified control-oriented dynamic model of a one-stage, one-load-demand cycle, which might be more suitable to model-based control strategies regarding computational load, but without excessive inaccuracy concerning the dominant dynamics. Moreover, the optimization of a one-stage, one-load-demand cycle is to be analysed. Given a certain cooling demand, the objective of the optimization stage is to calculate the steady-state cycle achieving the highest energy efficiency which satisfies the cooling demand while satisfying some operating constraints. This optimal cycle implies an upper bound on the achievable energy efficiency, whereas some of the variables defining the cycle are to be set as optimal references to the tracking controller. Nevertheless, similar approaches in the literature have been shown to fail in achieving the optimal cycle in the control stage. That is why the controllability of the one-stage, one-load-demand cycle is to be analysed using linear theory and a nonlinear pointwise analysis based on the phase portrait method. Concerning control, the objective is to make a contribution in the two main schemes in the literature about control of refrigeration systems: the conventional scheme and the energy-efficiency-aware one. In the first one, in addition to the reference on the outlet temperature of the evaporator secondary flux imposed by the cooling demand, a low but constant set point on the degree of superheating is applied, and the controller is merely designed to get these two variables to track their references. However, energy efficiency is explicitly considered when calculating the references for the low-level controller in the second scheme. Regarding the tracking control, a multivariable centralised H∞ controller, based on the 22 Chapter 2. Modelling Different modelling techniques are considered and the MB approach is selected to model the refrigerant behaviour at the heat exchangers, namely the evaporator and condenser. Steady-state models are developed and detailed for all elements, paying special attention to the heat exchangers. Furthermore, dynamic modelling is also addressed. The MB approach is selected due to its trade-off between accuracy and computational cost. However, although the complexity and computational load of the MB model is minor compared to the FV formulation while achieving comparable accuracy [ 24 ], up to nine state variables for the condenser and up to six variables for the evaporator are necessary to completely describe the system behaviour. It makes this model still too complex to be used within a model-based strategy. Some considerations are taken into account in order to reduce the model order and therefore the complexity of the MB formulation. Typical internal volumes of evaporator and condenser are taken into account, as well as their consequences for the system dynamics. Moreover, fast dynamics of some states of the condenser model are disregarded, giving rise to a simplified control-oriented nonlinear model of the whole cycle, characterised by only three state variables. The complete evaporator state is disregarded as well as some fast states of the condenser, which implies a novelty in modelling of refrigeration systems and reduces the model complexity and computational cost. Some simulation results are presented comparing the original MB model with the control-oriented MB when some step changes on manipulated variables are imposed. Inaccuracy of the control-oriented dynamic model with respect to the original MB model is analysed as well as the difference concerning computational load. The contents of this Chapter have been included in some of the publications presented in Chapter 1, Section 1.5 [51, 53, 60, 65]. 2.1 Steady-state modelling Before addressing the dynamic modelling, the steady-state behaviour of a typical one-stage, one-load-demand cycle is studied. Nonlinear steady-state models are presented below for all elements in the cycle. 2.1.1 Expansion valve The steady-state model of the expansion valve is based on the methodology presented by Schurt et al. [ 11 , 33 ]. This model receives the pressures Pc , Pe , and the inlet specific enthalpy hc,out as inputs, and the output is the expanded refrigerant mass flow ˙me , as shown in Equation Set (2.1) . ˙me is denoted this way in order to highlight that it matches the inlet mass flow of the evaporator in a one-stage, one-load-demand cycle, whereas Av is the expansion valve opening. ˙me=ceev Avq2ρc,out (Pc−Pe) he,in =hc,out (2.1) ρc,out =ρ(Pc,hc,out ) is a refrigerant-specific function to compute the corresponding thermodynamic property. The CoolProp tool is used to calculate all properties of the 2.1 Steady-state modelling 23 refrigerant and secondary fluxes [ 69 ]. The expansion through the valve is considered to be isenthalpic, as indicated in Equation Set (2.1) . The coefficient ceev is a characteristic parameter of the valve. 2.1.2 Compressor Similarly, the steady-state model of the compressor is based on the method described by Schurt et al. [ 11 , 33 ]. The model receives the pressures Pc and Pe and the suction specific enthalpy he,out as inputs, whereas the outputs are the refrigerant mass flow compressed and discharged to the condenser ˙mc , the discharge specific enthalpy hc,in , and the power consumption of the compressor ˙ Wcomp , as indicated in Equation Set (2.2) . Note again that the notation is adapted to the one-stage, one-load-demand configuration, since the compressor suction matches the evaporator outlet and the compressor discharge matches the condenser inlet. ˙mc="St−c Pc Pecv,e,g cp,e,g−1!# N ve,out ˙ Wcomp =b˙mc(hc,in,is −he,out ) Tc,in,is =Tc+hc,in,is −hc,g cp,c,g hc,in =he,out + ˙ Wcomp −UA(Tc,in,is −Tsurr) ˙mc (2.2) The parameters b , c , St , and UA are characteristics of the compressor. Some refrigerantspecific functions are used to estimate the thermodynamic properties; their prototypes and arguments are indicated in Table 2.1. Note that q refers to the vapour quality and q=1 involves that the corresponding thermodynamic property is calculated considering saturated vapour. Once again the CoolProp tool is used [69]. Table 2.1 Functions used to calculate the thermodynamic properties of the refrigerant involved in the compressor model. Variable Function Argument 1 Argument 2 ve,out v(P,h)Pehe,out cv,e,gcv(P,q)Pe1 cp,e,gcp(P,q)Pe1 se,out s(P,h)Pehe,out hc,in,is h(P,s)Pcse,out TcT(P,q)Pc1 hc,gh(P,q)Pc1 cp,c,gcp(P,q)Pc1 24 Chapter 2. Modelling 2.1.3 Evaporator The equations which describe the steady-state model of the evaporator depend on its specific configuration. In the case of study, a brazed-plate counter-current heat exchanger is considered and hereafter modelled. A diagram of heat transfer along the evaporator is represented in Figure 2.1(a), whereas Figure 2.1(b) shows a qualitative temperature diagram of both refrigerant and secondary flux as they pass through the evaporator. According to this configuration, the equations which describe heat transfer at the evaporator are expressed in Equation Set (2.3) , regarding the superheated vapour zone, and Equation Set (2.4), for the two-phase zone. The effectiveness-NTU method is used when calculating both partial thermal powers ˙ Qe,sh and ˙ Qe,t p [70, 71]. ˙ Qe,sh =εe,sh(˙mcp)min(Te,sec,in −Te) Ce,sh =(˙mcp)min (˙mcp)max NTUe,sh =UAe,sh (˙mcp)min UAe,sh =αe,sh (1−ζe,t p)Ae,trns f Le εe,sh =ghe(NTUe,sh,Ce,sh) ˙ Qe,sh =˙me(he,out −he,g) Te,sec,mid =Te,sec,in − ˙ Qe,sh ˙me,sec cp,e,sec (2.3) ˙ Qe,t p =εe,t p(˙mcp)min(Te,sec,mid −Te) Ce,t p =(˙mcp)min (˙mcp)max ={phase change}=0 NTUe,t p =UAe,t p (˙mcp)min UAe,t p =αe,t p ζe,t p Ae,trns f Le εe,t p =ghe(NTUe,t p,Ce,t p) ˙ Qe,t p =˙me(he,g−he,in) Te,sec,out =Te,sec,mid − ˙ Qe,t p ˙me,sec cp,e,sec (2.4) The function ε=ghe(NTU,C) in Equation Sets (2.3) and (2.4) refers to the heatexchanger-specific curve used within the effectiveness-NTU method [ 72 ], whereas Te is the saturation temperature at evaporator pressure Pe . The thermodynamic properties of the refrigerant at the evaporator outlet are calculated using Equation Set (2.5). Te,out =T(Pe,he,out )TSH =Te,out −Te(2.5) 2.1 Steady-state modelling 25 SuperheatedTwo-phase Plate Plate Plate Secondary flux ˙me he,in ˙me,sec Te,sec,out Te Pe ˙me he,out ˙me,sec Te,sec,in Te,sec,mid ζe,tp ζe,sh ˙ Qe,tp ˙ Qe,sh (a) Heat transfer diagram along the evaporator according to the brazed-plate counter-current configuration 0 1 Evaporator normalised length Temperature Temperature diagram at the evaporator Te,sec,out Te,sec,mid Te,sec,in Te,out Te Te ˙ Qe,sh ˙ Qe,tp TWO−PHASE ZONE ζe,tp Refrigerant Secondary flux SUPERHEATED VAPOUR ZONE TSH (b) Temperature diagram of both the refrigerant and secondary flux along the evaporator length Figure 2.1 Heat transfer and temperature diagram of both the refrigerant and secondary flux along the evaporator length according to the brazed-plate counter-current configuration. 26 Chapter 2. Modelling Once again Te,out =T(Pe,he,out ) is a thermodynamic function included in the CoolProp tool. The degree of superheating TSH is calculated as the difference between the refrigerant outlet temperature and the saturation temperature at evaporator pressure. αe,sh and αe,t p are the overall heat transfer coefficients in the superheated vapour zone and the two-phase zone, respectively. Note that in a one-stage, one-load-demand cycle, a certain degree of superheating is desired for the sake of safe operation of the compressor, then the evaporator mode 1 is the most probable in steady state. This is the reason why only the evaporator mode 1 has been described and modelled in this Subsection. Anyway, the steady-state model in mode 2 can be easily developed from Equation Set (2.4) imposing that ζe,t p =1. 2.1.4 Condenser The equations which describe the steady-state model of the condenser also depend on its particular configuration. In the case of study, an air-cooled cross-flow heat exchanger is modelled. A diagram of the heat transfer along the condenser is represented in Figure 2.2(a), whereas Figure 2.2(b) shows a qualitative temperature diagram of the refrigerant and secondary flux as they pass through the condenser when operating in mode 1. The same is represented in Figures 2.3(a) and 2.3(b) when the condenser works in mode 2. Note that the most plausible condenser modes in steady state are modes 1 and 2, being possible the remaining ones shown in Figure 1.8 in startup and shutdown processes, as detailed by Li [22]. Superheated Two-phase Tube Secondary flux Subcooled Secondary flux Tube ˙mc hc,out Tc,sec,sc,out Tc,sec,tp,out Tc,sec,sh,out Pc ˙mc hc,in Tc ˙mc,sec ˙mc,sec Tc,sec,in Tc,sec,out ˙ Qc,sc ˙ Qc,sh ˙ Qc,tp ζc,sh ζc,tp ζc,sc (a) Heat transfer diagram along the condenser according to the cross-flow configuration and considering mode 1 Figure 2.2 Heat transfer and temperature diagram of both the refrigerant and secondary flux along the condenser length according to the cross-flow configuration and considering mode 1. 2.1 Steady-state modelling 27 0 1 Condenser normalised length Temperature Temperature diagram at the condenser in mode 1 Tc,sec,in Tc,sec,sh,out Tc,sec,in Tc,sec,tp,out Tc,sec,in Tc,sec,sc,out Tc,in TcTc Tc,out ˙ Qc,sh ˙ Qc,tp ˙ Qc,sc TWO−PHASE ZONE ζc,sh ζc,sh +ζc,tp = 1 −ζc,sc TSC Refrigerant Secondary flux SUBCOOLED LIQUID ZONE SUPERHEATED VAPOUR ZONE (b) Temperature diagram along the condenser length considering mode 1 Figure 2.2 Heat transfer and temperature diagram of both the refrigerant and secondary flux along the condenser length according to the cross-flow configuration and considering mode 1. According to this configuration, the equations which describe heat transfer at the condenser are gathered in Equation Sets (2.6) – (2.8) , regarding the superheated vapour zone, the two-phase zone, and the subcooled liquid zone (only in mode 1). The effectiveness-NTU method is used again when calculating all partial thermal powers [70, 71]. ˙ Qc,sh =εc,sh(˙mcp)min(Tc,in −Tc,sec,in) Cc,sh =(˙mcp)min (˙mcp)max NTUc,sh =UAc,sh (˙mcp)min UAc,sh =αc,sh ζc,sh Ac,trns f Lc εc,sh =ghe(NTUc,sh,Cc,sh) ˙ Qc,sh =˙mc(hc,in −hc,g) (2.6) 28 Chapter 2. Modelling Superheated Tube Secondary flux Secondary flux Tube Two-phase ˙mc hc,out Tc,sec,tp,out Tc,sec,sh,out Pc ˙mc hc,in Tc ˙mc,sec ˙mc,sec Tc,sec,in Tc,sec,out ˙ Qc,sh ˙ Qc,tp ζc,sh ζc,tp (a) Heat transfer diagram along the condenser according to the cross-flow configuration and considering mode 2 0 1 Condenser normalised length Temperature Temperature diagram at the condenser in mode 2 Tc,sec,in Tc,sec,sh,out Tc,sec,in Tc,sec,tp,out Tc,in TcTc,out ˙ Qc,sh ˙ Qc,tp TWO−PHASE ZONE ζc,sh Refrigerant Secondary flux SUPERHEATED VAPOUR ZONE (b) Temperature diagram along the condenser length considering mode 2 Figure 2.3 Heat transfer and temperature diagram of both the refrigerant and secondary flux along the condenser length according to the cross-flow configuration and considering mode 2. 2.1 Steady-state modelling 29 ˙ Qc,t p =εc,t p(˙mcp)min(Tc−Tc,sec,in) Cc,t p =(˙mcp)min (˙mcp)max ={phase change}=0 NTUc,t p =UAc,t p (˙mcp)min UAc,t p =αc,t p ζc,t p Ac,trns f Lc εc,t p =ghe(NTUc,t p,Cc,t p) ˙ Qc,t p =˙mc(hc,g−hc,f) (2.7) ˙ Qc,sc =εc,sc(˙mcp)min(Tc−Tc,sec,in) Cc,sc =(˙mcp)min (˙mcp)max NTUc,sc =UAc,sc (˙mcp)min UAc,sc =αc,sc ζc,sc Ac,trns f Lc εc,sc =ghe(NTUc,sc,Cc,sc) ˙ Qc,sc =˙mc(hc,f−hc,out ) (2.8) When modelling the mode 2 the equations which describe the two-phase zone are slightly different, as indicated in Equation Set (2.9). ˙ Qc,t p =εc,t p(˙mcp)min(Tc−Tc,sec,in) Cc,t p =(˙mcp)min (˙mcp)max ={phase change}=0 NTUc,t p =UAc,t p (˙mcp)min UAc,t p =αc,t p ζc,t p Ac,trns f Lc εc,t p =ghe(NTUc,t p,Cc,t p) ˙ Qc,t p =˙mc(hc,g−hc,out ) (2.9) Once again the function ε=ghe(NTU,C) in Equation Sets (2.6) – (2.9) refers to the heat-exchanger-specific curve used within the effectiveness-NTU method [ 72 ]. αc,sh , αc,t p , and αc,sc are the overall heat transfer coefficients in the superheated vapour zone, the two-phase zone and the subcooled liquid zone, respectively. Tc refers to the saturation temperature at condenser pressure Pc . Eventually, Equation Sets (2.10) or (2.11) are used to calculate the refrigerant thermodynamic state at the condenser outlet, in mode 1 or mode 30 Chapter 2. Modelling 2, respectively. hc,out =hc,sc Tc,out =T(Pc,hc,out ) TSC =Tc−Tc,out (2.10) ¯qc=Thom−1(Pc,¯ γc) qc,out =2 ¯qc−1 hc,out = (1−qc,out )hc,f+qc,out hc,g Tc,out =Tc TSC =0K (2.11) Tc,out =T(Pc,hc,out ) is calculated using the CoolProp tool, and the degree of subcooling TSC in mode 1 is defined as the difference between the saturation temperature Tc and the refrigerant temperature at the condenser outlet Tc,out , whereas in mode 2 the function ¯q=Thom−1(P,¯ γ) refers to the void fraction correlation [ 73 ], in this case inversely applied, which allows to calculate the mean vapour quality at the two-phase zone ¯qc and consequently the outlet vapour quality qc,out and the outlet specific enthalpy hc,out . Finally, the secondary flux outlet temperature Tc,sec,out is calculated as a weighted average of the partial outlet temperatures Tc,sec,sh,out , Tc,sec,t p,out , and Tc,sec,sc,out , in accordance with the zone lengths ζc,sh,ζc,t p, and ζc,sc, as indicated in Equation Set (2.12). Tc,sec,sh,out =Tc,sec,in + ˙ Qc,sh ˙mc,sec ζc,sh cp,c,sec Tc,sec,t p,out =Tc,sec,in + ˙ Qc,t p ˙mc,sec ζc,t p cp,c,sec Tc,sec,sc,out =Tc,sec,in + ˙ Qc,sc ˙mc,sec ζc,sc cp,c,sec Tc,sec,out =ζc,sh Tc,sec,sh,out +ζc,t p Tc,sec,t p,out +ζc,sc Tc,sec,sc,out (2.12) 2.2 Dynamic modelling 2.2.1 Starting point As commented in Chapter 1, only the heat exchangers are usually dynamically modelled, since their dynamics are dominant over the remaining components of the cycle. The dynamic model of the heat exchangers is based on the SMB model [ 21 , 23 ]. As stated in Chapter 1, this methodology divides the heat exchanging volume into variablelength zones, namely superheated vapour, two-phase and/or subcooled liquid zones, giving rise to different modes, depending on the refrigerant specific enthalpy at the heat exchanger 2.2 Dynamic modelling 31 inlet and outlet, related to saturated liquid and vapour specific enthalpies. Hence, the evaporator modes differentiate according to the existence (mode 1) or not (mode 2) of superheated vapour at its outlet, whereas the five condenser modes are defined by the refrigerant enthalpy at its inlet and outlet, in relation to saturated vapour and liquid specific enthalpies (see Figures 1.7 and 1.8). For all modes, the state vectors in Equation Set (2.13) gather all relevant information concerning the condenser and the evaporator at each instant. The state vectors xc and xe include heterogeneous variables such as the condensation and evaporation pressures Pc and Pe , specific enthalpies describing the single-phase zones such as hc,sh , hc,sc , and he,sh , the zone lengths ζc,sh , ζc,t p , and ζe,t p , and the mean void fractions within both heat exchangers ¯ γcand ¯ γe. xc= [hc,sh Pchc,sc ζc,sh ζc,t p ¯ γc]T xe= [ζe,t p Pehe,sh ¯ γe]T xcycle = [xcxe]T (2.13) These state vectors have been obtained by removing the wall temperature states from the original modelling developed by Li and Alleyne [ 23 ], since the thermal powers are calculated using the effectiveness-NTU method [ 70 , 71 ]. The uniformity of the state vectors, regardless of the mode, allows to retain a constant dynamic structure formulated in the nonlinear descriptor form shown in Equation Set (2.14) , for the evaporator and condenser, where wc and we refer to the input vector to the corresponding heat exchanger, as explained by Li and Alleyne [23]. Zc(xc,wc)˙xc=fc(xc,wc) Ze(xe,we)˙xe=fe(xe,we)(2.14) Coefficient matrices Zi(xi,wi) and forcing functions fi(xi,wi)∀i=c,e can be derived from the governing equations detailed in the related literature for each mode, which include refrigerant thermodynamic properties and geometric parameters, in addition to mass flow and energy balance terms [ 21 , 22 ]. It is important to remark that the elements in coefficient matrices Zc and Ze depend only on refrigerant thermodynamic properties and geometric parameters, such as the internal volume, which is usually known. However, the elements in forcing functions fc and fe , which define the system statics, depend strongly on thermal powers transferred at each zone. Tracking references are imposed on the variables related to a specific zone when it is inactive. The dynamic models are simulated by calculating ˙xi(k) = Z−1 i(k)fi(k)∀i=c,e at each time step k and numerically integrating to obtain the state vector xi(k+1)∀i=c,e. The time-varying model inputs, which are provided by other component models, are gathered in Equation Set (2.15) . The outputs of each model depend on the interfaces with the other component models. wc= [ ˙mc,sec Tc,sec,in ˙mc,in ˙mc,out hc,in]T we= [ ˙me,sec Te,sec,in ˙me,in ˙me,out he,in]T(2.15) 38 Chapter 2. Modelling be calculated using exclusively the mode 2 state information, as presented in Equation Set (2.11). Once determined the boundary conditions at each instant, the valve expansion model (Equation Set (2.1) ), the compressor model (Equation Set (2.2) ), and the evaporator steadystate model (Equation Sets (2.3) – (2.5) ) constitute a nonlinear equation system with two unknown variables: Pe and he,out . An iterative procedure has been designed to solve the system, which is described in Figure 2.4. As a result, not only the unknown evaporator variables ( Pe and he,out ) are obtained, but also some other cycle variables are calculated and gathered in the output vector y , as shown in Figure 2.4 and Equation (2.28) . Specifically, the refrigerant mass flow ˙m and the condenser inlet specific enthalpy hc,in are included in y, therefore all information of the condenser input vector wcis now calculated. y=           hc,in he,in ˙ Wcomp ˙m TSH Te,sec,out ζe,t p           (2.28) 2.2.4 Comparison simulations Some simulation results are presented below comparing the original SMB model with the control-oriented one when some step changes on both manipulated variables ( N and Av ) are imposed. In Figure 2.5 a step change on N is applied, holding the expansion valve opening, whereas in Figure 2.6 a step change on Av is implemented. In this case both step changes are applied while the condenser works in mode 1. 2.2 Dynamic modelling 39 Condenser state information Manipulated variables Disturbances Supposing Expansion valve model Compressor model Evaporator steady-state model Pe Pe Pe Pe Pehe,out he,out he,out hc,out hc,out Pc Pc Pc Av Av N N Tsurr Tsurr ˙me,sec ˙me,sec Te,sec,in Te,sec,in he,in he,in ˙m ˙m ˙m h′ e,out h′ e,out hc,in ˙ Wcomp h′ e,in TSH Te,sec,out ζe,tp eit,he,in eit,he,out eit eit y Equation Set (2.1) Equation Set (2.2) Correction Equation Sets (2.3) – (2.5) Figure 2.4 Iterative procedure to solve evaporator variables. 40 Chapter 2. Modelling Time [min] 0 25 50 75 100 125 Compressor speed [Hz] 30 32.5 35 Manipulated variables 42 42.5 43 Valve opening [%] Time [min] 0 25 50 75 100 125 Mass flow [g s-1] 5.61 5.62 5.63 5.64 Refrigerant mass flow Original SMB, Compressor Original SMB, Valve Control-oriented SMB (a) Manipulated variables and refrigerant mass flow Time [min] 0 25 50 75 100 125 Pressure [bar] 17.14 17.16 17.18 17.2 Condenser pressure Original SMB Control-oriented SMB Time [min] 0 25 50 75 100 125 Pressure [bar] 1.06 1.1 1.14 1.18 Evaporator pressure Original SMB Control-oriented SMB (b) Condenser and evaporator pressures Figure 2.5 Comparison between the original SMB model and the control-oriented SMB model (mode 1) when applying a step change on the compressor speed N. 2.2 Dynamic modelling 41 Time [min] 0 25 50 75 100 125 Temperature [K] 0.52 0.56 0.6 0.64 Degree of subcooling at condenser outlet Original SMB Control-oriented SMB Time [min] 0 25 50 75 100 125 Temperature [K] 9 10 11 12 13 Degree of superheating at evaporator outlet Original SMB Control-oriented SMB (c) Degree of subcooling and superheating Time [min] 0 25 50 75 100 125 Temperature [ºC] -22.47 -22.46 -22.45 Outlet temperature of evaporator secondary flux Original SMB Control-oriented SMB Time [min] 0 25 50 75 100 125 COP 1.096 1.114 1.132 1.15 COP Original SMB Control-oriented SMB (d) Outlet temperature of evaporator secondary flux and Coefficient of Performance Figure 2.5 Comparison between the original SMB model and the control-oriented SMB model (mode 1) when applying a step change on the compressor speed N. 42 Chapter 2. Modelling Time [min] 0 25 50 75 100 125 Compressor speed [Hz] 29 30 31 Manipulated variables 42 46 50 54 Valve opening [%] Time [min] 0 25 50 75 100 125 Mass flow [g s-1] 5.6 6.1 6.6 7.1 7.6 Refrigerant mass flow Original SMB, Compressor Original SMB, Valve Control-oriented SMB (a) Manipulated variables and refrigerant mass flow Time [min] 0 25 50 75 100 125 Pressure [bar] 17.9 18.3 18.7 19.1 Condenser pressure Original SMB Control-oriented SMB Time [min] 0 25 50 75 100 125 Pressure [bar] 1.16 1.24 1.32 1.4 Evaporator pressure Original SMB Control-oriented SMB (b) Condenser and evaporator pressures Figure 2.6 Comparison between the original SMB model and the control-oriented SMB model (mode 1) when applying a step change on the valve opening Av. 2.2 Dynamic modelling 43 Time [min] 0 25 50 75 100 125 Temperature [K] 0 1 2 3 4Degree of subcooling at condenser outlet Original SMB Control-oriented SMB Time [min] 0 25 50 75 100 125 Temperature [K] 3 5 7 9 11 Degree of superheating at evaporator outlet Original SMB Control-oriented SMB (c) Degree of subcooling and superheating Time [min] 0 25 50 75 100 125 Temperature [ºC] -23.1 -22.9 -22.7 -22.5 -22.3 Outlet temperature of evaporator secondary flux Original SMB Control-oriented SMB Time [min] 0 25 50 75 100 125 COP 1.07 1.1 1.13 1.16 COP Original SMB Control-oriented SMB (d) Outlet temperature of evaporator secondary flux and Coefficient of Performance Figure 2.6 Comparison between the original SMB model and the control-oriented SMB model (mode 1) when applying a step change on the valve opening Av. 44 Chapter 2. Modelling As shown in Figures 2.5 and 2.6, minor differences are appreciated between the original SMB model and the control-oriented one. Only immediately after the step change, when there is not refrigerant mass flow equilibrium and the assumption of a unique refrigerant mass flow in the cycle does not exactly apply there are appreciable differences. As shown in Figures 2.5 and 2.6, the step response of all cycle variables, both in the original SMB and the control-oriented one, includes two dynamic components: one very fast, due to the mass flow imbalance and the evaporator intrinsic dynamics, and another one slower, due to the dominant condenser dynamics. As expected, since it does not consider the evaporator dynamics and the fast effects of mass flow imbalance, the control-oriented SMB model does not reproduce properly the fast response of the variables, especially those related to the evaporator, such as the degree of superheating and the output temperature of the evaporator secondary flux. However, the control-oriented SMB model does represent properly the slow dominant dynamics of the cycle. The transient differences cause the small steady-state errors, since it will be explained in Chapter 4 that the system has multiple equilibria given a couple of manipulated variables N and Av , and the steady-state equilibrium depends on the transient trajectory. Changes on the compressor speed cause greater mass flow imbalance, therefore the errors in Figure 2.5 are greater than in Figure 2.6. Anyway, the differences are small enough to be compensated by the controller robustness, and the computation time of the control-oriented model is about 2% of that of the original SMB model. Moreover, some simulation results when the condenser operates in mode 2 are presented. In Figure 2.7 a step change on N is applied, whereas in Figure 2.8 a step change on Av is implemented. In this case, since no degree of subcooling exists, the refrigerant vapour quality at the condenser outlet is represented, on account of that it is no longer zero in mode 2. The conclusions are similar to those derived when the condenser works in mode 1: there are more differences when the compressor speed changes, because the transient mass flow imbalance is greater, and only immediately after the step change there are noticeable differences. 2.3 Final remarks This Chapter has been devoted to modelling of vapour-compression refrigeration systems, focusing on the one-compression-stage, one-load-demand configuration. Firstly, steadystate modelling has been addressed and each element of the cycle has been separately modelled, paying special attention to the heat exchangers: evaporator and condenser. Regarding heat transfer, the effectiveness-NTU method has been used to model the different modes considered for each heat exchanger. Secondly, the dynamic modelling of such systems has been also addressed, also focusing on the heat exchangers, since they concentrate the dominant dynamics. Different approaches from the literature have been assessed and the switched moving boundary (SMB) model has been taken as starting point. However, the development of a controloriented simplified model seems interesting, since the still high order of the SMB model continues to limit its use in practical model-based control strategies. Some considerations have been taken into account to reduce the model order and thus the complexity of the original SMB model. Then, frequency differences between the dynamics of both 2.3 Final remarks 45 the evaporator and the condenser concerning their internal volumes and the refrigerant thermodynamic properties have been detailed and assessed. It has been concluded that the evaporator dynamics are usually much faster than those of the condenser, because of the lower internal volume of the first one, and also given that the refrigerant properties, especially the density, cause the evaporator dynamics to be intrinsically faster than those of the condenser. Moreover, the dynamics due to the refrigerant mass flow imbalance within the cycle considered in the original SMB model have been assessed and it has been concluded that the refrigerant mass flow equilibrium between the condenser inlet and outlet is very fast compared to the heat transfer dominant dynamics, therefore a unique refrigerant mass flow passing through the condenser might be assumed, disregarding the fast transient due to mass flow imbalance. A simplified control-oriented nonlinear model of the whole one-stage, one-load-demand cycle has been proposed, whose state vector is made up of only three variables, all related to the condenser, which turns out to generate the slowest and thus dominant dynamics. Differences between the control-oriented model and the original SMB model have been assessed in simulation, applying step changes on the manipulated variables. The differences show to be negligible regarding dominant dynamics, whereas the computation load is definitely reduced. All of this causes this simplified model to fit the requirements usually imposed when using model-based control strategies: reasonable accuracy and low computational load. The control-oriented dynamic model, which reduces the complexity of the original SMB model while retaining the dominant dynamics, is one of the main contributions of the Thesis. The steady-state models developed in Chapter 2 are to be used both in the parameter identification of the plant, performed in Chapter 3, and also when analysing the optimization stage of the model-based control strategy described in Chapter 4. Furthermore, the control-oriented dynamic model proposed in this Chapter is used in Chapter 4 to study the controllability of the one-stage, one-load-demand cycle, along with the control simulations included in Chapter 5. 46 Chapter 2. Modelling Time [min] 0 25 50 75 100 125 Compressor speed [Hz] 30 32.5 35 Manipulated variables 42 42.5 43 Valve opening [%] Time [min] 0 25 50 75 100 125 Mass flow [g s-1] 5 5.1 5.2 5.3 Refrigerant mass flow Original SMB, Compressor Original SMB, Valve Control-oriented SMB (a) Manipulated variables and refrigerant mass flow Time [min] 0 25 50 75 100 125 Pressure [bar] 16.55 16.6 16.65 Condenser pressure Original SMB Control-oriented SMB Time [min] 0 25 50 75 100 125 Pressure [bar] 1 1.07 1.14 Evaporator pressure Original SMB Control-oriented SMB (b) Condenser and evaporator pressures Figure 2.7 Comparison between the original SMB model and the control-oriented SMB model (mode 2) when applying a step change on the compressor speed N. 2.3 Final remarks 47 Time [min] 0 25 50 75 100 125 Vapour quality 0.005 0.015 0.025 Vapour quality at condenser outlet Original SMB Control-oriented SMB Time [min] 0 25 50 75 100 125 Temperature [K] 11.1 12.7 14.3 Degree of superheating at evaporator outlet Original SMB Control-oriented SMB (c) Vapour quality at condenser outlet and degree of superheating Time [min] 0 25 50 75 100 125 Temperature [ºC] -22.3 -22.28 -22.26 Output temperature of evaporator secondary flux Original SMB Control-oriented SMB Time [min] 0 25 50 75 100 125 COP 1.08 1.12 1.16 COP Original SMB Control-oriented SMB (d) Outlet temperature of evaporator secondary flux and Coefficient of Performance Figure 2.7 Comparison between the original SMB model and the control-oriented SMB model (mode 2) when applying a step change on the compressor speed N. 54 Chapter 3. Experimental plant Main compressor Booster compressor Air-cooled condenser Accumulator Evaporator 1 Pump 1 CV2 CV1 EEV1 Pressure upholding valve EEV2 Evaporator 2 Pump 2 TT PT PT PT TT TT TT PT TT PT TT PT TT TT TT TT TT FT FT Filter/Dryer Bubble viewer Liquid line Fan Water-cooled condenser Pump 3 TT FT CV3 CV4 Tank 3 Tank 1 R1 Tank 2 R2 PT Figure 3.1 Experimental plant layout. 3.1 Description of the experimental plant 55 Booster compressor Evaporator 1 Pump 1 CV2 EEV1 Pressure upholding valve Pump 2 TT PT PT PT TT TT TT PT TT PT TT PT TT TT TT TT TT FT FT Filter/Dryer Bubble viewer Liquid line Fan Water-cooled condenser Pump 3 TT FT CV4 Tank 3 Tank 1 R1 Tank 2 R2 PT Air-cooled condenser Accumulator Evaporator 2 CV1 CV3 Main compressor EEV2 Figure 3.2 One-stage, one-load-demand configuration highlighted on the experimental plant layout. 56 Chapter 3. Experimental plant Booster compressor Pump 1 CV2 Pump 2 TT PT PT PT TT TT TT PT TT PT TT PT TT TT TT TT TT FT FT Filter/Dryer Bubble viewer Liquid line Fan Water-cooled condenser Pump 3 TT FT CV4 Tank 3 Tank 1 R1 Tank 2 R2 PT Air-cooled condenser Accumulator Evaporator 2 CV1 CV3 Main compressor Evaporator 1 Pressure upholding valve EEV1 EEV2 Figure 3.3 One-stage, two-load-demand configuration highlighted on the experimental plant layout. 3.1 Description of the experimental plant 57 Booster compressor Pump 1 Pump 2 TT PT PT PT TT TT TT PT TT PT TT PT TT TT TT TT TT FT FT Filter/Dryer Bubble viewer Liquid line Fan Water-cooled condenser Pump 3 TT FT CV4 Tank 3 Tank 1 R1 Tank 2 R2 PT Air-cooled condenser Accumulator Evaporator 2 CV1 CV3 Main compressor Evaporator 1 Pressure upholding valve EEV1 EEV2 CV2 Figure 3.4 Two-stage, two-load-demand configuration highlighted on the experimental plant layout. 58 Chapter 3. Experimental plant Main compressor Booster compressor Tank 2 Tank 1 Pump 2 Pump 1 Pump 3 Watercooled condenser Accumulator Air-cooled condenser Figure 3.5 Compressors placed at the experimental facility. opening to hold a specified set point on the degree of superheating. In order to perform the calculation of the degree of superheating, each SVC has a pressure transducer and a thermocouple analogue input. Evaporators, expansion valves, and secondary flux circuits are represented in Figure 3.6. Finally, the accumulator depicted in Figure 3.1 is only used to store the refrigerant when the plant is not working; in no case it is used in operation as a device which ensures that the refrigerant is saturated liquid at the expansion valve inlet, as other experimental plants do, for instance that described by Rasmussen et al. [1]. 3.1.2 Sensors The experimental plant has full instrumentation for research purposes. Thermocouples and pressure transducers have been placed at selected points to calculate all the refrigerant enthalpies and the degree of superheating at both evaporator outlets, whereas the brine volumetric flow rates are measured by magnetic-inductive flow sensors. All sensors are indicated in Figure 3.1, following standard nomenclature, and some of them are also visible in Figures 3.5 and 3.6. Nevertheless, there is no refrigerant mass flow sensor at the facility. This is a common issue at commercial refrigeration systems, since such sensors are too expensive and the measurement accuracy usually shows not to be as high as required. A summary of the experimental setup of apparatus is presented in Table 3.1. 3.1 Description of the experimental plant 59 Evaporator 2 Evaporator 1 Tank 1 Tank 2 Pump 2 Pump 1 EEV2 EEV1 R2 R1 Filter/dryer Figure 3.6 Evaporators, expansion valves, and secondary flux circuits. 3.1.3 Programmable Logic Controller (PLC) A PLC is used for low-level control of the refrigeration plant. Along with the power module, a digital input module and a digital output module are used to manage digital signals. Furthermore, three analogue input/output mixed modules are set up, along with two specific thermocouple signal converters. In Appendix A the signals connected through the analogue mixed modules and the thermocouple signals wired through the converters are detailed. Moreover, all digital and analogue signals from VFDs and SVCs are connected to the PLC through a communication bus over RS-485. The ModBus protocol [ 76 ] is selected due to its simplicity and reliability. The PLC is the master of the communication, whereas all VFDs and SVCs are slaves. The bus access method is time division, setting a maximum slot time of 100 ms. Table 3.2 summarises the total number of variables managed by the PLC. The PLC also manages the control sample time, which basically depends on the communication bus cycle time. The data flow from and to the VFDs is always imperative, but the same is not suitable for the SVCs. If the EEVs are controlled through their SVCs, all information related to these controllers (which is read and written through the bus) is relevant. However, if the expansion valves are controlled straightly by the PLC through digital signals, most reading and writing operations concerning the SVCs are not required. This fact allows to reduce the bus cycle time, since the number of reading and writing 60 Chapter 3. Experimental plant Table 3.1 Feature summary of the experimental plant. Main Compressor Booster Compressor Power [kW] 6.4 Power [kW] 1.8 Number of cylinders 4 Number of cylinders 2 Speed range [Hz] [25 70] Speed range [Hz] [25 70] Variable frequency drive (VFD) Yes Variable frequency drive (VFD) Yes Air-cooled condenser Water-cooled condenser Type Cross-flow tubes Type Brazed-plate Secondary flux Air Secondary flux Water Current direction Cross flow Current direction Counter-current Fan power [kW] 1.05 Pump power [kW] 0.25 Evaporator 1 Evaporator 2 Type Brazed-plate Type Brazed-plate Secondary flux Brine of 28% mono-ethylene glycol aqueous solution Secondary flux Brine of 60% propylene glycol aqueous solution Current direction Counter-current Current direction Counter-current Resistance power [kW] 2 Resistance power [kW] 2 Pump power [kW] 0.25 Pump power [kW] 0.25 Electronic expansion valve 1 (EEV1) Electronic expansion valve 2 (EEV2) Method Pulse-width modulation (PWM) Method Pulse-width modulation (PWM) Thermocouples Pressure Transducers Range [ºC] [-40 180] High pressure range [bar] [0 30] Precision [ºC] 0.1 Low pressure range [bar] [0.5 7.5] Volumetric flow sensor Precision [bar] 0.01 Range [L min-1][0.1 25] Refrigerant Precision [L min-1]0.05 R404a operations is lower. Moreover, due to the unreliability of the SVCs regarding the bus data flow, it is possible to reduce the slot time up to 50 ms. Table 3.3 shows the total number of reading and writing operations per bus cycle and the timing for both valve control choices. Even once selected a valve control mode, there are some data that the PLC must read imperatively once a bus cycle, whereas some others are not so relevant and they can be read at lower frequency. Taking this into account, the bus cycle is organised in such a way that all the mandatory required data are read once a cycle and the rest of them are read by using a shift scheme. A detailed description of the reading and writing operations carried out in both valve control modes as well as the shift scheme used in each case are provided in Appendix A. 3.1 Description of the experimental plant 61 Table 3.2 Variables managed by the PLC. Type Number Digital inputs 15 Digital outputs 8 Analogue inputs 12 Analogue outputs 3 Thermocouple inputs 8 ModBus variables 60 Table 3.3 Bus timing for the two valve control modes. Valve control mode Reading and writing operations Slot time [ms] Bus cycle time [s] Through SVCs 30 100 3 Straight PLC control 20 50 1 3.1.4 Supervision software and OPC communication A Supervisory Control And Data Acquisition (SCADA) system is developed and works as human-machine interface. It allows to control the plant manually in open loop and monitor all cycle variables when a high-level controller is applied. Figure 3.7 shows the SCADA screen when a high-level controller, namely through OLE for Process Control (OPC) communication standard, is performed. OPC is a communication standard used to connect devices from different manufacturers [ 77 ]. Here, OPC is used for communication between the PLC and MATLAB ® software [ 78 ]. As explained previously, the PLC is used to data acquisition and low-level control, whereas the computing program is devoted to implementing high-level controllers. As seen in Figure 3.8, the PLC is connected to the computer by an Ethernet link. The SCADA system and the PLC are directly communicated as long as they are from the same manufacturer. MATLAB ® , however, needs to use an OPC server to access process variables. OPC is set up in MATLAB ® /Simulink environment where the number of read and written data is gathered in Table 3.4. 62 Chapter 3. Experimental plant Figure 3.7 SCADA screen under OPC control. PlantPLC SCADA OPC server Matlab® Simulink Ethernet Computer Figure 3.8 Communication diagram. Table 3.4 Process variables managed through OPC. Analogue Digital Read data 50 5 Written data 13 - 3.2 Identification 63 3.2 Identification As stated in Chapter 2, all elements of the cycle are separately modelled, creating dynamic submodels for the refrigerant flow along the heat exchangers and static models for the compressors, expansion valves, and secondary fluxes at the evaporators and condenser. The identification of all components according to experimental data is addressed below. 3.2.1 Expansion valves The static model considered for both expansion valves has been previously presented in Chapter 2, Subsection 2.1.1. Notwithstanding, the mass flow characteristic of the valve is included in Equation (3.1) as a remainder, once the notation has been adapted for the multi-load-demand case. ˙me=ceev Avq2ρc,out ∆Peev =ceevAvq2ρc,out (Pc,out −Pe)(3.1) The coefficient ceev is a characteristic parameter of the valve, which allows to predict the refrigerant mass flow when varying the valve opening and the cycle variables, such as pressures and enthalpies. Test data and mass flow models of EEVs are very limited in open literature, as stated by Park et al. [ 79 ]. The simplest model consisting of a unique value of the valve coefficient ceev , despite being in general valid for single-phase flows, produces non-negligible deviations when considering several operating points, because more complex phenomena occur when the refrigerant expands through the valve and becomes two-phase. The original works by Schurt et al. proposed a correlation for ceev where its value depends only on the degree of subcooling [ 11 , 33 ]. Considering the facility under study, the component selection and performance range of the manipulated variables allow only to achieve cycles with very low degree of subcooling, even zero in some cases. It causes this correlation to produce very similar values of ceev , and it is not valid when TSC reaches zero. A sensitivity analysis has been carried out and the conclusions have shown that inaccuracy when estimating the refrigerant mass flow through the valve model involves great differences between the experimental and simulated variables when validating the whole cycle. Therefore, a more complex model is better suited than that proposed by Schurt et al. Then, as proposed by Shanwei et al. [ 80 ], a polynomial correlation for ceev is considered, where its dependence on all input parameters of the expansion valve model is highlighted in Equation Set (3.2) . The degree of subcooling TSC is also included as recommended by Shanwei et al. [ 80 ], since that dependence becomes relevant when analysing experimental data with very low degree of subcooling. ceev =ceev(Av,∆Peev,ρc,out ,TSC) ceev =10−9 n0+n1 Av 100 +n2 ∆Peev Pc,out +n3 ρc,out ρc,f,out +n4 TSC Tc,sat,out !(3.2) In Equation Set (3.2) the temperature Tc,sat,out refers to the saturation temperature at the outlet condenser pressure Pc,out , whereas the density ρc,f,out matches the density of saturated liquid at the outlet condenser pressure. An identification method to calculate 70 Chapter 3. Experimental plant Superheated vapour zone As stated in Equation Set (3.12) , energy balance in steady state at the superheated vapour zone involves that ˙ Qe,sh =˙me(he,out −he,g) . The pressure Pe and the outlet temperature Te,out are measured, whereas ˙me is estimated applying the global energy balance indicated in Equation (3.3). Once calculated ˙ Qe,sh , Equation (3.15) shows how to estimate the efficiency at the superheated vapour zone εe,sh from available data. εe,sh = ˙ Qe,sh (˙mcp)min(Te,sec,in −Te)(3.15) If the heat-exchanger specific curve ε=ghe(NTU,C) is inversely applied, the values of NTUe,sh and then UAe,sh can be calculated, which allows to achieve an experimental relationship between αe,sh and ζe,t p . The latter is indicated in Equation (3.16) and depicted in Figure 3.10. Although many steady-state points are considered for identification, only one representative curve corresponding to a single steady-state experimental point is shown in Figure 3.10, since its purpose is only to describe qualitatively the relationship between αe,sh and ζe,t p . For the sake of brevity, only figures related to Evaporator 1 are henceforth shown, since those corresponding to Evaporator 2 are qualitatively equivalent due to the same reason given above. αe,sh =UAe,sh Ae,trns f Le 1 ζe,sh =UAe,sh Ae,trns f Le 1 1−ζe,t p (3.16) As seen in Figure 3.10, the greater ζe,t p is, the greater αe,sh is, since the superheated vapour zone is smaller and the overall coefficient must be as high as necessary to effectively transfer the calculated thermal power ˙ Qe,sh. Two-phase zone Energy balance at the two-phase zone implies that ˙ Qe,t p =˙me(he,g−he,in) . Considering isenthalpic expansion at the corresponding valve, he,in =hc,out can be measured, thus the value of ˙ Qe,t p is obtained. The efficiency at the two-phase zone εe,t p can be calculated as indicated in Equation (3.17) , where Te,sec,mid is previously calculated applying secondary flux energy balance at the superheated vapour zone. Te,sec,mid =Te,sec,in − ˙ Qe,sh ˙me,sec cp,e,sec εe,t p = ˙ Qe,t p (˙mcp)min(Te,sec,mid −Te) (3.17) Once again applying inversely the curve ε=ghe(NTU,C) , the values of NTUe,t p and UAe,t p can be computed. The experimental relationship between αe,t p and ζe,t p shown 3.2 Identification 71 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Two-phase zone normalised length 10 30 50 70 90 110 130 150 170 Heat transfer coefficient [W m-2 K-1] Qualitative relationship between αe5,sh and ζe5,tp Figure 3.10 Qualitative relationship between αe,sh and ζe,t p at Evaporator 1 (at 5 °C). in Equation (3.18) is achieved, which is also represented in Figure 3.11. Once again the purpose of this figure is only to show qualitatively the relationship between αe,t p and ζe,t p . αe,t p =UAe,t p Ae,trns f Le 1 ζe,t p (3.18) As expected, if the two-phase zone is greater, αe,t p must be smaller to effectively transfer the measured thermal power ˙ Qe,t p, as indicated in Figure 3.11. Hypothesis between αe,sh and αe,tp Combining Equations (3.16) and (3.18) , the ratio between αe,t p and αe,sh is a function of ζe,t p , for each experimental point. This ratio is represented in Figure 3.12 for a representative experimental point as depicted in Figures 3.10 and 3.11. The values close to limits ( ζe,t p =1 and ζe,t p =0 ) are not represented because it would entail considering no superheated vapour zone and no two-phase zone, respectively. As previously stated and represented by means of an electric analogy in Figure 3.9, the overall heat transfer coefficients αe,t p and αe,sh must be calculated considering a serial 72 Chapter 3. Experimental plant 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Two-phase zone normalised length 0 500 1000 1500 2000 2500 3000 3500 Heat transfer coefficient [W m-2 K-1] Qualitative relationship between αe5,tp and ζe5,tp Figure 3.11 Qualitative relationship between αe,t p and ζe,t p at Evaporator 1 (at 5 °C). combination of all thermal resistances, as shown in Equation Set (3.19). 1 αe,t p =1 λe,t p +Re,cond +1 λe,sec 1 αe,sh =1 λe,sh +Re,cond +1 λe,sec (3.19) Reference values or at least orders of magnitude of convective heat transfer coefficients in the literature [ 70 ] are collected in Table 3.11, as well as the conduction thermal resistance, which has been estimated considering the thickness and conductivity of the plate material. From these reference data on, the estimated ratio between αe,t p and αe,sh is in the order of 10. Considering this value and using Figure 3.12, a probable range for the length of the twophase zone ζe,t p can be obtained. Ranges of ζe5,t p ∈[0.58,0.7] and ζe20,t p ∈[0.81,0.93] are achieved. Considering such plausible ranges, it is possible to calculate average values of the overall heat transfer coefficients αe,sh and αe,t p using experimental data such as those depicted in Figures 3.10 and 3.11. The average values and standard deviations of all coefficients are gathered in Table 3.12, considering both evaporators. Note that the range for ζe20,t p causes the length of the superheated vapour zone ζe20,sh to be small. It causes the identification of αe20,sh to be more difficult and numerically sensitive than those of αe5,t p and αe5,sh : that is why the standard deviation of αe20,sh is 3.2 Identification 73 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 Two-phase zone normalised length 0 25 50 75 100 125 150 175 200 Heat transfer coefficient ratio Qualitative relationship between the ratio αe5,tp/αe5,sh and ζe5,tp Figure 3.12 Qualitative relationship between the ratio αe,t p/αe,sh and ζe,t p at Evaporator 1 (at 5 °C). Table 3.11 Orders of magnitude of the thermal resistances involved in heat transfer at the evaporator. Heat transfer coefficients and thermal resistances Two-phase (tp) Superheated vapour (sh) λe[W m-2 K-1]∼105∼5·101 Re,cond [m2K W-1]2.05 ·10−6 λe,sec [W m-2 K-1]∼5·102 αe[W m-2 K-1]∼499 ∼45.4 much greater than the others shown in Table 3.12. The main reason is the asymptotic feature of the relationship between αe,sh and ζe,t p , when the latter tends to a value close to one, as shown in Figure 3.10. Furthermore, the highly nonlinear features of the refrigeration system are known to cause the coefficients to vary for different operating points. Indeed, correlations usually produce different values for heat transfer coefficients when the flow conditions vary and specially when the temperature difference between the refrigerant and the secondary flux 74 Chapter 3. Experimental plant Table 3.12 Identified overall heat transfer coefficients at both evaporators. Evaporator Heat transfer coefficient Mean value [W m-2 K-1] Standard deviation [W m-2 K-1] Evaporator 1 b αe5,t p 520.5 40.11 b αe5,sh 51.9 8.85 Evaporator 2 b αe20,t p 144.8 13.8 b αe20,sh 14.87 7.26 is modified. However, the identification methodology has been applied for a wide (as wide as possible considering the features of the experimental facility concerning the component selection and the performance range of the manipulated variables) set of steady-state operating points. It causes the standard deviations of the heat transfer coefficients to be greater than wished. If a more accurate estimation of the heat transfer coefficients is desired, considering diverse values for a range of operating conditions, the identification method could be applied online using filtered measurable variables, such as the degree of superheating, pressures, and inlet and outlet temperatures of the heat exchangers, to estimate the working conditions and calculate the corresponding heat transfer coefficients, in addition to an estimation of the heat exchanging zone lengths. 3.2.4 Condenser Similarly, the MB approach is also applied when modelling the condenser. In this case the condenser length is divided into three zones, since three refrigerant states may simultaneously exist along the condenser: superheated vapour, two-phase fluid, and subcooled liquid. It has been studied in Chapter 2 that the system dynamics of a heat exchanger depend only on refrigerant thermodynamic properties and the heat exchanger internal volume, which can be easily calculated or measured. However, the partial thermal powers ˙ Qc,sh , ˙ Qc,t p , and ˙ Qc,sc , related to the three possible zones, define system statics, and their calculation does depend to a great extent on the specific configuration of the heat exchanger and its ability to transfer heat in all zones. Thus, a steady-state model of the condenser is also considered for identification purposes. As commented in Chapter 2, the identification procedure is focused on the air-cooled condenser, which has a cross-flow tube configuration. The equations which detail heat transfer at the condenser were in depth analysed in Chapter 2, Subsection 2.1.4, but they are included here in Equation Sets (3.20) – (3.24) to improve the comprehensibility of the identification procedure. Note that only condenser mode 1 is considered and pressure drop along the heat exchanger is also modelled, unlike the ideal model described in Subsection 2.1.4. Then, Equation Set (3.20) describes heat transfer at the superheated vapour zone, Equation Set (3.21) does the same for the two-phase zone, and finally Equation Set (3.22) refers 3.2 Identification 75 to the subcooled liquid zone. ˙ Qc,sh =εc,sh(˙mcp)min(Tc,in −Tc,sec,in) Cc,sh =(˙mcp)min (˙mcp)max NTUc,sh =UAc,sh (˙mcp)min UAc,sh =αc,sh ζc,sh Ac,trns f Lc εc,sh =ghe(NTUc,sh,Cc,sh) ˙ Qc,sh =˙mc(hc,in −hc,g) Pc,g=Pc,in −ζc,sh kdrop ˙m2 c (3.20) ˙ Qc,t p =εc,t p(˙mcp)min(Tc,g−Tc,sec,in) Cc,t p =(˙mcp)min (˙mcp)max ={phase change}=0 NTUc,t p =UAc,t p (˙mcp)min UAc,t p =αc,t p ζc,t p Ac,trns f Lc εc,t p =ghe(NTUc,t p,Cc,t p) ˙ Qc,t p =˙mc(hc,g−hc,f) Pc,f=Pc,g−ζc,t p kdrop ˙m2 c (3.21) ˙ Qc,sc =εc,sc(˙mcp)min(Tc,f−Tc,sec,in) Cc,sc =(˙mcp)min (˙mcp)max NTUc,sc =UAc,sc (˙mcp)min UAc,sc =αc,sc ζc,sc Ac,trns f Lc εc,sc =ghe(NTUc,sc,Cc,sc) ˙ Qc,sc =˙mc(hc,f−hc,out ) (3.22) Eventually, Equation Set (3.23) indicates how to calculate the refrigerant state at the condenser outlet, while in Equation Set (3.24) the secondary flux outlet temperature model 76 Chapter 3. Experimental plant is detailed. Pc,out =Pc,f−ζc,sc kdrop ˙m2 c Tc,out =T(Pc,out ,hc,out ) Tc,sat,out =T(Pc,out ,q=0) TSC =Tc,sat,out −Tc,out (3.23) Tc,sec,sh,out =Tc,sec,in + ˙ Qc,sh ˙mc,sec ζc,sh cp,c,sec Tc,sec,t p,out =Tc,sec,in + ˙ Qc,t p ˙mc,sec ζc,t p cp,c,sec Tc,sec,sc,out =Tc,sec,in + ˙ Qc,sc ˙mc,sec ζc,sc cp,c,sec Tc,sec,out =ζc,sh Tc,sec,sh,out +ζc,t p Tc,sec,t p,out +ζc,sc Tc,sec,sc,out (3.24) Once again the identification objective is to obtain average estimations of all overall heat transfer coefficients αc,sh , αc,t p , and αc,sc . In this case, to calculate the heat transfer coefficients it is necessary to know the corresponding zone lengths, namely ζc,sh and ζc,t p , but their measurements are not available from experimental data. Similarly, according to the dominance of heat transfer when the fluid is two-phase, it is guessed that αc,t p is higher than αc,sh and αc,sc [ 70 , 71 ]. Nevertheless, there is no previous assumption of the relationship between single-phase heat transfer coefficients. Moreover, just like at the evaporators, other thermal resistances must be taken into account, since the coefficients to be estimated are global according to the proposed formulation. A diagram which illustrates this issue by applying an electric analogy is shown in Figure 3.13. Superheated vapour zone As indicated in Equation Set (3.20) , energy balance in steady state at the superheated vapour zone involves that ˙ Qc,sh =˙mc(hc,in −hc,g) . The inlet pressure Pc,in and the inlet temperature Tc,in are measured, whereas the refrigerant mass flow ˙mc is estimated applying the global energy balance indicated in Equation (3.3) . Nevertheless, the refrigerant pressure distribution cannot be calculated without zone length estimations, as expressed in Equation Set (3.25), where kdrop concerns the pressure drop factor along the whole condenser. Pc,g=Pc,in −ζc,sh kdrop ˙m2 c Pc,f=Pc,g−ζc,t p kdrop ˙m2 c Pc,out =Pc,f−ζc,sc kdrop ˙m2 c (3.25) Linear modelling of pressure drop involves that partial pressure drops at each zone depend on zone lengths. Therefore, Pc,g and, as a result, hc,g and ˙ Qc,sh , cannot be computed without supposing ζc,sh. 3.2 Identification 77 SuperheatedTwo-phase Tube Secondary flux Subcooled Secondary flux Tube ˙mc,out hc,out Tc,sec,sc,out Tc,sec,tp,out Tc,sec,sh,out Tc,g Pc,g ˙mc,in hc,in Tc,f Pc,f ˙mc,sec ˙mc,sec Tc,sec,in Tc,sec,out 1/ λ c,sec 1/ λ c,sec 1/ λ c,sec 1/ λ c,sh 1/ λ c,tp 1/ λ c,sc Rc,cond Rc,cond Rc,cond Figure 3.13 Electric analogy of heat transfer at the condenser. For a certain value of ζc,sh , according to Equation Set (3.20) , the efficiency at the superheated vapour zone εc,sh is computed from the calculated value of ˙ Qc,sh as shown in Equation (3.26). εc,sh = ˙ Qc,sh (˙mcp)min(Tc,in −Tc,sec,in)(3.26) To calculate (˙mcp)min it is also necessary to suppose ζc,sh , since the secondary flux mass flow ˙mc,sec is supposed to be distributed among the zones according to their lengths. Just like at the evaporator, assuming that the heat-exchanger specific curve ε=ghe(NTU,C) is inversely applied, the values of NTUc,sh and then UAc,sh can be calculated. An experimental relationship between αc,sh and ζc,sh is obtained, which is indicated in Equation (3.27) and depicted in Figure 3.14. Many steady-state points are considered for identification, but only one curve corresponding to a single representative steady-state experimental point is shown in Figure 3.14 in order to represent qualitatively the relationship between αc,sh and ζc,sh. αc,sh =UAc,sh(ζc,sh) Ac,trns f Lc 1 ζc,sh (3.27) The dependence of UAc,sh on ζc,sh is explicitly indicated in Equation (3.27) . As expected, greater values of ζc,sh involve smaller values of αc,sh , since the superheated vapour zone is longer. 78 Chapter 3. Experimental plant 0.24 0.29 0.34 0.39 0.44 Superheated vapour zone normalised length 10 12 14 16 18 20 22 24 Heat transfer coefficient [W m-2 K-1] Qualitative relationship between αc,sh and ζc,sh Figure 3.14 Qualitative relationship between αc,sh and ζc,sh. Two-phase zone Energy balance at the two-phase zone involves that ˙ Qc,t p =˙mc(hc,g−hc,f) . Once estimated ˙mc , pressures Pc,f and Pc,g allow to completely determine ˙ Qc,t p . As indicated in Equation Set (3.25) , Pc,f depends explicitly on ζc,t p , whereas it is also influenced by ζc,sh through Pc,g. Therefore, estimations of both ζc,sh and ζc,t p are required to compute ˙ Qc,t p. Given reasonable values of ζc,sh and ζc,t p , the efficiency at the two-phase zone εc,t p can be calculated as shown in Equation (3.28), according to Equation Set (3.21). εc,t p = ˙ Qc,t p (˙mcp)min Tc,g−Tc,sec,in(3.28) By applying the specific curve ε=ghe(NTU,C) at the two-phase zone, the values of NTUc,t p and UAc,t p can be calculated. The experimental relationship between αc,t p , ζc,sh , and ζc,t p indicated in Equation (3.29) is obtained. αc,t p =UAc,t p(ζc,sh,ζc,t p) Ac,trns f Lc 1 ζc,t p (3.29) Figure 3.15 illustrates the double dependence of αc,t p . Once again, although many steadystate points are considered for identification, curves relative to a single experimental point are shown in Figure 3.15, since its purpose is only to represent qualitatively the relationship 3.2 Identification 79 between αc,t p , ζc,t p , and ζc,sh . On the one hand, in Figure 3.15(a) the experimental relationship between αc,t p and ζc,t p is represented, for fixed ζc,sh . The greater ζc,t p is, the smaller αc,t p is, since the two-phase zone is longer. On the other hand, Figure 3.15(b) shows the experimental relationship between αc,t p and ζc,sh , for fixed ζc,t p . This dependence is minor, since for fixed ζc,t p the calculation of αc,t p is only affected by ζc,sh through pressure Pc,g . Supposing a certain ζc,t p , if the superheated vapour zone is longer, it involves a greater pressure drop at this zone and thus Pc,g is lower. As stated in Equation Set (3.21) , the hottest point related to heat transfer at the two-phase zone is Tc,g , which is lower when Pc,g decreases. It entails lower temperature difference between the refrigerant and the secondary flux, and that is the reason why αc,t p is slightly greater in order to actually transfer ˙ Qc,t p. Subcooled liquid zone Energy balance at the subcooled liquid zone imposes that ˙ Qc,sc =˙mc(hc,f−hc,out ) . Despite measuring the outlet pressure Pc,out and the outlet temperature Tc,out , and estimating ˙mc , ˙ Qc,sc cannot be calculated without estimations of ζc,t p and ζc,sh , since Pc,f depends on both zone lengths. Supposing plausible values of ζc,sh and ζc,t p , the efficiency at the subcooled liquid zone εc,sc can be computed as shown in Equation (3.30), according to Equation Set (3.22). εc,sc = ˙ Qc,sc (˙mcp)min(Tc,f−Tc,sec,in)(3.30) By applying inversely the curve ε=ghe(NTU,C) , the values of NTUc,sc and then UAc,sc can be computed. The experimental relationship between αc,sc , ζc,sh , and ζc,t p indicated in Equation (3.31) and represented in Figure 3.16 arises. Similarly, for the sake of clarity, only the curves related to a single experimental steady-state point are shown in Figure 3.16. αc,sc =UAc,sc(ζc,sh,ζc,t p) Ac,trns f Lc 1 ζc,sc =UAc,sc(ζc,sh,ζc,t p) Ac,trns f Lc 1 1−ζc,sh −ζc,t p (3.31) Firstly, in Figure 3.16(a) the experimental relationship between αc,sc and ζc,t p is depicted, for fixed ζc,sh . As expected, greater values of ζc,t p involves shorter subcooled liquid zone, which causes αc,sc to be greater. Secondly, Figure 3.16(b) shows the experimental relationship between αc,sc and ζc,sh , for fixed ζc,t p ; the effect of ζc,sh on αc,sc is shown to be similar. An asymptotic behaviour is observed, since when ζc,sh or ζc,t p increases, the subcooled liquid zone length ζc,sc tends ultimately to zero. Hypothesis between αc,sh,αc,tp, and αc,sc According to Equations (3.27) and (3.29) , the ratio between αc,t p and αc,sh turns out to be a function of ζc,t p and ζc,sh , as represented in Figure 3.17(a). Similarly, the ratio between αc,sh and αc,sc is a function of the same variables, as shown in Equations (3.27) and (3.31) . This ratio is also shown in Figure 3.17(b). Note that the depicted surface refers to a single representative steady-state experimental point. 86 Chapter 3. Experimental plant Cooling chamber Evaporator 1 Condenser Main compressor Expansion valve 1 Secondary flux Secondary flux Freezing chamber Evaporator 2 Expansion valve 2 Secondary flux Pressure upholding valve Figure 3.19 One-stage, two-load-demand refrigeration cycle. assumed to be isenthalpic, whereas mass and energy balances are posed at the compressor intake, since at this point two refrigerant flows with different mass flows and specific enthalpies are mixed together: that circulating through the Evaporator 2 at -20 ° C, and the other expanded through the pressure upholding valve. In Figure 3.20 an ideal P-h diagram of such a cycle is depicted (without considering the pressure drop at the condenser), along with the refrigerant-specific saturation curves, where the function of such valve is clarified. In Figure 3.20 temperature levels of all secondary fluxes are also represented: • Freezing level: represents the refrigerant saturation pressure at the inlet temperature of the secondary flux at Evaporator 2 (-20 °C). • Cooling level: refers to the refrigerant saturation pressure at the inlet temperature of the secondary flux at Evaporator 1 (5 °C). • Ambient: refers to the refrigerant saturation pressure at the inlet temperature of the condenser secondary flux, which matches ambient temperature when using the air-cooled condenser. They have been included in the P-h diagram only to represent qualitatively the sign of the temperature difference between the refrigerant and the secondary flux at both heat exchangers. Since all components have been identified in steady state, the global validation is carried out considering equilibrium points. The pressure upholding valve has been identified using 3.3 Validation 87 Evaporator 1 Evaporator 2 Compressor Condenser EEV1 EEV2 Pressure upholding valve Specific enthalpy Pressure P−h diagram R404a saturated liquid R404a saturated vapour Cycle Freezing level Ambient Cooling level Te5,SH Te20,SH Pe20 Pe5 hc,out Pc he5,out hc,in TSC he20,out Figure 3.20 P-h diagram of a one-stage, two-load-demand cycle. the same methodology explained in Subsection 3.2.1 for the electronic expansion valves, giving rise to the model shown in Equation (3.34). ˙me5,out =cpuv q2ρe5,out (Pe5−Pe20)(3.34) No valve opening factor is considered in Equation (3.34) , since the manual open setting is not changed during the experiments and they are all carried out with the same open setting. This information is implicit in parameter cpuv , whose estimated value and standard deviation are stated in Equation (3.35) . In this case a unique value of the coefficient cpuv is valid to describe the mass flow characteristic of this valve, since the refrigerant at the pressure upholding valve is not a two-phase fluid, but superheated vapour, as appreciated in the P-h diagram shown in Figure 3.20. b cpuv =2.86 ·10−6m2σcpuv =5.5·10−8m2(3.35) In Figure 3.21 some experimental and simulated P-h diagrams are overlaid at various operating points. Table 3.15 gathers some relative errors of the estimation on measurable variables for different operating points. Note that in Figure 3.21 pressure drop at the condenser is indeed represented both in the experimental cycle and in the simulated one. Experimental data show that pressure drop at the evaporators is negligible, which is related to heat exchanger design, since plate heat exchangers usually cause lower pressure drop than tube ones. It can be observed in 88 Chapter 3. Experimental plant 225 275 325 375 425 475 1 2 5 10 20 30 40 Specific enthalpy [kJ kg−1] Pressure [bar] P−h diagram Condenser EEV1 EEV2 Evaporator 1 Evaporator 2 Compressor Experimental cycle Simulated cycle Freezing level Cooling level Ambient (a) Operating point 1 (OP1) 225 275 325 375 425 475 1 2 5 10 20 30 40 Specific enthalpy [kJ kg−1] Pressure [bar] P−h diagram Condenser EEV1 EEV2 Evaporator 2 Evaporator 1 Compressor Experimental cycle Simulated cycle Freezing level Cooling level Ambient (b) Operating point 2 (OP2) Figure 3.21 Experimental and simulated P-h diagrams at different operating points considering the one-stage, two-load-demand configuration. 3.3 Validation 89 225 275 325 375 425 475 1 2 5 10 20 30 40 Specific enthalpy [kJ kg−1] Pressure [bar] P−h diagram Condenser EEV1 EEV2 Evaporator 1 Evaporator 2 Compressor Experimental cycle Simulated cycle Freezing level Cooling level Ambient (c) Operating point 3 (OP3) 225 275 325 375 425 475 1 2 5 10 20 30 40 Specific enthalpy [kJ kg−1] Pressure [bar] P−h diagram Condenser EEV1 EEV2 Evaporator 1 Evaporator 2 Compressor Experimental cycle Simulated cycle Freezing level Cooling level Ambient (d) Operating point 4 (OP4) Figure 3.21 Experimental and simulated P-h diagrams at different operating points considering the one-stage, two-load-demand configuration. 90 Chapter 3. Experimental plant Table 3.15 Relative errors [%] on diverse measurable variables in the one-stage, two-load- demand cycle. Variable Operating point OP1OP2OP3OP4 Pe20 4.71 2.80 3.11 0.22 Pe515.13 9.31 18.39 12.21 Pc,in 10.06 8.77 12.35 10.30 he20,out 1.69 1.94 1.14 0.45 he5,out 0.40 0.76 1.06 0.64 hc,out 1.70 3.91 3.85 2.55 Figure 3.21 and Table 3.15 that the identified overall heat transfer coefficients at all heat exchangers, along with the parameters of the other elements of the cycle, allow the whole steady state cycle to be correctly estimated. Although some errors could be seen as too high, the identified steady-state model is intended to be used within model-based control strategies, and it is expected that the robustness of the controllers will be able to deal with such modelling errors. 3.3 Validation 91 3.3.2 Two-compression-stage, two-load-demand configuration This configuration was illustrated in Figure 1.5 and highlighted on the experimental plant layout in Figure 3.4. Once again for the sake of clarity the layout is replicated in Figure 3.22. Expansion valve 1 Secondary flux Expansion valve 2 Secondary flux Figure 3.22 Two-compression-stage, two-load-demand refrigeration cycle. In this configuration the booster compressor is responsible for increasing the refrigerant pressure from Evaporator 1 to Evaporator 2. Thus, the pressure upholding valve located at Evaporator 1 outlet is fully opened in this configuration, since the pressure difference between both evaporators is held by the booster compressor. A typical P-h diagram of the two-stage, two-load-demand cycle is represented in Figure 3.23, as well as the R404a saturation curves. Once again for the sake of simplicity, no pressure drop at the condenser is represented. Some experimental and simulated P-h diagrams are plotted together in Figure 3.24, considering different steady-state operating points. Moreover, Table 3.16 gathers some estimation errors, both on intensive and extensive variables. It is important to remark that the P-h diagrams in Figure 3.24 do include pressure drop at the condenser, both in the experimental cycle and in the simulated one. Minor deviations are appreciated when comparing simulated and experimental P-h diagrams, therefore the separately identified models allow to correctly estimate the whole steady-state cycle. As previously stated, it must be taken into account that the proposed steady-state model is intended to be included 92 Chapter 3. Experimental plant Specific enthalpy Pressure P−h diagram Condenser R404a saturated vapour R404a saturated liquid EEV1 EEV2 Evaporator 1 Evaporator 2 Main compressor Booster compressor Cycle Freezing level Cooling level Ambient Pc hboost,dis TSH,main TSH,boost Pe5 Pe20 hc,out TSC he5,out hc,in hmain,suc he20,out Figure 3.23 P-h diagram of a two-stage, two-load-demand cycle. in model-based control strategies. This is a control-oriented model whose main features are simplicity and high calculation speed, while it is expected that the robustness of controllers can cope with the modelling errors. 3.3 Validation 93 250 300 350 400 450 500 1 2 5 10 20 30 40 Specific enthalpy [kJ kg−1] Pressure [bar] P−h diagram Condenser EEV1 EEV2 Evaporator 1 Evaporator 2 Main compressor Booster compressor Experimental cycle Simulated cycle Freezing level Cooling level Ambient (a) Operating point 1 (OP1) 250 300 350 400 450 500 1 2 5 10 20 30 40 Specific enthalpy [kJ kg−1] Pressure [bar] P−h diagram Condenser EEV1 EEV2 Evaporator 1 Evaporator 2 Booster compressor Main compressor Experimental cycle Simulated cycle Freezing level Cooling level Ambient (b) Operating point 2 (OP2) Figure 3.24 Experimental and simulated P-h diagrams at different operating points considering the two-stage, two-load-demand configuration. 94 Chapter 3. Experimental plant 250 300 350 400 450 500 1 2 5 10 20 30 40 Specific enthalpy [kJ kg−1] Pressure [bar] P−h diagram Condenser EEV1 EEV2 Evaporator 1 Evaporator 2 Booster compressor Main compressor Experimental cycle Simulated cycle Freezing level Cooling level Ambient (c) Operating point 3 (OP3) 250 300 350 400 450 500 1 2 5 10 20 30 40 Specific enthalpy [kJ kg−1] Pressure [bar] P−h diagram Condenser EEV1 EEV2 Evaporator 1 Evaporator 2 Main compressor Booster compressor Experimental cycle Simulated cycle Freezing level Cooling level Ambient (d) Operating point 4 (OP4) Figure 3.24 Experimental and simulated P-h diagrams at different operating points considering the two-stage, two-load-demand configuration. 3.4 Final remarks 95 Table 3.16 Relative errors [%] on diverse measurable variables in the two-stage, two-load- demand cycle. Variable Operating point OP1OP2OP3OP4 Pe20 3.596 0.160 0.474 8.179 Pe50.007 0.003 0.011 0.005 Pc,in 2.745 0.294 0.400 2.263 Pc,out 2.716 0.440 0.298 2.431 he20,out 0.003 0.255 0.020 0.652 he5,out 1.269 1.411 1.025 1.329 hboost,dis 3.910 0.199 1.182 2.995 hmain,suc 0.422 0.007 0.143 0.960 hmain,dis 1.285 1.022 0.316 1.843 hc,out 1.386 0.111 0.469 2.168 3.4 Final remarks This Chapter has been devoted to the two-stage, two-load-demand experimental refrigeration plant mentioned in Chapter 1. The design, implementation, and automation of the plant have been detailed, whereas parameter identification and validation have been also addressed. The modelling conclusions obtained in Chapter 2 have been applied to the identification procedure, namely the steady-state models presented in Chapter 2 for the compressors and expansion valves have been identified using experimental data. Regarding the heat exchangers, parameter dependence of heat exchanger dynamics and statics has been studied, concluding that the heat transfer coefficients influence essentially system statics but they do not affect much system dynamics. Therefore, steady-state identification of the heat exchangers has been also addressed. A novel identification procedure focused on the heat exchangers has been presented. Diverse refrigerant phases along each heat exchanger, according to the MB modelling approach, have been considered, thus diverse zones are differentiated, whose lengths are inaccessible state variables. It has been assumed that a unique overall heat transfer coefficient can be identified for each zone. Coherent values have been calculated taking into account some orders of magnitude of convective heat transfer coefficients in the literature and their influence on overall coefficients. In spite of individually modelling all elements of the cycle, global validation considering the one-compression-stage, two-load-demand configuration and the two-compression- stage, two-load-demand one has been performed. Minor differences have been noticed when comparing simulated and experimental P-h diagrams, giving rise to estimation errors around 10% in the key variables. The identification procedure of the heat exchangers is 102 Chapter 4. Optimization and controllability analysis Cooling power [W] 375 425 475 525 575 625 675 725 775 825 875 COP 1.08 1.13 1.18 1.23 1.28 Optimal COP Cooling power [W] 375 425 475 525 575 625 675 725 775 825 875 Temperature [K] 1 7 13 19 25 Optimal degree of superheating (a) Optimal COP and degree of superheating Cooling power [W] 375 425 475 525 575 625 675 725 775 825 875 Pressure [bar] 14 15 16 17 18 Condenser pressure Cooling power [W] 375 425 475 525 575 625 675 725 775 825 875 Pressure [bar] 0.8 1 1.2 1.4 Evaporator pressure (b) Condenser and evaporator pressures Figure 4.2 Optimal cycle variables for a given cooling load range. 4.1 Global optimization 103 Cooling power [W] 375 425 475 525 575 625 675 725 775 825 875 Mass flow [g s-1] 3 4 5 6 7 8 9 10 Refrigerant mass flow (c) Refrigerant mass flow or active charge Cooling power [W] 375 425 475 525 575 625 675 725 775 825 875 Compressor speed [Hz] 30 40 50 Compressor speed Cooling power [W] 375 425 475 525 575 625 675 725 775 825 875 Valve opening [%] 25 35 45 55 65 75 Expansion valve opening (d) Steady-state manipulated variables Figure 4.2 Optimal cycle variables for a given cooling load range. 104 Chapter 4. Optimization and controllability analysis Firstly, and most importantly, it is highlighted that optimal cycles are not achieved with minimum degree of superheating for all the cooling load range, but for low and medium cooling load the optimal cycle requires a greater degree of superheating. It involves that the control strategy conventionally applied (setting a constant degree of superheating as reference) might drive the cycle far away from the optimum, specially when the cooling load is lower than the design value. Therefore, an improvement on energy efficiency could be achieved if the cycle was optimally operated, in particular for low cooling load. Secondly, it is observed that the COP is greatly reduced when, to satisfy the cooling demand, the cycle needs to increase the compressor speed above its minimum. This occurs when it is no longer possible to satisfy the cooling demand by increasing the expansion valve opening, which implies that a greater refrigerant mass flow circulates. Thirdly, concerning the pressures, for very low cooling demand the dominant constraint is that of minimum condenser pressure, which is determined by the inlet temperature of the condenser secondary flux and the minimum temperature difference between the refrigerant and the secondary flux imposed at the condenser for heat exchanger efficiency purposes. Obviously the results cannot be extrapolated to all refrigeration systems nor all cooling demands, since the objective function is as nonlinear as the system model. However, it may be worth calculating the optimal cycle for the facility under study to satisfy a given cooling demand using the steady-state model, because it might provide some guidelines on set point selection, as analysed in Chapter 5. 4.2 Controllability analysis As stated in Section 4.1, a one-compression-stage, one-load-demand cycle is completely defined by a set ψcycle made up of three variables, whether considering only the compressor speed N and the expansion valve opening Av as manipulated variables. The variable set may be selected as indicated in Equation (4.1) , but any set including three independent variables is valid to describe the cycle. It has been also stated that the simplified control-oriented model presented in Chapter 2 for the whole cycle is considered as starting point for the controllability analysis. It concentrates the dominant dynamics at the condenser and proposes a three-variable state vector, both in mode 1 and 2, as indicated in Equations (4.4) and (4.5). xc,m1= [Pcζc,sc hc,sc]T(4.4) xc,m2= [Pcζc,sh ¯ γc]T(4.5) Since the variables included in the state vector are independent, both in mode 1 and 2, the state vector can represent a valid variable set to completely describe the cycle. Therefore, it is suitable to study the controllability of the condenser dynamic model, since driving the state from a random initial point to the desired point corresponds to driving the whole cycle to the point defined by the desired state. The reduced model of the condenser when working in mode 1 has been detailed in Chapter 2, Subsection 2.2.3, but it is also included in Equation Sets (4.6) and (4.7) to 4.2 Controllability analysis 105 improve the comprehension of the controllability analysis. Moreover, the expressions of the elements fc,m1,i∀i=1,2,3 are also replicated in Equation Set (4.8) for comparison purposes. Zc,m1(xc,m1,wc)˙xc,m1=fc,m1(xc,m1,wc) xc,m1= [Pcζc,sc hc,sc]T wc= [ ˙mc,sec Tc,sec,in ˙m hc,in]T (4.6) Zc,m1(xc,m1,wc) =   zc,m1,11 −1zc,m1,13 zc,m1,21 0zc,m1,23 zc,m1,31 −1zc,m1,33   fc,m1(xc,m1,wc) = [fc,m1,1fc,m1,2fc,m1,3]T (4.7) fc,m1,1(xc,m1,wc) = 1 ρc,t p VR 1 hc,f−hc,sc ˙ Qc,sc +˙m(hc,f−hc,sc) fc,m1,2(xc,m1,wc) = 1 ρc,t p ζc,t p VR ˙ Qc,t p +˙m(hc,g−hc,t p)+ hc,f−hc,t p hc,f−hc,sc ˙ Qc,sc! fc,m1,3(xc,m1,wc) = 1 ρc,sc VR 1 hc,f−hc,sc ˙ Qc,sc +˙m(hc,f−hc,sc) (4.8) Similarly, the reduced model of the condenser when working in mode 2 is indicated in Equation Set (4.9) , where the structure of coefficient matrix Zc,m2 and force function fc,m2 are presented in Equation Set (4.10) . The expressions of the elements fc,m2,i∀i=1,2,3 are detailed in Equation Set (4.11). Zc,m2(xc,m2,wc)˙xc,m2=fc,m2(xc,m2,wc) xc,m2= [Pcζc,sh ¯ γc]T(4.9) Zc,m2(xc,m2,wc) =   zc,m2,11 0zc,m2,13 zc,m2,21 1 0 zc,m2,31 −1zc,m2,33   fc,m2(xc,m2,wc) = [fc,m2,1fc,m2,2fc,m2,3]T (4.10) 106 Chapter 4. Optimization and controllability analysis fc,m2,1(xc,m2,wc) = 1 ρc,t p ζc,t p VR˙ Qc,t p −˙m(hc,out −hc,t p)+ +hc,g−hc,t p hc,g−hc,sh ˙ Qc,sh +˙m(hc,in −hc,sh) fc,m2,2(xc,m2,wc) = 1 ρc,sh VR 1 hc,g−hc,sh −˙ Qc,sh +˙m(hc,g−hc,in) fc,m2,3(xc,m2,wc) = 1 ρc,t p VR 1 hc,g−hc,sh ˙ Qc,sh −˙m(hc,g−hc,in) (4.11) 4.2.1 Analysis based on linear theory A convenient model representation in mode 1 can be obtained applying algebraic manipulation, as shown in Equation Set (4.12) , where force function ˆ fc,m1 has a null element. The same strategy can be applied to obtain the model representation in mode 2 indicated in Equation Set (4.13), where force function ˆ fc,m2has also a null element. ˆ Zc,m1˙xc,m1=ˆ fc,m1 ˆ Zc,m1=    zc,m1,11 −1zc,m1,13 zc,m1,21 0zc,m1,23 zc,m1,31 −ρc,t p ρc,sc zc,m1,11 −1+ρc,t p ρc,sc zc,m1,33 −ρc,t p ρc,sc zc,m1,13     ˆ fc,m1= [fc,m1,1fc,m1,20]T (4.12) ˆ Zc,m2˙xc,m2=ˆ fc,m2 ˆ Zc,m2=    zc,m2,11 0zc,m2,13 zc,m2,21 1 0 zc,m2,31 +ρc,sh ρc,t p zc,m2,21 −1+ρc,sh ρc,t p zc,m2,33     ˆ fc,m2= [fc,m2,1fc,m2,20]T (4.13) Considering the representation shown in Equation Set (4.12) , elements fc,m1,1 and fc,m1,2 can be considered as virtual manipulated inputs, since for a given state and disturbances their values depend on the actual manipulated inputs ˙m and hc,in , whose values depend in turn on the actual control actions N and Av . In other words, given the state and disturbances, the actual manipulated inputs ˙m and hc,in can be calculated from the virtual manipulated inputs fc,m1,1 and fc,m1,2 by solving the nonlinear system generated by the two first equations of Equation Set (4.8) . The same can be applied to elements fc,m2,1 and fc,m2,2 in Equation Set (4.13) , where the two first equations of Equation Set (4.11) generate the nonlinear system from which ˙m and hc,in are calculated. Then, the underactuated features of the 4.2 Controllability analysis 107 whole cycle also emerge here, when studying only the condenser model, since the control objective is to regulate the three states by manipulating only two variables. Coefficient matrices ˆ Zc,m1 and ˆ Zc,m2 depend on the state, manipulated variables, and disturbances, in addition to thermodynamic properties of the refrigerant when evaluated at the intensive variables gathered in the state, specially the condenser pressure Pc . None of them are singular for a wide enough pressure range, reasonable values of disturbances, and achievable control actions, thus each one can be inverted giving rise to the model representation indicated in Equation Set (4.14) , which is valid for modes 1 and 2 if changing the corresponding subscript mi∀i=1,2. ˆ Zc,mi˙xc,mi=ˆ fc,mi ˆ Z−1 c,mi≡ˆ Zc,mi,inv ˆ Zc,mi,inv =  ˆzc,mi,inv,11 ˆzc,mi,inv,12 ˆzc,mi,inv,13 ˆzc,mi,inv,21 ˆzc,mi,inv,22 ˆzc,mi,inv,23 ˆzc,mi,inv,31 ˆzc,mi,inv,32 ˆzc,mi,inv,33   (4.14a) ˙xc,mi=ˆ Zc,mi,inv ˆ fc,mi=  ˆzc,mi,inv,11 ˆzc,mi,inv,12 ˆzc,mi,inv,13 ˆzc,mi,inv,21 ˆzc,mi,inv,22 ˆzc,mi,inv,23 ˆzc,mi,inv,31 ˆzc,mi,inv,32 ˆzc,mi,inv,33    fc,mi,1 fc,mi,2 0 = =  ˆzc,mi,inv,11 ˆzc,mi,inv,12 ˆzc,mi,inv,21 ˆzc,mi,inv,22 ˆzc,mi,inv,31 ˆzc,mi,inv,32  fc,mi,1 fc,mi,2(4.14b) The model representation shown in Equation (4.14b) has a linear structure if ˆ fc,mi is considered as the manipulated input vector, in spite of all elements in ˆ Zc,mi,inv and ˆ fc,mi being nonlinear functions of the state xc,mi , the actual manipulated inputs ˙m and hc,in , and the disturbances ˙mc,sec and Tc,sec,in included in wc . Given the linear structure of the model representation, and despite the system nonlinearities, the linear controllability matrix CCC gives a first overall assessment of the difficulties which emerge when controlling the system. As described in Equation Set (4.15) , its rank is shown to match at most the number of columns in matrix Bc,mi , namely the number of manipulated inputs. Then, the controllability analysis based on linear theory suggests that the system is not fully controllable. ˙xc,mi=Ac,mixc,mi+Bc,miuc,mi(4.15a) 108 Chapter 4. Optimization and controllability analysis Ac,mi=  000 000 000  Bc,mi=ˆ Zc,mi,inv =  ˆzc,mi,inv,11 ˆzc,mi,inv,12 ˆzc,mi,inv,21 ˆzc,mi,inv,22 ˆzc,mi,inv,31 ˆzc,mi,inv,32   uc,mi=ˆ fc,mi=fc,mi,1 fc,mi,2 (4.15b) CCC =Bc,miAc,miBc,miA2 c,miBc,mi rank(CCC) = rank(Bc,mi)⩽2(4.15c) 4.2.2 Nonlinear analysis based on the phase portrait method However, all elements in matrix Bc,mi∀i=1,2 depend on the disturbances and the state. Therefore, the conclusions provided by the controllability analysis on the simplified SMB model based on linear theory could be questionable, since it is a highly nonlinear system. That is the reason why the linear controllability study is extended to a nonlinear analysis based on the phase portrait method. The dynamic equations of the control-oriented SMB model are highly nonlinear and include indeed thermodynamic functions which are used to evaluate refrigerant properties. This causes an analytical study to be very difficult to carry out, therefore a graphic method is instead addressed to evaluate the controllability of the system. If the model representation shown in Equation (4.14b) is considered, for both condenser modes, vectors χc,m1 and χc,m2 can be defined as the deviations of the states with respect to the desired values, as indicated in Equation Set (4.16) . For the sake of simplicity, the desired state is not expected to vary with time, thus the error model shown in Equation (4.17) arises, valid for both modes. χc,m1≡xc,m1−xre f c,m1=  Pc−Pre f c ζc,sc −ζre f c,sc hc,sc −hre f c,sc   χc,m2≡xc,m2−xre f c,m2=  Pc−Pre f c ζc,sh −ζre f c,sh ¯ γc−¯ γre f c   (4.16) ˙χc,mi=˙xc,mi−˙xre f c,mi=˙xc,mi=  ˆzc,mi,inv,11 ˆzc,mi,inv,12 ˆzc,mi,inv,21 ˆzc,mi,inv,22 ˆzc,mi,inv,31 ˆzc,mi,inv,32  fc,mi,1 fc,mi,2(4.17) 4.2 Controllability analysis 109 For the sake of notation simplicity, let rename the elements of the matrix defining the dynamic model of χc,mias shown in Equation (4.18). ˙χc,mi=  ˙ χc,mi,1 ˙ χc,mi,2 ˙ χc,mi,3 =  ˆzc,mi,inv,11 ˆzc,mi,inv,12 ˆzc,mi,inv,21 ˆzc,mi,inv,22 ˆzc,mi,inv,31 ˆzc,mi,inv,32  fc,mi,1 fc,mi,2≡ ≡  bc,mi,11 bc,mi,12 bc,mi,21 bc,mi,22 bc,mi,31 bc,mi,32  fc,mi,1 fc,mi,2(4.18) In order to make the graphic results of the nonlinear controllability analysis more intelligible, the three-dimensional problem is simplified into a two-dimensional one. It is considered that one of the virtual manipulated inputs is devoted to controlling one of the state variables, in such a way that this state holds its reference value. Therefore, the control problem is reduced to driving the two remaining states to their references by using only one manipulated input. It allows the phase portrait of this reduced problem to be represented in a two-dimensional chart. Consider for instance mode 1. If a secondary controller is responsible for holding hc,sc at its reference value hre f c,sc at any instant, it can be assumed that ˙ χc,m1,3=0 . The dynamic model of the remaining states is simplified as shown in Equation Set (4.19) , where the secondary control law is explicitly calculated in Equation Set (4.19a). ˙ χc,m1,3=0⇒bc,m1,31 fc,m1,1+bc,m1,32 fc,m1,2=0⇒fc,m1,1=−bc,m1,32 bc,m1,31 fc,m1,2(4.19a) ˙ χc,m1,1 ˙ χc,m1,2=bc,m1,11 bc,m1,12 bc,m1,21 bc,m1,22  fc,m1,1 fc,m1,2= =    bc,m1,12 −bc,m1,32 bc,m1,31 bc,m1,11 bc,m1,22 −bc,m1,32 bc,m1,31 bc,m1,21    fc,m1,2≡ ≡dc,m1,1 dc,m1,2fc,m1,2= =Dc,m1fc,m1,2 (4.19b) The reduced two-dimensional problem shows not to be fully controllable if the elements dc,m1,1 and dc,m1,2 are constant. In particular, the one-dimensional controllable subspace SSSc,m1 can be analytically calculated as shown in Equation Set (4.20) and it corresponds to the line represented in Figure 4.3 considering random values of dc,m1,1 and dc,m1,2 . The arrows in Figure 4.3 represent that the movement due to the available control action is possible along the line in both directions, but in any case its magnitude is limited by the physically achievable control action. As Figure 4.3 illustrates, only initial points within the controllable subspace SSSc,m1 can be driven to the origin by any controller, and moreover it 110 Chapter 4. Optimization and controllability analysis will fail in achieving the origin from any initial point out of the controllable subspace SSSc,m1 , since controllability is a system property, not a controller-dependent one. Furthermore, it is shown in Equation Set (4.21) that the distance δc,m1 from any initial point to the controllable subspace SSSc,m1does not vary for any control action fc,m1,2. SSSc,m1=(χc,m1,1 χc,m1,2∈R2|χc,m1,1dc,m1,1 χc,m1,2dc,m1,2=0)= =(χc,m1,1 χc,m1,2∈R2|χc,m1,1dc,m1,2−χc,m1,2dc,m1,1=0)= =(χc,m1,1 χc,m1,2∈R2|χc,m1,2=dc,m1,2 dc,m1,1 χc,m1,1) (4.20) δc,m1≡χc,m1,1dc,m1,1 χc,m1,2dc,m1,2 ˙ δc,m1=˙ χc,m1,1dc,m1,2−˙ χc,m1,2dc,m1,1= = (dc,m1,1dc,m1,2−dc,m1,2dc,m1,1)fc,m1,2=0∀fc,m1,2 (4.21) -1 -0.5 0 0.5 1 -0.1 0 0.1 0.2 0.3 0.4 Controllable subspace SSSc,m1 δc,m1 IP χc,m1,1 χc,m1,2 Figure 4.3 Controllable subspace of the system described in Equation (4.19b) for constant values of parameters dc,m1,1and dc,m1,2. 4.2 Controllability analysis 111 Nevertheless, the parameters dc,m1,1 and dc,m1,2 of the system to be controlled indicated in Equation (4.19b) are not constant, but state-dependent, as well as the virtual manipulated input fc,m1,2 . Moreover, they depend on the actual manipulated variables ˙m and hc,in . Therefore a pointwise analysis is carried out, considering the nonlinear equations of the system. First of all, some equilibrium points have been considered on the phase portrait shown in Figure 4.4. All of them observe the limitation of hc,sc =hre f c,sc imposed previously, thus the two-dimensional control problem of driving χc,m1,1 and χc,m1,2 to the origin by using only fc,m1,2as manipulated input is addressed. -0.5 0 0.5 1 1.5 2 -0.2 -0.1 0 0.1 0.2 0.3 0.4 χc,m1,1=Pc−Pref c[bar] χc,m1,2=ζc,sc −ζref c,sc Phase plane for χc,m,=hc,sc −href c,sc =0 kJ kg-1 Origin Controllable subspace from origin Random equilibrium point χc,m1 Controllable subspace from χc,m1 Figure 4.4 Phase plane of the system described in Equation (4.19b). Then, the values of parameters dc,m1,1 and dc,m1,2 corresponding to each equilibrium point have been precisely calculated for a variety of actual manipulated inputs ˙m and hc,in effectively achievable by the compressor and the expansion valve. All pairs { ˙m,hc,in } must satisfy the secondary control law shown in Equation Set (4.19a) to ensure that the assumption of χc,m1,3=0 holds at any instant. The values of dc,m1,1 and dc,m1,2 define the slope of the line which represents the controllable subspace for each equilibrium point, as shown in Figure 4.4 for reasonable ranges of χc,m1,1 and χc,m1,2 around the origin. Note that the origin and its controllable subspace SSSc,m1 have been graphically emphasized. Naturally, the values of dc,m1,1 and dc,m1,2 achieved for each pair { ˙m,hc,in } are different and thus generate different slopes, but the standard deviation of such slopes with respect 118 Chapter 5. Control Co ntrol of one-stage, one-load-demand vapour-compression refrigeration systems is addressed in this Chapter. It was stated in Chapter 1 that there are mainly two schemes in the literature about control. The conventional scheme in industry is very simple: in addition to the reference on Te,sec,out imposed by the cooling demand, a low but constant set point on TSH is applied and the controller is designed to get these two variables to track their references as efficiently as possible in presence of disturbances by manipulating N and Av . Optimal steady-state operation is not explored through the references in the conventional control scheme, but the low degree of superheating set as reference is believed to lead the cycle to high energy efficiency. However, there are other control strategies in the literature where energy efficiency is considered when calculating the references for the low-level controller, which may be updated according to the dynamic behaviour of the system. The objective of this Chapter is to make a contribution in both control schemes: tracking control applied to the conventional scheme, and energy-efficiency-aware control. Regarding the tracking control, most developed techniques are linear controllers. However, as analysed in the literature, when modelling a refrigeration system linearly, some problems could arise and lead to closed-loop instability when applying a linear controller. This is due to the high thermal inertia, dead times, high coupling between system variables, and strong nonlinearities. If the controller was robust enough, it could cope with those issues even at operating points further from the design point, despite its linearity. Indeed, these uncertainties should be studied not only along frequency, as proposed by Larsen and Holm [ 38 ], but also at several operating points which could be reached. These operating points depend mainly on the cooling demand, which may vary according to the external conditions of the plant (for instance ambient temperature, inlet temperature of evaporator secondary flux, etc). Therefore, a controller which outperforms at the design point is not intended, but a controller which could work well at several and distant operating points, along with diverse frequency ranges. According to this, the development of a multivariable centralised H∞ controller, based on the S/KS/T Mixed Sensitivity Problem, is proposed in this Chapter. The objective of this controller is only to track the references on TSH and Te,sec,out , which are selected as indicated within the conventional control scheme: the cooling demand is expressed as a desired value for Te,sec,out , while a low but constant degree of superheating is set as reference for the tracking controller. First of all, multivariable linear models are proposed and identified at a number of operating points, in order to characterise the main system dynamics. These operating points are selected according to the desired working area of the system, consisting of positive degree of superheating TSH , and also bearing in mind that, for energy-efficiency purposes, some zones of the working area are more likely to be explored in steady state. The linear models have been identified using the step response method, and one of the operating points has been selected as the design one, seeking to minimise the uncertainty region. Secondly, a controllability analysis of the nominal model (that identified at the design point) is performed to find out possible dynamic constraints in its closed-loop performance, being imposed by the inherent features of the system. Finally, a methodology which easily allows to impose the desired performance of the closed-loop system on the calculation of the weighting matrices of the S/KS/T Problem is described. This methodology has already successfully applied to other thermodynamic systems, and it has evidenced to be effective in enhancing the closed-loop performance in spite of its 5.1 Robust H∞tracking controller 119 simplicity [86, 87]. Concerning the energy-efficiency-aware control, it is observed that all strategies in the literature show trouble in reaching the optimal cycle satisfying a certain cooling demand, as detailed in Chapter 1. Indeed, the mentioned control strategies described in Chapter 1 do not even aspire to achieve the global optimal cycle, but they propose suboptimal solutions manipulating the available control actions. Given the lack of full controllability suggested by the existing control strategies and shown by the controllability analysis performed in Chapter 4, a suboptimal model-based control architecture is proposed. It has been stated that a cycle is defined by three variables but only two can be effectively controlled, and one of them must be devoted to satisfying the cooling demand. Therefore, an online one-degree-of-freedom optimization is proposed to obtain the best reference for the remaining variable to be controlled according to the selected energy-efficiency metrics. The hierarchical control structure is similar to that used by Jain [ 12 ]; the difference lies that in this case the global optimal solution is not projected from the optimization space onto the control space, but a suboptimal solution is calculated online considering the dynamic behaviour of the uncontrolled variable. The optimizer considers the cooling demand constraint, as well as some limitations of the components and operating constraints, while a nonlinear steady-state model of the whole cycle is used to characterise the feasible variable sets. Furthermore, the low-level controller which realises the references provided by the optimizer is a novel technique in predictive control field, called Practical Nonlinear Model Predictive Control (PNMPC) [88, 89]. Regarding the online optimizer, the different suboptimal cycles for a given facility and a certain cooling demand are compared regarding the COP according to the value considered for the uncontrolled variable. These optimization results are compared to those of the global optimization presented in Chapter 4, Section 4.1. The performance of the SHC architecture is compared to the FB+FF strategy proposed by Jain [ 12 ] and the TV-ESC developed by Guay [ 48 ], regarding the COP achieved in steady state and the dynamic behaviour of the controlled variables. The contents of this Chapter have been included in some of the publications presented in Chapter 1, Section 1.5 [50, 51, 53, 55, 56, 60, 62, 64]. 5.1 Robust H∞tracking controller 5.1.1 Modelling Only the heat exchangers are dynamically modelled, since their dynamics are dominant over the remaining components of the cycle, as stated in Chapter 1. The compressor and expansion valve are statically modelled as detailed in Chapter 2, Subsections 2.1.1 and 2.1.2 respectively. However, the original dynamic modelling of the heat exchangers developed by Li and Alleyne is considered for the design of the robust controller [ 23 ]. This model was mostly described in Chapter 2, Subsection 2.2.1, but the latter was in turn a simplification of the original modelling by Li and Alleyne. Only the significant differences are remarked here, since the complete model equations can be found in the work cited for the keen reader. 120 Chapter 5. Control For all modes, the state vectors in Equation Set (5.1) gather all relevant information concerning the condenser and the evaporator at each instant. Note that in this case the wall temperatures at all zones are included as states, which increases the dimension of both state vectors. xc= [hc,sh Pchc,sc ζc,sh ζc,t p Tc,t p,wTc,sc,w¯ γc]T xe= [ζe,t p Pehe,sh Te,t p,wTe,sh,w¯ γe]T(5.1) The constant dynamic structure formulated in the nonlinear descriptor form shown in Equation Set (2.14) holds for all modes. In case of losing any zone, the related state variables hold their values, being only transformed the appropriate equations of the models into auxiliary ones which simply get the corresponding variable to track a reference value. Some technological requirements are taken into account for the design of the robust controller. As commented in Chapter 1, the refrigerant flow at the compressor intake should never be two-phase, leading the evaporator to work only in mode 1. Furthermore, due to energy-efficiency purposes, some combinations of the manipulated variables are more likely to be explored in steady state. Then, the COP is represented in the space of the manipulated variables for a given facility in Figure 5.1. Note that zero COP values correspond to infeasible combinations of Nand Av. 90 80 COP in the space of manipulated variables 70 60 Valve opening [%] 50 8 40 5 7 30 6 4 3 2 1 20 1800 2475 3150 3825 1.5 1 0.5 0 4500 Compressor speed [rpm] COP Figure 5.1 Achieved COP values in the space of manipulated variables. A maximum COP line from N=1800 rpm to N=4500 rpm and from Av=35% to Av=65% can be noticed in Figure 5.1. The working area, defined as the combinations of manipulated variables which lead to the desired operation of the evaporator (mode 1) while ensuring high COP, is highlighted in Figure 5.2. 5.1 Robust H∞tracking controller 121 90 80 70 60 Desired working area Valve opening [%] 50 8 5 40 7 6 4 30 3 2 1 20 1800 2475 3150 3825 0 1 4500 Compressor speed [rpm] Figure 5.2 Working area represented in the space of manipulated variables. A series of eight operating points have been selected based on the working area defined in Figure 5.2, and taking into account the maximum COP line depicted in Figure 5.1, which approximately matches the right border in Figure 5.2. These points are not uniformly arranged over the working area, but focused around the line which ensures higher COP values. Table 5.1 gathers the selected operating points. Table 5.1 Selected operating points. Operating point N[rpm] Av[%] OP11800 25 OP22138 30 OP32475 35 OP42813 40 OP53150 50 OP63488 45 OP73825 50 OP84163 60 The main dynamics of the system at the selected operating points are characterised by 122 Chapter 5. Control means of low-order MIMO linear models. These models have been identified using the well-known step response method. Figure 5.3 shows a typical step response (specifically at the operating point OP 4 ), being those corresponding to the other operating points qualitatively similar. Since the Te,sec,out step response when applying a change on N shows overshooting but not oscillations, and the rest of step responses exhibit non-minimum phase behaviour (although this effect is very fast in the case of the TSH step response when applying a change on N and it is not appreciable in Figure 5.3(a)), every model is expressed in the continuous transfer matrix form shown in Equation (5.2) , whose transfer functions have been modelled as a static gain ( Ki j ) and the time constants corresponding to a zero (τzi j ) and two poles, one fast (τf asti j ) and another one slower (τslowi j ). ∆TSH (s) ∆Te,sec,out (s)=G(s)∆N(s) ∆Av(s) G(s) =        K11(τz11 s+1) (τf ast11 s+1)(τslow11 s+1) K12(τz12 s+1) (τf ast12 s+1)(τslow12 s+1) K21(τz21 s+1) (τf ast21 s+1)(τslow21 s+1) K22(τz22 s+1) (τf ast22 s+1)(τslow22 s+1)        (5.2) The selected nominal operating point is OP 4 , since it minimises the uncertainty region. The numerical values of Ki j,τzi j ,τf asti j , and τslowi j are collected in Table 5.2. 5.1.2 Controllability analysis After identifying a linear model which describes the main dynamics of the system at the nominal operating point, a controllability analysis is carried out with the aim of finding out possible limitations on its closed-loop performance [ 90 ]. A study of the zeros and poles and their directions, singular values, and condition number of the system is performed. Prior to the controllability analysis, the system has to be adjusted in scale, according to the maximum deviation of each input and output, which are indicated in Equation Set (5.3). u   ∆Nmax =337 rpm ∆Amax v=10% y   ∆Tmax SH =4°C ∆Tmax e,sec,out =0.05°C (5.3) Using these values, the system can be now scaled according to Equation (5.4) , where Gnom(s) is the transfer matrix identified at the nominal operating point and ˆ Gnom(s) is the 5.1 Robust H∞tracking controller 123 Time [min] 19 21 23 25 27 29 Temperature [K] 5 7 9 11 13 15 Degree of superheating TSH N Compressor speed [rpm] 2400 2500 2600 2700 2800 2900 (a) Dynamic response of TSH when applying a step change on N Time [min] 19 21 23 25 27 29 Temperature [K] 10 12 14 16 18 20 Degree of superheating TSH Av Valve opening [%] 34.5 36.5 38.5 40.5 (b) Dynamic response of TSH when applying a step change on Av Figure 5.3 Step response at the nominal operating point (OP4). 124 Chapter 5. Control Time [min] 19 21 23 25 27 29 Temperature [ºC] -20.76 -20.74 -20.72 -20.7 -20.68 -20.66 Outlet temperature of evaporator secondary flux Te,sec,out N Compressor speed [rpm] 2400 2500 2600 2700 2800 2900 (c) Dynamic response of Te,sec,out when applying a step change on N Time [min] 19 21 23 25 27 29 Temperature [ºC] -20.8 -20.76 -20.72 -20.68 -20.64 -20.6 Outlet temperature of evaporator secondary flux Te,sec,out Av Valve opening [%] 34.5 36.5 38.5 40.5 (d) Dynamic response of Te,sec,out when applying a step change on Av Figure 5.3 Step response at the nominal operating point (OP4). 5.1 Robust H∞tracking controller 125 Table 5.2 Numerical values of the identified transfer functions. Operating point KK rpm K % K rpm K % τz[s] τfast [s] τslow [s] OP11.4·10−2−1.1 −6.3·10−5−1.6·10−2  −1−9.1 200 −3.6  0.3 6.7 18.2 3.6  28.6 20 105.3 28.6 OP21.8·10−2−1.5 −6.5·10−5−1.6·10−2  −1−5.5 133 −3.6 0.25 4 15.4 3.6  25 26.3 66.7 33  OP31.6·10−2−1.3 −6.7·10−5−1.6·10−2  −1−5.5 133 −5  0.3 4 15.4 6.7  25 25 66.7 30.3 OP41.4·10−2−1.1 −7.1·10−5−1.7·10−2  −1−9.1 142.9−6.7  0.3 6.7 14.3 10  28.6 20 66.7 25  OP51.3·10−2−1.14 −9.1·10−5−1.5·10−2 −13 −13.9 55 −12.5 11.1 9.1 10.8 20  12.5 22 20 20  OP61.1·10−2−1 −6.2·10−5−1.7·10−2  −1−10 105.3−5  0.3 4 12.5 7.1 13.2 20 58.8 25  OP79.9·10−3−1 −6.8·10−5−1.8·10−2  −1−11.8 105.3−5.5  0.3 5 12.5 10.5 13.2 18.2 58.8 25  OP81.7·10−2−1.1 −1.6·10−4−1.5·10−2 −1.25 −11.8 41.7−10   0.5 6.7 13.6 11.8 22.2 20 33.3 20  scaled transfer matrix at the same point. ˆ Gnom(s) = SM−1 err Gnom(s)SMu(5.4) Scaling matrices SMerr and SMu have been calculated as indicated in Equation Set (5.5). SMerr =∆Tmax SH 0 0∆Tmax e,sec,out =4 0 0 0.05  SMu=∆Nmax 0 0∆Amax v=337 0 0 10 (5.5) The first step of the analysis consists in studying the poles and zeros of the system, as explained by Skogestad and Postlethwaite [ 90 ]. Table 5.3 shows these values in the continuous time domain. 126 Chapter 5. Control Table 5.3 Poles and zeros of the nominal model. Roots Value Poles −7·10−2 −5·10−2 −4·10−2 −3.5·10−2 −1.5·10−2 −1.5·10−1 −1·10−1 −3.7 Zeros −4.51 ·10−2±9.45 ·10−3j −1.224 ·10−2 −1.224 ·10−2 −1.1·10−1 1.21 ·10−1 The system is stable, as observed in Table 5.3, since there are no unstable poles. However, there is a transmission zero placed in the right-half plane (RHP), which imposes some constraints considering the closed-loop performance [ 91 ]. The RHP-zero is located at zRHP =1.21 ·10−1 . As stated by Skogestad and Postlethwaite [ 90 ], this yields the upper bound to the attainable control bandwidth ωcr stated in Equation (5.6) , which leads to the approximate lower bound for the rise time of at least one output tr indicated in Equation (5.7). ωcr <|zRHP| 2.8=4.32 ·10−2rad/s (5.6) tr≈π 2ωcr >36.4s (5.7) Considering the output direction of the system for the RHP-zero is also relevant given the multivariable nature of the system. So, zRHP has the output direction shown in Equation (5.8). TSH Te,sec,out =0.97 0.245 (5.8) It can be seen in Equation (5.8) that the first component is quite bigger than the second one. Therefore, although the degrading effect of a RHP-zero can be moved to a given output [ 92 ], its natural effect is focused on the first variable. This way, the control of the degree of superheating at the evaporator outlet is more limited than that of the outlet temperature of the evaporator secondary flux. Therefore, if a tight control on TSH is imposed, the behaviour of Te,sec,out would be strongly degraded. 5.1 Robust H∞tracking controller 127 The singular values of the system along frequency are the next topic to be analysed, focusing on the minimum singular value because of its relevance as a controllability index. The singular values of the nominal model are shown in Figure 5.4, whereas Table 5.4 presents their directions at low frequency [ 90 ]. It can be observed that from ωmin,SV =10−1 rad/s onwards the minimum singular value becomes less than 0 dB, what means that until that frequency independent output changes can be applied. From ωmin,SV onwards, the dependence between both outputs becomes remarkable. Nevertheless, taking into account the limitations obtained from Equations (5.6) and (5.7) , both outputs are independent for all allowed frequencies. Frequency [rad s-1] 10-3 10-2 10-1 100101 Magnitude 10-3 10-2 10-1 100 101Singular values of the nominal model ωCN ωmin,SV ωcr σmax σmin Frequency [rad s-1] 10-3 10-2 10-1 100101 Condition number 0 2.5 5 7.5 10 Relationship between maximum and minimum singular values ωCN ωmin,SV ωcr Figure 5.4 Singular values and condition number of the nominal model along frequency. Table 5.4 Directions of the singular values at low frequency. Singular Input direction Output direction value [N Av]T[TSH Te,sec,out]T Maximum [0.0871 0.996]T[0.31 0.95]T Minimum [0.9962 −0.0871]T[0.95 −0.31]T Eventually, the condition number for all considered frequencies is also depicted in Figure 5.4, where the decrease in the minimum singular value is faster than that in the maximum value from frequencies roughly ωCN =6·10−1 rad/s. This implies a significant