Energy-efficiency control of two-stage two-load-demand refrigeration cycles.

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 Energy Efficiency Optimization and Control of Two-stage Two-load-demand Refrigeration Cycles Autor: José Enrique Alonso Alfaya 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, 2018 Tesis Doctoral Ingeniería Industrial Energy Efficiency Optimization and Control of Two-stage Two-load-demand Refrigeration Cycles Autor: José Enrique Alonso Alfaya 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 2018 Tesis Doctoral: Energy Efficiency Optimization and Control of Two-stage Two-load- demand Refrigeration Cycles Autor: José Enrique Alonso Alfaya 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: A mi familia 1 Acknowledgements This Thesis is not only consequence of my effort and dedication to it, but also the help, support, consideration, cooperation, even friendship of many people. Please, forgive me for writing the acknowledgments in my language. Muchas son las personas que han contribuido, influido o ayudado a que esta tesis doctoral pueda salir adelante. Seguro que me dejo alguna por el camino y espero que comprenda. En primer lugar, a mis directores de tesis, Manolo Gil y Paco Rodríguez Rubio, tengo agradecerles especialmente el honor de haber sido su alumno durante estos más de cuatro años, por la paciencia que han tenido, y por ofrecerme su tiempo y apoyo. Me han hecho crecer mucho más de lo que jamás hubiese pensado. Muchísimas gracias. A mis compañeros del proyecto. David, por su tranquila y constante ayuda y conversación que nos ha aportado tanto. También me quiero acordar de Jose Joaquín, que estuvo haciendo el trabajo final de máster con nosotros y nos aportó su conocimiento del campo. Imposible es olvidarme de Ramón, aunque no forme parte del proyecto ha sido un gran compañero de doctorado, siempre con una sonrisa y dispuesto a ayudarte; y Manu, también ajeno al proyecto pero siempre echando una mano cuando la necesitábamos. Por último y quizás más importante, tengo que agradecer a Guillermo toda su colaboración y ayuda durante estos cuatro años. No hay palabras para agradecerle toda la ayuda prestada, así como los tantos consejos y charlas sobre el proyecto que me ha hecho avanzar mucho más rápido de lo que pudiese imaginar. Quién me iba a decir que su ausencia durante las estancias de cada uno fuese a notarse tanto. También dentro de la universidad, no puedo olvidarme de Fernando Castaño, sin el cual la planta habría tardado muchísimo más en funcionar, y de Daniel Limón, mi tutor del ya tan lejano proyecto final de carrera, y que fue el primero que me sugirió la idea de hacer el doctorado. Además, me gustaría recordar a todos los profesores de las asignaturas en las que he participado y algunos más que he conocido por la escuela por su ayuda desinteresada para afrontar la impartición de estas y la confianza depositada en mí. De II Contents 3.6.4 Two-compressor, one-load-demand cycle 43 4 Experimental plant 47 4.1 Layouts 48 4.1.1 Plant configurations 49 4.2 Static elements 52 4.2.1 Air-cooled condenser 52 4.2.2 Water-cooled condenser 52 4.2.3 Heat exchanger for evaporators 53 4.2.4 Secondary tanks 54 4.2.5 Fluids 54 4.3 Sensors 55 4.3.1 Resistance thermometer 55 4.3.2 Pressure sensor 58 4.3.3 Flowmeter sensor 58 4.4 Actuators 59 4.4.1 Electronic expansion valves 59 4.4.2 Variable-speed compressors 60 4.4.3 Liquid pumps 61 4.4.4 Thermal resistance 62 4.4.5 Condenser fan 63 4.5 Programmable logic controller 64 4.5.1 Configuration 64 4.5.2 Input/Output modules 65 4.5.3 Modbus protocol 67 4.6 SCADA 70 4.7 OPC 71 4.8 Future work: Thermal energy storage extension 71 5 Identification 75 5.1 VCRS Identification 77 5.1.1 Steady-state parameter identification 77 5.1.2 Evaporators 81 5.1.3 Condenser 88 5.1.4 Validation 97 5.2 Conclusions for VCRS identification 98 5.3 Pump identification and control 100 5.3.1 Cooling secondary circuit pump 100 5.3.2 Pump at -20oC 103 5.4 Secondary circuits identification 107 5.4.1 Tank schemes 107 5.4.2 Cooling-demand tank identification 108 5.4.3 Freezing-demand tank identification 112 Contents III 5.4.4 Glycol substitution 114 5.4.5 Conclusions for tank identifications 114 6 Control 117 6.1 Input-output pairing 118 6.2 Robust control 119 6.2.1 Process model 119 6.2.2 Controllability analysis 122 6.2.3 H∞controller synthesis 125 6.2.4 Simulation results 128 6.2.5 Conclusions of robust control 132 6.3 Experimental control 135 6.3.1 Cascade control for cooling evaporator 135 6.3.2 Cascade control for freezing evaporator 140 6.3.3 PID control for freezing secondary circuit 143 6.3.4 One-compressor and two-load-demand configuration for control 144 6.3.5 Conclusions for experiments 147 7 Optimization 149 7.1 Optimization of one-compressor two-load-demand configuration 152 7.1.1 Optimization results for one-compressor and two-load-demand configuration 153 7.2 Optimization of two-compressor, two-load-demand configuration 157 7.2.1 Optimization results for two-compressor and two-load-demand configuration 158 7.3 Conclusions 163 8 Conclusions and future work 167 8.1 Thesis contributions and conclusions 167 8.2 Future work 169 List of Figures 171 List of Tables 177 Bibliography 179 Acronyms CV configuration valve COP Coefficient of Performance EEV electronic expansion valve GWP Global Warming Power HVAC Heat, Ventilation and Air Conditioning IC Initial Condition ISE Integral of Square Error ITAE Integral of Time multiplied by Absolute Error 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 PID proportional-integral-derivative PLC Programmable Logic Controller PO Porcentage of overshooting [%] PWM Pulse-Width Modulation RE Relative error RGA Relative Gain Array RHP Right-half plane RMS Root Mean Square SCADA Supervisory Control And Data Acquisition SISO single-input-single-output SMB Switched moving boundary V VI Acronyms SVC Superheat valve controller TES Thermal Energy Storage VFD Variable frequency drive VSC Variable speed compressor VCRS Vapor compression refrigeration system ε-NTU Effectiveness-NTU (method) Notation Latin symbols Atrns f Heat transfer perimeter [m] AvExpansion valve opening [%] bCompressor parameter CThermal capacity ratio C(s)Continuous controller transfer function 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 DExternal diameter of the condenser fan [m] dInternal diameter of the condenser fan [m] diElement {i} in vector D eError eError vector fForce function in the heat exchanger model ˆ G(s)Scaled continuous transfer matrix G(s)Continuous transfer function gHeat exchanger characteristic function hSpecific enthalpy [J kg-1] IElectric current [A] ˆ K(s)Scaled controller continuous transfer matrix KStatic gain kDiscrete time step ki j Static gain of the element i, j. VII VIII Notation kCi j Closed-loop static gain of the element i, j. KRResistance gain [W % -1] kdrop Pressure drop factor [kg-1 m-1] LHeat exchanger length [m] ˙mMass flow [kg s-1] (˙m·cp) Result of multiply the mass flow and its specific heat at constant pressure [J s-1 K-1] NCompressor speed [Hz] NTU Number of Transfer Units PPressure [Pa] PO Porcentage of overshooting [%] QCooling power [W] QRThermal power generated by the electric resistance [W] Qre f r Refrigerant power [W] Qloss Thermal power lost at the tank [W] RConduction thermal resistance [m2K W-1] RElectrical resistance [Ω] R(t)Time-variable thermal resistance [m2K W-1] RE Relative error [%] 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] Tamb Ambient temperature [K] TDDerivative time [s] TIIntegral time [s] Tre f Reference temperature for Coolprop tables [K] TSC Degree of subcooling [K] TSH Degree of superheating [K] tTime trRise time [s] UA Global heat transfer coefficient [W K-1] VVoltage [V] Vtank secondary circuit tank volume [m3] vSpecific volume [m3kg-1] vair air speed [m s-1] WKS(s)Control sensitivity weighting matrix (continuous) WS(s)Output sensitivity weighting matrix (continuous) WT(s) Output complementary sensitivity weighting matrix (continuous) Wcomp cons Mechanical power consumed by the compressor [W] wInput vector of the heat exchanger model Notation IX xState vector yOutput vector ZCoefficient matrix ZTotal number of experimental data zGeneralised plant error vector 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 β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 δDifference εHeat exchange effectiveness ζNormalised heat exchanger zone length θParameter vector in the identification process κBandwidth parameter within the H∞controller design λRelative Gain Array resulting matrix λConvective heat transfer coefficient [W m-2 K-1] λ Refrigerant power coefficient in the identification process ξResidue vector in the identification process ρDensity [kg m-3] σSingular value for the H∞control design σThermal loss coefficient [W K-1] τTime constant [s] υDifferential increment ϕExperimental input matrix in the identification process χVapor quality ψOptimization variable set ωExogenous signal vector at robust control design ω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] XNotation Subscripts and superscripts Avvalve opening acc air-cooled condenser air air amb ambient BC booster compressor ccondenser comp compressor cond conduction cycle vapor-compression cycle CL Closed loop diag diagonal desired desired element or value disc discharge e evaporators or related to secondary circuit associated to both evaporators e5 cooling evaporator (5 o C) or related to secondary circuit associated to this evaporator e20 evaporator at -20 o C or related to secondary circuit associated to this evaporator eev electronic expansion valve exp experimental est estimated fsaturated liquid gsaturated vapour IC initial condition IL inner loop in inlet, input in f er inferred int intake is isentropic it iterative m1condenser mode 1 m2condenser mode 2 MC main compressor max maximum mean mean value meas measured mid middle min minimum OL outer loop out outlet, output Notation XI pconstant pressure pred predicted puv pressure upholding valve pump pump or related to this pump RHP right-half plane re f reference re ft refrigerant sc subcooled liquid zone sec secondary flux sh superheated vapour zone tank secondary circuit tank t p two-phase zone uinput V1version 1 V2version 2 6Chapter 2. Introduction 2.2.2 Control and Optimization The literature provides several examples of basic refrigeration cycle where SISO and MIMO controls are successfully implemented typically for cooling purposes. Some SISO linear controllers can be found in the literature for experimental plant. Marcinichen et al. have studied a SISO dual strategy for N and Av simultaneous control, specifically through PI controllers, where the experimental plant consists in a basic vapor compression refrigeration cycle and the reference temperature is set around 3.2 o C [ 22 ]. Jiangjiang et al. implement an hybrid PID-Neural Network controller, where the neural network tunes PID parameters online [ 23 ]. Meanwhile, Underwood proposes a decoupling control strategy where PID controllers are jointly tuned through optimization techniques [ 24 ]. Moving forward to non-linear control, Baeyoung et al. developed an algorithm to design sliding-mode controllers and carried some simulations [ 25 ]. Pangborn proposes a LQR MIMO controller with successfully reference tracking [ 26 ]. Regarding the robust control strategy, [ 27 ] design a MIMO H∞ controller which solves the S/KS problem. Here the authors bear in mind that, for the system they consider, the coupling between the outlet temperature of evaporator secondary flux and the valve opening is weak. The reduced order controller is compared to a SISO H∞ controller, obtained through the easing of the closed-loop system bandwidth. However, regarding the S/KS problem, only uncertainties along the frequency are considered. Regarding experimental control, most authors applied their control techniques or achievement on a basic cycle with one compression stage and one load demand, which typically is studied for cooling purposes (around 5o C), using air or water at the evaporator secondary side [ 22 , 24 , 26 , 27 ]. Senger and Trofino [ 28 ], use and experimental fridge for applied their studies of switching rules, which consists in determining the compressor speed and the distribution valve, which decides the evaporator is satisfied. They work with a refrigerator of which reference temperature are −4o C for the fresh-food evaporator and −21o C for the freezing one. This reference temperatures are similar to the one we use, however, their secondary sides of the evaporator are the refrigerator chambers, which depend on physical elements (typically food) and the air within the chambers. 2.3 Objectives •Modeling : Basic configuration of vapor compression refrigeration systems are studied before by Bejarano [ 1 ]. Extended configurations are to be studied considering up to two compressors and two evaporators. •Identification : Each element is expected to be identified separately while the global validation will be done. A novel methodology for heat transfer coefficient regarding heat exchangers. •Control : Designed controllers for simulation and then in experimentation are intended to be done. •Optimization : Variable maps for extended configuration will be carried out using as efficient index the corresponding Coefficient of Performance. 2.3 Objectives 7 Regarding modeling, switching moving boundary methods are considered to characterize the heat exchangers. Based on Figure 2.3, different modes can be established depending on the presence or absence of any fluid states. Typically the condenser works in mode 1 in standard operation or mode 2 if the subcooled zone has disappeared. The rest are intended for start-up and shut down operation which are not studied in this Thesis. For the evaporator, both modes can arise while standard operation, however the second one is undesirable since, it leads the possibility to pass liquid droplet to the compressor and cause unpleasant breakdowns. Mode 1 Mode 3 Mode 4 Mode 5 Superheated Superheated Superheated CONDENSER Subcooled Two-phase Two-phase Mode 2 Subcooled Two-phase Two-phase EVAPORATOR Superheated Two-phase Two-phase Mode 1 Mode 2 Figure 2.3 Condenser and evaporator modes according to Bin Li et al. [50]. One-compression-stage, two-load-demand refrigeration system is represented in Figure 2.4. It can be considered as an extension of the basic configuration shown in Figure 2.1, where a new expansion valve and evaporator appear. The manual pressure upholding valve creates a pressure drop to allow the control of both evaporators at two different pressure levels. In this case both cold room are particular named according to the refrigerant aims: Cooling chamber for secondary inlet temperatures around 5o C and freezing chamber when the desired secondary inlet temperatures is around −20oC. P-h diagram shown in Figure 2.5 is similar to the one for the basic configuration shown in Figure 2.2, where the second evaporator pressure line is included, as well the corresponding enthalpy line associated to the second electronic expansion valve and the upholding pressure valve. Besides, another dotted pressure line associated wit the corresponding secondary inlet temperature is added to diagram. The degree of subcooling is defined as done before, as well as the degree of superheating 8Chapter 2. Introduction Cooling chamber Evaporator Condenser Main compressor Expansion valve Secondary flux Secondary flux Freezing chamber Evaporator Expansion valve Secondary flux Pressure upholding valve S S d dd f ff l ll Figure 2.4 One-stage, two-load-demand refrigeration cycle. 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 Pc Pe5 Pe20 hc,out hc,in Te20,SH TSC he20,out he5,out Figure 2.5 One-stage, two-load-demand P-h diagram. at outlet temperature of both evaporators. However, the important one is the one related to the compressor intake, since its refrigerant inlet consists in a mix of refrigerant mass flows of both evaporators. The main objective at this configuration is to control both refrigerant power at the same 2.3 Objectives 9 time as long as the degree of superheating at the compressor intake maintains a value over zero. The two-compressor, two-load-demand configuration, depicted in Figure 2.6, maintains the same scheme than in Figure 2.4, where an additional compressor, called booster, is included in at the outlet evaporator line associated to the freezing chamber. The upholding pressure valve remains in the figure because is part of the experimental plant, but it is fully opened for this configuration. So, in this case the booster compressor manages the pressure jump between both evaporators, while the main compressor is devoted to the jump between cooling evaporator and condenser. Cooling chamber Evaporator Condenser Main compressor Expansion valve Secondary flux Secondary flux Booster compressor Freezing chamber Evaporator Expansion valve Secondary flux Pressure upholding valve Figure 2.6 Two-stage, two-load-demand refrigeration cycle. As can be seen in Figure 2.7, the P-h diagram remains similar to the previously shown, where the compressor effects over pressure and specific enthalpy are included. In this case two degrees of superheating are defined, one for each compressor intake. Having a multi-configuration purpose-built experimental plant where advances controllers,in addition to optimization techniques, as well as simulation models can be tested, results of vital importance to deep in the understanding of the multi-stage, multi-load vapor-compression refrigeration systems. However, bad dimensioning, accuracy problems derive from sensors or actuators, or several technical problems, such as a bad placement of a sensor, can lead us to an unapproachable testbed with numerous source of problems very coupled among which them. 10 Chapter 2. Introduction 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 Pe5 Pe20 hc,out hc,in hboost,dis he5,out he20,out TSH,main TSH,boost TSC hmain,suc Figure 2.7 Two-stage, two-load-demand P-h diagram. For example, a bad temperature-sensor accuracy at reasonably low heat power absorb at evaporators suppose up to a variability of 130 W of a theoretical average power of 1000 W: around a 13%. If the uncertainties are higher, the error increases enormously. 2.4 Thesis overview Chapters are distributed as follow: •Chapter 3 presents the statics models for each element involved in the one-compressor one-load-demand configuration. For further configurations, explained before, static upholding pressure valve model is added as well as some mass and energy balances required for closing iterative procedures. Theses procedures are also explained and some simulations show tracking reference results. Secondary circuits are also modeled dynamically, considering the experimental design. •Chapter 4 include all information about the design, layout and automation of the experimental plant. Characterization and calibrations are obtained for some sensors and actuators. •Chapter 5 addresses the identification procedure of the experimental plant element according to the models proposed in Chapter 3. A novel methodology to identify heat transfer coefficients at condenser and evaporators is proposed, based on some relationships between and some order of magnitude found in the literature. Steadystate experimental data allow to obtain all parameter values. Global validation is carried out using the two-compressor two-load-demand configuration. Furthermore, secondary circuits are also studied, since they play an important role at the experimental plant. Identification methodology, results, and some problems are described. 2.4 Thesis overview 11 •Chapter 6 presents some contributions in control, from the point of view of the simulations, and later the experimentation. Decentralized multivariable H∞ controller based on S/KS/T Mixed Sensitivity Problem is designed for one-compressor one-load-demand, considering several operating points. Controllability is studied before looking for the best input-output pairing and any other plant constraints. Results comparing robust, PID and model predictive controllers are depicted. •Chapter 7 describes the optimization method, based on statics models commented before, for one-compression-stage refrigeration cycle and for the extended configuration with two compressors. Some maps showing the results are depicted for both configuration. •Chapter 8 summarizes the mains conclusions and contributions and proposes some future research lines. 12 Chapter 2. Introduction 2.5 List of publications The following papers have been published during the elaboration of this Thesis: Journal papers: 1) [ 29 ] D. Rodríguez, G. Bejarano, J. A. Alfaya, M. G. Ortega, and F. Castaño, “Parameter identification of a multi-stage, multi-load-demand experimental refrigeration plant,” Control Eng. Pract., vol. 60, pp. 133–147, 2017. 2) [ 30 ] G. Bejarano, J. A. Alfaya, M. G. Ortega, and M. Vargas, “On the difficulty of globally optimally controlling refrigeration systems,” Appl. Therm. Eng., vol. 111, pp. 1143–1157, 2017. 3) [ 31 ] G. Bejarano, D. Rodríguez, J. A. Alfaya, M. G. Ortega, and F. Castaño, “On identifying steady-state parameters of an experimental mechanical-compression refrigeration plant,” Appl. Therm. Eng., vol. 109, pp. 318–333, 2016. 4) [ 32 ] G. Bejarano, J. A. Alfaya, M. G. Ortega, and F. R. Rubio, “Multivariable analysis and H∞ control of a one-stage refrigeration cycle,” Appl. Therm. Eng., vol. 91, pp. 1156–1167, 2015. 5) [ 33 ] J. A. Alfaya, G. Bejarano, M. G. Ortega, and F. R. Rubio, “Controllability analysis and robust control of a one-stage refrigeration system,” Eur. J. of Control, vol. 26, pp. 53–62, 2015. Conference papers: 1) [ 34 ] D. Rodríguez, G. Bejarano, J. A. Alfaya, and M. G. Ortega, “Robust and decoupling approach to pid control of vapour-compression refrigeration systems,” in 3rd IFAC Conference on Advances in Proportional-Integral-Derivative Control, 2018. 2) [ 35 ] J. Garrido, M. L. Ruz, F. Vázquez, J. A. Alfaya, and F. Morilla, “Decentralized PID control with inverted decoupling and superheating reference generation for efficient operation: Application to the benchmark pid 2018,” in 3rd IFAC Conference on Advances in Proportional-Integral-Derivative Control, 2018. 3) [ 36 ] G. Bejarano, J. A. Alfaya, D. Rodríguez, F. Morilla, and M. G. Ortega, “Benchmark for pid control of refrigeration systems based on vapour compression,” in 3rd IFAC Conference on Advances in Proportional-Integral-Derivative Control, 2018. 4) [ 37 ] D. Rodríguez, J. A. Alfaya, G. Bejarano, and M. G. Ortega, “Optimización y control en cascada de temperatura de recinto mediante sistemas de refrigeración” in XXXVIII Jorn. de Autom., Gijón (Spain), 2017, pp. 146–153. 5) [ 38 ] D. Rodríguez, J. A. Alfaya, G. Bejarano, M. G. Ortega, and F. Castaño, “Identificación paramétrica del condensador de una planta experimental de refrigeración,” in XXXVII Jorn. de Autom., Madrid (Spain), 2016, pp. 27–34. 2.5 List of publications 13 6) [ 39 ] D. Rodríguez, J. A. Alfaya, G. Bejarano, M. G. Ortega, and F. Castaño, “Steadystate parameter estimation of an experimental vapour compression refrigeration plant,” in Eur. Control Conf. (ECC), Aalborg (Denmark). IEEE, 2016, pp. 43–48. 7) [ 40 ] D. Rodríguez, J. A. Alfaya, G. Bejarano, M. G. Ortega, and F. Castaño, “Estimación de parámetros de una planta experimental de refrigeración,” in XXXVI Jorn. de Autom., Bilbao (Spain), 2015, pp. 951–958. 8) [ 41 ] J. A. Alfaya, G. Bejarano, M. G. Ortega, and F. R. Rubio, “Multi-operating- point robust control of a one-stage refrigeration cycle,” in Eur. Control Conf. (ECC), Linz (Austria). IEEE, 2015, pp. 3490–3495. 9) [ 42 ] G. Bejarano, J. A. Alfaya, M. G. Ortega, and F. R. Rubio, “Design, automation and control of a two-stage, two-load-demand experimental refrigeration plant,” in 23rd Mediterranean Conf. on Control and Autom., Torremolinos (Spain), 2015, pp. 537–544. 10) [ 43 ] J. A. Alfaya, G. Bejarano, M. G. Ortega, and F. R. Rubio, “Control robusto multivariable de un ciclo de refrigeración,” in XXXV Jorn. de Autom., Valencia (Spain), 2014. 3 Modelling Have no fear of perfection, you’ll never reach it Salvador Dalí Contents 2.1 Context 1 2.2 State of the art 5 2.2.1 Modeling 5 2.2.2 Control and Optimization 6 2.3 Objectives 6 2.4 Thesis overview 10 2.5 List of publications 12 In this chapter the modeling of VCRS is addressed, from the development of each element (compressor, condenser, expansion valve and evaporator) to the different configurations considering up to two compression stage and/or two load demands. Furthermore, evaporator secondary circuit are also modeling. Finally, some simulations results are shown for three different VCRS configurations: one compressor and one load demand; one compressor and two load demands, and two-compressor and two load demands. These configurations are taking into account for the experimental plant design as shown in Chapter 4. Regarding plant configurations, some different types can be found in the literature. Basic configuration with one compression stage and one load demand are widely studied in the literature. Li tests the Switching Moving Boundary models for condenser and evaporator at an experimental plant for this configuration and also for one compressor and two load demands [44]. 15 22 Chapter 3. Modelling 0 1 Condenser normalised length Temperature Temperature diagram at the condenser Tc,sec,in Tc,sec,sh,out Tc,sec,in Tc,sec,tp,out Tc,sec,in Tc,sec,sc,out Tc,in Tc,g Tc,f Tc,out Tc,sat,out ˙ Qc,sh ˙ Qc,tp ˙ Qc,sc TWO−PHASE ZONE TSC 0 1 Condenser normalised length Pressure Refrigerant pressure diagram at the condenser Pc,in Pc,g Pc,f Pc,out TWO−PHASE ZONE ζc,sh ζc,sh +ζc,tp = 1 −ζc,sc Refrigerant Secondary flux SUPERHEATED VAPOUR ZONE SUPERHEATED VAPOUR ZONE SUBCOOLED LIQUID ZONE SUBCOOLED LIQUID ZONE Figure 3.4 Temperature and pressure diagram along the condenser. Qc t p =εc t p ·(˙mcp)min ·(Tc re f r,g−Tc sec,in) Cc t p =(˙mcp)min (˙mcp)max ={phasechange}=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 =g(NTUc t p,Cc t p) Qc t p =˙mc re f r ·(hc re f r,g−hc re f r,f) Pc re f r,f=Pc re f r,g−ζc t p ·kdrop ·(˙mc re f r)2 (3.8) 3.2 Dynamics model: Condenser 23 Qc sc =εc sc ·(˙mcp)min ·(Tc re f r,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 =g(NTUc sc,Cc sc) Qc sc =˙mc re f r ·(hc re f r,f−hc re f r,out ) (3.9) Again ε=g(NTU,C) is the corresponding curve for this condenser extracted from [ 49 ]. Besides, thermodynamics properties, such as cc p or hc re f r,g are obtained again from the CoolProp tool [ 45 ]. αc sh , αc t p , and αc sc are the overall heat transfer coefficients, while kdrop is the pressure drop coefficient. Finally, Equations (3.10) and 3.11 show the refrigerant and secondary output variables, respectively. The first ones defines the degree of subcooling as commented before. The output temperature of the condenser secondary flux is calculated as the sum of the corresponding output temperatures weighted by the dimensionless length at each zone. Pc re f r,out =Pc re f r,f−ζc sc ·kdrop ·(˙mc re f r)2 Tc re f r,out =T(Pc re f r,out ,hc re f r,out ) Tc re f r,out,g=T(Pc re f r,out ,χ=0) Tc SC =Tc re f r,out,g−Tc re f r,out (3.10) Tc sec,sh,out =Tc sec,in +Qc sh ˙mc sec ·ζc sh ·cc p,sec Tc sec,t p,out =Tc sec,in +Qc t p ˙mc sec ·ζc t p ·cc p,sec Tc sec,sc,out =Tc sec,in +Qc sc ˙mc sec ·ζc sc ·cc p,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.11) These equations solve the problem when the three state zones exist at the condenser. this is the typical mode during operations. However it is also possible that the subcooled zone disappear. In this case, ζc sc =0 and Equations (3.9) becomes unnecessary. Therefore, last two of Equations (3.8) change the variables hc re f r,f and Pc re f r,f into hc re f r,out and Pc re f r,out , respectively. 3.2 Dynamics model: Condenser The system dynamics models are well expressed by Bejarano in [ 1 ], where the evaporator dynamics is disregarded due to lower volume and refrigerant density compared to the 24 Chapter 3. Modelling condenser ones which entails that the condenser dynamics are one order of magnitude lower. Furthermore, the inlet and outlet mass flow dynamics at the condenser are faster than the heat transfer dominant dynamics, disregarding the fast transient due to mass flow imbalance. The first two modes are only considered for the condenser from the ones proposed by B. Li [ 50 ], since they are the most plausible modes in steady state. According to G. Bejarano, the state vector xc is reduced to those elements of the modes 1 and 2, as shown in Equations (3.12) . Matrix Zc and forcing function vector fc structures are derived from the governing equations detailed in Modelling section of [ 1 ]. It is important to remark that Zc coefficients depend on thermodynamics variables and geometric parameters without forgetting the condenser state vector xc and the time-varying model inputs wc . On the other hand, elements of the forcing function fc strongly depend on the thermal power transferred at each zone that define the system statics. xc= [hc re f r,sh Pchc re f r,sc ζc sh ζc t p ¯ γc]T xc m1= [Pcζc sc hc re f r,sc]T xc m2= [Pcζc sh ¯ γc]T wc= [ ˙mc sec Tc sec,in ˙mc re f r,in ˙mc re f r,out hc re f r,in]T Zc(xc,wc)˙xc=fc(xc,wc) (3.12) Active variables of the condenser state vector xc depends the mode. For mode 1 state variables that evolved are stated in Equations (3.12) as xc m1 . The same occurs in mode 2, where the state vector xc m2 is used. Correspondingly, the structures of Zc and fc switching depending on the modes. Note that dimensionless length for the condenser zones ζi∀i=sh,t p,sc are interchangeable since the sum of them must be equal to 1. The dynamic models are simulated by computing the next step as ˙xi(k)=(Zc)−1(k)fc(k) at each time step kand numerically integrating to obtain the state vector xc(k+1). 3.3 Solution strategy Solution strategy is well explained by Bejarano, [ 1 ]. Since it is important for the following steps, it is summarized here. The time evolution of the simulation consist in two steps. Having a initial condition the first step is to integrate condenser dynamics model summarized in Equation (3.12) . The input vector ωc must be calculated at each instant, and, therefore, hc re f r,in also, depending on the steady-state models of the remaining elements of the cycle. The condenser pressure and the specific enthalpy at the condenser outlet are boundary conditions of this calculation. Both variable variables can be obtained independently of the condenser mode. Once determined the boundary conditions at each instant, the expansion valve model (Equations (3.1) ), the compressor model (Equations (3.3) ), and the evaporator steadystate model (Equations (3.4) – (3.5) ) constitute a nonlinear equation system with two unknown variables: Pe re f r and he re f r,out . Thus, the best way to solve them is to apply an iterative procedure (see Figure 3.5) proposed by [ 1 ]. The iterative solution gives some extra variables described in Equation (3.13) , besides the evaporator pressure and the specific 3.4 Implementation in Simulink S-functions: the Benchmark 25 enthalpy at the evaporator outlet. Finally, the variables needed to complete condenser input vector ωcare obtained in y y y, so the iterative procedure can be continued. 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 3.5 Iterative procedure to solve basic cycle variables. y y yc=hhc re f r,in he in Wcomp cons ˙mre f r TSH Te sec,out ζe t pi(3.13) 3.4 Implementation in Simulink S-functions: the Benchmark Interpreted Matlab code is a great tool for development and analysis of bulky amount of the lines, having the possibility to do any modification while the simulation is on progress. The main inconvenience is the inefficient of the execution, since the code has to be translated into machine code each time the simulation is executed. Another drawback is the difficulty to manage a large amount of lines as well as the editing time cost of reprogramming simulations. Simulink S-function blocks provides a interest tool to execute compiled codes in an user-friendly environment. Each element can be codified independently and integrated in a Simulink page. Also, the integration procedure is managed by Simulink, so the user just avoid dealing with these issues. Similar program solution can be found in the literature. Thermosys is a commercial library for Matlab-Simulink that allow to arrange different refrigeration-element blocks to create a vapor compression refrigeration system [ 51 ]. By connecting input-output variables and initializing each block Simulink is able to run a simulation of the system according to the desired control input. It is a very versatile product but it does not allow us to implement precisely our model. 26 Chapter 3. Modelling Figure 3.6 shows the Simulink implementation of the vapor compression refrigeration system models. Each element implements the equations explained in this section. Typically, program tag connects the output of the elements with the corresponding input of others elements. Some of these tag are input signals of the system. Simulink solves the iterative procedure described in Section 3.3 for each time step, as well as it handles algebraic loops and the step size for a variable-step solver. domingo, 20 de mayo de 2018 1:11 Latex página 1 Figure 3.6 Simulink scheme for the basic refrigeration cycle simulator. Talking of execution times, the refrigeration system simulator had been developed first in interpreted Matlab code, and then this closed simulation based on compiled S- function blocks was built. Comparing both running the same simulations the results are considerable. For example, in the case of a 60-minute simulation, the time expended by the interpreted code simulation could be around 20 minutes of real time. On the contrary, in the case of compiled code by Simulink results speed up to 4 real-time minutes. The Simulink S-function model has been proposed as Benchmark Challenge for control test at the 3rd IFAC Conference on Advances in Proportional-Integral-Derivative Control with a successfully results. Up to 14 conference papers have been presented using this simulator. Furthermore, 3 of them have personal contribution for this thesis [34, 35, 36]. Furthermore, this simulator has been proposed as Benchmark for the engineering control contest at the Spanish Committee of Automation (http://intranet.ceautomatica.es/og/ingenieriade-control/benchmark-2017-2018). 3.5 Secondary circuits 27 3.5 Secondary circuits Secondary circuits are designed to simulate the chambers to be cooled or frozen. It consist of a tank within a thermal resistance, that simulated the thermal load in a chamber, and a pump to move the secondary mass flow and the evaporator. They are also explained in Section 4.2.4 and after their parameter are identified in Section 5.4.1. Figure 3.7 shows the scheme of the variables involved. ˙mC sec Tsec,C out Trefr,C out PR,out Tsec,C in Tsec,C tank T ank Pump Evaporator EEV R QR Qloss Qrefr Figure 3.7 Variables involved in energy balances at secondary circuits. 3.5.1 Dynamic model The model consists in the energy balance of the tank as shown in Equations (3.14) . Qloss are the thermal losses of the system, Qre f r is the refrigeration heat, while QR stands for the thermal resistance. It is important to remark that Qre f r should be negatives in order to compensate the term in the energy balance, since QRand Qloss are always positives. The model receives as input the secondary mass flow ˙msec , and the outlet temperature of the secondary flow at the evaporator Te sec,out , while the output variable is the inlet temperature at the evaporator secondary flow Te sec,in.R(t)is the resistance control signal. Vtank ·ρe sec ·ce p,sec d(Te sec,in(t)) dt =Qre f r(t)+QR(t)+Qloss(t) Qre f r(t) = λ·˙msec(t)·[ce p,sec,out (t)·(Te sec,out (t)−Tre f )−ce p,sec,in(t)·(Te sec,in(t)−Tre f )] QR(t) = KR·R(t) Qloss(t) = σ·(Tamb(t)−Te sec,in(t)) (3.14) Where K , Vtank , α , and β are parameters, while the remain variables are thermodynamic properties of the secondary fluids, which are explained in Section 4.2.5. For the cooling fluid, Zitrec MC thermodynamics properties are obtained from the CoolProp tool, since this substance is characterized as an incompressible fluid and can be compute with its concentration in water. On the other hand, freezing substance it is not characterized in 28 Chapter 3. Modelling CoolProp tool. So, the assumption that its properties does not vary around the operation point leads to the following values for the required thermodynamics properties: •cFluidosol p=3370 W kg-1 K-1 •ρFluidosol =1046 kg m-3 3.6 Configurations Basic cycle configurations with one compressor and one load demand are well studied in the literature, providing different models [ 11 , 10 , 44 ]. Others authors try increasing the numbers of evaporators [ 52 , 53 , 54 ] adding a pressure upholding valve or just having the same pressure at every evaporators. One step further are the configurations with two compressor and two load demands as shown in [ 55 , 56 ]. This configuration makes use of the additional compressor to achieve the pressure difference between both evaporator operation point. The drawback is the addition of another constrains in the form of degree of superheating. Here, three different plant configuration are proposed: • Firstly, the basic one-compressor one-load-demand cycle has been explained before. In this case the main different is the possibility to work with two different operating point: cooling temperature, around 5o C and freezing one, at −20o C approximately. So, two simulators with the different operation point have been developed. • Next, one-compressor two-load-demand configuration is also addressed in simulation. The elements are the same as in the basic configurations. The operation point is established for 5o C and −20o C. Some additional equations have to be added to close the cycle variables, and a new iterative algorithm to solve the steady state is proposed. • Similarly, for the two-compressor and two-load-demand configuration a simulator has been developed based on the previous work. The addition of another compressor changes closing-equation of the simulator as well as the resolution of the iterative algorithm for the steady-state part. These configurations are considered on the experimental plant design explained in Chapter 4. 3.6.1 Basic cycle: one-compressor, one-load-demand cycle This configuration is explained along this chapter. In this case the same configuration has been evaluated for the two different level of temperature. Table 3.2 summarized the input and out variable range, initial condition for both secondary circuits considered. The COP is reminded in Equation (3.15) , since it is an important efficiency indicator typically used. COP =Qe Wcomp cons = ˙me re f r ·(he re f r,out −he re f r,in) ˙me re f r ·(hc re f r,in −he re f r,out )=he re f r,out −he re f r,in hc re f r,in −he re f r,out (3.15) 3.6 Configurations 29 Table 3.2 Initial operating point for basic cycles. Variable Range Initial conditions for −20oC Initial conditions for 5oC Unit Manipulated variables Av[10 ,100]50 53 % N[30 ,50]40 30 Hz Disturbances Tc sec,in [27 ,33]30 28 oC ˙mc sec [ 125 ,175 ] 150 1950 g s-1 Pc sec,in -1 1 bar Te sec,in [ −22 ,−18 ] −20 5 oC ˙me sec [55 ,75]64.5 64.5g s-1 Pe sec,in -1 1 bar Tamb [20 ,30]25 28 oC Output variables Te sec,out -−22.15 1.82 oC TSH -14.65 31.2oC Internal initial condition he,IC out -351 367.71 kJ kg -1 K -1 Simulation results for freezing temperature configuration Figure 3.9 represents the output variables when an 5% step has been carried out in the valve opening at minute 2, and a 5 Hz on the compressor at minute 10 (Figure 3.8). Figures 3.10 and 3.11 show the COP and the condenser and evaporator pressures along the simulation. Steps on the secondary pressures have omitted because of lack of interest. 30 Chapter 3. Modelling 0 5 10 15 20 40 42 44 46 Speed [Hz] Compressor speed 0 5 10 15 20 50 52 54 56 Time [min] Opening [%] Expansion valve opening Figure 3.8 Simulation of the manipulated variables for steps at Nand Ae20 v. 0 5 10 15 20 −22.4 −22.2 −22 T [ºC] Evolution of Tsec,out e20 0 5 10 15 20 12 14 16 T [ºC] Evolution of TSH e20 Figure 3.9 Simulation of the output variables for steps at Nand Ae20 v. 3.6 Configurations 31 0 5 10 15 20 0.985 0.99 0.995 1 1.005 1.01 1.015 1.02 1.025 Time [min] COP Coefficient of Performance Figure 3.10 Simulation of the COP for steps at Nand Ae20 v. 0 5 10 15 20 16.6 16.7 16.8 Pressure [bar] Condenser pressure 0 5 10 15 20 0.9 0.95 1 Time [min] Pressure [bar] Evaporator pressure Figure 3.11 Simulation of the Pressures for steps at Nand Ae20 v. 38 Chapter 3. Modelling 0 10 20 30 40 50 60 22.5 23 P[bar] Condenser pressure 0 10 20 30 40 50 60 5.8 6 P[bar] Cooling Evaporator pressure 0 10 20 30 40 50 60 2.15 2.2 2.25 2.3 Time [min] P[bar] Freezing Evaporator pressure Figure 3.21 Simulation pressures at one-compressor and two-load-demand configuration when steps are done at Nand Ae20 vand Ae5 v. 3.6 Configurations 39 3.6.3 Two-compressor two-load-demand cycle This configuration makes use of a second compressor, called booster, to manage the pressure difference between both evaporators. Although the pressure upholding valve at the evaporator refrigerant mass outlet is considered, it is set at the maximum valve opening possible in order to avoid strong pressure loss at this element. Again this configuration uses the steady-state elements explained in Section 3.1, where the expansion valve model and evaporator model are considered for colling and freezing purposes, and compressor models are used for main and booster configurations. Figure 3.22 shows the element placement for this configuration. Cooling chamber Evaporator Condenser Main compressor Expansion valve Secondary flux Secondary flux Booster compressor Freezing chamber Evaporator Expansion valve Secondary flux Pressure upholding valve Figure 3.22 Scheme for the two-compressor two-load-demand configuration. The closing equations are slightly different in this configuration. Equations (3.18) show the mass and energy balance are applied, where the refrigerant mass flows of each evaporator are obtained by the pressure upholding valve model and the booster compressor model. hMC re f r,int =he5 re f r ·˙me5 re f r +hBC re f r,disc ·˙mBC re f r,disc ˙mre f r ˙mMC re f r,int =˙me5 re f r +˙mBC re f r,disc ˙me5 re f r =fpuv(Pe5,Pmid) ˙mBC re f r,disc =fBC(hBC re f r,int ,Pe20,Pmid) (3.18) 40 Chapter 3. Modelling The iterative procedure used for this configuration again is very similar compared to the rest of configurations. Newly, for each step the condenser information is available at the beginning, Pc and hc re f r,out as well as the Manipulated variables and Disturbances - extended for this configuration. Pe20 , Pe5 , Pmid , he20 re f r,out , and he5 re f r,out form the supposed variables. The procedure uses the expansion valve models to solve both evaporator valves and the pressure upholding one, obtaining the corresponding refrigerant mass flows and specific enthalpies outlets. Next step solves compressor models for both elements to obtain the inlet and outlet specific enthalpies and the power consumed, and the mass flow balance. Finally, evaporator model and the energy balance are applied to obtain the last corresponding specific enthalpies at the element inlet. Up to five variables are double obtained in this procedure. ˙me5re f r , he5 re f r,in , he20 re f r,in , hBC re f r,int , and hMC re f r,int and the corresponding duplicates are used for correction. If the comparison is suitable, the steady-state procedure finished for the present iteration and the dynamic part can be done. In the negative case, some corrections are done on the supposed variables and the iteration is repeated. The COP for this configuration is expressed in Equation (3.19) . Moreover, the degrees of superheating for both compressor intakes are also described. COP =Qe5+Qe20 ˙ WMC cons +˙ WBC cons =... ... = ˙me5 re f r ·(he5 re f r,out −he5 re f r,in)+ ˙me20 re f r ·(he20 re f r,out −he20 re f r,in) ˙me20 re f r ·(hBC re f r,disc −hBC re f r,int )+ ˙mre f r ·(hc re f r,in −hMC re f r,int ) TMC re f r,int =T(Pmid,hMC re f r,int ) TMC SH,int =TMC re f r,int −TMC re f r,int,g TBC re f r,int =T(Pe20,he20 out ) TBC SH,int =TBC re f r,int −TBC re f r,int,g (3.19) The initial conditions used for this configuration are described in Table 3.4. Simulation results Figures 3.24 – 3.27 show the simulation results of the two-compressor two-load-demand configuration for a step in each manipulated variable. In minute 2 a step is done at the freezing expansion valve opening, after, at minute 22 the turn is for the main compressor speed, and at minute 42 the step is in the cooling expansion valve opening. Finally, at minute 62 the step is done at booster compressor speed. 3.6 Configurations 41 Table 3.4 Initial operating point for two-compressor and two-load-demand configuration. Variable Initial condition Unit Manipulated variables Ae20 v34 % Ae5 v75 % NMC 30 Hz NBC 34.9Hz Disturbances Tc sec,in 25 oC ˙mc sec 140 g s-1 Te20 sec,in −20 oC Te5 sec,in 5oC ˙me20 sec 74 g s-1 ˙me5 sec 320 g s-1 Tamb 25 oC Output variables Te5 sec,out 4.2oC Te20 sec,out −21.5oC TMC SH,int 33.6oC TBC SH,int 2oC Internal initial condition hBC,IC re f r,int 351.3 kJ kg -1 K -1 Initial suppositions hMC re f r,int 380.3kJ kg-1 Pe20 2bar he5 re f r,out 369.1kJ kg-1 Pe55.76 bar hc re f r,out 267.7kJ kg-1 Pc re f r,out 20 bar 0 20 40 60 80 30 32 34 36 38 Speed [Hz] Compressor speed 0 20 40 60 80 30 40 50 60 70 Opening [%] Time [min] Expansion valve opening Main Booster Freezing Cooling Figure 3.24 Manipulated variables simulation at two-compressor and two-load-demand configuration for steps at Nand Ae20 vand Ae5 v. 42 Chapter 3. Modelling 0 10 20 30 40 50 60 70 80 4.18 4.2 4.22 4.24 T [ºC] Evolution of Tsec,out e5 0 10 20 30 40 50 60 70 80 −21.55 −21.5 −21.45 −21.4 −21.35 T [ºC] Evolution of Tsec,e20,out 0 10 20 30 40 50 60 70 80 33 34 35 T [ºC] Evolution TSH MC 0 10 20 30 40 50 60 70 80 2 4 6 T [ºC] Time [min] Evolution of TSH BC Figure 3.25 Output variables simulation at two-compressor and two-load-demand configuration for steps at Nand Ae20 vand Ae5 v. 0 10 20 30 40 50 60 70 80 0.97 0.98 0.99 1 1.01 1.02 1.03 1.04 Time [min] COP Coefficient of Performance Figure 3.26 COP simulation at two-compressor and two-load-demand configuration for steps at Nand Ae20 vand Ae5 v. 3.6 Configurations 43 0 20 40 60 80 19.9 20 20.1 P [bar] Condenser pressure 0 20 40 60 80 5.6 5.7 5.8 P [bar] Cooling evaporator pressure 0 20 40 60 80 2 2.2 Time [min] P [bar] Freezing evaporator pressure Freezing middle Figure 3.27 Simulation pressures at one-compressor and two-load-demand configuration for steps at Nand Ae20 vand Ae5 v. 3.6.4 Two-compressor, one-load-demand cycle This configuration uses both compressors to achieve the pressure difference between the evaporator and compressor. Theoretically, using two smaller compressors could be able to achieve the same compression ratio than a bigger one. Therefore, this configuration is not further studied since there are no possibilities of comparing between compressors. 44 Chapter 3. Modelling Condenser state information Manipulated variables Disturbances Supposing Upholding pressure valve model Compressor model Expansion Valve 5 model Expansion valve -20 model Mass Flow connection Evaporator -20 steady-state model Evaporator 5 steady-state model Pe20 Pe20 Pe20 Pe20 Pe20 he20 out he20 out he20 out he20 out he20 out Pe5 Pe5 Pe5 Pe5 Pe5 he5 out he5 out he5 out he5 out hc out hc out hc out hc out Pc Pc Pc Pc Ae20 v Ae20 v Ae5 v Ae5 v NMC NMC Tamb Tamb ˙msec20 ˙msec20 ˙msec5 ˙msec5 Te20 sec,in Te20 sec,in Te5 sec,in Te5 sec,in he20 in he20 in he5 in he5 in ˙me20 ˙me20 ˙me20 ˙me5 ˙me5 ˙me5 hc in hMC in hMC in hMC in ˙ WMC he20′ in he5′ in Te20 sec,out Te5 sec,out ζe20 tp ζe5 tp eit eit y ˙mALT e5 Correction ˙mR ˙mR Figure 3.17 Iterative procedure to solve evaporator variables for the one-compressor two-load-demand configuration. 3.6 Configurations 45 Condenser state information Manipulated variables Disturbances Supposing Upholding pressure valve model Expansion Valve 5 model Expansion valve -20 model Mass Flow connection Evaporator 5 steady-state model Main Compressor model Energy Balance equation Booster Compressor model Evaporator -20 steady-state model Pe20 Pe20 Pe20 Pe20 Pe20 Pe5 Pe5 Pe5 Pe5 Pe5 he20 out he20 out he20 out he20 out he5 out he5 out he5 out hcond out hcond out hcond out hcond out Pmiddle Pmiddle Pc Pc Pc Pc Pint Pint Ae20 v Ae20 v Ae5 v Ae5 v NMC NMC NBC NBC Tamb Tamb Tamb hBC out hBC out ˙msec20 ˙msec20 ˙msec5 ˙msec5 Tsec20 in Tsec20 in Tsec5 in Tsec5 in he20 in he5 in ˙me20 ˙me20 ˙me20 ˙me20 ˙me5 ˙me5 ˙me5 ˙me5 ˙me5 hc in hMC in hMC in ˙ WMC ˙ WBC he20′ in he5′ in Te20 sec,out Te5 sec,out ζe20 tp ζe5 tp hBC in eit eit ˙m′ e5 Correction ˙mR ˙mR ˙mR Figure 3.23 Iterative procedure to solve evaporator variables for the two-compressor two-load-demand configuration. 4 Experimental plant Engineering is the art of modelling materials we do not wholly understand, into shapes we cannot precisely analyse so as to withstand forces we cannot properly assess, in such a way that the public has no reason to suspect the extent of our ignorance. Dr. A.R. Dykes Contents 3.1 Statics elements 16 3.1.1 Electronic expansion valve 16 Pressure upholding valve 16 3.1.2 Compressor 16 3.1.3 Evaporator heat transfer 17 3.1.4 Condenser heat transfer 20 3.2 Dynamics model: Condenser 23 3.3 Solution strategy 24 3.4 Implementation in Simulink S-functions: the Benchmark 25 3.5 Secondary circuits 27 3.5.1 Dynamic model 27 3.6 Configurations 28 3.6.1 Basic cycle: one-compressor, one-load-demand cycle 28 Simulation results for freezing temperature configuration 29 Simulation results for cooling temperature configuration 32 3.6.2 One-compressor two-load-demand cycle 34 Simulation results 35 3.6.3 Two-compressor two-load-demand cycle 39 Simulation results 40 47 54 Chapter 4. Experimental plant 4.2.4 Secondary tanks Tanks simulate the desired chambers to be cooled or frozen. Figure 4.8 shows an example of them, where the secondary mass flow inlet at the upper left part of the tank, while the outlet is just at the bottom. The thermal resistance (explained in 4.4.4) is placed at one third of the bottom, where the secondary mass outlet flow influenced more. As technical aspect, for each secondary circuit, the total volume of each one is around 12.34 liters. Figure 4.8 Example of tank for simulated cooling or freezing chambers at the experimental plant. 4.2.5 Fluids The refrigerant is the most important element of the cycle, since the selection of it affects to the design of the other elements. Some important point are the evaporation pressure, since it is desirable to be higher that atmospheric pressure in order to avoid entering air in the system; condensing pressure, on which depend the robustness of the system; and high latent heat, which adjust the refrigerant mass flow, among others. R404a is selected because of is appropriate for refrigeration at low and average temperatures. Comparing this refrigerant to the ones that are used for the same purposes, for example R502, the Ozone Depletion Potential presented by R404a is almost none and the Global Warming Power (GWP) is very low, improving the respective data for the R502 [57, 58]. Tanks are filled with different anti-freezing fluid. In these cases the selected one is glycol, a liquid fluid at standard conditions with interesting characteristic at operating ranges, as high thermal conductivity, low viscosity and almost constant density. For the tank that simulate freezing conditions the selected one is propylene glycol called Fluidosol of the manufacturer Termicol at 60% of concentration volume [ 59 ]; and for the cooling tank, mono-ethylene glycol called Zitrec MC of the manufacturer Arteco at 40% of concentration volume [ 60 ]. Both glycols are mixed with distilled water. The 4.3 Sensors 55 concentrations are chosen considering the corresponding freezing point under the chamber operation points. 4.3 Sensors Sensors are important part to acquire information from the plant. Then this data can be studied to perform the system dynamics in terms of transfer functions. However a bad selection of sensors, bad placement of just malfunctioning would bring several problems that will be very difficult to find. In this section, different sensors are studied to have as high accuracy as possible. 4.3.1 Resistance thermometer Resistance thermometers measure the temperature at each point of interest of the plant. The most typical is the Pt100 due to the variability electrical resistance of this material, Platinum, when the surrounding temperature change. Moreover, resistance variations are usually related linearly with temperature, thus it makes easier to work with this sensors. Some different sensors are used for this plant. First, Pt100 are placed in the positions marked in Figure 4.1, except the outlet temperatures of each evaporator secondary fluxes and each tank temperature (between the pump and the tank). Then, the outlet temperatures of each evaporator secondary fluxes have been introduced to be able to measured the temperature differences between the inlet and the outlet and therefore the refrigerant power. All this sensors are tree-wire connection. After, all sensors have been duplicated to have more accurate temperature metrics at each point. This time the selected sensors have been the Pt1000 resistance thermometer. These ones are useful to avoid parasitic resistance from wires, and also improve the measurement thanks to the four-wire connections. Finally, up to 19 extra Pt100 resistance thermometers from a different manufacturer have been used at the same key point in the plant. Hereafter, this sensors are called by TTi where i is the numeration assigned. Table 4.3 summarized sensors’ names, types, manufactures, accuracies, and number of connection wires. All sensors are connected to the Programmable Logic Controller (PLC) via ART 0814 Input/Output Card (see 4.5.2) except for two of then: Te5,Pt100 sec,out and Te20,Pt100 sec,out . This sensors use transducers which are connected to AMM 0660 Input/Output Card. ART 0814 Input/Output Cards satisfy the electric generation and measurement, while AMM 0660 Input/Output Card only measures the electric current from the transducer, so the electric power is provided via a continuous power supply. Calibration Sensor calibration has been carried out to know the accuracy errors of each sensor in order to be able to understand the process. Since there exits different resistance thermometers from different manufacturers, the same test has been done to all the sensor at the same time. This test consists in immerse all temperature sensors in the same two-phase water solution and let them reach zero Celsius degrees. A portable cooler isolated contains the 56 Chapter 4. Experimental plant Table 4.3 Technical data of the resistance thermometer used at the experimental plant. Name Sensor type Manufacturer Accuracy [ o C] Wires TBC re f r,int Pt100 Unknown 0.3 3 Pt1000 TE connectivity 0.1 4 TMC re f r,int Pt100 Unknown 0.3 3 Pt1000 TE connectivity 0.1 4 TMC re f r,disc Pt100 Unknown 0.3 3 Pt1000 TE connectivity 0.1 4 Tliquid re f r Pt100 Unknown 0.3 3 Pt1000 TE connectivity 0.1 4 Te5 re f r,out Pt100 Unknown 0.3 3 Pt100,TT5 RSC 0.1 4 Te20 re f r,out Pt100 Unknown 0.3 3 Pt100,TT15 RSC 0.1 4 Tamb Pt100 Unknown 0.3 3 Pt1000 TE connectivity 0.1 4 Tacc sec,out Pt100 Unknown 0.3 3 Pt1000 TE connectivity 0.1 4 Te5 sec,out Pt100 Unknown 0.3 3 Pt1000 TE connectivity 0.1 4 Pt100,TT11 RSC 0.1 4 Te5 sec,in Pt100 Unknown 0.3 3 Pt1000 TE connectivity 0.1 4 Pt100,TT9 RSC 0.1 4 Te5 sec,tank Pt1000 TE connectivity 0.1 4 Pt100,TT10 RSC 0.1 4 Te20 sec,out Pt100 Unknown 0.3 3 Pt1000 TE connectivity 0.1 4 Pt100,TT3 RSC 0.1 4 Te20 sec,in Pt100 Unknown 0.3 3 Pt1000 TE connectivity 0.1 4 Pt100,TT4 RSC 0.1 4 Te20 sec,tank Pt1000 TE connectivity 0.1 4 Pt100,TT6 RSC 0.1 4 water solution where the sensor are immersed. The solution consist of 4kg of ice cubes and 1l of liquid water. Data acquisition was active before the sensor have been immersed in the solution, until the ice has melted down. The length is around 15 hours of which only one and a half hours are selected as valid data (time lapse from 2 to 3.5 hours). It is important to remind that the sampling time of the data acquisition system is 1 second. Therefore, enough data is recovered to validate the experiment. Table 4.4 shows the offset values obtained from experiments. These values have been 4.3 Sensors 57 applied at the PLC in order to correct these temperature deviations. The Pt1000 sensors associated with TMC re f r,disc has been omitted for technical reasons. Most accurate sensors named as TTi are shown in the table. Nevertheless, the worst precision written is far enough from the manufacturer’s specifications. As general assessment, practically all sensor shows a higher errors than the accuracy expected by the manufacturer. Table 4.4 Experimental results for the resistance thermometer offset measured. Variable Sensor type Offset [ o C] TBC re f r,int Pt100 0.40816 Pt1000 0.5 TMC re f r,int Pt100 0.23555 Pt1000 0.08978 TMC re f r,disc Pt100 0.83255 Pt1000 * Tliquid re f r Pt100 0.39883 Pt1000 0.22166 Te5 re f r,out Pt100,OSAKA 0.5 Pt100,TT5 0.21737 Te20 re f r,out Pt100,OSAKA 0.9 Pt1000,TT15 0.60019 Tamb Pt100 1.3697 Pt1000 0.413 Tacc sec,out Pt100 0.48943 Pt1000 0.56641 Te5 sec,out Pt100 0.86091 Pt1000 0.73457 Pt100,TT11 0.3 Te5 sec,in Pt100 0.85562 Pt1000 0.58358 Pt100,TT9 0.12937 Te5 sec,tank Pt1000 0.99435 Pt100,TT10 0.51918 Te20 sec,out Pt100 1.0753 Pt1000 0.14621 Pt100,TT3 0.21128 Te20 sec,in Pt100 0.074431 Pt1000 0.14644 Pt100,TT4 0.26273 Te20 sec,tank Pt1000 0.70391 Pt100,TT6 1.1925 58 Chapter 4. Experimental plant 4.3.2 Pressure sensor These sensors measure the relative pressure along primary circuit. Their data lead one to understand the state of the refrigerant mainly at the evaporator primary outlets but also at the condenser primary inlet and outlet. The manufacturer is OSAKA and the models are PP 8H for low pressures and PP 30H for high ones. It is a two-wire connection with an accuracy of 0.5% of the range ( [0.5,8] bar for PP 8H and [0,30] bar for PP 30H). This means a maximum error of 0.04 bar in the worst case, which is enough accuracy for this process. Figure 4.9 Example of pressure sensor. 4.3.3 Flowmeter sensor Flow measurement provides important extensive information for thermodynamics. These measurements allow to study the volume passing through heat exchanger and, therefore, the energy. Due to high cost of these sensors for two-phase high-pressure refrigerant, they are only placed in the secondary flows, where glycol is easier to be measured. The selected sensors are the magnetic-inductive flow meter SM6000 of the manufacturer ifm electronic. These sensors are connected with 4 wires that allows a range of [0.1,25] L min-1 with an accuracy of 0.05 L min-1. Figure 4.10 Example of flowmeter sensor installed at the experimental plant. 4.4 Actuators 59 4.4 Actuators 4.4.1 Electronic expansion valves Electronic Expansion Valves (EEV) control the valve opening and therefore the refrigerant flow that comes from the condenser, at high pressure, to the evaporator, at lower pressure. Secondarily, the both high and low pressure are affected by the amount of mass crossing the valve. As particular matter, the refrigerant change the state from liquid one at high pressure to two-phase state a low pressure. The selected elements are relay expansion valves of the manufacturer Danfoss, type AVK10 (Figure 4.11). The valve opening ranges depend on the duty cycle, which is 6 s by default, but we have reduced it to 1 s for experimental reasons. Thus, the valve opening is [0,100] % of the duty cycle. Typically the lower value is around 30 % due to the lack of pressure at downstream. Correspondingly, the heat capacity range permit is [0.6,14] kW at an evaporating temperature range of [−50 ,60]o C which are enough for our control purposes. Figure 4.11 Electronic expansion valve installed. Super Heat Controllers or PWM signals Superheat Valve Controllers (SVC) are devices which can self-control the expansion valve opening by calculating the degree of superheating using internal reference tables of the refrigerant and two measurements. They have integrated two measure input, both for pressure and temperature at the outlet of the primary flux. They use a Proportional-Integral controller which can be adjusted internally. Figure 4.12 shows an example of the VEX-200-PU model of the manufacturer OSAKA used for both chambers. 60 Chapter 4. Experimental plant Figure 4.12 Example of Super Heat Controller of the manufacturer OSAKA. This device is useful for industrial purposes. However, since it is used for experimental problems, this system turn out to be very stiff and inaccurate. For instance, the pressuresensor accuracy is one order less at this device’s inputs than at the PLC’s ones. This means a degree-of-superheating error around 4 o C which is unacceptable. Therefore, these devices have been replaced by solid-state relays controlled with by a Pulse-Width Modulation (PWM) signals from the PLC. The duty cycle is still the same, 1s. 4.4.2 Variable-speed compressors Variable-Speed Compressors (VSC), Figure 4.13 are the most electric-power consumers in the system.The semihermetic reciprocating compressors designed by Bitzer provides enough energy to afford a compression ratio around 12 at the Main compressor. Table 4.5 summaries the technical data of both VSC. Figure 4.13 Variable-Speed Compressors (VSC) installed at the vapor-compression refrigeration system. 4.4 Actuators 61 Table 4.5 Technical data of the Variable-Speed Compressor installed at the experimental plant. Compressor Main Booster Model 4FES-3Y 2KES-0.5Y Power [kW] 6.4 1.8 Number of cylinders 4 2 Experimental Speed range [Hz] [30 ,50] [30 ,50] Experimental Pressure Intake [10 -5 · Pa] [1,3] [1,3] Experimental Pressure Discharge [10 -5 · Pa] [17 ,24] [2,3] Variable frequency drives Variable Frequency Drivers (VFD) are installed to control compressor speed as well as pump speed. The VLT r HVAC Basic Drive FC 101 drives provide an optimal control for pumps, while VLT r HVAC Drive FC 102 drives do the same for compressors. The drives, shown in Figure 4.14 communicate with the PLC using Modbus protocol. Figure 4.14 Variable Frequency Drives (VFD) installed on the experimental plant. VLT r HVAC Basic Drive FC 101 on the left and VLT r HVAC Drive FC 102 on the right. 4.4.3 Liquid pumps Pumps are used to move secondary fluxes at the evaporators, and also at the water-cooled condenser. The selected ones are in-line pumps model Omega DN25/32 2-80-2 of the manufacturer Smedegard. (Figure 4.15). The pump can be use with water and glycol fluids. Properties and data of interest are summary in Table 4.6, where some information is empirically obtained. 62 Chapter 4. Experimental plant Table 4.6 Summary of data for Smedegard in-line pumps. Variable Value Working temperature range [ o C] [−30 ,50] Volumetric flow range [L/min] [0.1,15.8043] Speed range [Hz] [0.5,45] Maximum power consumption [kW] 0,25 Figure 4.15 In-line Smedegard pump for secondary fluxes at evaporators and water-cooled condenser. 4.4.4 Thermal resistance Thermal resistances are used to simulate the load demand on each chamber. The manufacturer is Electrifor, and the model is ED002C. It is a three-phase element connected in star which provided around 1800 W of calorific power. It is controlled by a PWM signal sent from the PLC to a solid-state relay. Calibration of the resistances For control purposes and to develop an accurate model it is necessary to verify the manufacturer’s products and specifications. So, the thermal power of the resistances has been measured and Figure 4.16 shows the results. For each chamber, Wmeas is the thermal power consumption directly measured with a wattmetter. Since line-current values I and the electric resistance at each phase, Rw90 Ω , is known for both chambers, the thermal power consumption can be calculated as Equation (4.1) shows, where V is the measured phase voltage. We,resistance in f er,V1=3·Vi·Ii We,resistance in f er,V2=3·Ri·I2 i (4.1) Measured values of thermal power consumption can be considered as linear within the 4.4 Actuators 63 PWM-signal range. Furthermore, the maximum values, for 100 % of the PWM signal, is 1640 W for both chambers’ resistances, since the wattmetter resolution is around 10 W. 0 10 20 30 40 50 60 70 80 90 100 0 0.5 1 1.5 2 Thermal power [kW] Static characteristic of the thermal resistance for cooling chamber Wmeas e,resistance Winfer,V1 e,resistance Winfer,V2 e,resistance Comparison line 0 10 20 30 40 50 60 70 80 90 100 PWM signal [%] 0 0.5 1 1.5 2 Thermal power [kW] Static characteristic of the thermal resistance for freezing chamber Wmeas e,resistance Winfer,V1 e,resistance Winfer,V2 e,resistance Comparison line Figure 4.16 Static characteristic for both resistances. 4.4.5 Condenser fan The fan is the main element to control the secondary mass flow at the air-cooled condenser. It is an threephase axial fan (shown in Figure4.17) manufactured by Ebmpapst with a exterior diameter of D=0.5 m and interior of d=0.166 m considering only the blade length. The input signal is controlled using a PWM sent from the PLC. Figure 4.17 Fan installed at the air-cooled condenser. Some data have been obtained by measuring the air speed passing through the fan with a 70 Chapter 4. Experimental plant 4.6 SCADA A SCADA is developed and works as human-machine interface. It allows to control the plant manually in open loop and to monitor all cycle variables when a high-level controller is applied. Figure 4.19 shows the SCADA screen when a high-level controller, namely through OPC communication standard, is performed. Figure 4.19 SCADA screen under OPC control. As generally statement here, gray boxes are devoted to pressure sensors; red ones to temperature sensors and green to flowmeter. Pale red is used to show inferred variables related to temperature. Talking of actuators, dark blue boxes shows the measured and reference speed of compressor and pumps; Light blue ones shows the valve openings, while the orange boxes are used to the resistance power. Finally, the fan speed have its own color for the boxes. On the right and left side, buttons open different plots where analogical variables are depicted in order to supervised the experiment behavior. The center-right column of buttons are used to start up and shut down the experimental plant, as well as to manage the OPC communication. Alarm management is implemented to shows the user any possible errors that come from the elements of the plant or from control. For example, low values of Main compressor intake pressure sensors force the compressor to stop and move to an alarm screen until the problem is solved; meanwhile the active valve/s open a little bit in order to increase the commented pressure. 4.7 OPC 71 4.7 OPC OPC (OLE (Object Linking and Embedding) for Process Control) is a communication standard used to connect devices from different manufacturers [ 62 ]. Here, OPC is used for communication between the PLC and MATLAB software [ 63 ]. As explained previously, the PLC is used to data acquisition and low-level control, whereas the numerical computing program is devoted to implement high-level controllers. As seen in Figure 4.20, the PLC is connected to the computer by an Ethernet link. SCADA and 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. PlantPLC SCADA OPC server Matlab® Simulink Ethernet Computer Figure 4.20 Communication diagram. OPC client is set up in MATLAB/Simulink environment, where the number of read and written data are collected in Table 4.13. This environment allows real-time communication and processing where the sample time is up to 1 s. Actually, only about 25% of this time is busy, which would let reduce the control sample time. Nevertheless, due to technological constraints concerning the communication bus, this time has not been improved yet. Table 4.13 Process variables managed through OPC. Analogue Digital Read data 82 5 Written data 13 0 4.8 Future work: Thermal energy storage extension As future work, a design of a Thermal Energy Storage which extend the use of the experimental plant is developed. Figure 4.21 shows the experimental plant extension respect to Figure 4.1, where refrigerant lines from compressor intakes to the liquid line through 72 Chapter 4. Experimental plant condensers remain similar. For both secondary circuits, the main difference resides in the addition of a second heat exchanger called TESi (i = 1,2) and the associated electronic expansion valve at the primary flux. Some pressure and temperature sensors have been added in order to characterize the plant at each point of interest. Furthermore, a second pair of pump and flowmeter per secondary circuit have been added to be able to control the secondary mass flow through the TES. Also, some check valves avoid secondary-flow reverse circulation. FT Glycol pump 1 CV2 CV1 EEV1 Pressure Upholding valve TT PT TT TT Liquid line TES1 TES pump 1 EEV for TES 1 Tank 1 R1 FT Glycol pump 2 EEV2 TT PT TT TT TES2 TES pump 2 EEV for TES 2 Tank 2 R2 Evaporator 1 Evaporator 2 FT TT TT TT FT TT PT PT PT PT PT PT PT PT PT Figure 4.21 Thermal Energy Storage extension for the experimental plant. Internally, the Thermal Energy Storage consists in several separate pipelines crossing from the top to the base of a cylindrical tank, as shown in Figure 4.22. The tank inside is filled with an appropriate glycol solution, while there are three different types of pipelines. Two of them are used to run the plant refrigerant and the secondary flux, which come from the old part of the experimental plant. The last type is a close one where a Phase-Change Material (PCM) is introduced. This material enable the storage or the release of energy, depending on the user decision. PCM is a novel application of element that can store bulky quantities of energy while they are phase changing. Thus, the right material for the present heat-exchange process expands the possibilities to satisfy constrains or any type of cost function. 4.8 Future work: Thermal energy storage extension 73 Figure 4.22 Appearance of the thermal energy storage tank. 5 Identification So through identification you have pleasure and pain Jiddu Krishnamurti Contents 4.1 Layouts 48 4.1.1 Plant configurations 49 4.2 Static elements 52 4.2.1 Air-cooled condenser 52 4.2.2 Water-cooled condenser 52 4.2.3 Heat exchanger for evaporators 53 4.2.4 Secondary tanks 54 4.2.5 Fluids 54 4.3 Sensors 55 4.3.1 Resistance thermometer 55 Calibration 55 4.3.2 Pressure sensor 58 4.3.3 Flowmeter sensor 58 4.4 Actuators 59 4.4.1 Electronic expansion valves 59 Super Heat Controllers or PWM signals 59 4.4.2 Variable-speed compressors 60 Variable frequency drives 61 4.4.3 Liquid pumps 61 4.4.4 Thermal resistance 62 Calibration of the resistances 62 4.4.5 Condenser fan 63 4.5 Programmable logic controller 64 4.5.1 Configuration 64 75 76 Chapter 5. Identification 4.5.2 Input/Output modules 65 4.5.3 Modbus protocol 67 4.6 SCADA 70 4.7 OPC 71 4.8 Future work: Thermal energy storage extension 71 In this chapter the system modeling and the experimental plant, explained respectively in Chapters 3 and 4, merged to approximate the model parameter of each element at the vapor compression refrigeration system and at the secondary tanks. The section is divided into two parts. The first one focuses on the VCRS to address the cycle parameters of the compressor, expansion valve, condenser and evaporators. Expansion valve and compressor parameters are obtained using the least square method, while a novel method is proposed for heat-transfer coefficient identification, supposing a relationship between them. The second part deals with the secondary circuits, since it is important to understand their behavior to have a controlled scenario for control purposes. Taking into account that the developed model is intended to be used within modelbased optimization and control strategies, the moving boundary approach is used for the heat exchangers, while all other elements are statically modeled. The moving boundary approach is a first-principle methodology, therefore the parameters of each heat exchanger model must be estimated. Some authors who make use of this approach propose to use correlations to calculate the heat-transfer-related parameters [ 19 , 50 ], using a white-box model approach. However, we consider that the steady-state parameter identification based on experimental data plays a key role in such systems to propose model-based optimization and control techniques, since calculating the heat transfer coefficients based only on correlations (when available for the refrigerant and the type of heat exchanger at the specific facility) might lead to modeling errors which may be greater than the generated when identifying these coefficients according to experimental data. Most correlations are used for heat exchanger design and not for real-time simulation, and not all flow conditions (laminar or turbulent regime) are usually considered for all fluids and types of heat exchanger. Moreover, having an estimation of heat transfer coefficients according to experimental data and therefore closer to real values than those provided by correlations enables to compute in simulation a more accurate estimation of the length of the different heat exchanger zones. Those non-measurable variables, among others, determine the cycle state and they are key variables to define optimal cycles corresponding to specified cooling demands [ 30 ]. Thus, an accurate estimator of such variables is required when trying to optimally control this kind of systems. The remainder of the chapter is organized as follows: the methodology applied to VCRS parameter identification in each case is explained in 5.1. Some global validation results are collected in Section 5.1.4, also where their implications for control and optimization are remarked. Finally, in Section 5.2 the main conclusions of this part to be derived are expressed. Moreover, secondary circuits parameter identification is explained in Section 5.4.1, where the energy balance approach is explained. Then, the problem formulation and solution for each secondary circuit is explained in 5.4.2 and 5.4.3 where the least square method is applied and identification and validation experiments are shown for each 5.1 VCRS Identification 77 secondary circuit. Finally, Section 5.4.4 and 5.4.5 show, respectively, some problems during the experimentation and the conclusion derived for this part. Some contents of this chapter have been included in the publications commented in Chapter 2, Section 2, [29, 31, 37, 38]. 5.1 VCRS Identification 5.1.1 Steady-state parameter identification As stated in Chapter 3, all components involved in the two-compression-stage, two-load- demand cycle are separately modeled. Dynamic submodels for the refrigerant flow along the condenser are generated, while static models are identified for the compressors, expansion valves, and secondary fluxes at the evaporators and condenser. Electronic expansion valves The static model of the electronic expansion valve is stated in Section 3.1.1. As shown in Equation 5.1, the hydraulics formula is used to describe the mass flow characteristic of the valve. Its inputs are the pressure at the condenser outlet, Pc re f r,out , and at the corresponding evaporator, Pe re f r , as well as the specific enthalpy at the expansion valve inlet, hc re f r,out . The expanded refrigerant mass flow is calculated as an output, which matches the refrigerant mass flow at the corresponding evaporator inlet ˙me re f r,in. ˙me re f r,in =cAv eev ·Av·q2·ρc re f r,out ·(Pc re f r,out −Pe re f r)(5.1) The refrigerant density ρc re f r,out =ρ(Pc re f r,out ,hc re f r,out ) is calculated using CoolProp tool. The valve coefficient cAv eev enables 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 electronic expansion valves are very limited in open literature, as stated by Park et al. [ 64 ]. The simplest model consisting of a unique value of the valve coefficient cAv eev , despite being in general valid for single-phase flows, produces nonnegligible 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 cAv eev where its value depends only on the degree of subcooling [ 12 , 46 ]. Considering the facility under study, the component selection and the performance range of the manipulated variables enable 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 cAv eev , and it is not valid when Tc SC achieves zero. A sensitivity analysis has been conducted and its conclusions have shown that inaccuracy when estimating the refrigerant mass flow through the models 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. [ 65 ], a polynomial correlation for cAv eev is considered, where its dependence on all input parameters of the expansion valve model is highlighted in Equation 5.2. The degree of subcooling Tc SC is also included as recommended by Shanwei 78 Chapter 5. Identification et al. [ 65 ], since this dependence becomes relevant when analyzing experimental data with very low degree of subcooling. cAv eev =cAv eev(Av,∆PAv,ρc re f r,out ,Tc SC) cAv eev =10−9 d0+d1 Av 100 +d2 ∆PAv Pc re f r,out +d3 ρc re f r,out ρc re f r,out,f +d4 Tc SC Tc re f r,out,f!(5.2) In Equation 5.2 the temperature Tc re f r,out,f refers to the saturation temperature at the outlet condenser pressure Pc re f r,out , whereas the density ρc re f r,out,f matches the density of saturated liquid at the outlet condenser pressure. An identification method to calculate the coefficient cAv eev and the parameters di∀i=0,1,...,4 which define its polynomial correlation from experimental data is intended. However, since there is no refrigerant mass flow sensor at the facility under study, it is not possible to directly apply Equation 5.1. 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. Nevertheless, an estimation of this refrigerant mass flow can be obtained in steady state by applying energy balance at the corresponding evaporator, as expressed in Equation 5.3. Qe=˙me sec ce p,sec (Te sec,in −Te sec,out ) = ˙me re f r (he re f r,out −he re f r,in)(5.3) In steady state ˙me re f r,in =˙me re f r,out ≡˙me re f r , thus ˙me re f r is calculated using Equation 5.3. As an example, Table 5.1 collects some experimental data set used to estimate the refrigerant mass flow ˙me5 , while Table 5.2 does the same regarding the refrigerant mass flow ˙me20 . The specific enthalpies required in the energy balance shown in Equation 5.3 are obtained using CoolProp tool from the pressure-temperature pairs Pe re f r – Te re f r,out and Pc re f r – Tc re f r,out collected in Tables 5.1 and 5.2, since the refrigerant is single-phase at the evaporator and condenser outlet. Table 5.1 Experimental data for the refrigerant mass flow estimation concerning the expansion valve for the evaporator at 5◦C. Pe refr [bar] Te refr,out [oC] Pc refr [bar] Tc out [oC]˙me sec [kg s-1]Te sec,in [oC] Te sec,out [oC] ˙me refr [g -1] 5.42 5.0 20.97 45.5 0.32 4.9 4.1 8.2 5.63 3.6 20.97 45.7 0.33 6.3 5.1 13.0 5.58 3.2 22.01 47.2 0.32 4.9 4.0 10.4 5.57 3.1 21.08 45.6 0.32 4.9 4.0 9.4 5.60 3.3 24.06 51.0 0.32 4.9 3.9 11.3 5.56 3.3 21.07 45.5 0.32 4.9 3.9 10.3 Once estimated ˙me re f r for a variety of steady-state experimental points, a least squares method is applied to minimize the errors between refrigerant mass flow calculated from Equations (5.1) and 5.2 and the values estimated by applying Equation 5.3 for all available steady-state experimental points. Identified values of all parameters included in Equation 5.2 for both expansion valves are shown in Table 5.3, while the mean value and root mean square (RMS) of mass flow relative errors are gathered in Table 5.4. Relative error and 5.1 VCRS Identification 79 Table 5.2 Experimental data for the refrigerant mass flow estimation concerning the expansion valve for the evaporator at −20◦C. Pe20 re f r [bar] Te20 re f r,out [oC] Pc re f r [bar] Tc re f r,out [oC] ˙me20 sec [g -1]Te20 sec,in [oC] Te20 sec,out [oC] ˙me20 re f r [g -1] 1.80 -29.0 20.97 45.5 74.1 -19.6 -21.1 4.8 1.87 -27.9 20.97 45.7 74.4 -19.9 -21.0 3.5 1.88 -23.9 22.01 47.2 66.2 -20.8 -22.4 4.5 2.11 -23.2 21.08 45.6 71.5 -20.2 -21.3 3.0 1.92 -26.8 24.06 51.0 72.0 -19.9 -21.3 4.7 1.84 -27.4 21.07 45.5 70.8 -19.6 -21.4 5.2 RMS are calculated using Equations (5.4) and (5.5) , respectively, where Z is the total number of experimental data. RE ˙me re f r =|˙me re f r,pred −˙me re f r,meas | ˙me re f r,meas (5.4) RMS =v u u tZ ∑ i=1 RE2 ˙me re f r Z(5.5) Table 5.3 Identified parameters of the cAv eev correlation for valves at evaporator at 5 o C and -20oC. Expansion valve ddddd Ae5 v328.53 0.88 -4.38 -323.52 960.07 Ae20 v805.43 0.39 98.15 -892.56 2524.11 Table 5.4 Mean value and RMS of the relative error between predicted and estimated refrigerant mass flow concerning valves at evaporators at 5oC and −20oC. Expansion valve RE [%] RMS [%] Ae5 v3.72 4.50 Ae20 v11.56 16.80 Compressors Likewise, the static model of the compressors, shown in Equations (5.6) , is based on that proposed by Schurt et al. [ 12 , 46 ]. Its inputs are the intake and discharge pressures, Pcomp re f r,int and Pcomp re f r,disc , and the intake specific enthalpy, hcomp re f r,int . The compressed refrigerant flow ˙mcomp re f r , the discharge specific enthalpy hcomp re f r,disc , and the electric power consumption ˙ Wcomp cons are the calculated outputs. ˙mcomp re f r = St−cPcomp re f r,disc Pcomp re f r,int ccomp v,re f r,int,g ccomp p,re f r,int,g−1 N υcomp re f r,int (5.6a) 86 Chapter 5. Identification 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 (5.18) 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 5.4 Qualitative relationship between αe t p and ζe t p at Evaporator at 5◦C. As expected, if the two-phase zone is longer, αe t p must be smaller to effectively transfer the measured thermal power Qe t p, as indicated in Figure 5.4. Hypothesis between αe sh and αe t p Combining Equations (5.16) and (5.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 5.5 for a representative experimental point as depicted in Figures 5.3 and 5.4. The values close to limits ( ζe t p =1 and ζe t p =0 ) are not represented because it would entail considering no superheated vapor zone and no two-phase zone, respectively. As previously stated and represented by means of an electric analogy in Figure 5.2, the overall heat transfer coefficients αe t p and αe sh must be calculated considering a serial combination of all thermal resistances, as shown in Equations (5.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 (5.19) Reference values or at least orders of magnitude of convective heat transfer coefficients in the literature [ 47 ] are collected in Table 5.7, as well as the conduction thermal resistance, which has been estimated considering the thickness and conductivity of the plate material. 5.1 VCRS Identification 87 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 5.5 Qualitative relationship between the ratio αe t p/αe sh and ζe t p at evaporator at 5◦C . Table 5.7 Orders of magnitude of the thermal resistances involved in heat transfer at evaporators. Zone Two-phase tp Superheated vapor 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 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 5.5, a probable range for the length of the two-phase 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 depicted in Figures 5.3 and 5.4. In Table 5.8 average values and standard deviations of all coefficients are gathered, considering both evaporators. Note that the range for ζe20 t p causes the length of the superheated vapor 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 much greater than the others shown in Table 5.8. 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 5.3. 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 88 Chapter 5. Identification Table 5.8 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 5oCb αe5 t p 520.5 40.11 b αe5 sh 51.9 8.85 Evaporator -20oCb αe20 t p 144.8 13.8 b αe20 sh 14.87 7.26 specially when the temperature difference between the refrigerant and the secondary flux 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) steady-state operating points, where diverse states are considered. 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 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. 5.1.3 Condenser Similarly, the moving boundary 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 vapor, two-phase fluid, and subcooled liquid. The dynamic model of the condenser is also based on the methodology proposed by Rasmussen and Larsen [ 66 ]. Condenser dynamics are mainly dependent on refrigerant thermodynamic properties and its internal volume, whereas the partial thermal powers related to the three zones determine system statics. These thermal powers depend on the heat exchanger configuration and the overall heat transfer coefficients, therefore the condenser identification is also carried out in steady state. The procedure is focused on the air-cooled condenser, which has a cross-flow tube configuration as extensively in Chapter 4. In Figure 5.6 a temperature and pressure scheme of both refrigerant and secondary flux is depicted. Refrigerant pressure drop along the condenser is taken into account, since experimental data show it not to be insignificant, unlike at the evaporators. Equations (5.20) – (5.22) detail heat transfer and pressure drop distribution at the superheated vapor zone, two-phase zone and subcooled liquid zone, respectively. Note that 5.1 VCRS Identification 89 0 1 Condenser normalised length Temperature Temperature diagram at the condenser Tc,sec,in Tc,sec,sh,out Tc,sec,in Tc,sec,tp,out Tc,sec,in Tc,sec,sc,out Tc,in Tc,g Tc,f Tc,out Tc,sat,out ˙ Qc,sh ˙ Qc,tp ˙ Qc,sc TWO−PHASE ZONE TSC 0 1 Condenser normalised length Pressure Refrigerant pressure diagram at the condenser Pc,in Pc,g Pc,f Pc,out TWO−PHASE ZONE ζc,sh ζc,sh +ζc,tp = 1 −ζc,sc Refrigerant Secondary flux SUPERHEATED VAPOUR ZONE SUPERHEATED VAPOUR ZONE SUBCOOLED LIQUID ZONE SUBCOOLED LIQUID ZONE Figure 5.6 Temperature and pressure diagram along the condenser. ε=g(NTU,C) refers to the heat-exchanger specific curve used within the effectiveness- NTU method, which is once again used to compute all partial thermal powers. Qc sh =εc sh ·(˙mcp)min ·(Tc re f r,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 re f r,sh =g(NTUc sh,Cc sh) Qc sh =˙mc re f r ·(hc re f r,in −hc re f r,g) Pc re f r,g=Pc re f r,in −ζc sh ·kdrop ·(˙mc re f r)2 (5.20) Qc t p =εc t p ·(˙mcp)min ·(Tc re f r,g−Tc sec,in) Cc t p =(˙mcp)min (˙mcp)max ={phasechange}=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 =g(NTUc t p,Cc t p) Qc t p =˙mc re f r ·(hc re f r,g−hc re f r,f) Pc re f r,f=Pc re f r,g−ζc t p ·kdrop ·(˙mc re f r)2 (5.21) 90 Chapter 5. Identification Qc sc =εc sc ·(˙mcp)min ·(Tc re f r,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 =g(NTUc sc,Cc sc) Qc sc =˙mc re f r ·(hc re f r,f−hc re f r,out ) (5.22) αc sh , αc t p , and αc sc are the overall heat transfer coefficients at the superheated vapor, twophase and subcooled liquid zone, respectively. Tc re f r,g and Tc re f r,f , which are the saturation temperatures at corresponding pressures Pc re f r,g and Pc re f r,f , vary due to refrigerant pressure drop. Eventually, Equations (5.23) indicates how to calculate the refrigerant state at the condenser outlet, while in Equations (5.24) the secondary flux outlet temperature model is detailed. The secondary mass flow is assumed to be divided based on the different zone lengths, whereas the secondary flux outlet temperature Tc sec,out is computed from partial output temperatures following considering the zone length as weights. Pc re f r,out =Pc re f r,f−ζc sckdrop ˙m2 c Tc re f r,out =T(Pc re f r,out ,hc re f r,out ) Tc re f r,out,f=T(Pc re f r,f,χ=0) Tc SC =Tc re f r,out,f−Tc re f r,out (5.23) Tc sec,sh,out =Tc sec,in +Qc sh ˙mc sec ·ζc sh ·cc p,sec Tc sec,t p,out =Tc sec,in +Qc t p ˙mc sec ·ζc t p ·cc p,sec Tc sec,sc,out =Tc sec,in +Qc sc ˙mc sec ·ζc sc ·cc p,sec Tc sec,out =ζc sh ·Tc sec,sh,out +ζc t p ·Tc sec,t p,out +ζc sc ·Tc sec,sc,out (5.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 [ 47 , 48 ]. Nevertheless, there is no previous assumption of the relationship between single-phase heat transfer coefficients. Moreover, 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 5.7. 5.1 VCRS Identification 91 SuperheatedTwo-phase Tube Secondary flux Subcooled Secondary flux Tube ˙mc out hc refr,out Tc sec,sc,out Tc sec,tp,out Tc sec,sh,out Tc ref r,g Pc ref r,g ˙mc refr,in hc refr,in Tc ref r,f Pc ref r,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 5.7 Electric analogy of heat transfer at the condenser. Superheated vapor zone As indicated in Equations (5.20), energy balance in steady state at the superheated vapor zone involves that Qc sh =˙mc re f r (hc re f r,in −hc re f r,g) . The inlet pressure Pc re f r,in and the inlet temperature Tc re f r,in are measured, and the refrigerant mass flow ˙mc re f r is estimated using the compressor model. Nevertheless, refrigerant pressure distribution cannot be calculated without zone length estimations, as expressed in Equations (5.25) , where kdrop concerns the pressure drop factor along the whole condenser. Pc re f r,g=Pc re f r,in −ζc sh ·kdrop ·˙m2 c Pc re f r,f=Pc re f r,g−ζc t p ·kdrop ·˙m2 c Pc re f r,out =Pc re f r,f−ζc sc ·kdrop ·˙m2 c (5.25) Linear modelling of pressure drop involves that partial pressure drops at each zone depend on zone lengths. Therefore, Pc re f r,g and, as a result, hc re f r,g and Qc sh , cannot be computed without supposing ζc sh. For a certain value of ζc sh , according to Equations (5.20) , the efficiency at the superheated vapor zone εc sh is computed from the calculated value of Qc sh as shown in Equation 5.26. εc sh =Qc sh (˙mcp)min ·(Tc re f r,in −Tc sec,in)(5.26) To calculate (˙mcp)min it is also necessary to suppose ζc sh , since the secondary flux mass flow ˙mc sec is distributed among the zones according to their lengths. Like at the evaporator, assuming that the heat-exchanger specific curve ε=g(NTU,C) is inversely applied, the values of NTUc sh and then UAc sh can be calculated. It enables to obtain an experimental relationship between αc sh and ζc sh , which is indicated in Equation 5.27 and depicted in Figure 5.8. As for the evaporators, many steady-state points are considered for identification, but only one curve corresponding to a single steady-state experimental 92 Chapter 5. Identification point is shown in Figure 5.8 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 (5.27) 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 5.8 Qualitative relationship between αc sh and ζc sh. The dependence of UAc sh on ζc sh is explicitly indicated in Equation 5.27. As expected, greater values of ζc sh involve smaller values of αc sh , since the superheated vapor zone is longer. Two-phase zone Energy balance at the two-phase zone involves that Qc t p =˙mc re f r (hc re f r,g−hc re f r,f) . Once estimated ˙mc re f r , pressures Pc re f r,f and Pc re f r,g allow to completely determine Qc t p . As indicated in Equations (5.25) , Pc re f r,f depends explicitly on ζc t p , whereas it is also influenced by ζc sh through Pc re f r,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 , according to Equations (5.21) , the efficiency at the two-phase zone εc t p can be calculated as shown in Equation 5.28. εc t p =Qc t p (˙mcp)min ·(Tc re f r,g−Tc sec,in)(5.28) By applying the specific curve ε=g(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 5.1 VCRS Identification 93 ζc t p indicated in Equation 5.29 is obtained. αc t p =UAc t p(ζc sh,ζc t p) Ac trns f Lc 1 ζc t p (5.29) Figure 5.9 illustrates the double dependence of αc t p . Once again, although many steadystate points are considered for identification, curves relative to a single experi mental point are shown in Figure 5.9, since its purpose is only to represent qualitatively the relationship between αc t p , ζc t p , and ζc sh . On one hand, in Figure 5.9 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 5.9 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 re f r,g . Supposing a certain ζc t p , if the superheated vapor zone is longer, it involves a bigger pressure drop at this zone and thus, Pc re f r,g is lower. As stated in Equations (5.21) , the hottest point related to heat transfer at the two-phase zone is Tc re f r,g , which is lower when Pc re f r,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. 0.45 0.47 0.49 0.51 0.53 0.55 Two-phase zone normalised length 40 50 60 70 80 90 100 110 120 Heat transfer coefficient [W m-2 K-1] Qualitative relationship between αc,tp and ζc,tp ζc,sh = 0.44 0.24 0.29 0.34 0.39 0.44 Superheated vapour zone normalised length 43 43.5 44 44.5 45 45.5 Heat transfer coefficient [W m-2 K-1] Qualitative relationship between αc,tp and ζc,sh ζc,tp = 0.55 Figure 5.9 Qualitative relationship between αc t p , ζc t p , and ζc sh . On the left: for fixed ζc sh = 0.44. On the right: for fixed ζc t p =0.55.. Subcooled liquid zone Energy balance at the subcooled liquid zone imposes that Qc sc =˙mc re f r (hc re f r,f−hc re f r,out ) . Despite measuring the outlet pressure Pc re f r,out and the outlet temperature Tc re f r,out , and estimating ˙mc re f r , Qc sc cannot be calculated without estimations of ζc t p and ζc sh , since Pc re f r,f depends on both zone lengths. Supposing plausible values of ζc sh and ζc t p , according to Equations (5.22) , the efficiency at the subcooled liquid zone εc sc can be computed as shown in Equation 5.30. εc sc =Qc sc (˙mcp)min(Tc re f r,f−Tc sec,in)(5.30) 94 Chapter 5. Identification 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 5.31 and represented in Figure 5.10 arises. Similarly, for the sake of clarity, only the curves related to a single experimental steady-state point are shown in Figure 5.31. αc sc =UAc sc(ζc sc) 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 (5.31) 0.45 0.47 0.49 0.51 0.53 0.55 Two-phase zone normalised length 0 5 10 15 20 25 Heat transfer coefficient [W m-2 K-1] Qualitative relationship between αc,sc and ζc,tp ζc,sh = 0.44 0.24 0.29 0.34 0.39 0.44 Superheated vapour zone normalised length 0 5 10 15 20 25 Heat transfer coefficient [W m-2 K-1] Qualitative relationship between αc,sc and ζc,sh ζc,tp = 0.55 Figure 5.10 Qualitative relationship between αc sc , ζc t p , and ζc sh . On the left: for fixed ζc sh =0.44. On the right: for fixed ζc t p =0.55.. Firstly, in Figure 5.10 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 5.10 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 t p, and αc sc According to Equations (5.27) and (5.29) , the ratio between αc t p and αc sh turns to be a function of ζc t p and ζc sh , as represented in Figure 5.11. Similarly, the ratio between αc sh and αc sc is a function of the same variables, as shown in Equations (5.27) and (5.31) . This ratio is also shown in Figure 5.11. Note that the depicted surface refers to a single representative steady-state experimental point. According to Figure 5.11, the initial guess about dominance of αc t p over αc sh and αc sc is confirmed. Regarding single-phase heat transfer coefficients, the ratio between αc sh and αc sc shows to be dependent on the zone lengths. When ζc sc =1−ζc sh −ζc t p gets close to zero, the ratio takes values less than one, while for greater ζc sc the ratio reaches values up to 15. In Figure 5.12 a detailed view of Figure 5.11 for small values of ζc sc is shown, where ratios αc sh/αc sc below one can be observed. Like at evaporators, the overall heat transfer coefficient at each condenser zone is a result of serial combination of all involved thermal resistances, as indicated in Equations 5.1 VCRS Identification 95 Qualitative relationship between the ratio αc,tp/αc,sh, ζc,tp, and ζc,sh ζc,tp ζc,sh 0 0.55 3 0.525 6 αc,tp/αc,sh 0.5 9 0.475 0.44 12 0.42 0.4 0.38 0.45 0.36 0.34 ζc,tp Qualitative relationship between the ratio αc,sh/αc,sc, ζc,tp, and ζc,sh ζc,sh 0 0.34 5 10 αc,sh/αc,sc 15 0.39 0.55 0.525 0.5 0.44 0.475 0.45 Figure 5.11 Qualitative relationships between the ratios αc t p/αc sh , on the left and αc sh/αc sc , on the right; and zone lengths ζc t p and ζc sh. αc,sh/αc,sc ζc,tp Qualitative relationship between the ratio αc,sh/αc,sc, ζc,tp, and ζc,sh ζc,sh 0.5 0.43 1 0.4325 1.5 0.435 2 0.55 0.5475 0.4375 0.545 0.5425 0.44 0.54 Figure 5.12 Detailed view of the qualitative relationship between the ratio αc sh/αc sc , ζc t p , and ζc sh, for small values of ζc sc. (5.32) and schematically represented in Figure 5.7. 1 αc sh =1 λc sh +Rc cond +1 λc sec 1 αc t p =1 λc t p +Rc cond +1 λc sec 1 αc sc =1 λc sc +Rc cond +1 λc sec (5.32) Reference values of convective heat transfer coefficients for the refrigerant (superheated vapor, two-phase and subcooled liquid zones) and secondary flux [ 47 , 48 ] are stated in Table 5.9, along with the calculated conduction thermal resistance. 102 Chapter 5. Identification 5 10 15 20 25 0 0.05 0.1 0.15 0.2 0.25 Npump e5 [Hz] Static characteristic of the cooling evaporator pump msec e5 [kg s−1] Ttank e5 = 5ºC Ttank e5 = 10ºC Ttank e5 = 15ºC Ttank e5 = 20ºC Straight line Figure 5.17 Static characteristic between the pump and the glycol in the tank at 5oC. Control Design Control has been developed to reject disturbances, since the aim is to maintain the mass flow while the experiment is carrying out. Equation (5.34) shows the first-order model identification, where Kmean is the mean value of the static gains shows in Figure 5.16, and τmean is the mean value of the time constant shows in Figure 5.16; sstands for the Laplace’s variable. Gpump(s) = Te5 tank(s) Ne5 pump(s)=Kmean τmean ·s+1=10.21 ·10−3 2.93 ·s+1(5.34) The proposed controller shown in Equation (5.35a) is an ideal-form PI. The control have been designed using the root-locus method, where the desired close-loop transfer function has been stablished as shown in 5.35b. Ccontroller(s) = Ne5 pump(s) eTsec,e5 tank (s)=Kc·1+1 TI·s(5.35a) Gdesired C−L(s) = Ne5 pump(s) Tsec,e5,re f tank (s)=1 3·s+1(5.35b) The main limitation here is the sampling time of which frequency is 1 [rad s -1 ], which has been considered in the design. Furthermore, the controller parameters have been readjusted experimentally. The selected controller’s parameter are in table 5.13. 5.3 Pump identification and control 103 Table 5.13 Pump controller parameters for secondary circuit at 5oC. Parameter Value Kc98 [kg s-1 Hz-1] TI3[s] Experimental Validation Figure 5.18 depicted the experimental validation for secondary 5-Celsius-degree mass-flow controller designed in the section before. This experiment shows how the tank temperature variation does not affect the mass-flow control and the pump speed has plenty of capacity to support stronger disturbances. The standard deviation is around 1.5% respect to the reference. 25 40 60 80 100 120 140 160 10 12 14 16 Npump e5 [Hz] Experimental data of mass flow controller at cooling evaporator 25 40 60 80 100 120 140 160 60 62 64 66 68 70 msec e5 [g s−1] 25 40 60 80 100 120 140 160 −3 0 3 6 9 12 15 Time [min] Ttank e5 [ºC] msec e5 msec e5,ref Figure 5.18 Experimental validation of Ne5 pump designed controller. 5.3.2 Pump at -20oC Similarly, the control design for the -20-Celsius-degree mass flow is designed considering almost the same details. There exist some difference, since fluids are different and the operating point is lower in this case. Experimental data We have carried out some experiment at the same secondary circuit with to obtain variations on the mass flow when a step at the pump speed is done. The intake pump temperature, Te20 tank is maintained at the same value during each step. This test has been done in step of 5 Hz from 5 to 25 Hz. Also, it is repeated for temperatures from −25o C to 5 , in steps of 5 . 104 Chapter 5. Identification The refrigeration cycle works during all the experiment to maintain the secondary flux outlet temperature and, therefore, the inlet one. Figure 5.19 shows an example of data obtained and processed for identification purposes, where temperature approximately remains at the same value; steps are done in the pump speed and the mass flow reacts as a first order model. 42.9 43 43.1 43.2 43.3 43.4 43.5 −10.3 −10.2 −10.1 −10 −9.9 Experimental data for freezing secondary pump identification Ttank e20 [ºC] 42.9 43 43.1 43.2 43.3 43.4 43.5 10 20 30 40 Npump e20 [Hz] 42.9 43 43.1 43.2 43.3 43.4 43.5 0 50 100 150 200 Time [min] msec e20 [kg/s] Figure 5.19 Example of experimental data used for pump identification at secondary circuit at -20oC. First order model identification The data have been identified as a first-order model with minimum delay. The pump has been study for the range [−25 ,5]o C at Te20 tank , however the lower-limit range has been modified to −20oC because the plant response is not appropriated. Figures 5.20 shows the trends for the static gain, K , and time constant, τ for the identification of a first order transfer function. The static gain increases linearly as long as the -20-Celsius-degree pump increases. Likewise, the static gain increases when the temperature does, which means that exists a decrement of the secondary flux density. 5.3 Pump identification and control 105 10 15 20 25 30 35 40 2 3 4 5 6 7 8 9 x 10−3 Npump e20 [Hz] Static gain at Npump e20 steps and different intake temperatures K [kg s−1 Hz−1] Ttank e20 = −20ºC Ttank e20 = −15ºC Ttank e20 = −10ºC Ttank e20 = −5ºC Mean value 10 15 20 25 30 35 40 1.6 1.8 2 2.2 2.4 2.6 2.8 3 Npump e20 [Hz] Time constant at Npump e20 steps and different intake temperatures τ [s] Ttank e20 = −20ºC Ttank e20 = −15ºC Ttank e20 = −10ºC Ttank e20 = −5ºC Mean value Figure 5.20 Static gain and time constant through Ne20 pump and Tsec,e20 tank . Time constant shows a more stable variation around 2.15 seconds. Again, delay represents a important fact in the control. Here it is caused by the PLC communication process (for more information see Section 4.5.3). Aforementioned table 5.12 shows the average delay time. In any case, the system is considered as a first-order model including the delay at the time constant. Static characteristic Figure 5.21 shows the static characteristic of the pump, where the the temperature has more influence on the mass flow than in the 5-Celsius-degree circuit. 10 15 20 25 30 35 40 0 0.05 0.1 0.15 0.2 0.25 Npump e20 [Hz] Static characteristic of the freezing evaporator pump msec e20 [kg s−1] Ttank e20 = −20ºC Ttank e20 = −15ºC Ttank e20 = −10ºC Ttank e20 = −5ºC Straight line Figure 5.21 Pump static characteristic of the secondary circuit at −20oC. 106 Chapter 5. Identification Control Design As done before, this control has been developed to reject disturbances, since the aim is to maintain the mass flow while the experiment is on progress. Equation (5.36) shows the first-order model identification, where Kmean is the mean value of the static gains shows in Figure 5.20, and τmean is the mean value of the time constants show in Figure 5.20; s stands for the Laplace’s variable. Gpump(s) = Te20 tank(s) Ne20 pump(s)=Kmean τmean ·s+1=5.7612 ·10−3 2.15 ·s+1(5.36) The proposed controller shown in Equation (5.37a) is a PI one in ideal form. Control has been designed using the root-locus method, where the desired close-loop transfer function has been established as shown in 5.37b. Ccontroller(s) = Ne20 pump(s) eTe20 tank (s)=Kc1+1 TI·s(5.37a) Gdesired CL (s) = Ne20 pump(s) Te20,re f tank (s)=1 (3·s+1)(5.37b) The main limitation here is the sampling time of which frequency is 1 [rad s -1 ], which has been considered in the design. Furthermore, the controller parameters have been readjusted experimentally. The selected controller’s parameter are in table 5.13. The main limitation here is the sampling time of which frequency is 1 [rad s -1 ], which has been considered in the design. Moreover, the controller parameters have been readjusted experimentally. The selected controller’s parameter are in table 5.14. Table 5.14 Pump controller parameter for secondary circuit at -20oC. Parameter Value Kc147 [kg s-1 Hz-1] TI3[s] Experimental Validation Figure 5.22 depicted the experimental validation for secondary -20-Celsius-degree massflow controller designed in the section before. This experiment shows how the tank temperature affects the compressor speed and the controlled mass flow. This disturbance is well rejected by the controller as seen at the compressor speed variations. The standard deviation is around 2% respect to the reference. 5.4 Secondary circuits identification 107 80 85 90 95 100 22 24 26 28 Npump e20 [Hz] Experimental validation of mass flow controller at freezing evaporator 80 85 90 95 100 48 50 52 54 56 58 msec e20 [g s−1] msec e20 msec e20,ref 80 85 90 95 100 −13 −12 −11 −10 −9 Time [min] Ttank e20 [°C] Figure 5.22 Experimental validation of Ne20 pump designed controller. 5.4 Secondary circuits identification Energy-balance identification at the secondary circuits are addressed here with the aim of characterize its dynamics. As a general assumption, the same procedure is applied to both circuit, highlighting the main differences in each one. 5.4.1 Tank schemes As explained before, tanks simulate the desired secondary conditions at the evaporator for cooling and freezing purposes. The secondary circuit is formed by a tank, an electric resistance, a pump, and the secondary circuit of the heat exchanger. The capacity of each tank is around 13 litres approximately where the inlet is in the upper part and the outlet lower part. Inside the tanks the electric resistances are placed which can heat up to 1.64 kW. The pump is a in-line glanded one with a speed range from 1 Hz to 45 Hz. The counter-current brazed-plate heat exchanger is used for transferring heat power. Temperature and flow sensors are placed along the secondary flux. Flow sensor is situated at the outlet of the secondary flux of the heat exchanger. Three thermal-resistance sensors are placed along the pipes of the secondary circuit: one at the outlet of the tank, one at the inlet and another at the outlet of the secondary flux of the heat exchanger. As general description, Figure 5.23 shows the variable names of data-acquisition measurements obtained from experimental data; the power balance on the tank scheme is plotted as well. Both tanks have the same sensors and actuators. 108 Chapter 5. Identification ˙me sec Te sec,out Te refr,out Pe refr,out Te refr,in Te tank T ank Pump Evaporator EEV R QR Qloss Qrefr Figure 5.23 Variables involved in energy balances at secondary circuits. Tank identification provides a simple model to understand its behaviour. The energy balance shown in Eq. (5.38) is proposed. Vtank ·ρe sec ·ce p,sec d(Te sec,in(t)) dt =Qre f r(t)+QR(t)+Qloss(t) Qre f r(t) = λ·˙msec(t)·[ce p,sec,out (t)·(Te sec,out (t)−Tre f )−ce p,sec,in(t)·(Te sec,in(t)−Tre f )] QR(t) = KR·R(t) Qloss(t) = σ·(Tamb(t)−Te sec,in(t)) (5.38) Where QR(t),Qloss(t),Qre f r(t) are the thermal powers involved in the tank, and which sign are depicted in Figure 5.23. Note that although Qre f r(t) is positive according to the scheme, its value must be negative (according to the temperature difference) in order to balance the equation, and taking into account that QR(t),Qloss(t) are always positive terms. The rest variables are measured by temperature, volumetric flow, and pressure sensors. R(t) is the resistance input signal in percentage. λ , and σ are the proposed identification parameters that should be always zero or positive, while K is a resistance parameter defined in Section 4. Finally, Tre f is the tabulated reference temperature at which the specific heats at constant pressure are obtained (According to CoolProps [45]). 5.4.2 Cooling-demand tank identification Experimental data from secondary circuit at 5oC Figure 5.24 shows the experimental data selected for the chamber at 5 o C involved in estimation and validation procedures, respectively. For both experiments, the mass flow controller designed in Section 5.3 is applied to maintain the reference mass flow ; the resistance reference is kept constant during the experiment. A cascade control has been implemented, according to the scheme 5.25. The main aim, which is to control tank temperature Te5 tank , is done by the external control which gives a cooling power reference. This reference is used by the internal controller to obtained the desired cooling power by attacking the valve opening Ae5 vrelated to the 5oC evaporator. Between both controllers there exist a change variable from Qe5 re f r to Te5 sec,out . 5.4 Secondary circuits identification 109 45 50 55 60 65 70 75 80 85 −5.5 −3.5 −1.5 0.5 2.5 4.5 T [ºC] Experimental data for cooling tank identification 45 50 55 60 65 70 75 80 85 55 65 75 85 Av e5 [%] 45 50 55 60 65 70 75 80 85 1.2 1.4 1.6 Q [kW] Time [min] Tsec,out e5 Ttank e5 Qrefr e5 QR e5 80 90 100 110 120 −8 −6 −4 −2 0 2 T [ºC] Experimental data for cooling tank validation Tsec,out e5 Ttank e5 80 90 100 110 120 60 80 100 Av e5 [%] 80 90 100 110 120 1.2 1.4 1.6 Q [kW] Time [min] Qrefr e5 QR e5 Figure 5.24 Experimental data to estimate the parameters. Both controllers are PID ones whose parameters have been empirically set aiming the best settling time and less overshooting possible. Te5,ref tank Te5 tank Te5 tank Te5 sec,out Te5 sec,out Qe5,ref refr Te5,ref sec,out Ae5 v Tank Experimental Internal V ariable External Control Control Plant Change Figure 5.25 Cascade control scheme implemented for controlling Ttank at both chambers. For both experiments, the reference value for ˙me sec is 0.0645 [kg s -1 ], while the resistance reference has been set to 77%. A two-Celsius-degree step change are done. Results As stated in Section 5.1.1, the least square methods is used here to identify the parameters of the secondary tank at 5 o C, given that the tank energy balance model is affine in the parameters to be estimated. In this case identification includes the dynamics since, all variables can be measured via sensors or obtained via CoolProps tool. Thus, a convenient formulation of the model is expressed in Equation (5.39) , where θ= [λ σ]T is the parameter vector to be estimated. According to Equations (5.38) , the experimental data generates vector Γ(i) and matrix ϕ(i) as detailed in Equations (5.40) , whereas ξ(i) is the residue vector for experimental point i∀i=1,2, ..., Z . Time delay, td=75 s, has been estimated to Tsec,e5 tank due to mass flow residence time along the secondary circuit. Γ=ϕ·θ Γ(i) = ϕ(i)·θ+ξ(i)∀i=1,2,...,Z(5.39) 110 Chapter 5. Identification Γ="ρe5 sec(i)·V·Ce5 p,tank(i)·Te5 tank(i)−Te5 tank(i−1) ts −KR·Re5(i)# ϕ= [ϕ(i)ϕ(i)] ϕ(i) = Tamb(i)−Te5 tank(i−td) ϕ(i) = me5 sec(i)·ce5 p,sec,out (i)·(Te5 sec,out (i)−Tre f )−ce5 p,sec(i)·(Te5 tank(i−td)−Tre f ) θ=λ σ i∀i=1,2, ..., Z (5.40) where ts is the sampling time. In this case the time dependence of both ρ and cp is neglected due to their small variations. Furthermore, here is important to evaluate each specific heat at constant pressure for each temperature, Te5 tank and Te5 sec,out , since the consideration of a cp evaluated at a mean temperature implies considerable errors around 80 W. The reference temperature, Tre f is set at 30o C according to thermodynamic tables used [45]. As indicated in Equation 5.41, the specific values of the parameters b θ are calculated by minimizing the quadratic norm of the residue vector ξ , where Γ∈R2Z×1 and ϕ∈R2Z×4 , being Z the total number of experimental data. Note that ϕ+ refers to the pseudoinverse of matrix ϕ. b θ=argmin θ||ξ||2→b θ=ϕ+·Γ(5.41) Table 5.15 shows parameter results after applying the least square method to an experimental data shown in Figure 5.24. The parameter α has been disregarded due to an inconsistency in the solution, since it should be always zero or greater. Table 5.15 Parameter estimation results for the secondary circuit at tank at 5ºC. Parameter Value Standard deviation [%] λ-[W K−1]- σ1.0087 1.532 Results are shown in Figures 5.26 and 5.27, where the first one cope with the identification problem, and the second one is devoted to validation. Figure 5.26 compares the experimental result used in this case to the simulation data obtained applying the least squares method. The estimated parameters gathers the dynamics information of the plant. 5.4 Secondary circuits identification 111 50 55 60 65 70 75 80 85 1 1.5 2 2.5 3 3.5 4 4.5 Time [min] T [ºC] 5ºC tank parameter identification Ttank sec,C5,est Ttank sec,C5,exp Figure 5.26 Identification results for the cooling-secondary-circuit dynamics. Validation results are shown in 5.27. Here the estimated temperature tracks the dynamic of the experimental one. 80 90 100 110 120 −1.5 −1 −0.5 0 0.5 1 1.5 2 Time [min] T [ºC] Cooling tank validation Ttank e5,est Ttank e5,exp Figure 5.27 Validation Results for freezing-demand secondary circuit. 118 Chapter 6. Control First order model identification 104 Static characteristic 105 Control Design 106 Experimental Validation 106 5.4 Secondary circuits identification 107 5.4.1 Tank schemes 107 5.4.2 Cooling-demand tank identification 108 Experimental data from secondary circuit at 5oC108 Results 109 5.4.3 Freezing-demand tank identification 112 Experimental data 112 Results 112 5.4.4 Glycol substitution 114 5.4.5 Conclusions for tank identifications 114 Th is chapter is devoted to study the control of the vapor compression refrigeration systems either in simulation and experimentation, depending on the existing testbed configurations: one stage and one load demand (cooling or freezing purposes); one stage and two load demands; and two compression stages and two load demands. The main objective regarding experiment plant available is focused in control experiments that improves the Coefficient of Performance (COP) of the plant. Some contents of this chapter have been included in the publications commented in Chapter 2, Section 2, [32, 33, 34, 37, 41, 43]. 6.1 Input-output pairing The Relative Gain Array (RGA) method is used find out the most suitable input-output matching. Considering a basic vapor-compression-cycle configuration, the system can be identified as linear model as shown in Equation (6.1) , where a static gain, ki j , a dominating zero, τzi j , and two poles (one fast, τfi j , and another one slower, τsi j ) are considered for every transfer function Gi j(s)∀i,j=1,2. ∆Te sec,out (s) ∆TSH (s)=G G G(s)∆Av(s) ∆N(s) G G G(s) =   k11(τz11 s+1) (τf11 s+1)(τs11 s+1) k12(τz12 s+1) (τf12 s+1)(τs12 s+1) k21(τz21 s+1) (τf21 s+1)(τs21 s+1) k22(τz22 s+1) (τf22 s+1)(τs22 s+1)  (6.1) The evaluation of the relationship between the system static gains and the closed-loop ones is proposed in Equation (6.2) , where kCi j ∀i,j=1,2 represent the closed-loop static 6.2 Robust control 119 gains and λis the relative gain matrix. k11 k12 k21 k22=−0.017 0.0016 −0.24 0.15  "1 kC11 1 kC12 1 kC21 1 kC22 #= k11 k12 k21 k22−1!T λ λ λ="k11 kC11 k12 kC12 k21 kC21 k22 kC22 #=1.1773 −0.1773 −0.1773 1.1773  (6.2) As stated by Bristol, the best input-output matching is that achieving the elements closest to 1 in matrix λ . Therefore, the RGA method confirms that the best input-output pairs are {Te sec,out ,Av} and {TSH ,N} . However, according to the static-gain elements obtained in 6.2, Av has a great influence on TSH , thus the partial decoupling strategy is expected to reduce this issue. This conclusion is extended for both secondary circuits. 6.2 Robust control Being aware of the complexity of the process, it could be appropriated to implement a controller which could cope with the coupling and uncertainties that may exist when modelling these systems linearly, due to the strong nonlinearities. These uncertainties should be studied not only along the frequency, but also at different operating points which could be reached. According to this, the development of a robust controller seems suitable. This paper proposes one in particular: a multivariable centralised H∞ controller, based on the S/KS/T Mixed Sensitivity Problem. Firstly, in order to describe the fundamental dynamic behaviour of the system, multivariable linear models are identified at several operating points by means of step response, selecting one of them as nominal model, which minimises the uncertainty region. Secondly, a controllability analysis of the nominal linear model is carried out to show some constraints in its performance. Finally, a systematic methodology for designing the weighting matrices of the S/KS/T Problem is detailed. This methodology has been successfully applied by the authors in other thermodynamic systems [ 67 ], proving to be simple and effective in improving the closed-loop dynamical performance. 6.2.1 Process model According to Section 3, the condenser dynamic models and evaporator static one are based on the Switching moving boundary method [ 50 ] and [ 68 ]; the other elements are statically modeled as well. This method splits the heat exchanger length in three zones, known as superheated vapor, two-phase fluid and subcooled liquid. Each part works as a heat exchanger with their particularities. The contribution to heat transfer and dynamic evolution of each zone depends on existence or absence of them, and they are called modes (See Section 2.2.1 for information). For example mode 1 at the condenser means that all zones are available, while mode 2 implies that the subcooled liquid zone has disappeared. 120 Chapter 6. Control Due to technological requirements, specially at the compressor, whose inlet flow should never be two-phase, it is desirable that both evaporator and condenser work in the mode that allows it [ 50 ]. A view in the space of manipulated variables of the combinations of inputs which lead to this combined behaviour in steady state is depicted in Figure 6.1. Figure 6.1 Working area in terms of control variables. On the working area defined by Figure 6.1 a few operating points have been selected, where the time response of the output variables with appropriate size steps in the manipulated variables is simulated. In Table 6.1 the operating points are detailed. Table 6.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 In order to describe the fundamental dynamic behaviour of the system, some linear models have been identified at the different operating points by means of step response. Figure 6.2 shows a typical step response (specifically at OP 4 ). Because the Te sec,out step 6.2 Robust control 121 response with a change in 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 with a change in N and it is not visible in upper-left subplot in Figure 6.2), every model is expressed in the continuous transfer matrix form shown in Equation (6.3). 19 21 23 25 27 29 5 7 9 11 13 15 TSH [oC] Degree of superheating at evaporator outlet 19 21 23 25 27 29 2400 2650 2900 N [rpm] 19 21 23 25 27 29 −20.76 −20.74 −20.72 −20.7 −20.68 −20.66 Te sec,out [oC] Outlet temperature of evaporator secondary flux Time [min] 19 21 23 25 27 29 2400 2650 2900 N [rpm] 19 21 23 25 27 29 10 12 14 16 18 20 TSH [oC] Degree of superheating at evaporator outlet 19 21 23 25 27 29 34.5 36.5 38.5 40.5 Av [%] 19 21 23 25 27 29 −20.8 −20.75 −20.7 −20.65 −20.6 T [oC] Outlet temperature of evaporator secondary flux Time [min] 19 21 23 25 27 29 34.5 36.5 38.5 40.5 Av [%] Figure 6.2 Identified nominal operating point (OP 4 ): (a) correspond to G11 , (b) to G12 , (c)to G21 and (d)to G22. ∆TSH (s) ∆Te sec,out (s)=G(s)∆N(s) ∆Av(s) G(s) =   k11(τz11 s+1) (τf11 s+1)(τs11 s+1) k12(τz12 s+1) (τf12 s+1)(τs12 s+1) k21(τz21 s+1) (τf21 s+1)(τs21 s+1) k22(τz22 s+1) (τf22 s+1)(τs22 s+1)  (6.3) The selected nominal operating point is OP 4 , since it minimises the uncertainty region. The numerical values of ki j,τzi j ,τfi j and τsi j are collected in Table 6.2. 122 Chapter 6. Control Table 6.2 Numerical values of identified transfer matrices. Point kK rpm-1 K%-1 K rpm-1 K%-1 τz[s] τf[s] τs[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  6.2.2 Controllability analysis In this section, a basic input-output controllability analysis [ 69 ] on the minimal linear model is made in order to check limitations of the closed-loop performance. The results are done for the nominal operating point and they are presented next. First of all, the system has to be adjusted in scale, in agreement with the maximum allowed deviation of each input and output. For this application, values shown in Equation (6.4) have been chosen as maximum variations. ∆N=337 rpm,∆Te sec,out =0.05 oC, ∆Av=10 % ,∆TSH =4oC.(6.4) Taking these values into account, the system can be scaled according to Equation (6.5) , where G(s) is the transfer matrix identified in the previous section and ˆ G(s) is the scaled system. The scaling matrices Du and Derr have been calculated as indicated in Equation 6.2 Robust control 123 (6.6). ˆ G(s) = D−1 err ·G(s)·Du(6.5) Derr =∆TSH,max 0 0∆Te,max sec,out =4 0 0 0.05  Du=∆Nmax 0 0∆Av,max =337 0 0 10 (6.6) The first step of the analysis consists in studying the poles and zeros of the system (as explained in [69]). Table 6.3 shows these values in the continuous time domain. Table 6.3 Poles and zeros of the nominal model. Roots Poles Zeros Value −7·10−2−4.51 ·10−2±9.45 ·10−3j −5·10−2−1.224 ·10−2 −4·10−2−1.224 ·10−2 −3.5·10−2−1.1·10−1 −1.5·10−21.21 ·10−1 −1.5·10−1 −1·10−1 −3.7 It can be seen that the system is stable, which does not imply any strong constraint in the closed-loop behaviour. Nevertheless, the system has a transmission zero placed in the right-half plane (RHP), which imposes some limitations in the system performance [ 70 ]. The RHP-zero is located at zRHP =1.21 ·10−1 . As stated in [ 69 ], this yields the upper bound to the achievable control bandwidth ωcr shown in Equation (6.7) , which gives the approximate lower bound for the rise time of at least one output specified in Equation (6.8). ωcr <|zRHP| 2.8=4.32 ·10−2rad s-1 (6.7) tr≈π 2ωcr >36.4s(6.8) Due to the multivariable system, it is also important to consider the output direction of the system for the RHP-zero. In this case, zRHP has the output direction shown in Equation (6.9). TSH Te sec,out =0.97 0.245 (6.9) 124 Chapter 6. Control It can be seen 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 [ 71 ], its natural effect is focused on the first variable. This way, the control of the degree of superheating at evaporator outlet is more limited than the control of the output temperature of evaporator secondary flux. Therefore, if a tight control on the degree of superheating is imposed, the behaviour of the output temperature of the secondary flux would be strongly degraded. The analysis of the singular values of the nominal model as a function of frequency is depicted in Figure 6.3, and Table 6.4 shows their directions at low frequency [ 69 ]. It can be observed that the minimum singular value, used as controllability index, is greater than 0dB for all frequencies below 10−1rad s-1 . Therefore independent output changes can be made until that frequency. From it, the system outputs become dependent. However, due to the limitation shown in Equation (6.7) and (6.8) , both outputs are independent for all allowed frequencies. 10−3 10−2 10−1 100101 10−3 10−2 10−1 100 101Singular values of the nominal model Magnitude [dB] σmax σmin 10−3 10−2 10−1 100101 0 5 10 Frequency [rad s−1] Condition number Relationship between maximum and minimum singular values Figure 6.3 Singular values (upper figure) and Condition number (lower figure) of the nominal model. Table 6.4 Directions of the singular values at low frequencies. Singular value Input direction [N Av]T Output direction [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 The condition number for all considered frequencies is also shown in Figure 6.3. It can be observed that from frequencies around 6·10−1rad s-1 , because the drop-off in the minimum value is more pronounced than the one in the maximum value, the condition number experiences a significant increase. Therefore, the system can be considered 6.2 Robust control 125 ill-conditioned from these frequencies on, making it strongly sensitive to uncertainties. This fact gives an additional approximate lower bound for the rise time. This bound is specified in Equation (6.10) , and it is shown to be less restrictive than the bound stated in Equation (6.8). tr≈π 2·6·10−1>2.634 s(6.10) 6.2.3 H∞controller synthesis A centralised multivariable robust controller is proposed, taking into account that the system must work properly at different operating points. It is difficult to find out the way variations in each model parameter affect the system outputs. Therefore, it is not reasonable to synthesise a controller based on parametric uncertainty, as usual in very complex systems. Thus, an H∞ controller, based on unstructured uncertainty, is designed, using the Mixed Sensitivity Problem approach. The feedback H∞ controller design can be formulated as an optimization problem, which can be posed under the configuration shown in Figure 6.4. where the optimal H∞ control problem consists of computing a controller in such a way that the energy ratio γ between z and ω is minimised. This optimal problem is still open, but a solution exists for the suboptimal case [ 72 , 73 ]. Therefore, the value of γ is decreased as much as possible through an iteration procedure. ω = r GP(s) y+ -e vu K(s) Ws(s) WKS(s) WT(s) G(s) Z1 Z3 ZZ2 Figure 6.4 S/KS/T Mixed Sensitivity configuration. The configuration used for building up the generalised plant is given by the S/KS/T Mixed Sensitivity Problem [ 74 ], which is shown in Figure 6.4. In this case, the expression of the resulting closed-loop transfer matrix, Tzω(s) , is shown in (6.11), where S0(s) is the output sensitivity transfer matrix, K0(s) the controller transfer matrix, and T0(s) output complementary sensitivity transfer matrix; and their expressions are collected in (6.12). The factors WS(s) , WT(s) and WKS(s) constitute the weighting matrices, which allow to 126 Chapter 6. Control specify the range of relevant frequencies for the corresponding closed-loop transfer matrix. Tzw(s) =   WS(s)S0(s) WKS(s)K(s)S0(s) WT(s)T0(s) (6.11) S0(s) = (I+G(s)K(s))−1 K(s)S0(s) = K(s)(I+G(s)K(s))−1 T0(s) = G(s)K(s)(I+G(s)K(s))−1 (6.12) As it is known, proper shaping of T0(s) is desirable for tracking problems, noise attenuation, and robust stability regarding multiplicative output uncertainties. Moreover, a convenient shaping of S0(s) allows to improve the performance of the system. In addition, matrix K(s)S0(s)helps avoid some numerical problems in the synthesis algorithms. Therefore, since the controller is obtained from the generalised plant, the synthesis problem with this configuration is reduced to the design of some appropriate weighting matrices which will impose the control specifications. Based on this, the generalised plant can be built up as shown in Figure 6.4, and consequently the controller can be calculated on a computer by using a synthesis algorithm. Design of weighting matrices The selection of the weighting matrices WS(s) and WT(s) has been accomplished following the design rules proposed in [ 75 ] and [ 67 ]. Thus, once the scaled nominal model is available (see Equation (6.5) ), the maximum singular values of the multiplicative output uncertainty can be estimated as shown in Equation (6.13) , where ˆ G(jω) is the scaled nominal model frequency response and ˆ GOPi(jω) stands for the system frequency response at the operating point OP i ( i=1,2,...,8 ). These maximum singular values are shown in Figure 6.5. σ(ˆ EOPi(jω)) = σ(( ˆ GOPi(jω)−ˆ G(jω)) ˆ G(jω)−1)(6.13) Next, WT(s) is designed as a square diagonal matrix with all its diagonal elements with the same transfer function, as shown in Equation (6.14) , where the transfer function WTdiag (s) must be stable, minimum phase, with a high gain at high frequencies, and with magnitude greater than the maximum singular value of the uncertainty previously computed, for each non-nominal model and frequency, as stated in Equation (6.15) . In this application, WTdiag (s) has been chosen as specified in 6.16. The magnitude of this transfer function is also plotted in Figure 6.5. It can be seen that it satisfies all the requirements imposed on it. WT(s) = WTdiag (s)·I2x2(6.14) |WTdiag (jω)|⩾σmax(ˆ EOPi(jω)) ∀ω,∀OPi(6.15) 6.2 Robust control 127 10−3 10−2 10−1 100 −20 −15 −10 −5 0 5 10 15 20 25 30 Frequency [rad s−1] |WT|, Maximum Singular Value [dB] |WT| as uncertainty upper bound σ max(ÊOP1) σ max(ÊOP2) σ max(ÊOP3) σ max(ÊOP5) σ max(ÊOP6) σ max(ÊOP7) |WTdiag(j ω)| Figure 6.5 WT(s) matrix as an upper bound of the maximum singular values of the multiplicative output uncertainties. WTdiag (s) = 31.773 ·s+0.79433 0.00565 ·s+1(6.16) WS(s) is taken as a square diagonal matrix of transfer functions, as stated in Equation (6.17) , where each diagonal element WSi(s) is designed according to the Equation (6.18) . WS(s) = "WSTSH (s)0 0WSTe sec,out (s)#(6.17) WSi(s) = αi·s+10(κi−1)·ωT s+βi·10(κi−1)·ωT ,i=TSH ,Te sec,out (6.18) ωT is the crossover frequency of WTdiag (s) and its value is about 1.92 ·10−2 rad s -1 . Parameters αi and βi are the transfer function gains at high and low frequencies, respectively. According to [ 75 ], the following values have been chosen: αTSH =αTe sec,out =0.5 and βTSH =βTe sec,out =10−4. Finally, dimensionless parameters κTSH and κTe sec,out determine the corresponding transfer function bandwidth, and their values must be higher as the desired response speed of the corresponding output increases. In this design values of κTSH =0.95 and κTe sec,out =1.05 are used. Then, the weighting design is completed after making weight WKS(s) equal to the identity matrix in order to avoid numerical problems in the synthesis algorithm. Once the weighting matrices have been designed, the generalised plant shown in Figure 6.4 can be built up, and therefore, the controller can be computed.