Dynamische thermodynamische modellering van koelmachines

Dynamische thermodynamische modellering van koelmachines Brecht Jacobs

Promotor: prof. dr. ir. Michel De Paepe Begeleiders: Bernd Ameel, Kathleen De Kerpel Masterproef ingediend tot het behalen van de academische graad van Master in de ingenieurswetenschappen: werktuigkunde-elektrotechniek

Vakgroep Mechanica van Stroming, Warmte en Verbranding Voorzitter: prof. dr. ir. Roger Sierens Faculteit Ingenieurswetenschappen en Architectuur Academiejaar 2010-2011

Dynamische thermodynamische modellering van koelmachines Brecht Jacobs

Promotor: prof. dr. ir. Michel De Paepe Begeleiders: Bernd Ameel, Kathleen De Kerpel Masterproef ingediend tot het behalen van de academische graad van Master in de ingenieurswetenschappen: werktuigkunde-elektrotechniek

Vakgroep Mechanica van Stroming, Warmte en Verbranding Voorzitter: prof. dr. ir. Roger Sierens Faculteit Ingenieurswetenschappen en Architectuur Academiejaar 2010-2011

Dankwoord Graag had ik enkele mensen bedankt die hebben bijgedragen tot het tot stand brengen van deze thesis. Dank aan mijn promotor prof. dr. ir. Michel De Paepe en mijn begeleiders ir. Bernd Ameel en ir. Kathleen De Kerpel, voor de tijd die ze voor mij hebben vrijgemaakt, de goede raad en de deskundige begeleiding. Ik kon steeds bij hen terecht voor vragen en problemen. Ook ben ik mijn familie erkentelijk voor hun steun en begrip gedurende het academiejaar en gedurende de hele opleiding. Ten slotte bedank ik mijn vrienden die van het universiteitsleven een mooie en onvergetelijke tijd hebben gemaakt.

De auteur geeft de toelating deze masterproef voor consultatie beschikbaar te stellen en delen van de masterproef te kopiëren voor persoonlijk gebruik. Elk ander gebruik valt onder de beperkingen van het auteursrecht, in het bijzonder met betrekking tot de verplichting de bron uitdrukkelijk te vermelden bij het aanhalen van resultaten uit deze masterproef.

Gent, 6 juni 2010 Brecht Jacobs

The author gives permission to make this master dissertation available for consultation and to copy parts of this master dissertation for personal use. In the case of any other use, the limitations of the copyright have to be respected, in particular with regard to the obligation to state expressly the source when quoting results from this master dissertation.

Gent, 6 June 2010 Brecht Jacobs

’Dynamische Thermodynamische Modellering van Koelmachines’ Brecht Jacobs

Promotor: prof. dr. ir. Michel De Paepe Begeleiders: ir. Bernd Ameel, ir. Kathleen De Kerpel Masterproef ingediend tot het behalen van de academische graad van Master in de Ingenieurswetenschappen: werktuigkunde-elektrotechniek Vakgroep Mechanica van Stroming, Warmte en Verbranding Voorzitter: prof. dr. ir. Roger Sierens Faculteit Ingenieurswetenschappen en Architectuur Academiejaar 2010-2011

Samenvatting Zowel in industriële als residentiële toepassingen worden koelmachines steeds vaker gebruikt. Er moet bijgevolg grote aandacht besteed worden aan goed ontwerp. Goed ontwerp omvat enerzijds de energie-efficiëntie en anderzijds het correct functioneren van de installatie. Door wijzigende koellast, wijzigende buitencondities of ten gevolge van overdimensionering bij het ontwerp, werkt een koelcyclus echter zeer zelden in vollast. Beide aspecten moeten aldus gegarandeerd worden zowel in vollast als in deellast. In de literatuur zijn al enkele cyclusmodellen ontwikkeld. Deze gaan er echter vaak van uit dat er voldoende gegevens beschikbaar zijn wat betreft de geometrie en de sturing van de componenten. Vaak is dit confidentiële informatie, zodat deze modellen niet zomaar gebruikt kunnen worden. In dit huidig werk wordt een theoretisch model ontwikkeld dat in staat is zowel het statisch als het dynamisch gedrag van een koelcyclus te voorspellen. Er zijn slechts weinig details gekend over de te simuleren koelcyclus, bijgevolg moeten er extra veronderstellingen en vereenvoudigingen gemaakt worden. Ook wordt er in beperkte mate gebruik gemaakt van data-fitting. Toch is het model in staat zowel het statisch als het dynamisch gedrag kwalitatief goed te voorspellen. De nauwkeurigheid van de simulatie is echter niet altijd even hoog. Er wordt onderzocht wat de oorzaak is van deze afwijkingen en er worden voorstellen gedaan om hieraan tegemoet te komen. Zo wordt er aangetoond waar de kritische punten bij cyclusmodellering liggen. Trefwoorden Koelcyclus, simulatie, theoretisch model, dynamisch, steady-state.

Dynamic Thermodynamic Modeling of a Refrigeration Cycle Brecht Jacobs Supervisors: Michel De Paepe, Bernd Ameel, Kathleen De Kerpel Abstract— A lumped-parameter theoretical model of a vaporcompression refrigeration cycle is developed. First, each individual component is modeled by simplifying the one-dimensional formulation of the conservation laws of mass, momentum and energy. Subsequently, these models are combined to form a simulation model of a full cycle. Given the lack of detailed information of the compressor and the valve, the modeling involved data-fitting. The results of the simulation are compared to the data obtained from a 90 ton (300 kW) centrifugal water chiller. Both steady-state as dynamic behavior are generally well predicted. However, the accuracy of the predicted COP is rather low. Also, the pressure in both the condenser and the evaporator are not very well predicted. Lastly, the transient behavior of the mass flow rate is not well captured during large disturbances. The cause of these discrepancies is analyzed and explained, resulting in suggestions for further research. Keywords—Refrigeration cycle, simulation, theoretical model, dynamic, steady-state.

I. I NTRODUCTION Experimental studies to examine the behavior of a refrigeration cycle are usually very costly and take a lot of time to perform. Also, the obtained results are not readily extendable to any other type of machine. Therefore, more and more computer models are being developed to predict the performance of chillers. These models are valuable tools to optimize the energy-performance of a system. Also, simulations serve well to analyze the design of a cycle. Third, designers can use it to implement and test control algorithms for the expansion valve and the compressor. Finally, this model will aid researchers to gain insight into the physical behavior of a refrigeration cycle under various operating conditions. The behavior of individual components is very well understood. However, combined as a system, some phenomena are not yet explained. II. T EST STAND DESCRIPTION Comstock and Braun [1] collected experimental data of a 90 ton (300 kW) centrifugal water chiller over a wide range of operating conditions. Here, a short overview is given of the test stand and its instrumentation. More information can be found in the above mentioned report. The chiller is a McQuay PEH048J, using R134a as refrigerant and water as secondary fluid. The evaporator and the condenser are both flooded shell-and-tube heat exchangers. Furthermore, the system is equipped with a pilot-driven thermostatic expansion valve and a centrifugal compressor. In both the condenser and the evaporator, water flows through the tubes, while the refrigerant is on the shell-side. At the outlet of the condenser, some of the refrigerant is led through a cooling line, first expanding across an fixed-size orifice to evaporator pressure and then used to cool the motor and its transmission. Afterwards, the refrigerant is mixed at the evaporator inlet with the stream that went through the thermostatic expansion valve. Capacity

control is achieved by varying the angle of the inlet guide vanes at the inlet of the compressor. Both transient as steady-state data were collected beginning with start-up, followed by 27 different operating conditions ending with shut-down. The different conditions were obtained by varying the evaporator water outlet set temperature and the condenser and evaporator water inlet temperature. III. M ODEL DEVELOPMENT A. Cooling line The cooling line is not modeled. After all, the amount of heat produced by losses in the transmission and the motor are several orders of magnitude smaller than the heat exchanged in evaporator and condenser. B. Thermostatic expansion valve Due to the lack of constructional details, the pilot-driven valve was simplified to a standard thermostatic expansion valve. As in [1], the flow through the valve is considered isenthalpic and stationary. Furthermore, the mass flow rate is expressed by the standard Bernoulli orifice equation. The flow area is derived from a force balance on the spindle, while neglecting its inertia. The dynamic response of the bulb is modeled by a 1st order differential equation, resulting from the conservation of energy, while treating the bulb and its contents as a single region, perfectly isolated from its surroundings. C. Centrifugal compressor The dynamics of the compressor are neglected. The mass flow rate through the compressor is computed as in [2]. Next, the compression is modeled as a polytropic process. Finally, the flow through the inlet guide vanes is modeled as an equivalent throttling process. Therefore, the pressure drop can be related to the mass flow rate by the Bernoulli orifice equation. D. Evaporator and condenser Both heat exchangers were modeled following the movingboundary approach as outlined in [3], [4] and [5]. The heat exchanger is subdivided into maximum three zones of which the boundaries coincide with the interfaces where the refrigerant transitions between single- and two-phase flow. Next, the conservation laws of mass, momentum and energy are simplified and integrated in each zone, yielding a set of 1st order linear differential equations.

E. System E.1 Dynamic model The dynamic model is obtained by combining the models of each individual component. This results in a system differential algebraic equations of index 1. E.2 Steady-state model The steady-state model is derived by setting all time derivatives in the dynamic model to zero, yielding a non-linear set of algebraic equations. The model is completed by adding the charge equation, which expresses the total amount of refrigerant charge contained in the system. IV. R ESULTS AND DISCUSSION

flow rate. As a consequence, the predicted mass flow rate is smaller than the measured values in order order to increase the pressure rise across the compressor. As the enthalpy at the inlet (hevap,in ) and the outlet (hevap,out ) of the evaporator are well predicted, this causes the underprediction of the cooling load at higher mass flow rate. The motor power is calculated by assuming a constant electro-mechanical efficiency (ηem ). However, this seems to be a bad approximation. From the measured data it was found that when the system works in part-load, the efficiency is approximately 62 %. At full-load however, the efficiency increases up to 72 %. As a higher COP corresponds to a higher part-load ratio, it is readily shown that this assumption also contributes to the underprediction of the high-range COPs.

A. Steady-state results

B. Dynamic results

The steady-state results show reasonable agreement with the measured data. However, as depicted in figure 1, the model fails to predict the higher-range COP. The lower-range COP on the other hand, is predicted within the ± 10 % bounds. This discrepancy is shown to be due to the underprediction of the cooling load (Qevap ) and the overprediction of the motor power (Pmotor ). These are defined as:

As can be seen on figure 2, the dynamic behavior of the system is well predicted. There is no overshoot nor undershoot and the duration of the simulated and the measured transient agree.

mass flow rate

predicted [−]

predicted [kg/s]

1.5 1 0.5 1

1.5 2 measured [kg/s]

2 120

8 7.5 7 6.5

3.1

3.2 time [s]

3.3

40

measured predicted 3.4 4

x 10

0.5

80 60

1

predicted measured

6 5.5 3

100 1.5 PLR [%]

mass flow rate [kg/s]

8.5

3

3.1

3.2 time [s]

3.3 4

x 10

20

3

3.1

3.2 time [s]

3.3 4

x 10

Fig. 2. Dynamic results compared to measured data

Figure 2 shows that the mass flow rate reacts too quick under large disturbances. This results from neglecting the inertia of the actuator of the inlet guide vanes and from assuming that the time-constant of the bulb is very small.

5 4 3 2

2.5

2

3 4 measured [−]

motor power

5

400 predicted [kW]

100 80 60 40 40

60 80 measured [kW] predicted

100 + 10%

V. C ONCLUSION A lumped parameter model for a vapor-compression refrigeration system has been developed. The predicted data is compared to the measured data. The cause of the discrepancies is thoroughly analyzed and explained, highlighting the critical points of a simulation model.

cooling load

120 predicted [kW]

x 10

6

2

20 20

(2)

COP = cooling load/motor power

2.5

0 0.5

(1)

pressure in the condenser [Pa]

Qevap = m ˙ (hevap,out − hevap,in ) 1 m ˙ (hcomp,out − hcomp,in ) Pmotor = ηem

5

9

300 200 100 0

R EFERENCES 0

100 200 300 measured [kW] measured

400

− 10%

Fig. 1. Steady-state results compared to measured data

The measured data show that a higher COP corresponds to a higher mass flow rate (m). ˙ At higher mass flow rate, the boiling process in the evaporator transitions from pool boiling to flow boiling. For evaporation on the outside of horizontal plain tube bundles under similar circumstances, Thome [6] has shown that, on average, the convective boiling heat transfer accounts for 40 % of the total heat transfer. The simulation on the contrary is carried out using a correlation which only predicts the nucleate boiling contribution, resulting in an underprediction of the heat transfer and thus the evaporator pressure at higher mass

[1] M. C. Comstock and J. E. Braun (1999). Experimental data from fault detection and diagnostic studies on a centrifugal chiller. Technical report, Ray W. Herrick Laboratories, Purdue University. [2] J. E. Bourdouxhe, M. Grodent en J. Lebrun (1996). HVAC1KIT - a toolkit for primary HVAC system energy caclulation. Universit de Lige. [3] N. B. O. L. Pettit, M. Willatzen en L. Ploug-Sørensen (1998). A general dynamic simulation model for evaporators and condensers in refrigeration. part ii: simulation and control of an evaporator. International Journal of Refrigeration, 21(5):404–414. [4] M. Willatzen, N. B. O. L. Pettit en L. Ploug-Sørensen (1998). A general dynamic simulation model for evaporators and condensers in refrigeration. part i: moving-boundary formulation of two-phase flows with heat exchange. International Journal of Refrigeration, 21(5):398–403. [5] W. J. Zhang en C. L. Zhang (2006). A generalized moving-boundary model for transient simulation of dry-expansion evaporators under larger disturbances. International Journal of Refrigeration, 29(7):1119–1127. [6] J. R. Thome (2010). Engineering Data Book III. Wolverine Tube, Inc.

Inhoudsopgave

i

Inhoudsopgave 1 Inleiding

1

2 De subkritische compressie-koelcyclus 2.1 Inleiding . . . . . . . . . . . . . . . . . 2.2 Algemene werking en begrippen . . . . 2.2.1 Werking . . . . . . . . . . . . . 2.2.2 Begrippen . . . . . . . . . . . . 2.3 De koelcyclus uit de proefstand . . . . 2.3.1 Opbouw . . . . . . . . . . . . . 2.3.2 Instrumentatie . . . . . . . . . 2.3.3 Testcondities . . . . . . . . . . 3 Literatuurstudie 3.1 Inleiding . . . . . . . . . . . . . . . . . 3.2 Thermostatische expansieklep . . . . . 3.2.1 Inleiding . . . . . . . . . . . . . 3.2.2 Klassieke methode . . . . . . . 3.2.3 Stromingsmechanische methode 3.3 Centrifugaalcompressor . . . . . . . . 3.3.1 Inleiding . . . . . . . . . . . . . 3.3.2 Steady-state modellen . . . . . 3.3.3 Dynamische modellen . . . . . 3.4 Verdamper en condensor . . . . . . . . 3.4.1 Inleiding . . . . . . . . . . . . . 3.4.2 Dynamische modellen . . . . . 3.5 Systeem . . . . . . . . . . . . . . . . . 4 Model koelmachine 4.1 Inleiding . . . . . . . . . . . . 4.2 Koellijn . . . . . . . . . . . . 4.3 Thermostatische expansieklep 4.3.1 Inleiding . . . . . . . . 4.3.2 Model . . . . . . . . . 4.4 Centrifugaalcompressor . . . 4.4.1 Inleiding . . . . . . . . 4.4.2 Model . . . . . . . . . 4.5 Verdamper en condensor . . . 4.5.1 Inleiding . . . . . . . . 4.5.2 Model . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . .

3 3 3 3 4 4 4 7 7

. . . . . . . . . . . . .

9 9 9 9 10 11 13 13 14 15 16 16 16 21

. . . . . . . . . . .

23 23 23 23 23 23 26 26 26 28 28 29

Inhoudsopgave 5 Steady-state model 5.1 Inleiding . . . . . . . . . . . . . . . . . . . . 5.2 Modelvergelijkingen . . . . . . . . . . . . . 5.3 Oplossingsmethodiek . . . . . . . . . . . . . 5.3.1 Systeemmodel . . . . . . . . . . . . 5.3.2 Componentmodellen . . . . . . . . . 5.4 Resultaten . . . . . . . . . . . . . . . . . . . 5.4.1 Nauwkeurigheid simulatieresultaten 5.4.2 Deellastgedrag . . . . . . . . . . . . 5.4.3 Opmerkingen . . . . . . . . . . . . . 5.5 Conclusie . . . . . . . . . . . . . . . . . . .

ii

. . . . . . . . . .

48 48 49 49 49 51 53 56 58 60 63

. . . . . . . .

65 65 65 66 66 69 69 71 74

7 Conclusie 7.1 Conclusie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2 Aanbevelingen voor verder onderzoek . . . . . . . . . . . . . . . . . . . . . . . .

75 75 75

A Proefstand

77

6 Dynamisch model 6.1 Inleiding . . . . . . . . . . . . . . . . . . . . 6.2 Modelvergelijkingen . . . . . . . . . . . . . 6.3 Oplossingsmethodiek . . . . . . . . . . . . . 6.4 Resultaten . . . . . . . . . . . . . . . . . . . 6.4.1 Uitvoeringssnelheid . . . . . . . . . . 6.4.2 Nauwkeurigheid simulatieresultaten 6.4.3 Schakelen tussen modes . . . . . . . 6.5 Conclusie . . . . . . . . . . . . . . . . . . .

. . . . . . . . . .

. . . . . . . .

. . . . . . . . . .

. . . . . . . .

. . . . . . . . . .

. . . . . . . .

. . . . . . . . . .

. . . . . . . .

. . . . . . . . . .

. . . . . . . .

. . . . . . . . . .

. . . . . . . .

. . . . . . . . . .

. . . . . . . .

. . . . . . . . . .

. . . . . . . .

. . . . . . . . . .

. . . . . . . .

. . . . . . . . . .

. . . . . . . .

. . . . . . . . . .

. . . . . . . .

. . . . . . . . . .

. . . . . . . .

. . . . . . . . . .

. . . . . . . .

. . . . . . . . . .

. . . . . . . .

. . . . . . . . . .

. . . . . . . .

. . . . . . . . . .

. . . . . . . .

. . . . . . . . . .

. . . . . . . .

. . . . . . . . . .

. . . . . . . .

. . . . . . . . . .

. . . . . . . .

B Dynamische zone-specifieke behoudsvergelijkingen verdamper en condensor 82 B.1 Verdamper . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 B.2 Condensor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 C Statische zone-specifieke behoudsvergelijkingen verdamper en condensor C.1 Verdamper . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . C.2 Condensor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

93 93 96

D Wiskundige technieken 100 D.1 De regel van Leibniz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 D.2 Het algoritme van Newton-Raphson [27] . . . . . . . . . . . . . . . . . . . . . . 100 E Werking simulatieprogramma 101 E.1 Werking steady-state simulatieprogramma . . . . . . . . . . . . . . . . . . . . . 101 E.2 Werking dynamisch simulatieprogramma . . . . . . . . . . . . . . . . . . . . . . 101 Bibliografie

106

Lijst van figuren

108

Lijst van tabellen

109

Nomenclatuur

iii

Nomenclatuur Afkortingen COP Coefficient Of Performance PLR Part Load Ratio

Symbolen A sectie Cp specifieke warmtecapaciteit bij constante druk Cv klepcoëfficiënt η rendement φ diameter F zone-grootte γ polytrope exponent of lokale void-fractie γ¯ zone-gemiddelde void-fractie g gravitationele versnelling h enthalpie k thermische conductiviteit L lengte µ dynamische viscositeit m ˙ massadebiet p druk P vermogen q warmteflux Q warmtevloed ρ densiteit τ tijdsconstante t tijd T temperatuur U tipsnelheid v absolute snelheid x dampfractie x toestandsvector y klepheffing

Subscript 2 secundair flu¨ıdum bulb reservoir cond condensor cwi inlaat koelwater evap verdamper ewi inlaat ijswater

m2 J/kg.K − − m − − − m2 /s J/kg W/m.K m P a.s kg/s Pa W W/m2 W kg/m3 s s K m/s m/s − m

Nomenclatuur f g in L out r set t TP TP − L TP − V V V − TP

iv gesatureerde vloeistof gesatureerde damp inlaat L-zone uitlaat koelmiddel instelpunt buis TP-zone interface tussen TP- en L-zone interface tussen TP- en V-zone V-zone interface tussen V- en TP-zone

Superscript e verdamper c condensor b reservoir

Hoofdstuk 1. Inleiding

1

Hoofdstuk 1

Inleiding Door de toename van de warmtedissipatie van kantoorapparatuur en de grotere bewustwording van thermisch comfort van gebouwgebruikers, is het bijna onvoorstelbaar dat een nieuw gebouw ontworpen wordt zonder koelmachines. Ook in veel industriële toepassingen is koeling onontbeerlijk. Er moet bijgevolg aandacht besteed worden aan goed ontwerp. Goed ontwerp omvat enerzijds de energie-efficiëntie en anderzijds het goed functioneren van de machine, dit zowel bij vollast als bij deellast. Een koelmachine werkt immers zo goed als nooit in vollast. Er is altijd een zekere afwijking van ontwerpomstandigheden. Dit kan te wijten zijn aan wijzigende koellast, wijzigende buitencondities, maar ook aan een zekere graad van overdimensionering bij het ontwerp.

Figuur 1.1: Koelmachine Een koelmachine kan zowel experimenteel als numeriek bestudeerd worden. Experimentele studies zijn kostelijk en vergen veel tijd en mankracht. Gestimuleerd door de sterk toenemende rekenkracht en de relatief lage kost van computers, wordt steeds vaker gebruik gemaakt van numerieke modellen. De componenten en de sturing kunnen op deze manier snel, eenvoudig en goedkoop geoptimaliseerd worden. Deze numerieke modellen zijn ook zeer algemeen en kunnen vrij gemakkelijk uitgebreid worden om ook het gedrag van andere types cycli te beschrijven. Er bestaan zowel dynamische als statische modellen. Een statisch model is voldoende om het gedrag van een koelmachine te bestuderen wanneer de koellast weinig variëert in de tijd. Dit is bijvoorbeeld het geval bij koeling van server rooms. De machine kan dan als quasistationair beschouwd worden. Dit is echter niet altijd het geval. Er zijn ook toepassingen waarbij de randvoorwaarden van de machine zeer snel wijzigen in de tijd, waardoor de inertie van de koelcyclus een belangrijke invloed heeft op het gedrag. Ook om het opstarten en uitzetten van de installatie te beschrijven, is een dynamisch model vereist. In de literatuur zijn

Hoofdstuk 1. Inleiding

2

al enkele modellen ontwikkeld om het dynamisch en het statisch gedrag van een koelcyclus te bepalen. Heel wat studies gaan er echter van uit dat er voldoende gegevens beschikbaar zijn wat betreft de geometrie en de sturing van de componenten. Vaak is dit confidentiële informatie zodat deze modellen niet zomaar gebruikt kunnen worden. Er moeten immers extra veronderstellingen en vereenvoudigingen gemaakt worden. De doelstelling van deze thesis omvat het ontwikkelen van een theoretisch model dat in staat is zowel het statisch als het dynamisch gedrag van de subkritische compressie-koelcyclus te voorspellen vetrekkend van een zeer beperkte hoeveelheid constructieve data. Vervolgens wordt dit ge¨ımplementeerd in MATLAB. Ten slotte wordt nagegaan in welke mate het bekomen model in staat is om de meetresultaten uit Comstock en Braun [11] te reproduceren. Er is tevens aandacht besteed aan de modulariteit en de gebruiksvriendelijkheid van het simulatieprogramma. Dit biedt de mogelijkheid om het cyclusmodel in een latere fase uit te breiden of om bepaalde componenten te vervangen of te verbeteren. In hoofdstuk 2 wordt kort de subkritische compressie-koelcyclus besproken en worden de belangrijkste begrippen nog eens opgefrist. In hetzelfde hoofdstuk wordt ook de koelmachine voorgesteld die gemodelleerd zal worden in de daarop volgende hoofdstukken. Vervolgens wordt er in hoofdstuk 3 een overzicht gegeven van de modellen die al in de literatuur ontwikkeld zijn. In hoofdstuk 4 worden de componenten van de cyclus 1 voor 1 gemodelleerd. Tot slot worden deze gecombineerd tot een steady-state en dynamisch cyclusmodel respectievelijk in hoofdstuk 5 en 6. In deze hoofdstukken worden de resultaten van de simulatie weergegeven en daarna vergeleken met de data opgemeten door Comstock en Braun [11]. Een conclusie en verdere aanbevelingen volgen in hoofdstuk 7.

Hoofdstuk 2. De subkritische compressie-koelcyclus

3

Hoofdstuk 2

De subkritische compressie-koelcyclus 2.1

Inleiding

Voor de volledigheid wordt in het eerste deel van dit hoofdstuk de algemene werking van een koelcyclus beschreven. Er bestaan zeer veel methodes om koeling te realiseren. Hier wordt enkel een subkritische compressie-koelmachine behandeld. Dit overzicht is gebaseerd op enkele referentiewerken [12, 24]. In het tweede deel wordt de koelmachine uit [11] beschreven. In deze sectie worden enkel de zaken besproken die van belang zijn voor dit huidig werk. Voor meer details over instrumentatie, dataverwerking en dergelijke wordt de lezer doorverwezen naar [11].

2.2 2.2.1

Algemene werking en begrippen Werking

In figuur 2.1 is de opbouw van een subkritische compressie-koelcyclus weergegeven. Deze cyclus bestaat uit 4 basiscomponenten: 1. de verdamper; 2. de compressor; 3. de condensor; 4. de klep. Het ge¨ıdealiseerd cyclusverloop is weergegeven door middel van een Ts- en ph-diagram in figuur 2.1. Het koelmiddel stroomt door de verdamper onder lage druk waar ze warmte onttrekt uit het secundair flu¨ıdum 1 en verdampt tot lichtjes oververhitte damp. Daarna wordt het koelmiddel door de compressor op condensordruk gebracht. Daar geeft ze warmte af aan het secundair flu¨ıdum waardoor ze condenseert tot onderkoelde vloeistof. Om de cyclus te vervolledigen, expandeert het koelmiddel in de klep opnieuw tot verdamperdruk. Het is belangrijk dat het koelmiddel bij intrede van de compressor gesatureerde damp is of oververhit is. Dit is noodzakelijk om vroegtijdig falen van de compressor te voorkomen. Immers, als het intredende flu¨ıdum een 2-fasig mengsel is, kunnen vloeistofdruppeltjes de compressor beschadigen. De inlaat van de compressor mag daarentegen niet te sterk oververhit zijn om de temperatuur in de compressor te beperken. Daarnaast moet ook de inlaat van de klep lichtjes onderkoeld zijn om instabiliteiten op systeemniveau te vermijden. 1

In het geval van water, heet het secundair flu¨ıdum in de verdamper ’het ijswater’, in de condensor wordt het ’het koelwater’ genoemd.


4

Figuur 2.1: Basis compressie-koelcyclus [12] Het werkelijk cyclusverloop ziet er anders uit. Ten gevolge van irreversibiliteiten gebeurt de compressie niet isentroop. Verder kan de uitlaat van de verdamper oververhit zijn en de uitlaat van de condensor onderkoeld. Ten slotte zijn er ladingsverliezen in de leidingen, maar ook in de verdamper en de condensor. Dit betekent dat de stroming doorheen de beide warmtewisselaars niet als isobaar kan beschouwd worden.

2.2.2

Begrippen

De ogenblikkelijke thermische prestatie van een koelmachine wordt typisch uitgedrukt door middel van de COP (Coefficient Of Perfomance). Deze wordt als volgt gedefinieerd: COP =

nuttig koeleffect netto energietoevoer

(2.2.1)

Een koelmachine kan in deellast werken. Dit betekent dat de huidige koellast afwijkt van de nominale. De mate van deellastwerking wordt aangegeven met de deellastverhouding of PLR (Part Load Ratio): nuttig koeleffect P LR = (2.2.2) nominaal nuttig koeleffect

2.3 2.3.1

De koelcyclus uit de proefstand Opbouw

In figuur 2.2 is de proefstand uit [11] schematisch weergegeven. De 4 basiscomponenten zijn duidelijk zichtbaar. Het systeem is echter uitgebreid met een koellijn en met regellussen. De koelcyclus is een McQuay koelmachine van het type PEH048J met een nominale koelcapaciteit van 300 kW . De verdamper en condensor zijn beiden 2-passige trommel-en-pijp warmtewisselaars. Het koelmiddel R134a loopt door de trommel (shell ). In tegenstelling tot de verdamper, loopt in de condensor het koelmiddel niet van onder naar boven, maar van boven naar onder. Door zijn grotere densiteit bevindt de vloeibare fase van het koelmiddel zich onderaan de trommel (liquid


5

water

condensor r134a

koellijn

motor expansie klep

compressor menger

verdamper

water

Figuur 2.2: Schematische weergave van de proefstand pool ). In de Engelstalige literatuur worden dergelijke warmtewisselaars aangeduid met de term flooded heat exchangers. Het secundair flu¨ıdum is water en loopt door de pijpen (tubes). Het stroomt binnen in de bovenste helft en verlaat de warmtewisselaar in de onderste helft. De constructieve details, zoals de lengte (L) en de binnen- en buitendiameter (respectievelijk φi en φo ) van de pijpen van beide warmtewisselaars zijn weergegeven in tabel 2.1. Dit zijn trouwens de enige gegevens die beschikbaar zijn over beide warmtewisselaars. Belangrijke gegevens die hier dus ontbreken zijn de eventuele aanwezigheid van tussenschotten, het al dan niet gevind zijn van de pijpen en hun schikking in de bundel.

verdamper condensor

φi cm 1.61 1.55

trommel φo L cm cm 1.96 243.84 1.91 243.84

aantal − 150 164

pijp φi cm 39.29 37.76

Tabel 2.1: Eigenschappen verdamper en condensor De compressor is een enkel-trappige centrifugaalcompressor. Het koelmiddel loopt eerst langs inlaatschoepen waardoor het expandeert. Een regellus stelt de hoek van de schoepen in om zo de ingestelde uitlaattemperatuur van het ijswater te bereiken (capacity control ). Op geregelde tijdstippen meet een controller de uitlaattemperatuur van het ijswater. Als ze de insteltemperatuur overschrijdt, stuurt de controller een signaal naar een hydraulisch systeem dat de inlaatschoepen opent. Bijgevolg neemt het koelmiddeldebiet toe en daalt de uitgangstemperatuur van het water. Er kan immers meer warmte overgedragen worden aan het koelmiddel zonder dat de uitlaat van de verdamper sterk oververhit wordt. Als de uitlaattemperatuur van het ijswater echter te laag is, moet de koelcapaciteit van de verdamper weer afnemen. Dit gebeurt door de schoepen opnieuw wat te sluiten. Van de compressor is enkel de diameter van de rotor en het toerental gekend, zie tabel 2.2. Van het stuuralgoritme is niets gekend, aangezien dit confidentiële informatie is. De expansieklep is een thermostatische expansieklep zoals in figuur 2.3. Deze bestaat uit een hoofdklep en een pilootklep. Boven het membraan van de pilootklep werkt de druk van een

Hoofdstuk 2. De subkritische compressie-koelcyclus diameter rotor toerental motor overbrengingsverhouding toerental rotor tipsnelheid rotor

6 0.122 3521 8.909 31368 200

m rpm − rpm m/s

Tabel 2.2: Eigenschappen compressor reservoir (bulb) dat gevuld is met het koelmiddel R500 dat zich altijd in gesatureerde toestand bevindt. Dit reservoir staat in thermisch contact met de uitgang van de verdamper. Dit betekent dat in steady-state de druk in dit reservoir gelijk is aan de saturatiedruk horend bij de uitgangstemperatuur van de verdamper. De kracht die door deze druk wordt uitgeoefend op het membraan wordt in de pilootklep in evenwicht gehouden door de veerkracht van de sluitveer en door de verdamperdruk. Strikt genomen werken er nog 2 extra krachten op de afsluiter: de traagheidskrachten werkend op de afsluiter en de kracht resulterend uit de drukval over de afsluiter. Deze laatste is gelijk aan de drukval vermenigvuldigd met zijn oppervlakte. In de meeste gevallen is de oppervlakte van het membraan groot genoeg ten opzichte van de werkzame oppervlakte van de afsluiter, zodat deze kracht kan verwaarloosd worden. De reservoirdruk stijgt wanneer de oververhitting aan de uitlaat van de verdamper toeneemt. Hierdoor zal de afsluiter in de pilootklep zakken zodat het massadebiet toeneemt. Een hoger massadebiet betekent dat de tussendruk toeneemt, zodat ook de afsluiter in de hoofdklep begint te zakken. Hierdoor vergroot het massadebiet doorheen de volledige klep en neemt de oververhitting af. De gewenste oververhitting kan ingesteld worden door middel van het voorspannen van de sluitveer. Buiten het koelmiddel in het reservoir is er van de klep geen enkel detail beschikbaar.

reservoir

membraan tussendruk

membraan

afsluiter

pilootklep

INGANG VERDAMPER

hoofdklep

UITGANG CONDENSOR

Figuur 2.3: De klep uit de proefstand (naar [11]) Slechts een deel van het koelmiddel dat uit de condensor komt, stroomt doorheen de expansieklep. Het andere deel wordt gebruikt om de motor en de transmissie te koelen. Dit koelmiddel loopt eerst door een nauwe opening waardoor het expandeert. Daarna neemt het een deel van de verlieswarmte op die in de motor en de transmissie wordt vrijgegeven. Aan de inlaat van de verdamper worden beide stromen weer gemengd. Ook van deze koellijn (cooling line) zijn er geen data gekend.


2.3.2

7

Instrumentatie

Het systeem is uitgerust met instrumentatie om debiet en temperatuur te meten, daarnaast wordt het ogenblikkelijk vermogen van de compressor geregistreerd. Deze sensoren sturen hun meetwaarden door naar een desktop PC waarop het programma VisSim draait. Dit programma is in staat is real-time datacollectie en -analyse te doen. Uit de meetdata berekent het waarden zoals de koellast, de COP en de compressor-efficiëntie. Een overzicht van de gemeten en berekende waarden en de bijhorende onzekerheden is gegeven in [11] en is terug te vinden in bijlage A.

2.3.3

Testcondities

Door de randvoorwaarden van verdamper en condensor te wijzigen, werd het systeem onderworpen aan 27 verschillende testcondities. Elke testconditie wordt aangeduid met een nummer zoals in tabel 2.3. De druk en het massadebiet van koelen ijswater worden constant gehouden. Het verloop van de inlaattemperatuur van het secundair flu¨ıdum in zowel condensor (Tcwi ) als verdamper (Tewi ) en het instelpunt van de ijswateruitlaattemperatuur (Tset ) zijn echter trapfuncties. Figuur 2.4 geeft hiervan het tijdsverloop. In deze figuur is ook het instelpunt van de koellast (P LRset ) uitgezet. Deze is gelijk aan P LRset =

m ˙ 2 Cp,2 (Tewi − Tset ) Qnom

(2.3.1)

temperatuur [K]

In deze vergelijking is m ˙ 2 het massadebiet dat door de verdamper loopt. Cp,2 is de specifieke warmtecapaciteit van het ijswater. Qnom is de nominale koelcapaciteit van de machine. Deze is gelijk aan 300 kW .

Tset

310

Tcwi

Tewi

300 290 280 0

1

2

3

4

tijd [s]

5 4

x 10

PLRset [%]

200 150 100 50 0 0

1

2

3 tijd [s]

4

5 4

x 10

Figuur 2.4: Testsequentie

Om de 10 seconden werd data geregistreerd. Om de 30 minuten werden de randvoorwaarden aangepast. Bij de gegeven testsequentie nam het overgangsverschijnsel van het systeem 10

Hoofdstuk 2. De subkritische compressie-koelcyclus Test 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27

Tset [K] 282.6 283.2 283.2 282.6 283.2 283.2 282.6 283.2 283.2 279.8 280.4 280.4 279.8 279.8 280.4 279.8 279.8 280.4 277.0 277.6 277.6 277.0 277.0 277.6 277.0 277.6 277.6

Tewi [K] 289.1 287.1 285.2 289.4 286.8 284.7 287.9 286.5 285.1 286.6 283.9 282.7 286.5 283.3 282.2 286.2 283.9 282.3 283.2 281.2 279.9 284.2 280.7 279.6 282.4 280.9 279.3

8 Tcwi [K] 303.0 302.8 302.7 297.3 297.0 297.1 294.3 291.7 290.1 302.8 302.8 302.6 297.2 297.2 297.1 294.4 291.5 291.4 299.9 299.9 299.9 294.4 293.9 293.3 290.6 289.4 289.2

Tabel 2.3: Testsequentie a` 15 minuten in beslag. Op deze manier worden dus zowel transiënt als steady-state data opgenomen.

Hoofdstuk 3. Literatuurstudie

9

Hoofdstuk 3

Literatuurstudie 3.1

Inleiding

Koelmachines kunnen theoretisch, maar ook op empirische wijze gemodelleerd worden. Empirische modellen beschrijven het gedrag van een koelcyclus aan de hand van correlaties die afgeleid zijn uit meetdata. Er kan een onderscheid gemaakt worden tussen black-box en graybox modellen. Het verschil tussen beiden is dat in een black-box model, in tegenstelling tot een gray-box model, de geschatte parameters geen fysische betekenis hebben. Een mooi overzicht van enkele empirische methodes is gegeven in [17]. In het huidig werk is er geopteerd voor een meer complex, maar ook veel krachtiger theoretisch model. Het voordeel van theoretische modellen ten opzichte van empirische modellen, is dat de bekomen modelvoorstelling niet zo sterk gebonden is aan het type koelmachine. Indien er bijvoorbeeld een theoretisch model ontwikkeld wordt voor een watergekoelde cyclus, dan zijn de bekomen vergelijkingen volledig analoog aan degene die men zou bekomen bij het modelleren van een luchtgekoelde machine. In secties 3.2, 3.3 en 3.4 wordt per component een overzicht gegeven van de modellen die in de literatuur terug te vinden zijn. Hierin worden voornamelijk de theoretische modellen behandeld. Er is ook aandacht besteed aan semi-theoretische modellen gezien de beperkte hoeveelheid beschikbare constructieve data. Uiteraard zijn er ook talloze empirische modellen beschikbaar, hiervoor wordt de lezer doorverwezen naar de literatuur. Dit hoofdstuk sluit af met sectie 3.5, waarin de systeemmodellering wordt besproken.

3.2

Thermostatische expansieklep

In de literatuur is niets terug te vinden over de modellering van een piloot-gestuurde expansieklep, vandaar dat dit overzicht zich beperkt tot standaard thermostatische expansiekleppen zoals in figuur 3.1

3.2.1

Inleiding

Een thermostatische expansieklep wordt gemodelleerd als isenthalp. Daarnaast wordt er verondersteld dat er geen massa-ophoping is in de klep. Er moet dus enkel nog een uitdrukking gevonden worden voor het massadebiet m ˙ r . In de literatuur komen 2 methodes voor. De eerste is een klassieke aanpak. Er wordt vertrokken van de Bernoulli-klepvergelijking die daarna gecorrigeerd wordt op basis van experimenten en theoretische resultaten. Daarnaast bestaat er ook een meer complexe methode gebaseerd op het stromingspatroon doorheen de klep. In de volgende 2 secties worden beide methodes meer in detail besproken.


10

membraan reservoir drukgelijkmaker uitlaat condensor

inlaat verdamper

klepafsluiter

sluitveer

Figuur 3.1: Standaard thermostatische expansieklep

3.2.2

Klassieke methode

De klassieke modellen steunen op de Bernoulli-klepvergelijking: p m ˙ r = Cv A 2ρin (pin − pout )

(3.2.1)

Het probleem is een geschikte uitdrukking te vinden voor de doorstroomsectie A enerzijds en de klepcoëfficiënt Cv anderzijds. Klepco¨ effici¨ ent Cv Bendapudi et al. [4] en Zhao en Zaheeruddin [48] gaan ervan uit dat Cv constant is. Strikt genomen is dit enkel geldig als de stroming 1-fasig en onsamendrukbaar is. In een koelcyclus is dit zeker niet het geval. Het koelmiddel expandeert immers en gaat meestal over van vloeistof naar 2-fasig (flashing). Verder treedt er choking op wanneer de uitgangsdruk voldoende laag is. Het massadebiet kan met andere woorden niet onbeperkt toenemen bij steeds afnemende uitgangsdruk, zoals zou blijken indien de klepcoëfficiënt constant is. Stearns et al. [34] proberen hieraan tegemoet te komen door Cv uit te drukken in functie van het Reynoldsgetal. Er zijn enkele functievoorschriften vooropgesteld waarbij de onbekende parameters bepaald worden door experimenten uit te voeren op de klep. Deze methode levert aanvaardbare resultaten op, maar heeft geen echte fysische achtergrond. Soms stelt de fabrikant gegevens ter beschikking onder de vorm van performantiekaarten [47]. Hierbij wordt Cv uitgedrukt in functie van bijvoorbeeld de drukval over de klep en de voorspanning van de veer. Deze curves zijn experimenteel bepaald. Doorstroomsectie A De doorstroomsectie is een functie van de klepheffing (y) en de geometrie van de klepafsluiter. In het geval van een conische afsluiter bijvoorbeeld [4, 40], geldt: A = a1 y + a2 y 2

(3.2.2)


11

De coëfficiënten a1 en a2 zijn zuiver geometrisch te bepalen. De klepheffing op zijn beurt wordt bepaald op basis van het krachtenevenwicht van de klepafsluiter [40, 23, 11]: Fbulb = Fref + Fspring + Fdyn + Fevap

(3.2.3)

Gezien de geringe inertie van de afsluiter, worden de dynamische krachten (Fdyn ) verwaarloosd. Deze vereenvoudiging heeft nauwelijks invloed op systeemniveau aangezien de dynamica van de cyclus gedomineerd wordt door de verdamper en condensor [4]. Er heerst dus een balans tussen de neerwaartse kracht ten gevolge van de druk in het reservoir (Fbulb ) en de opwaartse kracht ten gevolge van de verdamperdruk (Fevap ) en de sluitveer (Fspring ). Door de drukval over de afsluiter moet er nog een extra kracht in rekening gebracht worden (Fref ). De werkzame oppervlakte van de afsluiter is echter zeer klein ten opzichte van de oppervlakte van het membraan, zodat deze term kan verwaarloosd worden. De reservoirdruk is gelijk aan de saturatiedruk overeenkomend met de temperatuur van het reservoir (Tbulb ). Bij perfecte isolatie van het reservoir, is deze temperatuur in steady-state gelijk aan de uitgangstemperatuur van de verdamper (Tout ). In dynamische simulaties moet er ook rekening gehouden worden met de thermische capaciteit van het koelmiddel in het reservoir, de wanden van het reservoir en de pijpwand. Tahat et al. [35] houden met al deze factoren rekening. In [40] worden de dynamica en de thermische weerstand verwaarloosd. Bijgevolg zijn beide temperaturen gelijk aan elkaar. Dit is een vrij grove vereenvoudiging. De tijdsconstante van het reservoir heeft immers dezelfde grootte-orde als de dynamica van het systeem [4]. De meeste auteurs [4, 11, 23] bepalen de reservoirtemperatuur dan ook door middel van een vereenvoudigde energievergelijking: K

dTbulb = Tout − Tbulb dt

(3.2.4)

K is hierbij een constante en wordt meestal experimenteel [4] of (benaderend) theoretisch bepaald [23].

3.2.3

Stromingsmechanische methode

Figuur 3.2: Analogie stroming doorheen een klep en stroming doorheen een convergentdivergent gevolgd door dumpdiffusie [20] Een meer complexe aanpak is gebaseerd op de bevindingen van Simoes-Moreira en Bullard [32]. Er wordt een analogie gemaakt met de stroming doorheen een convergent-divergent gevolgd


12

door dumpdiffusie. De grootte van de keelsectie (throat) is afhankelijk van de positie van de klepafsluiter. Er wordt verondersteld dat het koelmiddel onderkoeld is aan de inlaat van de klep. Na een zeer snelle drukdaling in het convergent wordt het koelmiddel oververhit. Zoals aangetoond door Abuaf et al. [1] en Schrock et al. [29] begint onder deze omstandigheden de fase-overgang (flashing inception) in de keelsectie. Er ontstaat een 2-fasig mengsel dat initieel echter afwijkt van thermodynamisch evenwicht. Dit metastabiel mengsel wordt verder stroomafwaarts stabiel door het optreden van een verdampingsgolf (evaporation wave). Als de uitgangsdruk voldoende laag is, versnelt het flu¨ıdum tot supersoon. Daarna treedt een schokgolf op, waardoor de druk sterk toeneemt. Ten slotte wordt nog wat druk herwonnen door dumpdiffusie. Het stromingspatroon en het drukverloop zijn weergegeven in figuur 3.3.

Figuur 3.3: Het stromingspatroon en het drukverloop in een klep [20] De kritische druk voor choking wordt afgeleid in [20]. Wanneer choking optreedt in de klep, kan het massadebiet bepaald worden met de formules van Abuaf et al. [1]. Wanneer de uitlaatdruk voldoende hoog is, treedt er geen choking op. Er is wel nog een metastabiele zone die stroomafwaarts wordt beëindigd door een verdampingsgolf. In de literatuur is echter niet volledig beschreven hoe de stroming dan kan worden uitgerekend. De toestand voor en na een verdampingsgolf kan berekend worden met de formules afgeleid door Simoes-Moreira [31]. Hoe een metastabiele stroming moet uitgerekend worden, is echter niet bekend. Hoewel Liu et al. [20] een metastabiele stroming moeten uitrekenen om het chokingcriterium af te leiden, zijn ze hier zeer onduidelijk over. Deze kennis is echter noodzakelijk om een uitdrukking te vinden voor het massadebiet. Chen et al. [10] stellen een alternatieve methode voor om het massadebiet doorheen een klep te bepalen. Opnieuw wordt de Bernoulli-klepvergelijking uitgeschreven, niet tussen inen uitlaat, maar tussen inlaat en de keelsectie. Aangezien de stroming doorheen het convergent in de vloeistoffase blijft, is de klepcoëfficiënt constant. Typisch is deze gelijk aan 0.94 ± 0.04 [33]. De grootte van de keelsectie (Ath ) wordt bepaald uit de geometrie van de klep en de klepheffing. Het komt er nu nog op neer de druk in de keelsectie (pth ) te bepalen. In [10] wordt deze grootheid gefit in volgende correlatie: 2 pth = a0 + a1 pin + a2 ∆Tsc + a3 ∆Tsc + a4 Ath

Hierin is ∆Tsc de mate van onderkoeling van het flu¨ıdum aan de inlaat.

(3.2.5)


13

De stromingsmechanische modellen zijn fysisch sterk onderbouwd en leveren bijgevolg ook meer betrouwbare en nauwkeurige resultaten op dan het klassiek model. Dit vereist echter wel voldoende kennis van de constructieve details van de klep. Om confidentiële redenen zijn deze in de meeste gevallen echter niet beschikbaar.

3.3 3.3.1

Centrifugaalcompressor Inleiding

Sommige onderzoekers gaan er van uit dat alle dynamica verwaarloosbaar is. Zij houden dus geen rekening met: • de dynamica van het flu¨ıdum (massa, impuls en energie); • de dynamica van de compressor (inertie rotor en energie wanden); • de dynamica van het regelsysteem. Zo wordt een volledig steady-state model bekomen [5]. Deze zijn echter weinig bruikbaar om opstarten uitzetverschijnselen te beschrijven of in het geval van een frequentiegestuurde compressor. Om hieraan tegemoet te komen zijn er modellen ontwikkeld die wel dynamica in rekening brengen. Gravdahl en Egeland [13] en Jiang et al. [16] veronderstellen dat enkel de rotor (impeller ) dynamica vertoont. Morini et al. [25] daarentegen hebben een model opgesteld dat de volledige dynamica van de compressor, zowel van de rotor als de stroming, beschrijft.

(a) Technische tekening

(b) Foto

Figuur 3.4: Centrifugaalcompressor met inlaatschoepen [2] In wat volgt, wordt een overzicht gegeven van zowel de steady-state als de dynamische modellen. De nadruk wordt gelegd op de theoretische modellen. Experimentele modellen zijn echter vrij eenvoudig op te stellen indien voldoende experimentele data beschikbaar zijn [36, 38]. Deze methodiek wordt vaak aangewend wanneer onvoldoende constructiedetails beschikbaar zijn, zoals in [11, 4]. Het bekomen model heeft echter een beperkte toepasbaarheid. Het is enkel geschikt voor 1 bepaald type compressor en werkingsflu¨ıdum. Daarnaast zijn er meestal niet genoeg meetgegevens beschikbaar om ook het off-design gedrag van de compressor nauwkeurig te voorspellen.


3.3.2

14

Steady-state modellen

Bourdouxhe et al. [5] hebben een steady-state model ontwikkeld van een centrifugaalcompressor. In [44, 8] wordt hiervan gebruik gemaakt. Het model bestaat uit 2 vergelijkingen: Ay κ px py κ−1 κ 2 m ˙ = ρy tan(β) −1 −U (3.3.1) U κ − 1 ρx px py px = κ (3.3.2) κ ρx ρy Deze vergelijkingen volgen uit de Euler-vergelijking voor turbomachines en de energiebalans over compressor en rotor. Er wordt tevens verondersteld dat de compressie isentroop verloopt. Er wordt met andere woorden van uitgegaan dat alle verliezen geconcentreerd zijn in de motor en transmissie. Verder wordt er aangenomen dat de ingaande kinetische energie geen significant effect heeft op de energiebalans. Met de index x en y worden respectievelijk de inen uitlaat van de impeller aangeduid. De meeste auteurs [8, 44, 5] verwaarlozen het effect van de nageschakelde diffusor en collector, zodat de toestand y overeenkomt met de uitlaat van de compressor (out). Indien er inlaatschoepen aanwezig zijn, komt de toestand aan de inlaat van de rotor niet overeen met die aan de inlaat van de compressor (in). Deze inlaatschoepen zijn aangebracht om de stoot aan de inlaat van de rotor te reduceren. In [5] wordt de stroming doorheen deze schoepen gemodelleerd als een equivalent smoringproces. De smoring neemt toe wanneer het instelpunt van de PLR verkleint. Bijgevolg daalt de densiteit aan de inlaat van de rotor zodat het massadebiet koelmiddel dat door de cyclus loopt, afneemt en dus de koelcapaciteit daalt. Hieruit volgt dat: hx = hin

(3.3.3)

px = λpin

(3.3.4)

Figuur 3.5: Compressie-koelcyclus van een koelmachine met inlaatschoepen [44] De mate van smoring wordt uitgedrukt met de factor λ. In [5] wordt geen uitdrukking gegeven min voor λ. Browne en Bansal [8] veronderstellen echter dat λ = y + z Tset − Tset . z en y zijn 2 constanten die experimenteel bepaald worden. Deze vergelijking drukt uit dat de hoek van inlaatschoepen zich instelt afhankelijk van het verschil tussen de huidige en de minimaal mogelijke instelwaarde van de uitlaattemperatuur van het ijswater. Als beide temperaturen gelijk zijn, werkt de machine in vollast en is λ = y. De schoepen zijn volledig open, maar er kan toch nog een zeker mate van smoring optreden. De smoring neemt toe wanneer Tset toeneemt. De machine werkt dan in deellast. Yu en Chan [44] bepalen het massadebiet m ˙ als αm ˙ max . Het maximale massadebiet m ˙ max


15

wordt berekend op basis van formule (3.3.1), waarbij de inlaatschoepen volledig open staan. De factor α wordt bepaald aan de hand van een polynomisch verband met de deellastverhouding. Het vermogen geleverd door de compressor aan het gecomprimeerd flu¨ıdum wordt in [8, 44] bepaald als: ˙ ˙ comp = Win W (3.3.5) ηpol Hierbij is: ˙ in W

κ px = m ˙ κ − 1 ρx

"

pout px

κ−1

#

κ

−1

(3.3.6)

Dit is het vermogen dat wordt overgedragen aan het flu¨ıdum indien de stroming isentroop is. Door verliezen ten gevolge van onder andere wrijving is de stroming niet isentroop, dus zijn beide vermogens verschillend. Het verlies wordt begroot aan de hand van een polytroop rendement ηpol . Browne en Bansal [8] bepalen deze aan de hand van de methode van Braun et al. [6]. Volgens Yu en Chan [44] daarentegen is dit rendement een kwadratische functie van de deellastverhouding. Door rekening te houden met de elektromechanische verliezen in de motor en tranmissie, kan het totale elektrische vermogen berekend worden als: ˙ tot = W

˙ in W ηpol ηem

(3.3.7)

In [8, 44] wordt het elektromechanisch rendement begroot door een constante ηem . Om tot vergelijking (3.3.1) te komen, werd door Bourdouxhe et al. [5] verondersteld dat de compressie isentroop verloopt. Ten gevolge van irreversibiliteiten, verloopt de compressie in werkelijkheid echter nooit isentroop. Toch zijn de resultaten van bovenstaande studies bruikbaar, aangezien het compressieproces met goede benadering als polytroop kan benaderd worden. Dit levert immers dezelfde vergelijkingen op, alleen moet de isentrope exponent κ vervangen worden door een polytrope exponent γ, die kleiner is en zo de verliezen in de stroming reflecteert.

3.3.3

Dynamische modellen

Jiang et al. [16] vertrekken van de karakteristieke vergelijkingen van een compressor: 1 [τmotor − τcomp ] J ˙ comp = m∆h ωτcomp = W ˙ γ ηi ∆h γ−1 pout = pin 1 + Tin Cp ω˙ =

(3.3.8) (3.3.9) (3.3.10)

Vergelijking (3.3.8) is de bewegingsvergelijking voor starre lichamen en drukt de dynamica van de rotor uit. Het verschil tussen het aandrijvende (τmotor ) en tegenwerkende koppel (τcomp ) doet de rotor met traagheidsmoment J accelereren (ω˙ > 0) of decelereren (ω˙ < 0). Door de dynamica van de stroming te verwaarlozen, reduceert het energiebehoud zich tot vergelijking (3.3.9). ∆h is de enthalpietoename resulterend uit arbeidsoverdracht van de rotor naar het koelmiddel. Ten slotte geeft vergelijking (3.3.10) het verband tussen druk aan inen uitlaat, met ηi het isentrope rendement van de compressor. In een tweede stap begroten Jiang et al. [16] ηi . Dit gebeurt door de verliezen door wrijving, backflow, verdringing en incidentie in de rotor te bepalen en de verliezen in het slakkenhuis


16

te modelleren. Door begroting van ηi en eliminatie van τcomp en ∆h wordt een stelsel bekomen. Deze vergelijkingen laten toe het massadebiet, de uitgangsdruk, de rotorsnelheid en de efficiëntie te berekenen in de tijd onder uiteenlopende werkingsomstandigheden. Gravdahl en Egeland [13] gaan op een gelijkaardige manier tewerk. Hier wordt een compressormodel opgesteld om daarna een surge- en snelheidscontrolewet op te stellen. Ook Morini et al. [25] vertrekken van de behoudswetten: • massa-, impuls- en energiebehoud van het koelmiddel; • bewegingsvergelijking van de rotor; • energiebehoud van het materiaal. Deze dynamische vergelijkingen worden uitgeschreven in verschillende controlevolumes. Daarna wordt het stelsel numeriek opgelost aan de hand van het 5de orde Runge-Kutta algoritme.

3.4 3.4.1

Verdamper en condensor Inleiding

In de literatuur zijn er zowel statische als dynamische modellen van warmtewisselaars terug te vinden. Stationair vollast-, maar ook deellastgedrag van een koelmachine kan perfect voorspeld worden met een statisch verdamperen condensormodel. Opstart- of uitzetverschijnselen zijn echter dynamische effecten. Om dit te simuleren, zijn statische modellen uiteraard niet voldoende. Ook wanneer er nagegaan moet worden hoe een systeem reageert bij vrij snel veranderende randvoorwaarden, is op zijn minst een dynamisch warmtewisselaarmodel vereist. Aangezien een dynamisch model eenvoudig kan gereduceerd worden tot een statisch model, worden in wat volgt enkel de dynamische behandeld.

3.4.2

Dynamische modellen

In de literatuur zijn er 2 gangbare modelleringsmethodes: finite-volume en moving-boundary. In beide gevallen start de modellering met het opdelen van de warmtewisselaar in controlevolumes. In deze controlevolumes worden de behoudswetten uitgeschreven en vereenvoudigd. Onderstaande lijst geeft kort weer welke vereenvoudigingen het meest in de literatuur voorkomen en hoe deze verrechtvaardigd worden [40, 48, 19, 47, 22, 18, 4, 45, 42, 26]. 1. De stroming van koelmiddel en secundair flu¨ıdum is 1-dimensionaal. In werkelijkheid is de stroming 3-dimensionaal, maar in de meeste warmtewisselaars is er een dominante stromingsrichting tussen inen uitlaat. Een meer-dimensionale behandeling van de stroming zou trouwens de rekentijd zeer sterk doen toenemen en niet zeer veel bijdragen op systeemniveau. 2. De drukval in het koelmiddel en het secundair flu¨ıdum is verwaarloosbaar. Op deze manier valt de impulsvergelijking weg uit het stelsel. In realiteit is de drukval bij goed onderhouden warmtewisselaars ook vrij beperkt en is dit een aanvaardbare veronderstelling. Het secundair flu¨ıdum is meestal water. Het is gekend dat de eigenschappen van water weinig afhankelijk zijn van de druk, ook dit rechtvaardigt deze aanname. Er zijn echter koelmiddelen waarbij de temperatuur merkbaar verandert ook al is de drukval klein. Dit effect kan een vrij grote invloed hebben op de prestaties van de warmtewisselaar. Daarom zijn er modellen ontwikkeld die toch rekening houden met de drukval.


17

In sommige werken wordt de impulsvergelijking vervangen door een algebra¨ısche vergelijking. Er wordt verondersteld dat de drukval voornamelijk te wijten is aan wrijving en dat het effect van inertie en impulsverandering verwaarloosbaar is. Een dergelijke aanpak is voorgesteld door [43]. Deze vereenvoudiging ontkoppelt de berekening van de drukval van het oplossen van de energie- en massavergelijkingen. Zhang et al. [46] daarentegen behouden de impulsvergelijking in zijn originele vorm als differentiaalvergelijking. Deze onderzoekers hebben de resultaten van een transiënte, een steady-state en een nul-drukval impulsvergelijking in een moving-boundary modellering met elkaar vergeleken. Ze kwamen tot de conclusie dat de transiënt en steady-state modellering de meest nauwkeurige resultaten opleverden. De nul-drukval methode is echter veel sneller. Ook in [15, 28] blijft de impulsvergelijking behouden. Er wordt echter gemodelleerd volgens de finite-volume methode. 3. De axiale conductie in het materiaal, koelmiddel en secundair flu¨ıdum is verwaarloosbaar. Zowel de temperatuur van het koelmiddel als van het secundair flu¨ıdum wijzigen tijdens hun passage door de warmtewisselaar. Hierdoor ontstaan er temperatuursgradiënten die axiale warmtevloed induceren. Bendapudi et al. [4] merken op dat, door het relatief hoge massadebiet, conductieve effecten in zowel het koelmiddel als het secundair flu¨ıdum verwaarloosd kunnen worden. Het relatieve belang van convectieve ten opzichte van diffusieve warmte-overdracht in een flu¨ıdum wordt uitgedrukt door het Péclet-getal: P eD =

Dv α

(3.4.1)

In deze vergelijking is D een karakteristieke lengte, v de snelheid en α de thermische diffusiviteit van het flu¨ıdum in kwestie. Inderdaad, wanneer het massadebiet groot is, domineert de convectieve warmte-overdracht. Hoewel de axiale conductie in het materiaal veel groter is dan in de flu¨ıda, wordt ook dit verschijnsel verwaarloosd. Dit heeft fysische, maar ook wiskundige redenen. Fysisch is de dwarssectie voor axiale conductie immers klein ten opzichte van het manteloppervlak waarlangs de convectieve warmte-overdracht plaatsvindt. In een moving-boundary model van een mee- of tegenstroom warmtewisselaar, vormt deze term wiskundig een probleem. Ze verhindert dat er een stelsel 1ste orde differentiaalvergelijkingen wordt bekomen na integratie over de zones. Dit is omdat de conductieterm een afgeleide van 2de orde is. Het relatief belang van deze term wordt echter steeds groter naarmate de lengte van een zone naar nul gaat [26]. Vandaar dat Pettit et al. [26] en Willatzen et al. [42] een algebra¨ısche correctieterm invoeren in de zone-specifieke vergelijkingen om met het conductie-effect rekening te houden. Wanneer de warmtewisselaar in dwarsstroom werkt, komt dit probleem niet voor [47]. 4. De thermische weerstand van het materiaal is verwaarloosbaar. De weerstand tegen radiale warmtetransfer is klein ten gevolge van het grote warmtewisselende oppervlak en de hoge conductiviteit. In vele werken wordt ook de geometrie van de warmtewisselaar vereenvoudigd. Bendapudi et al. [4] bijvoorbeeld modelleren een 2-passige trommel-en-pijp warmtewisselaar als een 1-passige. Verder worden aannames gemaakt wat betreft het verloop van de toestandsgrootheden. Specifiek voor moving-boundary modellen, wordt er van uitgegaan dat de temperatuur van het secundair flu¨ıdum lineair verloopt [4, 48]. Ook het enthalpieprofiel wordt lineair verondersteld in iedere zone [45, 47, 4, 18, 22]. In het geval van luchtgekoelde machines, wordt de dynamica van het secundair flu¨ıdum verwaarloosd [22, 18]. In watergekoelde cycli, wordt hier wel rekening mee gehouden. In dit geval worden er echter andere veronderstellingen gemaakt. Er wordt meestal van uitgegaan dat de warmtewisselaar perfect ge¨ısoleerd is [4, 48]. Dit betekent dat de warmtevloed naar de omgeving niet in rekening gebracht wordt. In werkelijkheid is de isolatie echter nooit perfect, maar altijd voldoende om het warmteverlies te kunnen verwaarlozen ten


18

opzichte van de warmte-overdracht tussen het koelmiddel en secundair flu¨ıdum. Daarnaast wordt ook meestal de capaciteit van de trommelwand niet in rekening gebracht [4, 48]. Hoewel deze in werkelijkheid vrij groot is, blijkt het dat er slechts een kleine hoeveelheid van de totale warmte naar de trommel overgedragen wordt [4]. Dit is te wijten aan de relatief kleine warmtewisselende oppervlakte tussen het flu¨ıdum in de trommel en de trommelwand ten opzichte van het oppervlak van de pijpen. De moving-boundary en finite-volume modellen onderscheiden zich op basis van de onderverdeling in controlevolumes en het daaruit resulterende stelsel vergelijkingen. Finite-volume methode In een finite-volume model wordt de warmtewisselaar opgedeeld in een vast aantal volumes met vaste grootte (zie figuur 3.6). Vervolgens worden in elk volume de behoudswetten uitgeschreven, vereenvoudigd en gediscretiseerd. Dit resulteert in een algebra¨ısch stelsel dat met de gepaste numerieke technieken kan opgelost worden. Bendapudi et al. [4] hebben aangetoond dat een 15-tal controlevolumes voldoende zijn om het gedrag van een warmtewisselaar nauwkeurig te beschrijven. Een goed uitgewerkt finite-volume model is terug te vinden in [40].

Figuur 3.6: Finite-volume methode [23]

Moving-boundary methode Volgens de moving-boundary methode wordt de warmtewisselaar onderverdeeld in maximum 3 zones. De rand van deze zones valt samen met de fasegrenzen van het koelmiddel die optreden in de warmtewisselaar. De inlaat van een condensor is typisch oververhitte damp, de uitlaat is typisch onderkoelde vloeistof. Zo ontstaan er 3 zones in de condensor. In een verdamper daarentegen zijn er meestal slechts 2 zones, aangezien het koelmiddel aan de inlaat 2-fasig is en aan de uitlaat normaal gezien licht oververhitte damp. In ieder controlevolume worden de behoudswetten uitgeschreven en vereenvoudigd. Op deze manier wordt een stelsel partiële differentiaalvergelijkingen bekomen. In een volgende stap worden de vergelijkingen ge¨ıntegreerd over iedere zone. Dit resulteert in een stelsel 1ste orde differentiaalvergelijkingen: ˙ A(x)x˙ = f(x; u; u) (3.4.2) Hierin is x de toestandsvector terwijl u de vector met randvoorwaarden voorstelt. De toestandsvector bevat de toestand aan inlaat en uitlaat, maar ook interne variabelen, zoals zonelengtes en wandtemperaturen. Door integratie komen er echter nieuwe onbekenden in het stelsel: de interface-temperaturen van de wand. In de literatuur wordt soms gebruik gemaakt van een discontinue stuksgewijze functie [22, 18]. De wandtemperatuur tussen de vloeistof- en


19

(a) Verdamper

(b) Condensor

Figuur 3.7: Optredende modes in verdamper en condensor volgens [18] 2-fasenzone (Tt,L−T P ) wordt dan bijvoorbeeld bepaald als: ( L Tt,T P als dd dt > 0 Tt,L−T P = L Tt,L als dd dt < 0

(3.4.3)

Hierin is dL de lengte van de vloeistofzone. Andere auteurs veronderstellen dat de interfacewaarden gelijk zijn aan de zone-waarden van 1 van de aanliggende zones [42, 4, 19]. In [45] wordt aangetoond dat de tweede aanpak fysisch onrealistische waarden kan opleveren voor de zone-temperaturen van de pijp. Wanneer een zone verdwijnt, kan de zone-lengte sneller naar nul gaan dan zijn tijdsafgeleide. Bijgevolg zal de tijdsafgeleide van de zone-temperatuur zeer grote waarden aannemen. In hetzelfde artikel wordt ook de eerste methode bekritiseerd. Het voorkomt abnormale waarden van de zone-temperatuur, maar het is een discontinue functie. Bijgevolg is er geen vlotte overgang tussen dddtL > 0 en dddtL < 0. Zhang en Zhang [45] benaderen de interface-waarden door een continue algebra¨ısche functie: Tt,L−T P = Tt,L

dT P dL + Tt,T P dL + dT P dL + dT P

(3.4.4)

Het bekomen stelsel (vergelijking (3.4.2)) is afhankelijk van welke zones zich voordoen in de warmtewisselaar. Het kan zijn dat door verstoring in de randvoorwaarden 1 van de huidige zones verdwijnt of 1 nieuwe zone tevoorschijn komt. De warmtewisselaar kan zich met andere woorden in verschillende modes bevinden. Als gevolg hiervan verdwijnt 1 set vergelijkingen of komt er 1 bij. Het probleem dat zich dus stelt is 3-voudig: 1. Hoe wordt een overgang tussen 2 modes gedetecteerd? 2. Hoe wordt bepaald welke overgang plaatsvindt? 3. Hoe wordt de overgang numeriek aangepakt? De overgang tussen 2 modes wordt gedetecteerd aan de hand van een schakelcriterium (switching criterion). Dit criterium geeft ook aan naar welke mode moet worden geschakeld. In figuur 3.8 is een schakelcriterium weergegeven op basis van enthalpie. Dit is fysisch het meest


20

Figuur 3.8: Schakelcriterium op basis van enthalpie [26] voor de hand liggende criterium en wordt vaak gebruikt [19, 26, 42, 47, 45]. Daarnaast is enthalpie een continue grootheid, zodat er bij het schakelen tussen modes geen discontinu¨ıteiten optreden. Een andere mogelijkheid is schakelen op void-fractie of op zone-lengte [22, 18]. Wanneer geschakeld wordt naar een andere mode, moeten de interface-waarden aangepast worden aan de randvoorwaarden. Verder moet er ook voor gezorgd worden dat wanneer een zone weer actief wordt, de parameters en stofeigenschappen in die zone bepaald worden op basis van fysisch realistische zone-waarden. Deze zone-waarden moeten met andere woorden goed ge¨ınitialiseerd worden. Meestal blijft de algemene vorm (de toestandsvector en het aantal vergelijkingen) van het stelsel behouden en worden de vergelijkingen van een inactieve zone vervangen door pseudo-vergelijkingen (pseudo equations) [42, 26, 18, 22, 45, 47]. Deze zorgen ervoor dat inactieve zone-waarden de corresponderende waarden van de naburige zone volgen. Wanneer bijvoorbeeld de vloeistof-zone zou wegvallen, wordt de pseudo-vergelijking voor de buiswand: dTt,L = k [Tt,T P − Tt,L ] (3.4.5) dt Op deze manier volgt Tt,L de temperatuur van de naburige zone Tt,T P . Een speciaal geval is echter de zone-lengte. Wanneer een zone opnieuw intreedt, mag de zone-lengte niet nul zijn, maar moet ze een kleine positieve waarde aannemen (), anders worden de vergelijkingen singulier. Om dit te vermijden wordt de pseudo-vergelijking voor de zone-lengte: ddL = dL − dt

(3.4.6)

Met deze pseudo-vergelijkingen moet echter omzichtig omgesprongen worden, aangezien deze mogelijks instabiliteiten kunnen veroorzaken. Bendapudi et al. [4] maken hier bijvoorbeeld geen gebruik van. Wanneer een zone wegvalt, wordt het stelsel eenvoudigweg gereduceerd.


21

De eerste moving-boundary modellen [48] waren niet in staat om te gaan met modeveranderingen. Willatzen et al. [42] en Pettit et al. [26] brachten daar verandering in. Als eerste hebben deze onderzoekers een schakelcriterium voorgesteld en op het idee gekomen om pseudo-vergelijkingen te gebruiken. Zhang en Zhang [45] hebben dit model verbeterd door een tijdsafhankelijke void-fractie in te voeren. Zoals eerder vermeld, hebben ze ook verbeteringen aangebracht aan de schatting van de interface-temperaturen van de buiswand. Deze ingrepen hebben de robuustheid van het model verbeterd zodat het ook kan fungeren onder grote verstoringen van de randvoorwaarden. Ook wordt hierdoor de continu¨ıteit van de toestandsgrootheden verzekerd bij overgang naar een andere mode. Zhang et al. [47] hebben het moving-boundary model van Zhang en Zhang [45] aangepast voor flooded trommel-en-pijp warmtewisselaars. Ten slotte hebben McKinley en Alleyne [22] een ander schakelcriterium voorgesteld op basis van void-fractie en zone-lengtes, daarnaast is de void-fractie opgenomen in de toestandsvector. Deze aanpassingen zijn aangebracht om de stabiliteit en robuustheid van het model van Willatzen et al. [42] en Pettit et al. [26] te verbeteren. Het model van Li en Alleyne [18] is een uitbreiding op het werk van McKinley en Alleyne [22]. Er wordt een warmtewisselaarmodel ontwikkeld dat ook in staat is opstarten uitzetverschijnselen te beschrijven. Hiervoor zijn er extra modes ingevoerd en is het schakelcriterium uitgebreid. Vergelijking finite-volume en moving-boundary Het grote voordeel van een moving-boundary ten opzichte van een finite-volume simulatie is zijn snelheid. Dit is omdat er in een moving-boundary model maximaal 3 zones optreden in een warmtewisselaar, terwijl bij modellering volgens de finite-volume methode, toch minstens 15 controlevolumes vereist zijn om een vergelijkbare nauwkeurigheid te bekomen [4]. Er moet immers op gelet worden dat het mesh fijn genoeg is om de optredende zones te detecteren. Bendapudi et al. [4] hebben aangetoond dat in dit geval een finite-volume simulatie 3 keer langer duurt, hoewel er 5 keer meer vergelijkingen op te lossen zijn. Een nadeel van het lage aantal controlevolumes is dat dit iets minder nauwkeurige resultaten oplevert, aangezien de zone-karakteristieken vrij sterk gelumpt worden. Het moving-boundary model heeft ook moeilijkheden om sterk dynamische verschijnselen te beschrijven, zoals bij opstarten of uitzetten van een koelmachine. Met het werk van Willatzen et al. [42], Pettit et al. [26], Zhang en Zhang [45], Li en Alleyne [18] en McKinley en Alleyne [22] is hier evenwel verbetering in gekomen door invoering van het schakelcriterium en de pseudo-vergelijkingen.

3.5

Systeem

Eenmaal elk onderdeel van de cyclus is gemodelleerd, kan een cyclusmodel opgesteld worden. Dit gebeurt door ieder deelmodel te interconnecteren. Strikt genomen kan het dynamisch gedrag van een koelmachine enkel beschreven worden als ook alle componenten dynamisch zijn gemodelleerd. In de meeste gevallen volstaat het echter enkel de verdamperen condensordynamica in rekening te brengen. Uit experimenten blijkt immers dat de tijdsconstante van het systeem van dezelfde orde is als de tijdsconstante van de verdamper en condensor, die op zijn beurt een orde groter is dan de tijdsconstante van klep en compressor [4]. In de literatuur over cyclusmodellering ligt de klemtoon dan ook voornamelijk op de ontwikkeling van robuuste en snelle, maar voldoende nauwkeurige dynamische warmtewisselaarmodellen. Er wordt nauwelijks aandacht besteed aan een grondigere, eventueel zelfs dynamische modellering van klep en compressor [4, 48, 40, 47, 19]. De klep wordt gewoonlijk gemodelleerd volgens de klassieke methode, hoewel er betere methoden beschikbaar zijn, zoals toegelicht in sectie 3.2. De compressor wordt statisch gemodelleerd, zoals in sectie 3.3, hoewel voor opstarten uitzetverschijnselen dit onvoldoende is. Hieruit blijkt dat de huidige cyclusmodellen nog op vele vlakken kunnen verbeterd worden. Echter, hoe nauwkeuriger een model wordt, hoe meer


22

constructieve data er nodig zijn. Dit is vaak een probleem aangezien de interne opbouw van een klep of compressor confidentiële informatie is.

Hoofdstuk 4. Model koelmachine

23

Hoofdstuk 4

Model koelmachine 4.1

Inleiding

In hoofdstuk 3 zijn de verschillende modelleringstechnieken besproken. In dit hoofdstuk worden de componenten 1 voor 1 gemodelleerd. Eerst wordt duidelijk gemaakt waarom er voor welk model is gekozen. Daarna worden de modelvergelijkingen uitgewerkt.

4.2

Koellijn

Het koelmiddel dat door de koellijn loopt, neemt een deel van de verlieswarmte op die geproduceerd wordt in de motor en transmissie van de compressor. De hoeveelheid warmte die hier wordt geproduceerd is echter enkele grootte-ordes kleiner dan de warmte overgedragen in condensor en verdamper. Bijgevolg kan het effect van de koellijn op cyclusniveau verwaarloosd worden.

4.3 4.3.1

Thermostatische expansieklep Inleiding

In deze sectie wordt het model van de thermostatische expansieklep opgesteld. Er is gekozen voor het klassiek model. De stromingsmechanische methode geeft echter meer nauwkeurige resultaten en is fysisch ook meer gegrond. Dit eerste voordeel valt eigenlijk weg, aangezien de piloot-gestuurde klep vereenvoudigd wordt tot een gewone thermostatische expansieklep. De winst in nauwkeurigheid weegt dus niet op tegen de fout die door deze vereenvoudiging gemaakt wordt. Daarnaast is het, in tegenstelling tot het klassiek model, een zeer complexe methode. In het klassiek model daarentegen, wordt het gedrag van de klep gevat in slechts 2 algebra¨ısche vergelijkingen, eventueel aangevuld met 1 eerste orde differentiaalvergelijking voor de reservoirtemperatuur. De modellering is overgenomen uit het werk van Comstock en Braun [11], Wang et al. [40] en Mithraratne et al. [23]. Buiten het koelmiddel in het reservoir, zijn er geen details gekend over de klep.

4.3.2

Model

De klep die moet gemodelleerd worden is een thermostatische piloot-gestuurde expansieklep (figuur 2.3). Gezien het gebrek aan constructieve data, wordt deze vereenvoudigd tot een gewone thermostatische klep (figuur 3.1).


24

D d F2

y

stroming koelmiddel

θ

F1

Fs

Figuur 4.1: Krachtenevenwicht van de klepafsluiter

De stroming doorheen de klep wordt gemodelleerd als isenthalp, zodat: hr,in = hr,out

(4.3.1)

Verder wordt er verondersteld dat er geen massa-ophoping is in de klep. Het massadebiet wordt gegeven door de Bernoulli-klepvergelijking: q m ˙ r = Cv A 2ρr,in (pin − pout ) (4.3.2) pin en pout staan respectievelijk voor de druk aan inlaat en uitlaat van de klep. Deze drukken komen dus overeen met condensor- en verdamperdruk. • Doorstroomsectie A Er wordt uitgegaan van een conische afsluiter en circulaire zitting. Bijgevolg is: π D2 − d2 4

(4.3.3)

d = D − 2y tan(θ)

(4.3.4)

A= Met:

Zodat vergelijking (4.3.3) kan herschreven worden als: A = a1 y + a2 y 2

(4.3.5)

De klepheffing y wordt bepaald door het krachtenevenwicht van de klepafsluiter uit te drukken: Fs = F2 − F1 = (p2 − p1 )Am

(4.3.6)

In deze vergelijking zijn de dynamische krachten en de kracht ten gevolge van de drukval over de klepafsluiter verwaarloosd. In werkelijkheid heeft de klepafsluiter een zekere inertie, waardoor er dynamische krachten op inwerken. Deze inertie is echter zo klein dat de tijdsconstante van de klep minstens een orde kleiner is dan die van de warmtewisselaars. Op systeemniveau kan de dynamica van de klep dus verwaarloosd worden [4]. Ook de


25

kracht geassocieerd met de drukval over de klep wordt niet in rekening gebracht. Deze drukval is groot, maar de werkzame oppervlakte van de klepafsluiter is veel kleiner dan de oppervlakte van het membraan. Verder is F1 de kracht uitgeoefend door de verdamperdruk (p1 ) die op het membraan met oppervlakte Am inwerkt. Deze kracht wordt tegengewerkt door F2 en Fs . F2 resulteert uit de druk (p2 ) die het koelmiddel in het reservoir uitoefent op het membraan. Fs is de kracht uitgeoefend door de veer. In het geval van een lineaire veer, kan deze bepaald worden als: Fs = F0 + ks y

(4.3.7)

Waarin F0 de voorspanning is van de veer en ks de veerconstante. Met een schroefknop kan de voorspanning aangepast worden. Op deze manier wordt de oververhitting aan de uitlaat van de verdamper geregeld. Tijdens de metingen werd F0 niet gewijzigd. Vergelijking (4.3.6) kan omgevormd worden tot: y = k [p2 − (p1 + p0 )]

(4.3.8)

In deze vergelijking is k = Akm en p0 een equivalente druk zodat F0 = p0 Am . s Ten slotte wordt de druk in het reservoir bepaald als de saturatiedruk overeenkomend met de temperatuur van het koelmiddel in het reservoir (T2 ). Dit reservoir staat in thermisch contact met de uitlaat van de verdamper (Tevap ). Door het reservoir en zijn inhoud als 1 perfect ge¨ısoleerd geheel te beschouwen, kan deze temperatuur bepaald worden uit het energiebehoud: dT2 k2 A2 W2 = (Tevap − T2 ) (4.3.9) dt y2 In deze vergelijking is W2 de thermische capaciteit van het reservoir en zijn inhoud. k2 en y2 zijn respectievelijk de thermische conductiviteit en de dikte van de wand tussen de verdamperuitlaat en het reservoir. A2 is het beschikbare warmtewisselend oppervlak. Deze laatste 3 waarden worden constant verondersteld. Bijgevolg is: τb

dT2 = Tevap − T2 dt

(4.3.10)

τb is de tijdsconstante van het reservoir en zijn inhoud. In steady-state is de tijdsafgeleide nul, zodat de reservoirtemperatuur gelijk is aan de uitlaattemperatuur van de verdamper. • Klepcoëfficiënt Cv Deze parameter wordt constant verondersteld.

a1 a2 ymax Amax Amin k p0 Cv

0.0204 -0.41618 0.025 250e-06 0.5e-06 42e-08 80e+03 0.4

m2 /m m2 /m2 m m2 m2 m/P a Pa −

Tabel 4.1: Model van de klep De 5 vergelijkingen (4.3.1), (4.3.2), (4.3.5), (4.3.8) en (4.3.10) vormen samen het klepmodel. De factoren in vergelijking (4.3.5) werden bepaald door te veronderstellen dat θ = 20◦ en D = 17.8 mm [4]. Aangezien A een kwadratische functie is van y, moet y altijd kleiner


26

zijn dan een maximale waarde ymax , zodat steeds dA dy > 0. Fysisch betekent dit dat als de klepheffing toeneemt, ook de beschikbare doorstroomsectie stijgt. ymax is eenvoudigweg bepaald als de klepheffing waarbij de doorstroomsectie maximaal is (Amax ). Verder kan de klepheffing niet negatief zijn. Dit betekent dat y naar beneden moet begrensd worden. Deze minimum klepheffing resulteert in een minimum doorstroomsectie (Amin ). De waarden p0 , k en Cv zijn gefit aan de meetdata. De keuze van Cv is eigenlijk volledig arbitrair, aangezien de doorstroomsectie bepaald is aan de hand van een vooropgestelde geometrie. Alle berekende en veronderstelde waarden zijn weergegeven in tabel 4.1.

4.4 4.4.1

Centrifugaalcompressor Inleiding

In deze sectie wordt het model van de centrifugaalcompressor opgesteld. Er is gekozen voor een steady-state model. Om in te spelen op een wijzigende koellast kan er ingegrepen worden door het toerental van de rotor te wijzigen, ook kan de hoek van de inlaatschoepen aangepast worden. Soms worden ook beide technieken toegepast. In dit huidig werk wordt een centrifugaalcompressor gemodelleerd met inlaatschoepen, zonder frequentiesturing. Uitgezonderd bij opstarten en uitzetten, draait de compressor dus op constant toerental onafhankelijk van de koellast. Aangezien hier het opstarten uitzetverschijnsel niet bestudeerd wordt, moet de dynamica van de rotor dus niet in rekening gebracht worden. De hoek van de inlaatschoepen is echter wel afhankelijk van de PLR. Er wordt verondersteld dat hun stand quasi-onmiddellijk aangepast wordt wanneer het instelpunt van de PLR wijzigt. Bijgevolg kan ook de dynamica van het regelsysteem verwaarloosd worden. Verder is de hoeveelheid massa die zich in de compressor bevindt, beperkt. Dit betekent dat er nauwelijks massa-ophoping is, zodat het massabehoud zich reduceert tot een statische vergelijking die uitdrukt dat het ingaande massadebiet gelijk is aan het uitgaande massadebiet. Ook de impuls- en energievergelijking reduceren zich tot hun statische vorm. Er wordt immers van uitgegaan dat de verandering in de tijd van impuls en energie van het flu¨ıdum, ten gevolge van wijzigende randvoorwaarden, een orde sneller is in de compressor dan in de warmtewisselaars. Vervolgens wordt aangenomen dat de compressor perfect is ge¨ısoleerd, zodat er geen warmte-overdracht is naar de omgeving. In werkelijkheid is er altijd warmte-overdracht naar de omgeving. De warmte die op deze wijze de cyclus verlaat, is echter enkele grootte-ordes kleiner dan de warmte die wordt uitgewisseld in de verdamper en de condensor. Ten slotte wordt ook de thermische capaciteit van de compressorwand en de rotor verwaarloosd. Deze thermische massa is immers verwaarloosbaar ten opzichte van die van de warmtewisselaars. De modellering steunt op het werk van Bourdouxhe et al. [5]. De correlatie voor de polytrope exponent en de drukval over de inlaatschoepen zijn zelf opgesteld. Van de compressor is enkel het toerental en de diameter van de rotor gekend (zie tabel 2.2).

4.4.2

Model

Een compressor kan opgedeeld worden in 4 delen: • de inlaatschoepen; • de rotor; • de diffusor; • het slakkenhuis.


27

In een compressie-koelcyclus zijn de rotorschoepen typisch matig tot sterk achteruitgeheld, zodat de druktoename groot is, terwijl de snelheidstoename beperkt blijft. Dit betekent dat er in de diffusor slechts een kleine hoeveelheid druk herwonnen wordt. Bijgevolg verschilt de toestand aan de uitlaat van de rotor slechts weinig van de toestand aan de uitlaat van de diffusor. Verder wordt ook de invloed van het slakkenhuis verwaarloosd aangezien deze het gecomprimeerde gas collecteert en zo efficiënt mogelijk naar de uitlaat van de compressor leidt. Er treedt met andere woorden geen merkbare verandering op van de toestandsgrootheden. De toestand aan de uitlaat van de rotor is bijgevolg benaderend gelijk aan de toestand aan de uitlaat van de compressor. Deze toestand wordt aangeduid met het subscript out. De index in verwijst naar de inlaat van de compressor. Dit komt dus overeen met de inlaat van de inlaatschoepen. Ten slotte wordt de inlaat van de rotor aangeduid met het subscript x. De compressie in de rotor wordt gemodelleerd als een polytroop proces, zodat: pr,x pr,out = γ ργr,x ρr,out

(4.4.1)

De polytrope exponent is γ. Er wordt verondersteld dat deze grootheid enkel afhankelijk is van het massadebiet (m ˙ r ): (4.4.2) γ = a1 m ˙ ar 2 + a3 γ is dimensieloos en m ˙ r wordt uitgedrukt in [kg/s]. Deze betrekking, alsook de datapunten, zijn voorgesteld in figuur 4.3. De coëfficiënten a1 , a2 en a3 volgen uit een curve-fitting analyse en zijn opgelijst in tabel 4.2. Uit het werk van Bourdouxhe et al. [5] volgt een uitdrukking voor het massadebiet doorheen de compressor. Er is verondersteld dat de kinetische energie aan inlaat verwaarloosbaar is in de energiebalans over de volledige compressor. Bourdouxhe et al. [5] veronderstellen eveneens dat de compressie isentroop verloopt. Deze laatste vereenvoudiging is echter niet realistisch. Ten gevolge van irreversibiliteiten in de rotor is het beter de compressie als polytroop de modelleren. De resulterende formule blijft dezelfde, alleen moet de isentrope exponent vervangen worden door de polytrope exponent: A γ pr,x pr,out γ−1 γ 2 m ˙ r = ρr,out tan(β) −U (4.4.3) U γ − 1 ρr,x pr,x In deze vergelijking is A de doorstroomoppervlakte aan de uitlaat van de rotor. Deze is gelijk aan πDh. Daarnaast is U de tipsnelheid en is β de schoephoek aan de uitlaat van de rotor, zoals ge¨ıllustreerd in figuur 4.2.

Figuur 4.2: Schematische weergave van de rotor [5]


28 A β a1 a2 a3 b1 b2

0.00035 120 0.007768 -2.416 1.058 6.267e+07 -0.05462

m2 ◦

m−3 −

Tabel 4.2: Model van de compressor De toestand aan de inlaat van de rotor kan gevonden worden door de stroming doorheen de inlaatschoepen te modelleren als een equivalent smoringproces [5]. Dit betekent dat de stroming isenthalp kan beschouwd worden: hr,in = hr,x

(4.4.4)

De drukval is een functie van het debiet en wordt uitgedrukt door de Bernoulli-klepvergelijking. Na kwadratering en herschikking van beide leden, volgt: pr,x − pr,in = K

m ˙ 2r ρr,in

(4.4.5)

Wanneer de hoek van de inlaatschoepen vastligt, is K ongeveer constant, omdat de stroming 1-fasig is en de densiteit weinig verandert. De hoek van de inlaatschoepen is echter niet constant, deze worden aangestuurd door een onbekend regelsysteem. De taak van dit systeem is daarentegen wel gekend. De hoek van de schoepen wordt aangepast met als doel de ingestelde koellast (P LRset ) te realiseren. Het is dus aannemelijk te veronderstellen dat K een unieke functie is van P LRset . Uit de meetdata blijkt dat dit inderdaad het geval is. Curve-fitting levert volgend verband op: K = b1 eb2 P LRset (4.4.6) De dimensie van K is [m−3 ], terwijl P LRset in [%] is uitgedrukt. Er is verondersteld dat β = 120◦ (matige achteruithelling), A is begroot op 0.00035 m2 . In de meetdata zijn echter enkel de inen uitlaattoestand van de compressor en het massadebiet gegeven. Dit betekent dat eerst de drukval over de inlaatschoepen moet bepaald woorden aan de hand van vergelijking (4.4.3). Deze vergelijking wordt samen opgelost met (4.4.1). Daarna wordt K bepaald uit betrekking (4.4.5). Vervolgens zijn a1 , a2 , a3 , b1 en b2 bepaald door een datafit in de 27 verschillende steady-state meetpunten (zie figuur 4.3). De resultaten zijn opgelijst in tabel 4.2. In figuur 4.4 is de resulterende compressorkarakteristiek weergegeven. Deze is bepaald op inlaatvoorwaarden hr,in = 400 kJ/kg en pr,in = 1.5 bar. Het is duidelijk dat dit model niet in staat is om fenomenen zoals surge te voorspellen. Dit vormt in normale omstandigheden geen probleem omdat het werkingspunt dan voldoende ver ligt van surge. Tijdens opstarten en uitzetten of in abnormale omstandigheden kan het zijn dat surge wel optreedt, in dit geval is bovenstaand compressormodel onvoldoende.

4.5 4.5.1

Verdamper en condensor Inleiding

In deze sectie wordt het model van de verdamper en de condensor opgesteld. Er is gekozen voor een dynamisch moving-boundary model. Een steady-state aanpak is onvoldoende om het


29

1.1

datapunten a*x^b + c

γ [−]

1.09 1.08 1.07 1.06 0.6

0.8

1

1.2

1.4 mr [kg/s]

1.6

1.8

2

2.2

6

x 10

K [m−3]

15

a*exp(b*x) datapunten

10 5 0

30

40

50

60

70 80 PLRset [%]

90

100

110

120

Figuur 4.3: Datafit compressor dynamische gedrag van een cyclus te beschrijven. De tijdsconstante van beide warmtewisselaars is immers van dezelfde grootte-orde als van het systeem. Daarnaast is geopteerd voor het moving-boundary model omwille van de hogere uitvoeringssnelheid ten opzichte van een finite-volume model dat ongeveer dezelfde nauwkeurigheid heeft. De modellering is gebaseerd op het werk van Bendapudi et al. [4], Pettit et al. [26], Willatzen et al. [42], Zhang en Zhang [45] en Zhang et al. [47]. Zowel verdamper als condensor zijn 2-passige trommel-en-pijp warmtewisselaars. Het secundair flu¨ıdum is water en stroomt door de pijpen, het koelmiddel is R134a en loopt door de trommel. Het aantal buizen en hun afmetingen zijn gekend, ook de diameter van de trommel is gegeven (zie tabel 2.1).

4.5.2 4.5.2.1

Model Modelleringsmethodiek

Om het overzicht te behouden, wordt hier de modelleringsmethodiek kort toegelicht. Eerst worden de algemene behoudsvergelijkingen van flu¨ıda en vaste stoffen herhaald. Daarna worden deze vereenvoudigd voor de buiswand, de trommelwand, het koelmiddel en het secundair flu¨ıdum. Dit resulteert in een stelsel partiële differentiaalvergelijkingen, die de vereenvoudigde behoudsvergelijkingen genoemd worden. In de volgende stap wordt dit omgevormd tot een stelsel 1ste orde differentiaalvergelijkingen door integratie over de verschillende zones die in de warmtewisselaar optreden. Dit zijn de zone-specifieke behoudsvergelijkingen. Ten slotte wordt een schakelcriterium opgesteld voor verdamper en condensor. De vereenvoudigde behoudsvergelijkingen zijn identiek aan die in [4, 42, 26, 45], aangezien dezelfde veronderstellingen gemaakt worden. Echter, omdat de warmtewisselaars in de te modelleren proefstand in dwarsstroom werken en niet in mee- of tegenstroom, zoals verondersteld in bovengenoemde werken, zien de zone-specifieke vergelijkingen opgesteld in 4.5.2.4 er anders


30

5

10

x 10

U = 200 m/s U = 180 m/s U = 220 m/s

9

pr,out [Pa]

8

7

6

5

4

3 0.4

0.5

0.6

0.7 mr [kg/s]

0.8

0.9

1

Figuur 4.4: Compressorkarakteristiek uit. In de literatuur is er slechts 1 werk terug te vinden waarin een moving-boundary model opgesteld wordt voor dwarsstroom-warmtewisselaars [47]. Zhang et al. [47] hebben bij het opstellen van de behoudsvergelijkingen van de wand echter geen rekening gehouden met tijdsafhankelijke zone-grenzen. Dit zondigt evenwel aan het energiebehoud. Daarnaast hebben ze, in tegenstelling tot dit huidig werk, verondersteld dat de fasen volledig gescheiden zijn, zodat er slechts 2 zones zijn in de condensor. 4.5.2.2

Algemene behoudsvergelijkingen

Flu¨ıda De behoudswetten voor flu¨ıda, uitgedrukt in een cartesiaans assenstelsel (x,y,z) zijn [39]: X ∂ ∂ ρ+ (ρvi ) = 0 (4.5.1) ∂t ∂i i=x,y,z

X ∂ X ∂ ∂ ∂ (ρvj ) + (vi vj ρ) + p = ρgj + τij ∂t ∂i ∂j ∂i i=x,y,z

(j = x, y, z)   X X X ∂ ∂ ∂ ∂ ∂ ρ h + qi v i h −  p + vi p = φ + ∂t ∂i ∂t ∂i ∂i 

(4.5.2)

i=x,y,z



i=x,y,z

i=x,y,z

(4.5.3)

i=x,y,z

Deze vergelijkingen zijn in de literatuur beter bekend als de Navier-Stokes vergelijkingen. Het massabehoud (4.5.1) is uitgedrukt in [ mkg3 s ], het impulsbehoud (4.5.2) in [ mN3 ] en het energiebeW houd (4.5.3) in [ m ıdum op de i-as, wordt 3 ]. De projectie van de absolute snelheid van het flu¨ aangeduid met vi . ρ en h zijn respectievelijk densiteit en enthalpie. gi is de component van de gravitationele versnelling volgens de i-as. qi is de warmteflux in de richting i. φ is de zogenaamde dissipatiefunctie en is gelijk aan de vervormingsarbeid van de wrijvingskrachten. Dit is echter irreversibele arbeid en resulteert bijgevolg in entropiegeneratie. τ is de viscositeitsspanningstensor. In het geval van lineair-isotrope flu¨ıda is volgens de hypothese van Stokes τij = 2µγij met µ de dynamische viscositeitscoëfficiënt en 1 ∂ ∂ 1 X ∂ γij = vj + vi − δij vk (4.5.4) 2 ∂i ∂j 3 ∂k k=x,y,z

Hoofdstuk 4. Model koelmachine ( 1 waarin δij = 0

31

als i = j . als i = 6 j

Vaste stoffen Het energiebehoud voor vaste stoffen, uitgedrukt in een cartesiaans assenstelsel (x,y,z) is [14]: X ∂ ∂ (ρCT ) + qi = 0 (4.5.5) ∂t ∂i i=x,y,z

Er is verondersteld dat er geen interne warmtegeneratie plaatsvindt in het materiaal. Deze verW gelijking heeft dimensie [ m 3 ]. ρ is de massadichtheid van de stof, terwijl C de warmtecapaciteit voorstelt. Volgens de wet van Fourier, is: ∂ ∂T ∂ k (4.5.6) qi = − ∂i ∂i ∂i Hierin is k de thermische conductiviteit van het materiaal. 4.5.2.3

Vereenvoudigde behoudsvergelijkingen

Trommel en pijp Zowel in het materiaal van de trommel als van de pijp, worden volgende veronderstellingen gemaakt [4, 42, 26, 45]. 1. De materiaaleigenschappen zijn constant. De temperatuursafhankelijkheid van de eigenschappen van metalen is immers beperkt, zeker in het temperatuursgebied dat optreedt in een warmtewisselaar. 2. De radiale conductieve thermische weerstand is nul. Gezien de hoge thermische conductiviteit van metalen is deze weerstand verwaarloosbaar ten opzichte van de weerstand tussen wand en het koelmiddel of secundair flu¨ıdum. Dit impliceert dat de temperatuur van de wand constant is in de radiale richting. 3. De warmte-overdracht naar de trommel is nul. Bendapudi et al. [4] hebben aangetoond dat dit een aanvaardbare veronderstelling is. Het warmtewisselend oppervlak van de trommel is immers zeer klein ten opzichte van het totale oppervlak van de pijpen.

Figuur 4.5: Infinitesimaal volume dVt in een pijp Door het verwaarlozen van de warmtewisseling met de trommel, valt de energievergelijking van de trommel weg uit het stelsel. In een infinitesimaal volume dVt van een pijp (figuur 4.5), geldt: ∂Tt ∂ 2 Tt ρ t Ct dVt = dQt + kt 2 dVt (4.5.7) ∂t ∂z


32

Deze vergelijking is uitgedrukt in [W ]. dQt is de warmte-overdracht naar het infinitesimaal volume. Deze term is bijgevolg de som van de warmte-overdracht door convectie naar de pijp vanuit het koelmiddel en vanuit het secundair flu¨ıdum. Isolatie Er wordt verondersteld dat de isolatie perfect is. Dit wilt zeggen dat er geen warmtewisseling is met de omgeving. In werkelijkheid is de isolatie niet perfect, maar toch voldoende zodat het warmteverlies naar de omgeving enkele grootte-ordes kleiner is dan de warmte-transfer tussen het secundair flu¨ıdum en het koelmiddel. Koelmiddel In het koelmiddel worden volgende veronderstellingen gemaakt [4, 42, 26, 45]. 1. De stroming is 1-dimensionaal. Door de aanwezigheid van een pijpenbundel en tussenschotten is de stroming in de trommel sterk 3-dimensionaal. Daarnaast kan er een 2-fasige zone optreden in de warmtewisselaar. Toch is er een dominante stromingsrichting tussen inen uitlaat. 2. De drukval is nul. De doorstroomsectie van het koelmiddel is vrij groot en de afgelegde weg van in- naar uitlaat is klein, zodat de drukval beperkt blijft. Verder zijn de eigenschappen van het koelmiddel R134a weinig gevoelig voor de kleine drukval die in werkelijkheid optreedt in de warmtewisselaar. De drukval heeft dus weinig effect op de prestaties van de warmtewisselaar. Bijgevolg loont het niet de moeite om de drukval in rekening te brengen. Ten slotte hebben Zhang et al. [46] aangetoond dat deze vereenvoudiging op systeemniveau nauwelijks verschil maakt. Deze veronderstelling impliceert verder dat ook alle visceuze effecten kunnen verwaarloosd worden, aangezien deze altijd drukval met zich meebrengen. 3. De diffusieve warmtetransfer is nul. Het massadebiet van het koelmiddel is immers voldoende groot zodat de diffusieve warmtevloed in het koelmiddel verwaarloosbaar is ten opzichte van de convectieve warmtetransfer tussen het koelmiddel en de buiswand. Met andere woorden: het Péclet-getal is voldoende groot om het effect van diffusieve warmte-overdracht te kunnen verwaarlozen ten op zichte van de convectieve warmteoverdracht. 4. Het koelmiddel is zuiver of een azeotroop mengsel. Dit is inderdaad voldaan in het geval van R134a. Dit impliceert dat het toestandsprincipe kan toegepast worden: elke toestandsgrootheid is uit te drukken in functie van 2 andere onafhankelijke toestandsgrootheden. Deze eigenschap gaat niet op voor zeotrope mengsels. Bij verdamping of condensatie, veranderen immers de volumefracties van de componenten. Dit betekent dat, naast de 2 onafhankelijke toestandsgrootheden, ook nog de mengverhouding moet opgegeven worden. Door het verwaarlozen van de drukval is de impulsvergelijking weggevallen uit het stelsel. De massa- en energievergelijking voor een infinitesimaal volume dVr in het koelmiddel (figuur 4.6), reduceren zich respectievelijk tot: ∂ρr ∂ρr vr dVr + dVr = 0 ∂t ∂z ∂ρr hr ∂ρr vr hr ∂pr dVr + dVr = dVr + dQr ∂t ∂z ∂t

(4.5.8) (4.5.9)


33

Figuur 4.6: Infinitesimaal volume dVr in het koelmiddel De massavergelijking is uitgedrukt in [kg/s], terwijl de dimensie van de energievergelijking [W ] is. dQt is de warmte-overdracht naar het infinitesimaal volume. Deze term is bijgevolg de warmte-overdracht door convectie naar het koelmiddel vanuit de pijp. Secundair flu¨ıdum In het secundair flu¨ıdum worden volgende veronderstellingen gemaakt [4, 42, 26, 45]. 1. De stroming is 1-dimensionaal. Aangezien het secundair flu¨ıdum in de pijpen loopt en zijn fase behoudt, is dit een aannemelijke veronderstelling. 2. De drukval is nul. Ook hier loont het niet om de weinige drukval die in werkelijkheid optreedt in rekening te brengen. De eigenschappen van water zijn immers zeer zwak afhankelijk van de druk, waardoor het effect van de druk op de prestaties van de warmtewisselaar zeer klein is. Ook hier impliceert deze aanname dat de visceuze effecten kunnen verwaarloosd worden. 3. De diffusieve warmtetransfer is nul. Het massadebiet van het secundair flu¨ıdum is immers voldoende groot om axiale diffusieve warmtevloed te kunnen verwaarlozen ten opzichte van radiale convectieve warmtetransfer. 4. De stofeigenschappen zijn constant. Dit betekent dat de specifieke warmtecapaciteit constant is en het flu¨ıdum als onsamendrukbaar kan beschouwd worden. De drukval die het secundair flu¨ıdum ondergaat is typisch zeer klein. Daarnaast is het temperatuursverschil tussen inen uitlaat maximaal enkele graden. Onder deze omstandigheden blijven de stofeigenschappen van water ongeveer constant, wat deze veronderstelling verrechtvaardigt.

Figuur 4.7: Infinitesimaal volume dV2 in het secundair flu¨ıdum


34

Door het verwaarlozen van de drukval is de impulsvergelijking weggevallen uit het stelsel. De massa- en energievergelijking voor een infinitesimaal volume dV2 in het secundair flu¨ıdum (figuur 4.7), reduceren zich respectievelijk tot: ∂v2 dV2 = 0 ∂z ∂T2 ∂v2 T2 ∂p2 ρ2 Cp,2 dV2 + ρ2 Cp,2 dV2 = dV2 + dQ2 ∂t ∂z ∂t ρ2

(4.5.10) (4.5.11)

De massavergelijking is uitgedrukt in [kg/s], terwijl de dimensie van de energievergelijking [W ] is. dQ2 is de warmte-overdracht naar het infinitesimaal volume. Deze term is bijgevolg de warmte-overdracht door convectie naar het secundair flu¨ıdum vanuit de pijp. 4.5.2.4

Zone-specifieke behoudsvergelijkingen

In een warmtewisselaars kunnen 3 zones onderscheiden worden: 1. vloeistofzone (L) • in deze zone bevindt het koelmiddel zich in vloeibare toestand; • het aantal pijpen dat zich in deze zone bevindt, is een fractie FL van het totaal aantal pijpen in de warmtewisselaar. 2. 2-fasenzone (TP ) • in deze zone bevindt het koelmiddel zich in 2-fasige toestand; • het aantal pijpen dat zich in deze zone bevindt, is een fractie FT P van het totaal aantal pijpen in de warmtewisselaar. 3. dampzone (V ) • in deze zone bevindt het koelmiddel zich in dampvormige toestand; • het aantal pijpen dat zich in deze zone bevindt, is een fractie FV van het totaal aantal pijpen in de warmtewisselaar. Deze 3 zones zijn disjunct, maar nemen samen de volledige warmtewisselaar in. Dit betekent dat FL + FT P + FV = 1. Het gemeenschappelijk oppervlak tussen 2 controlevolumes wordt de interface genoemd. Zo is er een TP-V en een L-TP interface in de verdamper. De interfaces in de condensor worden aangeduid met V-TP en TP-L. Het massadebiet dat door de interface stroomt, wordt genoteerd als m ˙ int , waarbij int een interface is. Als gevolg van de zone-opdeling is de toestand van het koelmiddel op een interface altijd een gesatureerde toestand. Op de TP-V en de V-TP interfaces is het koelmiddel gesatureerde damp, terwijl de toestand op de L-TP en de TP-L interfaces gesatureerde vloeistof is. Iedere zone wordt op zijn beurt onderverdeeld in 3 disjuncte controlevolumes: 1. rond het secundair flu¨ıdum (2 ); 2. rond de buiswand (t); 3. rond het koelmiddel (r ). De behoudswetten, zoals opgesteld in de vorige paragraaf, worden ge¨ıntegreerd over elk controlevolume. Zo worden de zone-specifieke vergelijkingen bekomen. Dit is een stelsel 1ste orde differentiaalvergelijkingen: ˙ A(x)x˙ = f(x; u; u) (4.5.12)


35

x is de toestandsvector en bevat de toestand van het koelmiddel en secundair flu¨ıdum aan inen uitlaat van de warmtewisselaar, maar ook de interne variabelen. Met de interne variabelen wordt bedoeld de grootte van iedere zone en de wandtemperatuur van de pijpen in iedere zone. u en u˙ bevatten de randvoorwaarden. De integratiegrenzen komen overeen met de zone-grenzen. Deze zijn tijdsafhankelijk zodat bij integratie de regel van Leibniz moet toegepast worden, zoals beschreven in bijlage D. De integratie wordt sterk vereenvoudigd door de dubbele pas te modelleren als een enkele pas met dubbel zoveel pijpen. Zonder deze vereenvoudiging is het trouwens noodzakelijk dat de schikking van de buizen in de bundel gekend is. Anders is het onmogelijk om de randvoorwaarden van het secundair flu¨ıdum van de ene zone in verband te brengen met die van een andere zone. Immers, als er 2 zones zijn die een gemeenschappelijke pijp hebben, komt de inlaattemperatuur van het secundair flu¨ıdum in de ene zone overeen met de uitlaattemperatuur in de andere zone. Echter, om na te gaan of er dergelijke zones zijn, moet de schikking van de pijpen in de bundel gekend zijn. Verdamper In wat volgt, worden de behoudsvergelijkingen van de verdamper afgeleid. Er wordt verondersteld dat de warmtewisselaar zich in de TP-V-mode bevindt, zoals in figuur 4.8. Er wordt met andere woorden van uitgegaan dat enkel de TP- en de V-zone optreden. Bijgevolg is FL = 0 en FT P + FV = 1. Deze afleiding is gebaseerd op het werk van [4, 42, 26, 45]. Enkel de behoudsvergelijkingen voor secundair flu¨ıdum en de buiswand verschillen. R,out

R,out

V 2,in

2,out TP

R,in

R,in

Figuur 4.8: De verdamper in de TP-V-mode

Buiswand Integratie van het energiebehoud in de buiswand, respectievelijk in de TP- en de V-zone,


36

geeft:

dTt,T P dFT P FT P + FT P (Tt,T P − Tt,V ) dt dt dTt,V dFV FV + FV (Tt,V − Tt,T P ) = ρt Ct Vt dt dt

Q2t,T P + Qrt,T P = ρt Ct Vt Q2t,V + Qrt,V

(4.5.13) (4.5.14)

De totale warmte-overdracht naar de buiswand is de som van de warmte-overdracht uit het koelmiddel (Qrt ) en uit het secundair flu¨ıdum (Q2t ). Vt is het totale volume buiswand in de warmtewisselaar. Gezien de definitie van FT P en FV , ligt hiervan een fractie FT P in de TPzone en een fractie FV in de V-zone. De temperatuur op de interface tussen 2 zones is bepaald als een gewogen gemiddelde tussen beide zone-gemiddelde temperaturen [45]: Tt,T P −V = FV Tt,T P + FT P Tt,V

(4.5.15)

Wanneer bijvoorbeeld FV vergroot, komt de interface-temperatuur dichter bij Tt,T P te liggen. Dit kan ook fysisch verklaard worden omdat de interface nu ligt waar voordien de TP-zone lag. Tt,T P en Tt,V zijn de zone-gemiddelde wandtemperaturen en zijn gedefinieerd als: Z 1 Tt,T P = Tt dV VT P VT P Z 1 Tt,V = Tt dV VV VV

(4.5.16) (4.5.17)

Secundair flu¨ıdum Integratie van het energiebehoud in het secundair flu¨ıdum, respectievelijk in de TP- en de V-zone, geeft: dT2,T P dFT P −Q2t,T P = ρ2 Cp,2 V2 FT P + FT P (T2,T P − T2,V ) (4.5.18) dt dt +m ˙ 2,T P Cp,2 [(T2,out )T P − T2,in ] dT2,V dFV −Q2t,V = ρ2 Cp,2 V2 FV + FV (T2,V − T2,T P ) (4.5.19) dt dt +m ˙ 2,V Cp,2 [(T2,out )V − T2,in ] De warmte-overdracht naar het secundair flu¨ıdum vanuit de buiswand, is aangeduid met −Q2t . V2 is het totale volume secundair flu¨ıdum in de warmtewisselaar. Gezien de definitie van FT P en FV en het feit dat het secundair flu¨ıdum door de pijpen loopt, ligt hiervan een fractie FT P in de TP-zone en een fractie FV in de V-zone. Naar analogie met de buiswand, wordt de temperatuur op de interface tussen 2 zones bepaald als een gewogen gemiddelde tussen beide zone-gemiddelde temperaturen: T2,T P −V = FV T2,T P + FT P T2,V

(4.5.20)

Verder is er verondersteld dat het totale massadebiet (m ˙ 2 ) gelijk verdeeld wordt over alle pijpen, bijgevolg is: m ˙ 2,T P = FT P m ˙2

(4.5.21)

m ˙ 2,V = FV m ˙2

(4.5.22)


37

T2,T P en T2,V zijn de zone-gemiddelde temperaturen van het secundair flu¨ıdum en zijn gedefinieerd als: Z 1 T2,T P = T2 dV (4.5.23) VT P VT P Z 1 T2,V = T2 dV (4.5.24) VV VV Het temperatuursverloop tussen inen uitlaat wordt lineair beschouwd, zodat: 1 [T2,in + (T2,out )T P ] 2 1 = [T2,in + (T2,out )V ] 2

T2,T P = T2,V

(4.5.25) (4.5.26)

Dit resulteert in: dT2,T P 1 dT2,in d(T2,out )T P = + dt 2 dt dt dT2,V 1 dT2,in d(T2,out )V = + dt 2 dt dt

(4.5.27) (4.5.28)

Ten slotte wordt de uitlaattemperatuur van het secundair flu¨ıdum bepaald als: T2,out = FV (T2,out )V + FT P (T2,out )T P

(4.5.29)

Dit volgt uit het energiebehoud bij menging en de aanname dat Cp,2 aan de uitlaat van iedere zone gelijk is. Koelmiddel Integratie van het massabehoud in het koelmiddel, respectievelijk in de TP- en de V-zone, geeft: m ˙ r,in − m ˙ r,T P −V dFT P = [(¯ γ − 1)(ρg,r − ρf,r )] Vr dt dρf,r dρg,r dpr ∂¯ γ + γ¯ + (1 − γ¯ ) + (ρg,r − ρf,r ) FT P dt dpr dpr ∂pr dhr,in ∂¯ γ + (ρg,r − ρf,r ) FT P dt ∂hr,in m ˙ r,T P −V − m ˙ r,out dFV = [ρr,V − ρg,r ] Vr dt ∂ρr,V dhg,r ∂ρr,V dpr 1/2 + FV + dt ∂hr,V dpr ∂pr ∂ρr,V dhr,out + 1/2 FV dt ∂hr,V

(4.5.30)

(4.5.31)


38

Optellen van beide vergelijkingen resulteert in: m ˙ r,in − m ˙ r,out dFT P dFV = [(¯ γ − 1)(ρg,r − ρf,r )] + [ρr,V − ρg,r ] Vr dt dt ∂ρr,V dhg,r ∂ρr,V dpr FV 1/2 + + dt ∂hr,V dpr ∂pr dρf,r dρg,r ∂¯ γ + (1 − γ¯ ) + (ρg,r − ρf,r ) +FT P γ¯ dpr dpr ∂pr dhr,in ∂¯ γ + (ρg,r − ρf,r ) FT P dt ∂hr,in ∂ρr,V dhr,out + 1/2 FV dt ∂hr,V (4.5.32) Integratie van het energiebehoud van het koelmiddel, respectievelijk in de TP- en de V-zone, geeft: m ˙ r,in hr,in − m ˙ r,T P −V hg,r Qrt,T P dFT P − = [(¯ γ − 1) ((ρh)g,r − (ρh)f,r )] (4.5.33) Vr Vr dt d(ρh)f,r d(ρh)g,r dpr + γ¯ + (1 − γ¯ ) dt dpr dpr ∂¯ γ ((ρh)g,r − (ρh)f,r ) − 1 FT P + ∂pr dhr,in ∂¯ γ ((ρh)g,r − (ρh)f,r ) FT P + dt ∂hr,in m ˙ r,T P −V hg,r − m ˙ r,out hr,out Qrt,V dFV − = [(ρh)r,V − (ρh)g,r ] (4.5.34) Vr Vr dt ∂(ρh)r,V dhg,r ∂(ρh)r,V dpr + 1/2 + − 1 FV dt ∂hr,V dpr ∂pr ∂(ρh)r,V dhr,out 1/2 + FV dt ∂hr,V De warmte-overdracht naar het koelmiddel vanuit de buiswand door convectie, is aangeduid met −Qrt . Vr is het totale volume koelmiddel in de warmtewisselaar. Hoewel Fzone gedefinieerd is als het aantal buizen dat zich in een zone bevindt, is verondersteld dat deze grootheid een goede benadering is voor de fractie van Vr die het koelmiddel in deze zone inneemt. Strikt genomen is dit enkel geldig als de buizen uniform verdeeld zijn over de dwarssectie van de warmtewisselaar. Dit is in realiteit echter niet het geval. Deze vereenvoudiging is gemaakt omdat het exacte bundelpatroon niet gekend is. ρr,V , ρr,T P , (ρh)r,V en (ρh)r,T P zijn allen zone-gemiddelde waarden: Z 1 ρr dV ρr,T P = VT P VT P Z 1 ρr,V = ρr dV VV VV Z 1 (ρh)r,T P = (ρh)r dV VT P VT P Z 1 (ρh)r,V = (ρh)r dV VV VV

(4.5.35) (4.5.36) (4.5.37) (4.5.38)


39

ρr,V en (ρh)r,V worden berekend op de zone-gemiddelde enthalpie in de V-zone (hr,V ) en verdamperdruk (pr ). ρr,T P en (ρh)r,T P daarentegen, worden bepaald in functie van de zonegemiddelde void-fractie γ¯ : 1

ρr,T P =

Z

ρr dV

VT P V Z T P 1 ρf,r (1 − γ) + ρg,r γdV = VT P VT P 1 [ρf,r (1 − γ¯ ) + ρg,r γ¯ ] = VT P 1

(ρh)r,T P =

Z

(ρh)r dV VT P VT P Z 1 = (ρh)f,r (1 − γ) + (ρh)g,r γdV VT P VT P 1 [(ρh)f,r (1 − γ¯ ) + (ρh)g,r γ¯ ] = VT P

Met: γ¯ =

(4.5.39)

1 VT P

(4.5.40)

Z γdV

(4.5.41)

VT P

Deze waarde is enkel afhankelijk van de inlaatenthalpie en de druk van de TP-zone, zodat: d¯ γ ∂¯ γ dp¯r ∂¯ γ dhr,in = + dt ∂pr dt ∂hr,in dt

(4.5.42)

Er wordt verder verondersteld dat het enthalpieverloop in iedere zone lineair is (figuur 4.9): 1 [hr,in + hg,r ] (4.5.43) 2 1 hr,V = [hg,r + hr,out ] (4.5.44) 2 In steady-state omstandigheden is dit een goede benadering, zoals aangetoond door Wedekind et al. [41]. Of dit ook het geval is onder dynamische voorwaarden, blijft een open vraag. Deze veronderstelling komt echter zeer vaak terug in de literatuur en geeft goede resultaten [4, 22, 18, 42, 26, 45]. hr,T P =

h

TP

V

hr,out hg,r hr,in

V 0

VTP VTP + VV

Figuur 4.9: Lineair verloop enthalpie in de verdamper Een lineair enthalpieverloop in de TP-zone impliceert eveneens een lineair verloop van de dampfractie. De zone-gemiddelde void-fractie kan dus berekend worden als: 1 γ¯ = xr,out − xr,in

xZ r,out

γdx xr,in

(4.5.45)


40

Deze veronderstelling heeft ook als gevolg dat: ∂ρr,V dpr ∂ρr,V 1 dhr,out dhg,r dpr dρr,V = + + dt ∂pr dt ∂hr,V 2 dt dpr dt d(ρh)r,V ∂(ρh)r,V dpr ∂(ρh)r,V 1 dhr,out dhg,r dpr = + + dt ∂pr dt ∂hr,V 2 dt dpr dt

(4.5.46) (4.5.47)

De modelvergelijkingen zijn afhankelijk van de mode die optreedt. Bovenstaande vergelijkingen zijn afgeleid in de veronderstelling dat de verdamper zich steeds in de TP-V-mode bevindt. Dit is typisch de meest voorkomende werkingsmode van een verdamper. Door grote verstoringen in de randvoorwaarden, of gedurende een overgangsverschijnsel, kan de V-zone echter verdwijnen. De verdamper werkt dan in de TP-mode. In dit geval verloopt de afleiding volledig analoog. Het resultaat hiervan is terug te vinden in bijlage B. Om na te gaan wanneer een zone verdwijnt of bijkomt, moet de uitlaatenthalpie van de verdamper (hr,out ) gecontroleerd worden. Indien de verdamper in de TP-V-mode werkt, moet hr,out groter zijn dan de damp saturatie-enthalpie (hg,r ) bepaald op verdamperdruk pr . Indien dit niet het geval is, moet geschakeld worden naar de TP-mode. In deze mode is hr,out altijd kleiner dan hg,r . Dit schakelcriterium is schematisch weergegeven in figuur 4.10.

TP

TP-V

Figuur 4.10: Het schakelcriterium in de verdamper Wanneer de V-zone wegvalt, verdwijnen er 3 differentiaalvergelijkingen uit het stelsel. Om de vorm van het stelsel te behouden, worden deze vervangen door een set pseudo-vergelijkingen: kF

dFV = − FV dt

dTt,V = Tt,T P − Tt,V dt d(T2,out )V k2 = (T2,out )T P − (T2,out )V dt kt

(4.5.48) (4.5.49) (4.5.50)

De wandtemperatuur en de uitlaattemperatuur van het secundair flu¨ıdum in de inactieve Vzone, volgen als het ware de overeenkomstige waarden van de naburige actieve TP-zone. Wanneer teruggeschakeld wordt naar TP-V, worden zo de stofeigenschappen in de V-zone bepaald op fysisch realistische waarden. De volumefractie FV , mag echter niet op nul ge¨ınitialiseerd worden. Dit zou immers singulariteiten met zich meebrengen bij herintrede van de V-zone. Vandaar dat ze asymptotisch naar een positieve waarde nadert. moet echter voldoende klein zijn opdat haar invloed op de resultaten beperkt zou blijven. Condensor In wat volgt, worden de behoudsvergelijkingen van de condensor afgeleid. Er wordt verondersteld dat de warmtewisselaar zich in de V-TP-L-mode bevindt, zoals in figuur 4.11. Er wordt met andere woorden van uitgegaan dat zowel de V-, de TP- als de L-zone optreden. De bekomen vergelijkingen zijn analoog aan die in de verdamper. Deze afleiding is gebaseerd op het werk van [4, 42, 26, 45]. Enkel de behoudsvergelijkingen voor secundair flu¨ıdum en de buiswand verschillen.


41

R,in

R,in

V 2,in

2,out TP

L

R,out

R,out

Figuur 4.11: De condensor in de V-TP-L-mode

Buiswand Integratie van het energiebehoud in de buiswand, respectievelijk in de V-, TP- en de L-zone, geeft: dTt,V dFV FV Q2t,V + Qrt,V = ρt Ct Vt FV + (Tt,V − Tt,T P ) (4.5.51) dt dt FV + FT P dTt,T P dFV FV Tt,T P + FT P Tt,V Q2t,T P + Qrt,T P = ρt Ct Vt FT P + (4.5.52) dt dt FV + FT P dFT P dFL FL Tt,T P + FT P Tt,L + Tt,T P + dt dt FL + FT P dTt,L dFL FL Q2t,L + Qrt,L = ρt Ct Vt FL + (Tt,L − Tt,T P ) (4.5.53) dt dt FL + FT P De totale warmte-overdracht naar de buiswand is de som van de warmte-overdracht uit het koelmiddel (Qrt ) en uit het secundair flu¨ıdum (Q2t ). Vt is het totale volume buiswand in de warmtewisselaar. Gezien de definitie van FV , FT P en FL , ligt hiervan een fractie FV in de V-zone, een fractie FT P in de TP-zone en een fractie FL in de L-zone. De temperatuur op de interface tussen 2 zones is bepaald als een gewogen gemiddelde tussen beide zone-gemiddelde temperaturen: FV Tt,T P + FT P Tt,V FV + FT P FT P Tt,L + FL Tt,T P = FT P + FL

Tt,V −T P =

(4.5.54)

Tt,T P −L

(4.5.55)

Tt,V , Tt,T P en Tt,L zijn de zone-gemiddelde wandtemperaturen. Secundair flu¨ıdum Integratie van het energiebehoud in het secundair flu¨ıdum, respectievelijk in de V-, TP- en


42

de L-zone, geeft:

dT2,V FV dFV (T2,V − T2,T P ) (4.5.56) −Q2t,V = ρ2 Cp,2 V2 FV + dt dt FV + FT P + FV m ˙ 2 Cp,2 [(T2,out )V − T2,in ] dT2,T P dFV FV T2,T P + FT P T2,V −Q2t,T P = ρ2 Cp,2 V2 (4.5.57) FT P + dt dt FV + FT P dFT P dFL FL T2,T P + FT P T2,L + T2,T P + dt dt FL + FT P + FT P m ˙ 2 Cp,2 [(T2,out )T P − T2,in ] dT2,L FL dFL (T2,L − T2,T P ) (4.5.58) FL + −Q2t,L = ρ2 Cp,2 V2 dt dt FL + FT P + FL m ˙ 2 Cp,2 [(T2,out )L − T2,in ] De warmte-overdracht naar het secundair flu¨ıdum vanuit de buiswand door convectie, is aangeduid met −Q2t . V2 is het totale volume secundair flu¨ıdum in de warmtewisselaar. Gezien de definitie van FV , FT P en FL en het feit dat het secundair flu¨ıdum door de pijpen loopt, ligt hiervan een fractie FV in de V-zone, een fractie FT P in de TP-zone en een fractie FL in de L-zone. Naar analogie met de buiswand, wordt de temperatuur op de interface tussen 2 zones bepaald als een gewogen gemiddelde tussen beide zone-gemiddelde temperaturen: FV T2,T P + FT P T2,V FV + FT P FT P T2,L + FL T2,T P = FT P + FL

T2,V −T P =

(4.5.59)

T2,T P −L

(4.5.60)

Zoals in de verdamper, is er verondersteld dat het totale massadebiet (m ˙ 2 ) gelijk verdeeld wordt over alle pijpen: m ˙ 2,V = FV m ˙2

(4.5.61)

m ˙ 2,T P = FT P m ˙2

(4.5.62)

m ˙ 2,L = FL m ˙2

(4.5.63)

T2,V , T2,T P en T2,L zijn de zone-gemiddelde temperaturen van het secundair flu¨ıdum. Het temperatuursverloop tussen inen uitlaat wordt opnieuw lineair beschouwd, zodat: 1 [T2,in + (T2,out )V ] 2 1 = [T2,in + (T2,out )T P ] 2 1 = [T2,in + (T2,out )L ] 2

T2,V = T2,T P T2,L

(4.5.64) (4.5.65) (4.5.66)

Dit resulteert in: dT2,V 1 dT2,in d(T2,out )V = + dt 2 dt dt dT2,T P 1 dT2,in d(T2,out )T P = + dt 2 dt dt dT2,L 1 dT2,in d(T2,out )L = + dt 2 dt dt

(4.5.67) (4.5.68) (4.5.69)


43

Ten slotte wordt de uitlaattemperatuur van het secundair flu¨ıdum op dezelfde manier als in de verdamper: T2,out = FV (T2,out )V + FT P (T2,out )T P + FL (T2,out )L

(4.5.70)

Dit volgt uit het energiebehoud bij menging en de aanname dat Cp,2 aan de uitlaat van iedere zone gelijk is. Koelmiddel Integratie van het massabehoud in het koelmiddel, respectievelijk in de V-, TP- en de L-zone, geeft: m ˙ r,in − m ˙ r,V −T P dFV = [ρr,V − ρg,r ] Vr dt ∂ρr,V dhg,r ∂ρr,V dpr 1/2 + FV + dt ∂hr,V dpr ∂pr ∂ρr,V dhr,in + 1/2 FV dt ∂hr,V m ˙ r,V −T P − m ˙ r,T P −L dFL dFT P [(¯ γ − 1)(ρg,r − ρf,r )] + [ρf,r − ρg,r ] = Vr dt dt dρf,r dρg,r dpr d¯ γ + γ¯ + (1 − γ¯ ) + (ρg,r − ρf,r ) FT P dt dpr dpr dpr m ˙ r,T P −L − m ˙ r,out dFL = [ρr,L − ρf,r ] Vr dt ∂ρr,L dhf,r ∂ρr,L dpr + 1/2 + FL dt ∂hr,L dpr ∂pr ∂ρr,L dhr,out 1/2 FL + dt ∂hr,L

(4.5.71)

(4.5.72)

(4.5.73)

Optellen van deze 3 vergelijkingen resulteert in: m ˙ r,in − m ˙ r,out dFV dFT P dFL = [ρr,V − ρg,r ] + [(¯ γ − 1)(ρg,r − ρf,r )] + [ρr,L − ρg,r ] Vr dt dt dt ∂ρr,V dhg,r ∂ρr,V dpr + FV 1/2 + dt ∂hr,V dpr ∂pr dρf,r dρg,r d¯ γ +FT P γ¯ + (1 − γ¯ ) + (ρg,r − ρf,r ) dpr dpr dpr (4.5.74) ∂ρr,L dhf,r ∂ρr,L +FL 1/2 + ∂hr,L dpr ∂pr ∂ρr,V dhr,in + 1/2 FV dt ∂hr,V ∂ρr,L dhr,out + 1/2 FL dt ∂hr,L

Integratie van het energiebehoud van het koelmiddel, respectievelijk in de V-, TP- en de L-zone, geeft:


44

Qrt,V m ˙ r,in hr,in − m ˙ r,V −T P hg,r dFV − = [(ρh)r,V − (ρh)g,r ] (4.5.75) Vr Vr dt ∂(ρh)r,V dhg,r ∂(ρh)r,V dpr 1/2 + − 1 FV + dt ∂hr,V dpr ∂pr ∂(ρh)r,V dhr,in + 1/2 FV dt ∂hr,V m ˙ r,V −T P hg,r − m ˙ r,T P −L hf,r Qrt,T P dFT P − = [(¯ γ − 1)((ρh)g,r − (ρh)f,r )] (4.5.76) Vr Vr dt dFL + [(ρh)f,r − (ρh)g,r ] dt d(ρh)f,r d(ρh)g,r dpr + γ¯ + (1 − γ¯ ) dt dpr dpr d¯ γ ((ρh)g,r − (ρh)f,r ) − 1 FT P + dpr m ˙ r,T P −L hf,r − m ˙ r,out hr,out Qrt,L dFL − = [(ρh)r,L − (ρh)f,r ] (4.5.77) Vr Vr dt ∂(ρh)r,L ∂(ρh)r,L dhf,r dpr + − 1 FL 1/2 + dt ∂hr,L dpr ∂pr ∂(ρh)r,L dhr,out 1/2 FL + dt ∂hr,L De warmte-overdracht naar het koelmiddel vanuit de buiswand door convectie is aangeduid met −Qrt . Vr is het toale volume koelmiddel in de warmtewisselaar. Zoals in de verdamper, wordt verondersteld dat Fzone een goede benadering is voor de fractie van Vr die het koelmiddel in deze zone inneemt. ρr,V , ρr,T P , (ρh)r,V en (ρh)r,T P worden op dezelfde manier bepaald als in de verdamper. De zone-gemiddelde void-fractie is hier echter niet afhankelijk van de inlaatenthalpie, zodat: d¯ γ d¯ γ dpr = dt dpr dt

h hr,in hg,r

V

TP

(4.5.78)

L

hf,r hr,out 0

V Vv

Vv+VTP Vv+VTP+VL

Figuur 4.12: Lineair verloop enthalpie in de condensor


45

Opnieuw wordt het enthalpieverloop lineair verondersteld in iedere zone (figuur 4.12): 1 [hr,in + hg,r ] (4.5.79) 2 1 hr,T P = [hg,r + hf,r ] (4.5.80) 2 1 hr,L = [hf,r + hr,out ] (4.5.81) 2 Dit impliceert een lineair verloop van de dampfractie, zodat de volume-gemiddelde void-fractie op dezelfde manier wordt berekend als in de verdamper. Ook de tijdsafgeleiden van ρr,V , ρr,L , (ρh)r,V en (ρh)r,L nemen dezelfde vorm aan. hr,V =

De modelvergelijkingen zijn afhankelijk van de mode die optreedt. Bovenstaande vergelijkingen zijn afgeleid in de veronderstelling dat de condensor zich steeds in de V-TP-L-mode bevindt. Dit is typisch de meest voorkomende werkingsmode van een condensor. Door grote verstoring in de randvoorwaarden, of gedurende een overgangsverschijnsel, kan de L-zone echter verdwijnen. De condensor werkt dan in de V-TP-mode. In dit geval verloopt de afleiding volledig analoog. Het resultaat hiervan is terug te vinden in bijlage B. Verder moet er nagegaan worden wanneer een zone verdwijnt of bijkomt. Dit gebeurt door de uitlaatenthalpie van de condensor (hr,out ) te controleren. Indien de condensor in de V-TP-L-mode werkt, moet hr,out kleiner zijn dan de vloeistof saturatie-enthalpie (hf,r ) bepaald op condensordruk pr . Indien dit niet het geval is, moet geschakeld worden naar de V-TP-mode. In deze mode moet hr,out altijd groter zijn dan hf,r . Dit schakelcriterium is schematisch weergegeven in figuur 4.13

V-TP

V-TP-L

Figuur 4.13: Het schakelcriterium in de condensor Wanneer de L-zone wegvalt, verdwijnen er 3 differentiaalvergelijkingen uit het stelsel. Net zoals in de verdamper, worden deze vervangen door een set pseudo-vergelijkingen om de vorm van het stelsel te behouden: dFL kF = − FL (4.5.82) dt dTt,L kt = Tt,T P − Tt,L (4.5.83) dt d(T2,out )L k2 = (T2,out )T P − (T2,out )L (4.5.84) dt 4.5.2.5

Correlaties voor void-fractie en warmte-overdracht

De lokale void-fractie wordt bepaald met volgende experimentele vergelijking: γ= 1 + 0.79

1 0.58 1−x 0.78 ρg,r x

(4.5.85)

ρf,r

Dit is de void-fractie correlatie van Smith [12]. De warmte-overdrachtstermen zijn: Q2t,zone = Fzone πφi,t Nt Lα2t,zone [T2,zone − Tt,zone ]

(4.5.86)

Qrt,zone = Fzone πφo,t Nt Lαrt,zone [Tr,zone − Tt,zone ]

(4.5.87)


46

φi,t en φo,t zijn respectievelijk de binnen- en buitendiameter van de pijpen. De lengte van de pijpen wordt aangegeven met L. Nt is het aantal pijpen in de pijpenbundel. Gezien een meer-passige warmtewisselaar hier wordt vereenvoudigd tot een enkel-passige is deze waarde dus gelijk aan het aantal passen vermenigvuldigd met het aantal pijpen per pas. Tr,zone is de zone-gemiddelde temperatuur van het koelmiddel in de desbetreffende zone. Deze wordt geëvalueerd op de zone-gemiddelde druk en enthalpie. De convectiecoëfficiënten worden berekend met behulp van experimentele correlaties: 1. tussen koelmiddel en buiswand (αrt ) • zone L en V P rb 0.25 voor Reb = 1 − 102 P rt 0.25 0.5 0.36 P rb N urt = 0.52cn Reb P rb voor Reb = 102 − 103 P rt 0.25 0.63 0.36 P rb voor Reb = 103 − 2 · 105 N urt = 0.27cn Reb P rb P rt 0.25 0.8 0.4 P rb N urt = 0.033cn Reb P rb voor Reb = 2 · 105 − 2 · 106 P rt 0.36 N urt = 0.9cn Re0.4 b P rb

(4.5.88)

Deze correlatie is opgesteld door Zukauskas [12] en is geldig voor dwarsstroom over gladde buizen in bundels met parallelle schikking. De coëfficiënt cn is nodig bij korte bundels. Voor het aantal rijen n > 16 is het effect van het aantal buizen echter verwaarloosbaar. Het subscript b verwijst naar de hoofdstroming, terwijl het subscript t de buiswand aanduidt. Het Reynoldsgetal wordt bepaald met de uitwendige buisdiameter en de maximale snelheid die optreedt in de bundel. • zone TP (condensatie) " #0.25 kf3 (hg,r − hf,r )ρf,r (ρf,r − ρg,r )g αrt = C fgeom (4.5.89) µf,r (Tr,sat − Tt ) Deze correlatie is opgesteld door Beatty en Katz [3] en is geldig voor condensatie op enkele, horizontale, laag gevinde en gladde buizen. Voor gladde buizen is C = 0.725 en fgeom = 1. De correctiefactoren voor andere types buizen zijn terug te vinden in de literatuur [37]. • zone TP (koken) #0.25 " kg3 (hg,r − hf,r )ρg,r (ρf,r − ρg,r )g αrt = C fgeom (4.5.90) µg,r (Tt − Tr,sat ) Deze correlatie is opgesteld door Bromley [7]. Ze is geldig voor pool-boiling op gladde oppervlakken. fgeom en C zijn 2 constanten die afhankelijk zijn van de precieze geometrie van het kookoppervlak. Deze waarden zijn terug te vinden in de literatuur. 2. tussen secundair flu¨ıdum en buiswand (α2t ) • zone L, TP en V N u2t =

1+

f 2 (Red − 1000)P r 12.7( f2 )1/2 (P r2/3 −

"

2/3 # d 1+ L 1)

f = (1.58 ln(Red ) − 3.28)−2

(4.5.91) (4.5.92)


47

Deze correlatie is opgesteld door Gnielinksi [12]. Ze is geldig voor turbulente, gedwongen convectie in buizen, waarbij Red > 2300 en 0.5 < P r < 2000. d is de binnendiameter van de buis.

Hoofdstuk 5. Steady-state model

48

Hoofdstuk 5

Steady-state model 5.1

Inleiding

In dit hoofdstuk wordt het steady-state gedrag van de koelcyclus uit [11] gesimuleerd en vergeleken met de meetdata. Het model is ge¨ımplementeerd in MATLAB R2010a. De stofeigenschappen worden geëvalueerd met behulp van de XProps plug-in voor MATLAB. Het hoofddoel van deze thesis is het opstellen van een dynamisch model. Toch wordt er in dit hoofdstuk een statisch model ontwikkeld door alle tijdsafgeleiden in het dynamisch model gelijk te stellen aan nul. Hiermee kan het stationair deellastgedrag van de cyclus bestudeerd worden. Er zijn echter betere steady-state modellen beschikbaar in de literatuur [9]. Daarnaast kan het stationair gedrag van de cyclus ook voorspeld worden met het dynamisch model, eenvoudigweg door de randvoorwaarden constant te houden en de simulatie lang genoeg te laten duren. Er zijn echter andere belangrijke voordelen verbonden aan een dergelijk statisch model. Het kan immers ook gebruikt worden om: 1. de modellen van compressor en expansieklep te testen; 2. beginvoorwaarden te bepalen voor de dynamische simulatie. Dit laatste puntje verdient wat extra uitleg. Naast beginvoorwaarden van de inen uitlaattoestand van koelmiddel en secundair flu¨ıdum moeten er ook beginvoorwaarden opgegeven worden voor de interne variabelen van beide warmtewisselaars. Deze interne variabelen omvatten enerzijds de buiswandtemperaturen en anderzijds de zone-groottes. Meestal zijn de initiële toestandsgrootheden gekend, maar de interne variabelen niet. De initiële temperatuur van de buiswand heeft enkel invloed op het overgangsverschijnsel. De initiële zone-groottes daarentegen bepalen de lading koelmiddel (charge) die in het systeem zit en be¨ınvloeden bijgevolg het statisch gedrag van de cyclus. Het is dus van belang dat de initiële waarden van de zone-groottes consistent zijn met de charge. Indien het dynamisch model in staat is het opstartverschijnsel te beschrijven, vormt dit in principe geen probleem. Bij stilstand staan condensor en verdamper op gelijke druk. In de condensor bevindt zich oververhitte damp op koelwatertemperatuur. Verder is het koelmiddel in de verdamper een 2-fasig mengsel in thermisch evenwicht met het ijswater. Zo kunnen de initiële waarden van de zone-groottes en buiswandtemperaturen eenvoudig bepaald worden. Het model opgesteld in vorig hoofdstuk is echter niet geschikt om het opstartverschijnsel te beschrijven. Er moet dus een manier gevonden worden om consistente beginvoorwaarden te bepalen. Dit kan met behulp van een steady-state model dat bekomen is door de tijdsafgeleiden in het dynamisch model gelijk te stellen aan nul. In tegenstelling tot de dynamische simulatie, fungeert de charge hier als input. De resulterende zone-groottes zijn dus altijd consistent. De oplossing van de steady-state simulatie wordt dan gebruikt als beginvoorwaarde voor de dynamische simulatie.


5.2

49

Modelvergelijkingen

Het model van verdamper en condensor is bekomen door de tijdsafgeleiden in de dynamische vergelijkingen uit hoofdstuk 4 gelijk te stellen aan nul. Het resultaat hiervan bevindt zich in bijlage C. Dit resulteert in 2 keer 10 algebra¨ısche vergelijkingen en 2 keer 13 onbekenden. Deze onbekenden zijn xe en xc , met: e e e e e e e e xe = [her,in her,out per (T2,out )L (T2,out )T P (T2,out )V Tt,L Tt,T ˙ er ] P Tt,V FL FT P FV m c

x =

[hcr,in

hcr,out

pcr

c (T2,out )L (T2,out )T P

c (T2,out )V

c Tt,L

c Tt,T P

c Tt,V

FLc

FTc P

FVc

m ˙ cr ]

(5.2.1) (5.2.2)

Dit stelsel wordt aangevuld met het model van compressor en klep. Deze zijn zonder wijziging overgenomen uit hoofdstuk 4. Het klepmodel bestaat uit het isenthalp verband tussen inen uitlaat en de uitdrukking van Bernoulli voor het massadebiet. De vergelijkingen in het compressormodel omvatten enerzijds de uitdrukking van Bourdouxhe et al. [5] voor het massadebiet en anderzijds het polytroop verband tussen inen uitlaat. Dit resulteert in 2 keer 2 extra vergelijkingen. Er komen geen extra onbekenden in het stelsel voor, deze vergelijkingen brengen immers de inen uitlaattoestand van beide warmtewisselaars in verband. In steadystate vindt er geen massa-ophoping plaats in het systeem zodat het massadebiet doorheen de verdamper (m ˙ er ) en de condensor (m ˙ cr ) gelijk zijn. Om het stelsel te sluiten wordt deze set vergelijkingen aangevuld met de ladingsvergelijking. Deze drukt uit hoeveel lading koelmiddel er in het systeem vervat zit. Uiteindelijk wordt een niet-lineair algebra¨ısch stelsel bekomen van de vorm:

Dit moet opgelost worden naar x = xe

5.3

F(x) = 0 0 xc .

(5.2.3)

Oplossingsmethodiek

Het stelsel (5.2.3) wordt opgelost aan de hand van een sequentieel oplossingsalgoritme. Dit is schematisch weergegeven in figuur 5.1. In de volgende 2 secties wordt dit algoritme meer in detail besproken. Gelet op de doelstellingen van het steady-state model, biedt deze aanpak het voordeel dat de oorzaak van eventuele divergentie relatief eenvoudig kan opgespoord worden.

5.3.1

Systeemmodel

Het systeemmodel omvat alle componentmodellen en wordt iteratief opgelost. Deze iteratielus wordt de outer loop genoemd en is weergegeven in figuur 5.1. De inputgegevens van het systeemmodel zijn: 1. de randvoorwaarden van condensor en verdamper (het massadebiet, de inlaattemperatuur en de druk van het secundair flu¨ıdum); 2. de instelwaarde van de deellastverhouding (PLR); 3. de constructieve data van de componenten; 4. de totale lading koelmiddel in het systeem; 5. het controle-algoritme dat de inlaatschoepen aanstuurt.


50

START

Verdamper: Bereken · Interne variabelen · Uitlaatenthalpie · Uitlaattemperatuur ijswater

INITIALISEER inlaat verdamper · druk · enthalpie · massadebiet

UPDATE inlaat verdamper · druk · enthalpie · massadebiet

nee

EINDE

ja

Is oplossing voldoende nauwkeurig?

Compressor: Bereken · uitlaatenthalpie · condensordruk

Condensor: Bereken · Interne variabelen · Uitlaatenthalpie · Uitlaattemperatuur koelwater

Randvoorwaarden verdamper

Instelwaarde koellast

Randvoorwaarden condensor

Klep: Bereken · Massadebiet · uitlaatenthalpie

Lading: Bereken · verdamperdruk

Totale lading koelmiddel

Figuur 5.1: Schematische weergave van het oplossingsalgoritme van het systeemmodel Daarnaast moet een set initiële waarden opgegeven worden voor x. Deze zijn in principe volledig willekeurig, maar moeten toch voldoende dicht bij de oplossing liggen om het algoritme te doen convergeren. In de outer loop wordt elke component 1 voor 1 afgelopen. Eenmaal de uitlaattoestand van een component gekend is, kan de uitlaattoestand van de volgende component berekend worden. Zo wordt de volledige cyclus doorlopen. Dit is een iteratief algoritme, waarbij 1 iteratiestap overeenkomt met het doorlopen van 1 volledige cyclus. Een iteratiestap start met een set initiële waarden ψi die worden ge¨ updatet naar ψi+1 na het oplossen van iedere deelcomponent. In de eerste iteratiestap van het algoritme komen deze initiële waarden overeen met de opgegeven initiële waarden. Eenmaal 1 iteratie voltooid is, zijn dit de waarden berekend in de vorige stap. Telkens de volledige cyclus is doorlopen, wordt de massa-, energie en drukbalans gecontroleerd. Op basis hiervan wordt beslist of de oplossing voldoende nauwkeurig is of niet. Het algoritme blijft zich herhalen tot het verschil tussen de update van de inlaattoestand van de verdamper en de waarden berekend in de vorige iteratiestap voldoende klein is.


51

De oplossing van iedere deelcomponent wordt echter niet volledig doorgegeven naar systeemniveau. Er wordt relaxatie toegepast. Dit betekent dat: ψi+1 = (1 − β)ψi + βψnew

(5.3.1)

ψnew is de oplossing van een deelcomponent, bepaald op componentniveau. Verder is β de relaxatieparameter. Dit is een positieve waarde, kleiner dan 1. Relaxatie heeft als gevolg dat de oscillaties in de oplossing afgezwakt worden. Dit is nodig zodat de toestandsgrootheden in de outer loop binnen de werkingsgrenzen van XProps blijven.

5.3.2

Componentmodellen

Ook de componentmodellen worden iteratief opgelost. Deze iteratielus wordt een inner loop genoemd en is weergegeven in figuur 5.2.

START

INITIALISEER inlaatvoorwaarden

Component: Bereken uitlaatvoorwaarden

Is oplossing voldoende nauwkeurig?

UPDATE inlaatvoorwaarden

ja

EINDE

nee

Figuur 5.2: Schematische weergave van het oplossingsalgoritme van een componentmodel • De inputgegevens van het verdamperen condensormodel zijn: 1. het massadebiet, de inlaattemperatuur en de druk van het secundair flu¨ıdum; 2. het massadebiet, de inlaatenthalpie en de druk van het koelmiddel; 3. de constructieve data. Daarnaast moet een set initiële waarden opgegeven worden voor de uitlaattoestand, de interne variabelen en het massadebiet doorheen de warmtewisselaar. Er is verondersteld dat alle 6 mogelijke modes kunnen optreden in zowel verdamper als condensor. In werkelijkheid is dit weliswaar niet het geval. Deze modes zijn echter ingevoerd om oscillaties in de toestandsgrootheden van de outer loop op te vangen op componentniveau. Het invoeren van extra modes vereist echter dat ook het schakelcriterium wordt uitgebreid. In de eerste iteratiestap van de outer loop wordt de mode meegegeven met de initiële waarden. De voorgestelde mode in iedere volgende stap, is de berekende mode uit de vorige stap. Eerst wordt de inlaatenthalpie gecontroleerd. Op basis daarvan wordt beslist of de mode behouden wordt of de inlaatzone wegvalt. Daarna wordt de component uitgerekend. Zoals in het dynamisch model wordt vervolgens de oplossing gecontroleerd en indien nodig wordt geschakeld naar een andere mode op basis van de uitlaatenthalpie. Tijdens convergentie kunnen bepaalde zone-groottes


52

negatief worden. Wanneer dit voorvalt, wordt het algoritme afgebroken en wordt de warmtewisselaar in een andere mode uitgerekend. Dit is noodzakelijk omdat anders de warmte-overdracht omkeert en bijgevolg de wandtemperaturen onfysisch groot worden. Het feit dat een zone-grootte negatief wordt, is dus eigenlijk al een vroege waarschuwing dat de huidige mode incorrect is. De beslissingslogica is weergegeven in figuur 5.3.

L – TP – V

TP – V

L – TP FTP of FL < 0

TP V

L

Figuur 5.3: Schakelcriterium voor verdamper in steady-state model • De inputgegevens van het compressormodel zijn: 1. het massadebiet, de inlaatenthalpie en de inlaatdruk van het koelmiddel; 2. het controle-algoritme dat de inlaatschoepen aanstuurt; 3. de instelwaarde van de koellast. Deze laatste 2 bepalen immers de stand van de inlaatschoepen. Ook hier moet een set initiële waarden opgegeven worden voor de uitlaattoestand. • De inputgegevens van het klepmodel zijn: 1. het massadebiet, de inlaatenthalpie en de inlaatdruk van het koelmiddel; 2. de temperatuur van het koelmiddel aan de uitlaat van de verdamper. De klepheffing wordt immers bepaald aan de hand van de reservoirdruk. In steady-state is deze druk gelijk aan de saturatiedruk bepaald op de temperatuur van het koelmiddel aan de uitlaat van de verdamper. Net als in de compressor, vereist het klepmodel een set initiële waarden voor de uitlaattoestand. • De inputgegevens van het ’ladingmodel’ zijn: 1. de inen uitlaattoestand van het koelmiddel in verdamper en condensor; 2. de zone-groottes in verdamper en condensor. Beiden zijn nodig om de totale lading in het systeem te kunnen berekenen. Hier zijn geen initiële waarden vereist aangezien de ladingsvergelijking op een andere manier wordt opgelost dan de andere modellen.


53

Een componentmodel wordt iteratief opgelost met behulp van het Newton-Raphson algoritme, zoals beschreven in bijlage D. In elke iteratiestap wordt de Jacobiaan berekend aan de hand van een eindige-differenzen benadering. Hiermee wordt de oplossing uit de vorige stap φi ge¨ updatet naar φi+1 . Indien het verschil tussen beide waarden voldoende klein is, wordt het algoritme afgebroken en is de component opgelost. Oorspronkelijk was er geopteerd om ook de ladingsvergelijking op te lossen naar de verdamperdruk aan de hand van het Newton-Raphson algoritme. Wanneer de charge slechts enkele procenten afwijkt van de werkelijke charge, resulteert dit echter in een grote toename van de verdamperdruk. Het grootste deel van de koelmiddelmassa bevindt zich immers in de condensor, aangezien er daar een onderkoelde vloeistofzone aanwezig is. Het koelmiddel in de verdamper daarentegen bevindt zich voornamelijk in 2-fasige toestand. Dit betekent dus dat de verdamperdruk zeer groot zou moeten worden om een klein verschil in lading weg te werken. Een kleine verandering van verdamperdruk resulteert echter al in vrij grote veranderingen in de rest van de cyclus. Bijgevolg induceert dit sterke oscillaties op systeemniveau. Vandaar is er geopteerd om de verdamperdruk te wijzigen met een percentage dat afhankelijk is van het teken en de grootte van het verschil tussen berekende charge en werkelijke charge. Net zoals in de outer loop, wordt in de inner loops relaxatie toegepast. Dit betekent dat φi+1 een gewogen gemiddelde is van φi en de update die in de huidige stap is berekend. Elke component heeft zijn eigen relaxatieparameter. Relaxatie wordt toegepast omwille van 2 redenen. Opnieuw is het noodzakelijk dat de toestandsgrootheden tijdens het itereren in de inner loop binnen de werkingsgrenzen van XProps blijven. Daarnaast mag het schakelcriterium niet valselijk in werking treden wanneer een zone-grootte negatief wordt door oscillaties in de oplossing.

5.4

Resultaten

Er zijn 27 steady-state simulaties uitgevoerd, waarvan de randvoorwaarden elk corresponderen met 1 van de testgevallen die opgelijst zijn in tabel 2.3. De simulatieresultaten kunnen zowel kwantitatief als kwalitatief vergeleken worden met de metingen. Beide aspecten komen in deze sectie aan bod. Eerst wordt de nauwkeurigheid van de simulatie besproken, daarna het gesimuleerde deellastgedrag van de koelmachine. De sectie sluit af met enkele opmerkingen over de relaxatieparameters, de lading in het systeem en de gebruikte correlaties voor warmte-overdracht en void-fractie.


per

5

4

54

x 10

Te2,out 300

290

3.6 simulatie [K]

simulatie [Pa]

3.8

3.4 3.2

280

3

270

2.8 2.6

3

3.2

3.6

260 276

3.8

x 10

4.3

280 meting [K]

282

284

her,out

5

x 10

4.2 simulatie [J/kg]


278

5

x 10

her,in

5

2.6

3.4 meting [Pa]

2.4 2.3 2.2

4.1 4 3.9

2.1 2.2

2.25

2.3 2.35 2.4 meting [J/kg]

2.45

3.8

2.5

4

4.01

5

x 10

simulatie

+ 5%

meting

4.02 4.03 4.04 meting [J/kg]

4.05

4.06 5

x 10

− 5%

Figuur 5.4: Simulatieresultaten inen uitlaattoestand verdamper met ± 5 % foutgrenzen pcr

5

x 10

Tc2,out 330

10

320

9

310

simulatie [K]

simulatie [Pa]

11

8 7 6 5

5

6

7 8 meting [Pa]

9

270 290

10 x 10

hcr,in

x 10

300 meting [K]

305

310

hcr,out

5

2.6

x 10


4.4

4.2

4 4.25

295

5


290 280

5

4.8

300

2.4 2.3 2.2

4.3

4.35 4.4 meting [J/kg]

4.45

2.1 2.2

4.5

2.25

5

x 10 simulatie

+ 5%

meting

2.3 2.35 2.4 meting [J/kg]

2.45

2.5 5

x 10

− 5%

Figuur 5.5: Simulatieresultaten inen uitlaattoestand condensor met ± 5 % foutgrenzen


55

mr

COP

2.5

5.5 5 4.5 simulatie [−]

simulatie [kg/s]

2 1.5 1

4 3.5 3

0.5 2.5 0 0.5

1

1.5 meting [kg/s]

2

2

2.5

2

2.5

3

3.5 meting [−]

Pmotor

4

4.5

5

Qevap

120

400 350 simulatie [kW]

simulatie [kW]

100 80 60

300 250 200 150

40 100 20 30

40

50

60 70 meting [kW]

80

simulatie

50

90

0

100

+ 10%

meting

200 meting [kW]

300

400

− 10%

Figuur 5.6: Simulatieresultaten systeemparameters met ± 10 % foutgrenzen mr

COP

2.5

5 4.5

2

COP [−]

mr [kg/s]

4 1.5 1

3.5 3

0.5 0 20

2.5

40

60

80 100 PLR [%]

120

2 20

140

40

60

Pmotor

80 100 PLR [%]

120

140

120

140

Qevap

120

400 350

100 Qevap [kW]

Pmotor [kW]

300 80 60

250 200 150

40 20 20

100 40

60

80 100 PLR [%]

120

140 meting

50 20

40

60

80 100 PLR [%]

simulatie

Figuur 5.7: Simulatieresultaten systeemparameters uitgezet ten opzichte van de PLR


5.4.1

56

Nauwkeurigheid simulatieresultaten

In- en uitlaattoestand verdamper In figuur 5.4 worden de simulatieresultaten van de inen uitlaattoestand van de verdamper met de metingen vergeleken. In 82 % van de simulaties wordt de verdamperdruk (per ) voorspeld binnen de grenzen ± 10 % van de opgemeten waarden. De uitlaattemperatuur van e het ijswater (T2,out ) en de inen uitlaatenthalpie van het koelmiddel (her,in en her,out ) worden in 98 % van de gevallen voorspeld met een nauwkeurigheid van ± 1 %. In- en uitlaattoestand condensor In figuur 5.5 worden de simulatieresultaten van de inen uitlaattoestand van de condensor met de metingen vergeleken. In 100 % van de simulaties wordt de condensordruk (pcr ) voorspeld binnen de grenzen ± 10 % van de opgemeten waarden. De uitlaattemperatuur van het c koelwater (T2,out ) en de inen uitlaatenthalpie van het koelmiddel (hcr,in en hcr,out ) worden in 98 % van de gevallen voorspeld met een nauwkeurigheid van ± 1 %. Systeemparameters In figuur 5.6 worden de simulatieresultaten van de systeemparameters met de metingen vergeleken. In 82 % van de simulaties wordt het massadebiet koelmiddel (mr ) voorspeld binnen de grenzen ± 15 %. Het motor- en koelvermogen (Pmotor en Qevap ) worden respectievelijk in 71 % en 82 % van de gevallen voorspeld met een nauwkeurigheid van 15 %. Ten slotte wordt slechts in 66 % van de simulaties de COP voorspeld binnen ± 15 %. Het valt op dat de COP sterk onderschat wordt bij hogere PLR. Dit gedrag kan verklaard worden door het verloop van het compressorvermogen en het koelvermogen te bestuderen. De COP is immers berekend als: Qevap COP = (5.4.1) Pmotor met: Qevap = m ˙ r her,out − her,in (5.4.2) 1 Pmotor = m ˙ r hcr,in − her,out (5.4.3) ηem Pmotor is het totaal vermogen geleverd door de motor die de compressor aandrijft. De verliezen die optreden in de motor en transmissie worden uitgedrukt aan de hand van het elektromechanisch rendement ηem . Dit rendement wordt constant verondersteld, gelijk aan 62 %. Bij lage PLR wordt het motorvermogen goed voorspeld, bij hogere PLR wordt Pmotor echter overschat. Het verschil neemt toe naarmate de PLR stijgt. Deze discrepantie valt te verklaren aan de hand van het elektromechanisch rendement. Tijdens de simulatie is deze constant verondersteld, gelijk aan 62 %. Uit de meetwaarden uitgezet in 5.8 blijkt dat dit enkel een goede waarde is bij lage PLR. Bij hogere PLR daarentegen neemt ηem toe tot 72 %. Bijgevolg wordt het motorvermogen bij lage PLR goed voorspeld, terwijl er bij hogere PLR een overschatting optreedt, die toeneemt naarmate de PLR stijgt. Het verband tussen PLR en ηem kan ook fysisch verklaard worden. De aandrijving van de compressor is zo gedimensioneerd dat ze zich in haar nominale werkingstoestand bevindt wanneer de cyclus in vollast draait. Dit betekent dat, als het systeem in deellast werkt, dus bij lagere PLR, ook de aandrijving van de compressor in deellast werkt. Het is gekend dat in een elektromotor de elektrische verliezen kleiner zijn in vollast dan in deellast. Bijgevolg is ηem bij hoge PLR groter dan bij lage PLR.


57

72

70

68

ηem [%]

66

64

62

60

58 datapunten a*x^2 + b*x + c 56 20

40

60

80 PLR [%]

100

120

140

Figuur 5.8: Het verloop van het elektromechanisch rendement in functie van de PLR 5

11

5

x 10

4 3.8

10

3.6 per [Pa]

pcr [Pa]

9 8 7

3.4 3.2 3

6 5 20

2.8 40

60

80 100 PLR [%]

120

140

2.6 20

meting 1

x 10

40

60

80 100 PLR [%]

120

140

simulatie 1

Figuur 5.9: Simulatieresultaten verdamperen condensordruk in functie van de PLR 0.8

0.8

0.6

0.6

Ook het koelvermogen Qevap wordt goed voorspeld bij lage PLR. Bij hogere PLR is het gesimuleerde koelvermogen echter kleiner dan het opgemeten vermogen. Aangezien de enthalpie0.4 0.4 toename in de verdamper voor alle deellastverhoudingen goed wordt voorspeld, is deze trend 0.2 0.2 te verklaren aan de hand van het verloop van het massadebiet koelmiddel m ˙ r . Qevap is immers gelijk aan het product van beide grootheden. Er wordt geacht dat dit gedrag te wijten is aan 0 0 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 een onderschatting van de warmte-overdrachtscoëfficiënt α in de verdamper bij 1hogere PLR. Immers, wanneer α daalt, moet het temperatuursverschil tussen secundair flu¨ıdum en het koelmiddel vergroten om dezelfde warmte-overdracht te realiseren. In de verdamper resulteert dit in een drukdaling, terwijl de condensordruk stijgt. Uit de resultaten van de condensordruk, uitgezet ten opzichte van de PLR in figuur 5.9, blijkt dat de warmte-overdrachtscoëfficiënt in de condensor altijd kleiner wordt voorspeld dan in werkelijkheid, onafhankelijk van de PLR. In deze figuur is ook de verdamperdruk uitgezet. Bij toenemende PLR wordt de verdamperdruk steeds sterker onderschat. Uit figuur 5.7 is af te leiden dat een hoge PLR overeenkomt met een hoog massadebiet. Dit betekent dus dat bij toenemend massadebiet, de voorspelling van de warmte-overdrachtscoëfficiënt te laag uitvalt. Naarmate het massadebiet toeneemt, worden convectieve effecten tijdens het koken steeds groter. Het kookverschijnsel gaat over van pool boiling naar flow boiling. Tijdens flow boiling kunnen 2 modes van warmte-overdracht beschouwd worden: bellenkoken (nucleate boiling) en convectief koken (convective boiling). Bellenkoken is onder deze omstandigheden gelijkaardig aan pool boiling. Enkel het effect van de stroming op de groei en het loslaten van de bel-


58

len komt erbij. Convectief koken verwijst naar de convectieve warmte-overdracht die optreedt tussen de verwarmde wand en de vloeistoffase [12]. De verschillende flow boiling modellen combineren deze 2 mechanismen door middel van een machtsfunctie: αtp = [(αnb )n + (αcb )n ]1/n

(5.4.4)

αnb is onafhankelijk van de massasnelheid, terwijl αcb meestal onafhankelijk is van de warmteflux. Thome [37] heeft aangetoond dat bij flow boiling over bundels gladde buizen, het aandeel van het convectief koken in de totale warmte-overdracht gemiddeld gelijk is aan 40 %. Deze conclusie is ook geldig indien de buizen laag gevind zijn. De omstandigheden waaronder deze studie is uitgevoerd komen overeen met de gesimuleerde omstandigheden, zoals weergegeven in tabel 5.1.

massaflux G warmteflux q dampfractie x

[37] 5 tot 41 2 tot 35 10 tot 87

Simulatie 5 tot 22 4 tot 14 10 tot 100

kg/(m2 s) kW/m2 %

Tabel 5.1: Overeenkomst met de testomstandigheden uit [37] De correlatie die tijdens de simulatie gebruikt is om de warmte-overdrachtscoëfficiënt in de verdamper te voorspellen (vergelijking (4.5.90)), is enkel geldig voor pool-boiling. Bijgevolg wordt bij hogere massadebieten, en dus bij hogere PLR, de warmte-overdrachtscoëfficiënt in de verdamper steeds sterker onderschat. Het drukverschil tussen verdamper en condensor is dus hoger dan in realiteit en de overschatting neemt steeds toe. Om dit drukverschil te realiseren, moet het massadebiet koelmiddel verlagen. Zo verkleint de drukval door de inlaatschoepen en neemt de opvoerhoogte van de compressor toe, zodat de totale druktoename over de compressor stijgt.

5.4.2 5.4.2.1

Deellastgedrag Invloed PLR op toestandsgrootheden en massadebiet

In figuur 5.10 zijn 3 verschillende testgevallen weergegeven. Deze komen overeen met testcondities 1, 2 en 3 uit tabel 2.3. De inlaattemperatuur van het koelwater en de insteltemperatuur van het ijswater zijn in alle 3 gevallen gelijk. De PLR daalt echter omdat de inlaattemperatuur van het ijswater afneemt. De inlaattemperatuur van het ijswater neemt af, terwijl het instelpunt constant blijft. Gezien het massadebiet ijswater en koelwater niet veranderen, betekent dit dat de warmte die in condensor en verdamper moet gewisseld worden, verkleint. Bijgevolg daalt de druk in de condensor, terwijl de druk in de verdamper toeneemt. Zo daalt het temperatuursverschil en neemt bijgevolg de warmte-overdracht tussen het secundair flu¨ıdum en het koelmiddel af. Dit heeft als gevolg dat het massadebiet verkleint om de oververhitting aan de uitlaat van de verdamper te behouden. Immers, als de oververhitting te laag is, daalt de druk in het reservoir van de expansieklep, waardoor de klepheffing kleiner wordt en het debiet afneemt. Ook beginnen de inlaatschoepen van de compressor te sluiten als respons op de lagere PLR. Als gevolg van het lager massadebiet, neemt de polytrope exponent toe. Dit betekent dat het compressieverloop in het ph-diagram steiler verloopt. Een stijging van de polytrope exponent betekent ook dat de verliezen tijdens compressie toenemen. Uit het Ts-diagram in figuur 5.11 blijkt immers dat de entropietoename over de compressor groter is in deellast dan in vollast. Dit is ook zichtbaar in het dalende verloop van de COP.


59

7

10

hg,r hf,r

p [Pa]

testgeval 1 testgeval 2 testgeval 3

6

10

5

10 1.5

2

2.5

290

3.5

4

4.5 5

3.2

289

x 10

2

287

mr [kg/s]

COP [−]

3

288

2.8

1.5

1

2.6

286 285 1

2 testgeval

3

1

2 testgeval

3

0.5

1

2 testgeval

Figuur 5.10: Invloed PLR op toestandsgrootheden en massadebiet 380 sf,r sg,r

360


340

320 T [K]

Tewi [K]

3 h [J/kg]

300

280

260

240 1100

1200

1300

1400 1500 s [J/kg.K]

1600

1700

Figuur 5.11: Entropietoename tijdens compressie

1800

3

Hoofdstuk 5. Steady-state model 5.4.2.2

60

Invloed PLR op COP

De algemene trend, waarbij de COP stijgt bij toenemende PLR is weergegeven in figuur 5.7. De componenten in het systeem worden ontworpen om in nominale omstandigheden te werken. Bijgevolg is de COP het grootst wanneer de PLR ongeveer 100 % is. Naast PLR hebben ook de condensor- en verdampertemperatuur een invloed op de COP van een koelmachine. De COP van een ideale Carnot-cyclus wordt gegeven door [24]: COP =

TC TH − TC

(5.4.5)

Hierin is TC de temperatuur van het koud reservoir waarvan de warmte wordt onttrokken. Deze warmte wordt afgegeven aan een warm reservoir met temperatuur TH . Het warm en koud reservoir komen dus overeen met het secundair flu¨ıdum in respectievelijk de condensor en de verdamper. Wanneer de koelwatertemperatuur stijgt, neemt de COP af. De COP neemt echter toe bij toenemende ijswatertemperatuur. Dit gedrag is ook terug te vinden in de simulatieresulaten, zoals weergegeven in figuuur 5.12. Hierin is Tcwi de inlaattemperatuur van het koelwater, terwijl Tewi de inlaattemperatuur van de ijswater voorstelt. PLR ≈ 60 % Tewi ≈ 287 K

PLR ≈ 110 % Tewi ≈ 286.5 K

4

4

3.5

3.5

3.5

3

278

280

282 284 Tewi [K] testgeval 9 testgeval 18 testgeval 27

286

COP [−]

4

COP [−]

COP [−]

PLR ≈ 30 % Tcwi ≈ 290 K

3

290

295 300 Tcwi [K] testgeval 2 testgeval 5 testgeval 8

305

3

290

295 300 Tcwi [K]

305


Figuur 5.12: Invloed temperatuur koelwater en ijswater op COP

5.4.3

Opmerkingen

Lading koelmiddel en correcties op de warmte-overdrachtsco¨ effici¨ enten De totale lading in het systeem is 137 kg [11]. In de simulatie resulteert deze charge echter in een te hoge condensordruk en een sterke onderschatting van de COP. Bij lagere charge, neemt de COP meer realistische waarden aan, maar wordt de condensordruk te klein. Dit is weergegeven in figuur 5.13. De te hoge condensordruk kan gecompenseerd worden door de voorspelde convectiecoëfficiënten met een correctiefactor te vermenigvuldigen. Met een charge van 137 kg vallen de correcties echter onrealistisch groot uit (grootte-orde 10 à 20). Wanneer er geopteerd wordt voor lagere, meer realistische correcties (grootte-orde 1 à 2), blijkt uit figuur 5.14 dat de druk in de verdamper te laag uitvalt, terwijl de condensordruk te hoog is. Dit kan evenwel gecompenseerd worden door de lading op een lagere waarde in te stellen. Dit zou echter resulteren in een onrealistisch lage charge (grootte-orde 30 à 40 kg). Daarom is tijdens de simulaties de charge ingesteld op 85 kg. De correcties op de warmte-overdrachtscoëfficiënten zijn van grootte-orde 1 à 5, zoals weergegeven in tabel 5.2.


61

7

10

hg,r hf,r

p [Pa]

charge = 119 kg charge = 85 kg charge = 63.75 kg

6

10

5

2

2.5

3 h [J/kg]

3.5 2.1

3

2.05

2.9 2.8 2.7 2.6 60

4

4.5 5

3.1

mr [kg/s]

COP [−]

10 1.5

x 10

2 1.95 1.9

70

80

90 100 lading [kg]

110

120

1.85 60

70

80

90 100 lading [kg]

110

120

Figuur 5.13: Invloed lading koelmiddel op testgeval 1

verdamper condensor

koelmiddel secundair flu¨ıdum koelmiddel secundair flu¨ıdum

L 1 1 1.5 1.5

TP 1 5 1.5 5

V 1 1 1.5 1.5

Tabel 5.2: Correctiefactoren op de convectiecoëfficiënt in iedere zone Deze correctiefactoren zijn ook nodig om bij hogere massadebieten de invloed van flow boiling in de verdamper in rekening te brengen, zoals eerder toegelicht. De factoren zijn bepaald aan de hand van enkele simulaties bij middelgrote PLR. Dit zijn echter constante waarden, onafhankelijk van het massadebiet. Bijgevolg wordt de warmte-overdracht bij lage PLR overschat omdat dan zuivere pool boiling optreedt. Bij hogere PLR daarentegen wordt de warmte-overdracht onderschat omdat de invloed van het massadebiet op de warmte-overdracht hier onvoldoende in rekening gebracht wordt. Deze trend is bijgevolg ook terug te vinden in de verdamperdruk zoals te zien is in figuur 5.9. Bij lage PLR wordt de druk overschat, omdat de gecorrigeerde convectiecoëfficiënt te groot uitvalt. Bij hoge PLR daarentegen wordt de druk sterk onderschat door een te kleine correctie. Er moet echter voorzichtig omgesprongen worden met deze correcties. Het is van belang dat de convectiecoëfficiënt van de V-zone voldoende kleiner is dan die in de TP-zone. Neem bijvoorbeeld de condensor. Bij voldoende hoge druk werkt deze in de V-TP-L-mode. Stel dat de inlaatvoorwaarden constant blijven, terwijl de druk daalt, zoals in figuur 5.15. De TP- en V-zone vergroten, terwijl de L-zone verkleint. Op een bepaald moment is de druk voldoende


62

7

10

hg,r hf,r

p [Pa]

testgeval 1 (correcties = 1 à 5) testgeval 1 (correcties = 1 à 2)

6

10

5

10 1.5

2

2.5

3 h [J/kg]

3.5

4

4.5 5

x 10

Figuur 5.14: Invloed correctiefactoren op de convectiecoëfficiënt in testgeval 1 laag zodat de condensor overgaat van V-TP-L naar V-TP. Vanaf dan begint de V-zone te vergroten en de TP-zone te verkleinen. Als de convectiecoëfficiënt in de V-zone te groot is, groeit deze zone niet snel genoeg aangezien de damp dan veel te snel condenseert. De oplossing bij lagere druk wordt bijgevolg onfysisch. De uitlaatenthalpie is groter dan de damp saturatie-enthalpie, terwijl de TP-zone nog vrij groot is. Fysisch gezien, moet de grootte van de TP-zone naar nul gaan naarmate de uitlaatenthalpie de damp saturatie-enthalpie nadert. Echter, het enthalpieverloop van het koelmiddel doorheen de warmtewisselaar is stuksgewijs lineair verondersteld. Hieruit volgt ook een ongeveer lineair temperatuursverloop in de V-zone. In werkelijkheid is het temperatuursverloop echter exponentieel onder steady-state voorwaarden [12]. Dit betekent dat in de simulatie het temperatuursverschil tussen secundair flu¨ıdum en het koelmiddel in de V-zone overschat wordt. Bijgevolg wordt ook de warmte-overdracht in deze zone overschat en neemt de grootte van de V-zone niet snel genoeg toe. Het is dus van belang dat de correctiefactoren voor de warmte-overdracht in de V-zone voldoende klein zijn om te compenseren voor de overschatting van het temperatuursverschil. Ten slotte wordt de totale lading in het systeem onderschat omwille van 3 redenen. Ten eerste is de 2-passige warmtewisselaar vereenvoudigd tot een 1-passige met dubbel zoveel pijpen. Door deze vereenvoudiging verkleint het volume koelmiddel in de warmtewisselaar. Daarnaast zit er nog een zekere hoeveelheid lading in de leidingen van het systeem. Ten slotte is de void-fractie correlatie eigenlijk niet geschikt voor stroming over buizenbundels. Waarschijnlijk wordt de void-fractie overschat zodat de voorspelling van lading in beide warmtewisselaars te laag uitvalt. Relaxatieparameters Zoals uitgelegd in sectie 5.3, wordt er relaxatie toegepast op de outer loop. Dit is om te voorkomen dat de toestandsgrootheden buiten de werkingsgrenzen van XProps komen te liggen. Om de convergentie te laten slagen, of op zijn minst sneller te laten lopen, is er geen relaxatie toegepast op de zone-groottes. Stel bijvoorbeeld dat de L-zone in de condensor wegvalt ten gevolge van een te lage druk. Aangezien een groot deel van de lading in de L-zone van de condensor zit, valt de berekende charge te laag uit, bijgevolg doet het algoritme de verdamperdruk toenemen. De klep zorgt er echter voor dat de oververhitting aan de uitlaat van de verdamper constant blijft, zodat het massadebiet afneemt. Bijgevolg neemt de condensordruk weer toe en wordt de L-zone opnieuw ge¨ıntroduceerd. Door te sterke relaxatie van de zone-


5

druk

x 10

4.5 4

h [J/kg]

pr [Pa]

6.5

6

5.5

5

enthalpie

x 10

massadebiet 2

in

1.5

3.5

mr [kg/s]

5

7

63

out 3

0.5

2.5

5

10 test

2

15

1

5

10 test

0

15

5

10 test

15

1 0.8 0.6 0.4 0.2 0

2

4

6

8

10

12

14

16

18

test FL

FTP

FV

Figuur 5.15: Verloop zone-groottes in condensor bij dalende druk en constante inlaatenthalpie en massadebiet groottes blijft deze reactie echter te lang uit. De ladingsvergelijking ’ziet’ niet dat de L-zone is weggevallen en zal dus ook geen grote correcties aanbrengen aan de verdamperdruk. De druk in de condensor blijft dalen, tot de grootte van de L-zone in de outer loop voldoende klein is. Daarop wordt gereageerd door een sprong in de verdamperdruk. Deze geeft echter aanleiding tot een te grote druk in de condensor, waardoor die weer moet afnemen. Het risico is bijgevolg dat hetzelfde fenomeen opnieuw optreedt en de convergentie aldus slechts zeer traag vordert of niet bereikt wordt door de grote oscillaties tijdens het itereren.

5.5

Conclusie

In dit hoofdstuk is een model ontwikkeld dat in staat is het steady-state gedrag van de koelmachine beschreven in [11] te reproduceren. Het gesimuleerde deellastgedrag komt kwalitatief overeen met het opgemeten gedrag. Verder worden alle toestandsgrootheden, uitgezonderd verdamperen condensordruk goed voorspeld. De gesimuleerde COP vertoont dezelfde trend als de opgemeten COP. Het valt echter op dat er een steeds sterkere onderschatting optreedt naarmate de PLR stijgt. Er wordt aangetoond dat de afwijking op de COP en de verdamperen condensordruk te herleiden is naar de gekozen correlaties voor de warmte-overdracht. De correlatie die de warmte-overdracht in de verdamper voorspelt is ongeschikt voor flow boiling, zodat bij hogere PLR de warmte-overdracht onderschat wordt, wat resulteert in een lagere verdamperdruk en aldus in een lagere COP. Verder wordt het effect van het elektromechanisch rendement toegelicht. Het volstaat niet deze waarde als constant te beschouwen, aangezien dit resulteert in een sterke overschatting van het motorvermogen bij hoge PLR. Ten slotte


64

wordt aan de hand van de void-fractie correlatie verklaard waarom de charge in de simulatie kleiner is dan die in werkelijkheid. Hieraan is tegemoet gekomen door de warmte-overdracht te corrigeren.

Hoofdstuk 6. Dynamisch model

65

Hoofdstuk 6

Dynamisch model 6.1

Inleiding

In dit hoofdstuk wordt het dynamisch gedrag van de koelcylus uit [11] gesimuleerd en vergeleken met de meetdata. Het model is ge¨ımplementeerd in MATLAB R2010a. De stofeigenschappen worden geëvalueerd met behulp van de XProps plug-in voor MATLAB.

6.2

Modelvergelijkingen

Het model van verdamper en condensor is opgesteld in hoofdstuk 4. Het gedrag van de verdamper wordt uitgedrukt aan de hand van 7 differentiaalvergelijkingen, aangevuld met 1 algebra¨ısche vergelijking. Deze laatste drukt uit dat FT P + FV = 1. Wanneer de verdamper zich in de TP-mode bevindt, zijn de behoudswetten opgesteld in de V-zone niet meer geldig en worden vervangen door 3 pseudo-vergelijkingen. De onbekenden die in deze set vergelijkingen voorkomen, zijn de interne variabelen van de verdamper en de toestandsgrootheden en het massadebiet aan inen uitlaat: e e e e e e ye = [her,in her,out per (T2,out )T P (T2,out )V Tt,T ˙ er,in m ˙ er,out ] P Tt,V FT P FV m

(6.2.1)

Het condensormodel bestaat uit 10 differentiaalvergelijkingen. Ook hier wordt de algebra¨ısche vergelijking FL + FT P + FV = 1 aan het stelsel toegevoegd. In de V-TP-mode, zijn de betrekkingen afgeleid in de L-zone niet meer geldig. Zoals in de verdamper, worden deze vervangen door 3 pseudo-vergelijkingen. De onbekenden die in deze set vergelijkingen voorkomen, zijn de interne variabelen van de condensor en de toestandsgrootheden en het massadebiet aan inen uitlaat: c c c c c c c c yc = hcr,in hcr,out pcr (T2,out )L (T2,out )T P (T2,out )V Tt,L Tt,T ˙ cr,in m ˙ cr,out (6.2.2) P Tt,V FL FT P FV m

Het stelsel wordt verder aangevuld met het model van compressor en expansieklep. Beide modellen zijn statisch, de modelvergelijkingen zijn bijgevolg algebra¨ısch. Het compressormodel omvat 2 vergelijkingen: het massadebiet volgens Bourdouxhe et al. [5] en de polytrope betrekking tussen inen uitlaat. Het klepmodel bestaat ook uit 2 vergelijkingen: de Bernoulliklepvergelijking en het isenthalp verband tussen inen uitlaat. Er worden geen extra onbekenden in het stelsel ge¨ıntroduceerd. Deze vergelijkingen brengen immers inen uitlaattoestand van verdamper en condensor in verband. Het reservoir van de klep wordt echter dynamisch gemodelleerd. Zo wordt een extra onbekende aan het stelsel toegevoegd: yb = Tbulb

(6.2.3)


66

Uiteindelijk wordt een niet-lineair stelsel differentiaal-algebra¨ısche vergelijkingen bekomen van de vorm: dy M(t, y) (t) = F(t, y) (6.2.4) dt 0 Dit moet opgelost worden naar het tijdsverloop van y = ye yc yb . M is de massamatrix, deze is singulier. In steady-state reduceert dit stelsel zich tot: F(t, y) = 0

(6.2.5)

Dit is de steady-state vorm van het stelsel, gebruikt in hoofdstuk 5.

6.3

Oplossingsmethodiek

Het stelsel (6.2.4) wordt opgelost aan de hand van de ingebouwde MATLAB-solver ode15s. Deze solver wordt aangeraden voor het oplossen van differentiaal-algebra¨ısche problemen met index 1 en hoge stijfheid. Voor meer informatie over dit specifiek algoritme, wordt de lezer doorverwezen naar [30] en [21]. Er is gekozen voor een voorgeprogrammeerde solver aangezien deze al is geoptimaliseerd om dergelijke problemen op te lossen. Dit heeft echter als nadeel dat de oorzaak van eventuele divergentie veel moeilijker op te sporen is. De individuele modellen zijn echter al vooraf grondig getest. Zowel het compressor- als klepmodel zijn ge¨ımplementeerd in de steady-state simulatie. Daarnaast zijn ook het dynamische condensor- en verdampermodel als aparate entiteiten gecontroleerd.

6.4

Resultaten

Er is 1 dynamische simulatie uitgevoerd. De randvoorwaarden komen overeen met het tijdsverloop zoals afgebeeld op figuur 2.4. Het model is echter niet in staat het opstartverschijnsel te beschrijven. Daarom is er voor gekozen het starttijdstip van de simulatie te verleggen naar t = 2700 s. Dit tijdstip wordt verder aangeduid met t = 0 s. Op dit moment bevindt het systeem zich in stationaire toestand, zodat de beginvoorwaarden van de simulatie overeenkomen met de steady-state simulatieresultaten van testgeval 1. Hoewel de volledige testcyclus nog verder doorloopt, eindigt de simulatie na 47000 s. Het model is immers ook niet geschikt om het uitzetverschijnsel van de installatie te voorspellen. In wat volgt, wordt nagegaan of het transiënt gedrag goed gecapteerd wordt. Het steadystate gedrag is immers al behandeld in hoofdstuk 5. Ook wordt de uitvoeringssnelheid van de simulatie besproken. Dit is een belangrijk aspect aangezien er gekozen is voor een movingboundary model omwille van zijn hoge uitvoeringssnelheid. Kiezen voor de moving-boundary aanpak brengt wel extra moeilijkheden met zich mee. De modelvergelijkingen zijn immers afhankelijk van de mode die optreedt in de warmtewisselaar. Wanneer een zone verdwijnt, blijft de vorm van het stelsel behouden door de behoudsvergelijkingen van die zone te vervangen door pseudo-vergelijkingen. Dit garandeert dat het schakelen tussen modes continu verloopt. In het laatste deel van deze sectie wordt dit gedemonstreerd.


67

5

x 10

4

verdamperdruk [Pa]

3.8 3.6 3.4 3.2 3 2.8 2.6

0

0.5

1

1.5

2

2.5 tijd [s]

simulatie

5

3

3.5

4

4.5 4

x 10 meting

x 10

11

condensordruk [Pa]

10 9 8 7 6 5 4

0

0.5

1

1.5

2

2.5 tijd [s]

3

3.5

4

4.5

5 4

x 10

Figuur 6.1: Simulatieresultaten verdamperen condensordruk Condensor

5

x 10

hr,out [J/kg]

2.5 2.4 2.3 2.2

0 5 x 10

0.5

1

1.5

2

hr,in [J/kg]

4.6

2.5 tijd [s]

3

2.5 tijd [s]

3

4

4.5 4

x 10

4.4

4.2

0

0.5

1

1.5

2

meting 4.1

3.5

4

4.5 4

x 10 simulatie

Verdamper

5

hr,out [J/kg]

3.5

x 10

4.05 4 3.95

hr,in [J/kg]

2.5

0 5 x 10

0.5

0

0.5

1

1.5

2

2.5 tijd [s]

3

2.5 tijd [s]

3

3.5

4

4.5 4

x 10

2.4 2.3 2.2

1

1.5

2

3.5

4

4.5 4

x 10

Figuur 6.2: Simulatieresultaten inen uitlaatenthalpie van verdamper en condensor


68

uitlaattemperatuur ijswater [K]

286 284 282 280 278 276 0

0.5

1

1.5

2

2.5 tijd [s]

simulatie

3

3.5

4

4.5 4

x 10 meting

uitlaattemperatuur koelwater [K]

310

305

300

295

290 0

0.5

1

1.5

2

2.5 tijd [s]

3

3.5

4

4.5 4

x 10

Figuur 6.3: Simulatieresultaten uitlaattemperatuur koelen ijswater 2.2 2

massadebiet [kg/s]

1.8 1.6 1.4 1.2 1 0.8 0.6 0.4 0

0.5

1

1.5 meting

2

2.5 tijd [s]

3

3.5

4

4.5 4

x 10

simulatie (comp)

simulatie (tev)

140

PLRset [%]

120 100 80 60 40 20

0

0.5

1

1.5

2

2.5 tijd [s]

3

3.5

4

4.5 4

x 10

Figuur 6.4: Simulatieresultaten massadebiet compressor en expansieklep en tijdsverloop PLR


6.4.1

69

Uitvoeringssnelheid

De simulatie werd uitgevoerd op een PC met 4 GB RAM-geheugen en een dual-core processor van 2.4 GHz. De volledige run nam ongeveer 19500 s in beslag. Dit betekent dus dat de simulatie ongeveer 2.4 keer sneller verliep dan real-time.

6.4.2

Nauwkeurigheid simulatieresultaten

Toestandsgrootheden In figuren 6.1, 6.2 en 6.3 worden de simulatieresulaten van de toestandsgrootheden vergeleken met de metingen. Het tijdsverloop van de inen uitlaatenthalpie in zowel verdamper als condensor komt goed overeen met de meetresultaten. Zoals aangekaart in hoofdstuk 5, wordt de verdamperdruk niet goed voorspeld. Ook de nauwkeurigheid van de voorspelde verdamperdruk is laag. Toch wordt het transiënt gedrag van beide grootheden goed voorspeld. Figuur 6.5 geeft de respons van de condensordruk weer op een grote verstoring van de randvoorwaarden. Zelfs hier komen de opgemeten en gesimuleerde transiënt kwalitatief goed overeen. Beide overgangsverschijnselen hebben immers dezelfde vorm. Er treedt met andere woorden geen overof undershoot op en de voorspelde tijdsconstante komt overeen met de opgemeten waarde. 5

10

5

x 10

x 10

9 8.5 condensordruk [Pa]

8 7 6 5 1.4

simulatie meting 1.45

1.5

1.55 tijd [s]

1.6

1.65

8 7.5 7 6.5 simulatie meting

6 5.5

1.7

3

310

200

3.3

3.4 4

x 10

200 Tewi Tcwi

100 290

PLR [%]

150

PLR

temperatuur [K]

Tcwi temperatuur [K]

3.2 tijd [s]

310

Tewi

300

3.1

4

x 10

300

150

PLR 100

PLR [%]

condensordruk [Pa]

9

290 50

50 280 1.4

1.45

1.5

1.55 tijd [s]

1.6

1.65

1.7 4

x 10

280 3

3.1

3.2 tijd [s]

3.3

0 3.4 4

x 10

Figuur 6.5: Overeenkomst tussen opgemeten en gesimuleerd transiënt gedrag condensordruk Ten slotte wordt ook het verloop van de uitlaattemperatuur van het koelwater goed voorspeld. Het gedrag van de ijswatertemperatuur komt echter minder goed overeen met de werkelijkheid. Massadebiet Het transiënt verloop van het massadebiet vertoont een goede gelijkenis met de meetresultaten (zie figuur 6.4). Echter, bij sterke verandering van de randvoorwaarden komt het gesimuleerde gedrag niet meer overeen met het opgemeten gedrag. In figuur 6.6 is links de overgang tussen


70

testgeval 9 en 10 weergegeven en rechts de overgang tussen 18 en 19. Beide komen overeen met een daling van de insteltemperatuur van het uitgaand ijswater. Daarnaast stijgt in beide gevallen de inlaattemperatuur van zowel het ijswater als het koelwater. Dit komt dus overeen met een sterke stijging in de PLR. 2

massadebiet [kg/s]

massadebiet [kg/s]

2

1.5

1

0.5 1.4

1.45

1.5 tijd [s]

1.55

1

0.5

1.6

3

3.1

x 10

simulatie (comp) 120

100

100 PLRset [%]

120

80 60 40 20 1.4

3.05

4

meting

PLRset [%]

1.5

3.15 tijd [s]

3.2

3.25

3.3 4

x 10

simulatie (tev)

80 60 40

1.45

1.5 tijd [s]

1.55

1.6 4

x 10

20

3

3.05

3.1

3.15 tijd [s]

3.2

3.25

3.3 4

x 10

Figuur 6.6: Afwijking tussen opgemeten en gesimuleerd transiënt gedrag massadebiet Er is verondersteld dat de tijdsconstante van het reservoir van de thermostatische expansieklep gelijk is aan 5 s. Deze waarde is van dezelfde grootte-orde als de kleppen gemodelleerd door Mithraratne et al. [23]. Zoals weergegeven in figuur 6.7 resulteert een te kleine tijdsconstante in sterke oscillaties in het massadebiet doorheen de klep. De druk in het reservoir reageert dan immers zeer hevig op kleine veranderingen in de oververhitting aan de uitlaat van de verdamper. De tijdsconstante mag echter niet te groot zijn. Figuur 6.8 toont de simulatieresultaten bij de overgang van testgeval 3 naar 4 met 2 verschillende tijdsconstanten. Tijdens de overgang daalt de inlaattemperatuur van het koelwater, maar neemt die van het ijswater toe. De insteltemperatuur van het uitgaande ijswater blijft echter constant (zie tabel 2.3). Omdat de koellast groter wordt, terwijl de insteltemperatuur gelijk blijft, daalt de druk in de verdamper. Deze drukdaling heeft als gevolg dat het massadebiet doorheen de klep stijgt. Ook het massadebiet door de compressor neemt toe als gevolg van het openen van de inlaatschoepen door de toename in PLR. Door de stijging van het massadebiet, neemt de oververhitting van het koelmiddel echter af. Wanneer de uitlaatenthalpie sneller daalt dan de damp saturatie-enthalpie, zal deze op een bepaald moment lager uitvallen. Dit betekent dat de dampzone aan de uitlaat van de verdamper zal verdwijnen. De klep reageert evenwel op de dalende oververhitting door de klepheffing te verkleinen. Door deze corrigerende actie, zakt op een bepaald moment het inlaatdebiet onder het uitlaatdebiet. Dit betekent dus dat de verdamperdruk vanaf dan sneller begint te dalen. Als deze reactie optreedt nog voor de uitlaatenthalpie kleiner wordt dan damp saturatie-enthalpie, verdwijnt de dampzone niet. Immers, de damp saturatie-enthalpie begint


71

samen met de druk sneller te dalen, zodat het teken van het verschil met de uitlaatenthalpie behouden blijft. Wanneer de tijdsconstante te groot is, blijft deze reactie te lang uit, zodat de dampzone wegvalt. Uit de metingen blijkt dat dit nooit voorvalt. Daarom is er gekozen voor een lagere tijdsconstante. Zoals te zien in figuur 6.8, daalt de druk hier sneller, zodat de damp saturatie-enthalpie altijd onder de uitlaatenthalpie blijft liggen en de dampzone dus behouden blijft. Eigenlijk hangt de bovengrens af van de veronderstelde tijdsconstante van de sturing van de inlaatschoepen. Er is verondersteld dat de dynamica van de aansturing verwaarloosbaar is, zodat deze tijdsconstante gelijk is aan nul. Echter, hoe groter deze waarde verondersteld wordt, hoe groter de tijdsconstante van het reservoir kan gekozen worden. Immers, hoe trager de inlaatschoepen reageren op een stijging in de PLR, hoe trager de uitlaatenthalpie daalt, want het massadebiet aan de uitlaat neemt niet zo snel toe. Dit compenseert dus de late reactie van de klep. τ = 10 s

τ=5s

τ=0s

1.32

1.32

1.32

1.3

1.3

1.3

1.28

1.28

1.28 [kg/s]

1.34

[kg/s]

1.34

[kg/s]

1.34

1.26

1.26

1.26

1.24

1.24

1.24

1.22

1.22

1.22

1.2 500

1000 1500 t [s]

2000

1.2 500

1000 1500 t [s] mtev

2000

1.2 500

1000 1500 t [s]

2000

mcomp

Figuur 6.7: Invloed te kleine tijdsconstante van het reservoir Uit figuur 6.6 blijkt dat het massadebiet in de simulatie veel sneller reageert op een stijging in PLR dan in realiteit. Het verschil met de meetresultaten suggereert dus dat beide tijdsconstanten te klein zijn gemodelleerd. Dit volgt ook uit een vergelijking met het werk van Comstock en Braun [11] die dezelfde cyclus gemodelleerd hebben. Deze onderzoekers veronderstellen dat de tijdsconstante van het reservoir gelijk is aan 100 s. Comstock en Braun [11] nemen ook de traagheid van de inlaatschoepen in rekening. Deze traagheid is volgens hen te wijten aan het hydraulisch systeem dat de schoepen aanstuurt. Echter, nergens in het werk wordt de grootte-orde van deze traagheid vermeld.

6.4.3

Schakelen tussen modes

Het schakelen wordt getest door de inlaatvoorwaarden uit tabel 6.1 constant te houden terwijl een sprong aangelegd wordt in het instelpunt van de uitlaattemperatuur van het ijswater. Dit komt dus neer op een sprong in de PLR, zoals weergegeven in figuur 6.9. Het resulterend tijdsverloop van de toestandsgrootheden, het massadebiet en de interne variabelen van de verdamper is weergegeven in figuren 6.10 en 6.11.


72 τ = 25 s

5

3.96

x 10

4.038 mr,in

3.94

0.85

3.88

4.036 her,out [J/kg]

mer [kg/s]

per [Pa]

3.9

0.8 0.75

3.86

4.034 hg,r

4.032

responsietijd

3.82 3740

3760

3780 tijd [s]

0.65 3740

3800

3760

x 10

4.031 3740

3800

4.038 mr,in

0.85

3.88

0.8 0.75 responsietijd

3.86

3760

3780 tijd [s]

3800

0.65 3740

her,out

x 10

4.035 4.034 4.033

0.7

3.84

3800

4.036 her,out [J/kg]

mer [kg/s]

3.9

3780 tijd [s]

4.037

mr,out

3.92

3760

5

0.9

3.94

3.82 3740

3780 tijd [s]

τ=5s

5

per [Pa]

4.035

4.033 0.7

3.84

3.96

x 10

4.037

mr,out

3.92

her,out

5

0.9

hg,r

4.032 3760

3780 tijd [s]

4.031 3740

3800

3760

3780 tijd [s]

3800

Figuur 6.8: Invloed te grote tijdsconstante van het reservoir 100 90 80 70

PLR [%]

60 50 40 30 20 10 0 600

800

1000

1200

1400

1600

tijd [s]

Figuur 6.9: Tijdsverloop PLR om schakelen te testen

verdamper condensor

T2,in [K] 284.65 297.14

m ˙ 2 [kg/s] 13.57 16.8

p2 [P a] 101356.5 101356.6

Tabel 6.1: Inlaatvoorwaarden om schakelen te testen De tijdsconstante van het reservoir is moedwillig ingesteld op 75 s. Zoals eerder aangetoond, verdwijnt hierdoor de dampzone aan de uitlaat van de verdamper. Deze zone verdwijnt ongeveer op t = 1035 s en treedt opnieuw in op t = 1385 s. In figuur 6.10 is het effect van de pseudo-vergelijkingen weergegeven. Zowel de uitlaattemperatuur van het secundair flu¨ıdum als de buiswandtemperatuur van de inactieve V-zone naderen asymptotisch naar de corres-


73

Fezone

[−]

1

0.5

0 600

700

800

900

1000

1100 tijd [s]

1200

1300

1400

1500

1600

1200

1300

1400

1500

1600

1200

1300

1400

1500

1600

(Te2,out)zone 290

[K]

285 280 275 600

700

800

900

1000

1100 tijd [s] Tet

290

[K]

285 280 275 600

700

800

900

1000

1100 tijd [s] zone TP

zone V

Figuur 6.10: Tijdsverloop zone-groottes en temperaturen in de verdamper tijdens schakelen per

5

4

x 10

[J/kg]

[Pa]

x 10

4

3.8

3.6

3.4

3.2 600

her

5

4.5

3.5

in out

3 2.5

800

1000 1200 tijd [s]

1400

2 600

1600

800

mer

1000 1200 tijd [s]

1400

1600

1400

1600

Te2,out

1.4

283.5 283

1.2

[K]

[kg/s]

282.5 1 282

0.8 281.5 0.6 0.4 600

in out 800

1000 1200 tijd [s]

1400

1600

281 280.5 600

800

1000 1200 tijd [s]

Figuur 6.11: Tijdsverloop toestandsgrootheden en massadebiet verdamper tijdens schakelen


74

ponderende waarden van de actieve TP-zone. Ook de zone-grootte nadert naar een kleine, maar positieve waarde . Bij herintrede, verdwijnen de pseudo-vergelijkingen uit het stelsel en worden vervangen door de behoudswetten opgesteld in de V-zone. Uit figuur 6.11 blijkt dat er zowel bij het heen- als terugschakelen geen discontinu¨ıteiten optreden in de toestandsgrootheden. Dit versterkt aldus het vertrouwen in een moving-boundary modellering.

6.5

Conclusie

In dit hoofdstuk is een model ontwikkeld dat in staat is het dynamisch gedrag van de koelmachine beschreven in [11] te reproduceren. Het gesimuleerde transiënt gedrag komt kwalitatief goed overeen met het opgemeten gedrag. Er wordt geduid dat de traagheid van zowel het reservoir van de expansieklep als de sturing van de inlaatschoepen van de compressor een invloed hebben op het verloop van het massadebiet in de cyclus. De tijdsconstante van het reservoir mag niet te klein zijn omdat dit oscillaties in het klepdebiet induceert. Deze waarde mag echter ook niet te groot zijn omdat anders zones verdwijnen. Uit de meetdata blijkt echter dat de uitlaat van de verdamper altijd oververhit is, terwijl de uitlaat van de condensor altijd onderkoeld is. Verder is aangetoond dat de traagheid van de inlaatschoepen de limietwaarden van de reservoirtijdsconstante be¨ınvloedt. Ten slotte is gedemonstreerd dat het verdwijnen van een zone op zich geen probleem stelt. Bij het schakelen blijft de continu¨ıteit van de toestandsgrootheden behouden dankzij het invoeren van de pseudo-vergelijkingen.

Hoofdstuk 7. Conclusie

75

Hoofdstuk 7

Conclusie 7.1

Conclusie

In dit huidig werk is een theoretisch model ontwikkeld dat in staat is zowel het statisch als het dynamisch gedrag van een subkritische compressie-koelcyclus te voorspellen. De beschikbare hoeveelheid constructieve data is echter zeer beperkt, daarom is er bij de modellering van klep en compressor in beperkte mate gebruik gemaakt van data-fitting. De warmtewisselaars zijn gemodelleerd volgens de moving-boundary methode. Vervolgens zijn de simulatieresultaten vergeleken met de meetdata. Het gesimuleerde steady-state deellastgedrag komt kwalitatief overeen met het opgemeten gedrag. Verder worden ook de toestandsgrootheden, uitgezonderd de verdamperen condensordruk, goed voorspeld. Hoewel de gesimuleerde en opgemeten trends overeenkomen, zit er ook een vrij grote afwijking op de voorspelling van de COP. Er wordt aangetoond dat beide problemen te herleiden zijn tot de gekozen correlaties voor de warmteoverdracht. Verder wordt geduid dat ook het elektromechanisch rendement van de compressor een rol speelt in de voorspelling van de COP. Ook wordt de invloed van de void-fractie correlatie aangekaart en het effect van de polytrope exponent op het deellastgedrag toegelicht. Daarnaast wordt ook het dynamisch gedrag van de koelcyclus kwalitatief goed beschreven. Zowel de duur als het verloop van de overgangsverschijnselen worden accuraat gereproduceerd. Hieruit blijkt dat een steady-state datafit voldoende is om het dynamisch gedrag te capteren. Ook wordt het effect ge¨ıllustreerd van de tijdsconstante van zowel het reservoir van de expansieklep als van de sturing van de inlaatschoepen van de compressor. Ten slotte wordt aangetoond hoe het warmtewisselaarmodel zich gedraagt in abnormale omstandigheden, waarin de uitlaat van de verdamper niet meer oververhit is of waarin de uitlaat van de condensor 2-fasig wordt. Er is ook grote aandacht besteed aan de modulariteit en de gebruiksvriendelijkheid van het simulatiemodel. Dit behoorde niet zozeer tot de hoofdopdracht van de thesis. Toch is dit een belangrijk aspect van een simulatiemodel. Het biedt immers de mogelijkheid om het model in een latere fase op een snelle en eenvoudige wijze uit te breiden en te gebruiken om ook andere koelcycli te simuleren en te bestuderen. Modulariteit is immers 1 van de grote voordelen van numerieke modellen.

7.2

Aanbevelingen voor verder onderzoek

Er moet nog verder onderzoek gebeuren naar de invloed van de veronderstellingen die gedaan zijn bij het modelleren. Ten eerste is er verondersteld dat de drukval in beide warmtewisselaars nul is. Er moet echter nagegaan worden welke invloed dit heeft, zowel statisch als dynamisch. Verder is de dynamica van klep en compressor verwaarloosd. Ook het gevolg van choking in de klep wordt niet in rekening gebracht. Er zijn evenwel modellen ontwikkeld die hier wel rekening mee houden. Gezien de modulariteit van het simulatiemodel kunnen deze eenvoudig

Hoofdstuk 7. Conclusie

76

ge¨ımplementeerd worden om het effect op cyclusniveau te bestuderen. Ten slotte is geponeerd dat het enthalpieverloop in de verdamper stuksgewijze lineair is. Het wordt aanbevolen om ook de weerslag van deze veronderstelling te analyseren. Een ander onderwerp van verder onderzoek omvat het uitbreiden van het model, zodat zijn toepasbaarheid vergroot. Het compressor- en klepmodel kunnen bijvoorbeeld uitgebreid worden zodat het cyclusmodel ook in staat is het start-stop gedrag te voorspellen. Daarnaast kan het verdamperen condensormodel aangepast worden, zodat ze ook bruikbaar zijn voor andere types warmtewisselaars.

Bijlage A. Proefstand

77

Bijlage A

Proefstand Tabel A.3 geeft een overzicht van alle gemeten en berekende waarden. Ook wordt de absolute fout van iedere meting weergegeven, zowel in deellast- als in vollastomstandigheden. De relatieve fout wordt bepaald door de absolute fout te delen door de meetwaarde. De gebruikte nomenclatuur wordt verklaard in tabel A.1. Deze tabellen komen rechtstreeks uit [11]. Voor meer informatie wat betreft instrumentatie en de opbouw van de proefstand wordt de lezer doorverwezen naar [11].

37

Bijlage A. Proefstand

78

Table 3.2: Exported data from experimental test runs. Designation

Source

Description

Units

Time TWE_set TEI

VisSim MicroTech JCI AHU (RTD) MicroTech (Thermistor) JCI AHU (RTD) MicroTech (Thermistor) JCI AHU (RTD) MicroTech (Thermistor) JCI AHU (RTD) MicroTech (Thermistor) JCI AHU (RTD) JCI AHU (RTD) JCI AHU (RTD) JCI AHU (RTD) VisSim VisSim VisSim

Real time counter Chilled water setpoint—control variable Temperature of Evaporator Water In

Seconds F F

Temperature of Evaporator Water In

F

Temperature of Evaporator Water Out

F

Temperature of Evaporator Water Out

F

Temperature of Condenser Water In

F

Temperature of Condenser Water In

F

Temperature of Condenser Water Out

F

Temperature of Condenser Water Out

F

Temperature of Shared HX Water In (in Condenser Water Loop) Temperature of Shared HX Water Out (in Condenser Water Loop) Temperature of Building Water In (in Evaporator Water Loop) Temperature of Building Water Out (in Evaporator Water Loop) Calculated Condenser Heat Rejection Rate Calculated City Water Cooling Rate Calculated Shared HX Heat Transfer (only valid with no water bypass) Calculated 1st Law Energy Balance for Condenser Water Loop (only valid with no water bypass) Calculated Evaporator Cooling Rate Calculated Shared HX Heat Transfer (should equal Shared Cond Tons with no water bypass) Calculated Steam Heating Load Calculated 1st Law Energy Balance for Evaporator Water Loop Watt Transducer Measuring Instantaneous Compressor Power Calculated Coefficient of Performance Calculated Compressor Efficiency Flow Rate of Condenser Water Flow Rate of Evaporator Water Evaporator Approach Temperature (TWEO-TRE) Condenser Approach Temperature (TRC-TWCO) Saturated Refrigerant Temperature in Evaporator

F

TWEI TEO TWEO TCI TWCI TCO TWCO TSI TSO TBI TBO Cond Tons Cooling Tons Shared Cond Tons Cond Energy Balance Evap Tons Shared Evap Tons Building Tons Evap Energy Balance kW COP kW/ton FWC FWE TEA TCA TRE

VisSim VisSim VisSim VisSim VisSim JCI AHU VisSim VisSim JCI AHU JCI AHU MicroTech MicroTech MicroTech

Tabel A.1: Nomenclatuur

F F F Tons Tons Tons Tons Tons Tons Tons Tons kW -kW/ton GPM GPM F F F

38 Bijlage A. Proefstand

79

Table 3.2: Continued. Designation

Source

Description

Units

PRE TRC PRC TRC_sub T_suc Tsh_suc TR_dis Tsh_dis P_lift Amps RLA% Heat Balance (kW) Heat Balance% Tolerance%

MicroTech MicroTech MicroTech MicroTech MicroTech MicroTech MicroTech MicroTech MicroTech MicroTech MicroTech VisSim

Pressure of Refrigerant in Evaporator Saturated Refrigerant Temperature in Condenser Pressure of Refrigerant in Condenser Liquid-line Refrigerant Subcooling from Condenser Refrigerant Suction Temperature Refrigerant Suction Superheat Temperature Refrigerant Discharge Temperature Refrigerant Discharge Superheat Temperature Pressure Lift Across Compressor Current Draw Across One Leg of Motor Input Percent of Maximum Rated Load Amps Calculated 1st Law Energy Balance for Chiller

PSIG F PSIG F F F F F PSI Amps % kW

VisSim VisSim

% %

Unit Status Active Fault TO_sump TO_feed PO_feed PO_net TWCD TWED VSS VSL VH VM VC VE VW TWI

MicroTech MicroTech MicroTech MicroTech MicroTech MicroTech MicroTech MicroTech JCI AHU JCI AHU JCI AHU JCI AHU JCI AHU JCI AHU JCI AHU JCI AHU (RTD) JCI AHU (RTD) JCI AHU (RTD) JCI AHU (RTD) VisSim VisSim VisSim

Calculated 1st Law Energy Balance for Chiller Calculated Heat Balance Tolerance (ARI 550 defined as allowable test tolerance on heat balance) Consult Table B.4 in Appendix Consult Table B.3 in Appendix Temperature of Oil in Sump Temperature of Oil Feed Pressure of Oil Feed Oil Feed minus Oil Vent Pressure Condenser Water Temperature Difference (TWCO-TWCI) Evaporator Water Temperature Difference (TWEI-TWEO) Small Steam Valve Position Large Steam Valve Position Hot Water Valve Position 3-way Mixing Valve Position Condenser Valve Position Evaporator Valve Position City Water Valve Position Temperature of City Water In

0 – 27 0 – 44 F F PSIG PSI F F % Open % Open % Open % Mix % Open % Open % Open F

Temperature of City Water Out

F

Temperature of Hot Water In

F

Temperature of Hot Water Out

F

Calculated City Water Flow Rate Calculated Hot Water Flow Rate Calculated Condenser Water Bypass Flow Rate

GPM GPM GPM

TWO THI THO FWW FWH FWB

Tabel A.2: Vervolg tabel A.1


Table 3.4: Complete listing of measured and calculated variables with corresponding80 absolute uncertainties. Designation

Source

Full Load Uncertainty

Low Load Uncertainty

Units

TEI TWEI

JCI AHU (RTD) MicroTech (Thermistor) JCI AHU (RTD) MicroTech (Thermistor) JCI AHU (RTD) MicroTech (Thermistor) JCI AHU (RTD) MicroTech (Thermistor) JCI AHU (RTD) JCI AHU (RTD) JCI AHU (RTD) JCI AHU (RTD) VisSim VisSim VisSim VisSim VisSim VisSim VisSim VisSim JCI AHU VisSim VisSim JCI AHU JCI AHU MicroTech MicroTech MicroTech MicroTech MicroTech MicroTech MicroTech MicroTech MicroTech MicroTech MicroTech MicroTech MicroTech MicroTech VisSim VisSim

±0.05 ±0.2

±0.05 ±0.2

F F

±0.05 ±0.2

±0.05 ±0.2

F F

±0.05 ±0.2

±0.05 ±0.2

F F

±0.05 ±0.2

±0.05 ±0.2

F F

±0.05 ±0.05 ±0.05 ±0.05 ±1.32 ±0.9 ±1.32* ±2.07* ±1.08 ±1.08 ±0.63 ±1.65 ±1.8 ±0.117 ±0.026 ±2.8 ±2.2 ±0.3 ±0.5 ±0.3 ±0.3 ±0.3 ±0.5 ±0.54 ±0.2 ±0.36 ±0.2 ±0.54 ±0.6 ±5 ±4 ±12.4 ±3.1

±0.05 ±0.05 ±0.05 ±0.05 ±0.7 ±0.85 ±0.7* ±1.3* ±0.75 ±0.75 ±1.05 ±1.49 ±1.8 ±0.153 ±0.071 ±2.8 ±2.2 ±0.3 ±0.5 ±0.3 ±0.3 ±0.3 ±0.5 ±0.54 ±0.2 ±0.36 ±0.2 ±0.54 ±0.6 ±5 ±4 ±6.72 ±5.6

F F F F Tons Tons Tons Tons Tons Tons Tons Tons kW --GPM GPM F F F PSIG F PSIG F F F F F PSI Amps % kW %

TEO TWEO TCI TWCI TCO TWCO TSI TSO TBI TBO Cond Tons Cooling Tons Shared Cond Tons Cond Energy Balance Evap Tons Shared Evap Tons Building Tons Evap Energy Balance kW COP kW/ton FWC FWE TEA TCA TRE PRE TRC PRC TRC_sub T_suc Tsh_suc TR_dis Tsh_dis P_lift Amps RLA% Heat Balance (kW) Heat Balance%

Tabel A.3: Absolute fouten


81

Table 3.4: Continued. Designation

Source

Full Load Uncertainty

Low Load Uncertainty

Units

Tolerance% TO_sump TO_feed PO_feed PO_net TWCD TWED TWI TWO THI THO FWW FWH FWB

VisSim MicroTech MicroTech MicroTech MicroTech MicroTech MicroTech JCI AHU (RTD) JCI AHU (RTD) JCI AHU (RTD) JCI AHU (RTD) VisSim VisSim VisSim

±5.7 ±0.2 ±0.2 ±0.5 ±0.5 ±0.6 ±0.4 ±0.05 ±0.05 ±0.05 ±0.05 ±4.5 ±3.5 ±4.9

±5.7 ±0.2 ±0.2 ±0.5 ±0.5 ±0.6 ±0.4 ±0.05 ±0.05 ±0.05 ±0.05 ±7.0 ±25 ±7.3

% F F PSIG PSI F F F F F F GPM GPM GPM

*

denotes a calculation whose results are unreliable if any condenser water is bypassed around shared heat exchanger

Tabel A.4: Vervolg tabel A.3

The values given in Table 3.4 are the absolute uncertainties for each measurement, to find the uncertainty it is necessary to divide by the measurement. Important measurements with the highest uncertainties are: •

Suction superheat with an uncertainty of 5 to 25%

•

Liquid subcooling with an uncertainty of 5 to 15%

•

Evaporator approach temperature with an uncertainty of 4 to 10%

•

Condenser approach temperature with an uncertainty of 8 to 25%

•

kW/ton with an uncertainty of 4 to 5%

Differences in uncertainties occur because of significant changes in the measurements at various operating conditions.

3.7

Test Matrix

Soon after commissioning the chiller test stand, tests were run at various temperature and loading extremes to determine the operating envelope of the chiller. ARI Standard 550 recommends testing at a chilled water temperature of 44°F (TEO) and a condenser water entering temperature of 85°F (TCI) for standard rating conditions.

Bijlage B. Dynamische zone-specifieke behoudsvergelijkingen verdamper en condensor

82

Bijlage B

Dynamische zone-specifieke behoudsvergelijkingen verdamper en condensor

Verdamper

Deze sectie geeft een overzicht van de dynamische zone-specifieke behoudsvergelijkingen van de verdamper. Het TP-V-model is afgeleid in hoofdstuk 4. De vergelijkingen die het gedrag in de TP-mode beschrijven, zijn ook weergegeven. Deze zijn op dezelfde manier afgeleid. In beide modes geldt dat: FT P + FV = 1

B.1.1

(B.1.1)

TP-V-mode

Buiswand dTt,T P dFT P = ρt Ct Vt FT P + FT P (Tt,T P − Tt,V ) dt dt dTt,V dFV = ρt Ct Vt FV + FV (Tt,V − Tt,T P ) dt dt

Q2t,T P + Qrt,T P Q2t,V + Qrt,V

(B.1.2) (B.1.3)

Secundair flu¨ıdum d(T2,out )T P dT2,in dFT P ρ2 Cp,2 V2 FT P = ρ2 Cp,2 V2 1/2 FT P + FT P (T2,T P − T2,V ) −Q2t,T P − FT P m ˙ 2 Cp,2 [(T2,out )T P − T2,in ] − 1/2 dt dt dt dT2,in d(T2,out )V dFV −Q2t,V − FV m ˙ 2 Cp,2 [(T2,out )V − T2,in ] − 1/2 ρ2 Cp,2 V2 FV = ρ2 Cp,2 V2 1/2 FV + FV (T2,V − T2,T P ) dt dt dt

(B.1.4) (B.1.5)

Koelmiddel • massabehoud:

(B.1.6)

83

dFT P [(¯ γ − 1)(ρg,r − ρf,r )] Vr dt dρf,r dρg,r dpr ∂¯ γ + γ¯ + (1 − γ¯ ) + (ρg,r − ρf,r ) FT P Vr dt dpr dpr ∂pr dhr,in ∂¯ γ + (ρg,r − ρf,r ) FT P Vr dt ∂hr,in

m ˙ r,in − m ˙ r,T P −V =


B.1

(B.1.7)

Optellen van beide vergelijkingen resulteert in: dFT P dFV [(¯ γ − 1)(ρg,r − ρf,r )] Vr + [ρr,V − ρg,r ] Vr dt dt dρf,r ∂ρr,V dhg,r ∂ρr,V dρg,r dpr ∂¯ γ + FV 1/2 + + FT P γ¯ + (1 − γ¯ ) + (ρg,r − ρf,r ) Vr dt ∂hr,V dpr ∂pr dpr dpr ∂pr dhr,in ∂¯ γ (ρg,r − ρf,r ) FT P Vr + dt ∂hr,in ∂ρr,V dhr,out 1/2 FV Vr + dt ∂hr,V

m ˙ r,in − m ˙ r,out =

(B.1.8)

• energiebehoud:

dFT P [(¯ γ − 1) ((ρh)g,r − (ρh)f,r )] Vr dt d(ρh)f,r d(ρh)g,r dpr ∂¯ γ + γ¯ + (1 − γ¯ ) + ((ρh)g,r − (ρh)f,r ) − 1 FT P Vr dt dpr dpr ∂pr dhr,in ∂¯ γ + ((ρh)g,r − (ρh)f,r ) FT P Vr dt ∂hr,in

m ˙ r,in [hr,in − hg,r ] − Qrt,T P =

− hg,r [m ˙ r,in − m ˙ r,T P −V ]

(B.1.9)


dFV [ρr,V − ρg,r ] Vr dt ∂ρr,V dhg,r ∂ρr,V dpr + 1/2 + FV Vr dt ∂hr,V dpr ∂pr ∂ρr,V dhr,out 1/2 + FV Vr dt ∂hr,V

m ˙ r,T P −V − m ˙ r,out =

84

(B.1.10)

− hg,r [m ˙ r,T P −V − m ˙ rout ]

B.1.2

TP-mode

Buiswand Q2t,T P + Qrt,T P = ρt Ct Vt

dTt,T P FT P dt

(B.1.11)

Secundair flu¨ıdum −Q2t,T P − FT P m ˙ 2 Cp,2 [(T2,out )T P

d(T2,out )T P dT2,in ρ2 Cp,2 V2 FT P = ρ2 Cp,2 V2 1/2 FT P − T2,in ] − 1/2 dt dt

(B.1.12) (B.1.13)


(B.1.14)

85

dρf,r dρg,r dpr ∂¯ γ m ˙ r,in − m ˙ r,out = γ¯ + (1 − γ¯ ) + (ρg,r − ρf,r ) FT P Vr dt dpr dpr ∂pr dhr,in ∂¯ γ + (ρg,r − ρf,r ) FT P Vr dt ∂hr,in dhr,out ∂¯ γ + (ρg,r − ρf,r ) FT P Vr dt ∂hr,out


dFV [(ρh)r,V − (ρh)g,r ] Vr dt ∂(ρh)r,V dhg,r ∂(ρh)r,V dpr + 1/2 + − 1 FV Vr dt ∂hr,V dpr ∂pr ∂(ρh)r,V dhr,out + 1/2 FV Vr dt ∂hr,V

m ˙ r,out [hg,r − hr,out ] − Qrt,V =

m ˙ r,in hr in − m ˙ r,out hr,out − Qrt,T P

d(ρh)f,r d(ρh)g,r dpr ∂¯ γ = γ¯ + (1 − γ¯ ) + ((ρh)g,r − (ρh)f,r ) − 1 FT P Vr dt dpr dpr ∂pr dhr,in ∂¯ γ ((ρh)g,r − (ρh)f,r ) FT P Vr + dt ∂hr,in dhr,out ∂¯ γ + ((ρh)g,r − (ρh)f,r ) FT P Vr dt ∂hr,out

(B.1.15)

Pseudo-vergelijkingen

kF

dFV = − FV dt

dTt,V = Tt,T P − Tt,V dt d(T2,out )V k2 = (T2,out )T P − (T2,out )V dt kt

(B.1.16) (B.1.17) (B.1.18)

werd tijdens de simulaties gelijk gesteld aan 0.001. Voor kF , kt en k2 zijn respectievelijk de waarden 1000, 15 en 15 s gekozen.

B.2

Condensor

Deze sectie geeft een overzicht van de dynamische zone-specifieke behoudsvergelijkingen van de condensor. Het V-TP-L-model is afgeleid in hoofdstuk 4. De vergelijkingen die het gedrag in de V-TP-mode beschrijven, zijn ook weergegeven. Deze zijn op dezelfde manier afgeleid. In beide modes geldt dat: FL + FT P + FV = 1

(B.2.1)


• energiebehoud:

86

V-TP-L-mode

Buiswand dTt,V dFV FV FV + (Tt,V − Tt,T P ) dt dt FV + FT P dTt,T P dFV FV Tt,T P + FT P Tt,V dFT P dFL FL Tt,T P + FT P Tt,L = ρt Ct Vt FT P + + Tt,T P + dt dt FV + FT P dt dt FL + FT P dTt,L dFL FL FL + (Tt,L − Tt,T P ) = ρt Ct Vt dt dt FL + FT P

Q2t,V + Qrt,V = ρt Ct Vt Q2t,T P + Qrt,T P Q2t,L + Qrt,L

(B.2.2) (B.2.3) (B.2.4)

Secundair flu¨ıdum d(T2,out )V dT2,in dFV FV ρ2 Cp,2 V2 FV = ρ2 Cp,2 V2 1/2 FV + (T2,V − T2,T P ) −Q2t,V − FV m ˙ 2 Cp,2 [(T2,out )V − T2,in ] − 1/2 dt dt dt FV + FT P d(T2,out )T P dT2,in dFV FV T2,T P + FT P T2,V ρ2 Cp,2 V2 FT P = ρ2 Cp,2 V2 1/2 FT P + −Q2t,T P − FT P m ˙ 2 Cp,2 [(T2,out )T P − T2,in ] − 1/2 dt dt dt FV + FT P dFT P dFL FL T2,T P + FT P T2,L + T2,T P + dt dt FL + FT P dT2,in d(T2,out )L dFL FL ρ2 Cp,2 V2 FL = ρ2 Cp,2 V2 1/2 FL + (T2,L − T2,T P ) −Q2t,L − FL m ˙ 2 Cp,2 [(T2,out )L − T2,in ] − 1/2 dt dt dt FL + FT P

(B.2.5) (B.2.6) (B.2.7) (B.2.8)


(B.2.9)

87

dFV [ρr,V − ρg,r ] Vr dt ∂ρr,V ∂ρr,V dhg,r dpr + 1/2 + FV Vr dt ∂hr,V dpr ∂pr ∂ρr,V dhr,in + 1/2 FV Vr dt ∂hr,V

m ˙ r,in − m ˙ r,V −T P =


B.2.1

dFL [ρr,L − ρf,r ] Vr dt ∂ρr,L dhf,r ∂ρr,L dpr + 1/2 + FL Vr dt ∂hr,L dpr ∂pr ∂ρr,L dhr,out 1/2 + FL Vr dt ∂hr,L

(B.2.10)

m ˙ r,T P −L − m ˙ r,out =

(B.2.11)

Optellen van deze 3 vergelijkingen resulteert in: dFV dFT P dFL [ρr,V − ρg,r ] Vr + [(¯ γ − 1)(ρg,r − ρf,r )] Vr + [ρr,L − ρg,r ] Vr dt dt dt dρf,r ∂ρr,V dhg,r ∂ρr,V ∂ρr,L dhf,r ∂ρr,L dρg,r d¯ γ dpr + + (1 − γ¯ ) + (ρg,r − ρf,r ) + FL 1/2 + FV 1/2 + FT P γ¯ Vr + dt ∂hr,V dpr ∂pr dpr dpr dpr ∂hr,L dpr ∂pr ∂ρr,V dhr,in 1/2 FV V r + dt ∂hr,V ∂ρr,L dhr,out + 1/2 FL Vr dt ∂hr,L (B.2.12)


• energiebehoud:


dFT P dFL [(¯ γ − 1)(ρg,r − ρf,r )] Vr + [ρf,r − ρg,r ] Vr dt dt dρf,r dρg,r dpr d¯ γ + γ¯ + (1 − γ¯ ) + (ρg,r − ρf,r ) FT P Vr dt dpr dpr dpr

m ˙ r,V −T P − m ˙ r,T P −L =

88

(B.2.13)

− hg,r [m ˙ r,in − m ˙ r,V −T P ]

dFT P [(¯ γ − 1)((ρh)g,r − (ρh)f,r )] Vr dt dFL + [(ρh)f,r − (ρh)g,r ] Vr dt d(ρh)f,r d(ρh)g,r dpr d¯ γ γ¯ + + (1 − γ¯ ) + ((ρh)g,r − (ρh)f,r ) − 1 FT P Vr dt dpr dpr dpr + hg,r [m ˙ r,in − m ˙ r,V −T P ] + hf,r [m ˙ r,T P −L − m ˙ r,out ]

m ˙ r,in hg,r − m ˙ r,out hf,r − Qrt,T P =

dFL [(ρh)r,L − (ρh)f,r ] Vr dt ∂(ρh)r,L dhf,r ∂(ρh)r,L dpr + 1/2 + − 1 FL Vr dt ∂hr,L dpr ∂pr ∂(ρh)r,L dhr,out + 1/2 FL Vr dt ∂hr,L

(B.2.14)

m ˙ r,out [hf,r − hr,out ] − Qrt,L =

− hf,r [m ˙ r,T P −L − m ˙ r,out ]

(B.2.15)


dFV [(ρh)r,V − (ρh)g,r ] Vr dt ∂(ρh)r,V dhg,r ∂(ρh)r,V dpr + 1/2 + − 1 FV V r dt ∂hr,V dpr ∂pr ∂(ρh)r,V dhr,in + 1/2 FV Vr dt ∂hr,V

m ˙ r,in [hr,in − hg,r ] − Qrt,V =

89

V-TP-mode

Buiswand dTt,V dFV FV FV + (Tt,V − Tt,T P ) = ρt Ct Vt dt dt FV + FT P dTt,T P dFT P FT P = ρt Ct Vt FT P + (Tt,T P − Tt,V ) dt dt FT P + FV

Q2t,V + Qrt,V Q2t,T P + Qrt,T P

(B.2.16) (B.2.17)

Secundair flu¨ıdum dT2,in d(T2,out )V dFV FV −Q2t,V − FV m ˙ 2 Cp,2 [(T2,out )V − T2,in ] − 1/2 ρ2 Cp,2 V2 FV = ρ2 Cp,2 V2 1/2 FV + (T2,V − T2,T P ) (B.2.18) dt dt dt FV + FT P dT2,in d(T2,out )T P dFT P FT P (T2,T P − T2,V ) (B.2.19) −Q2t,T P − FT P m ˙ 2 Cp,2 [(T2,out )T P − T2,in ] − 1/2 ρ2 Cp,2 V2 FT P = ρ2 Cp,2 V2 1/2 FT P + dt dt dt FV + FT P Koelmiddel • massabehoud:

(B.2.20)

dFT P [(¯ γ − 1)(ρg,r − ρf,r )] Vr dt dρf,r dρg,r dpr ∂¯ γ + γ¯ + (1 − γ¯ ) + (ρg,r − ρf,r ) FT P Vr dt dpr dpr ∂pr dhr,out ∂¯ γ + (ρg,r − ρf,r ) FT P Vr dt ∂hr,out

(B.2.21)

m ˙ r,V −T P − m ˙ r,out =

90

dFV [ρr,V − ρg,r ] Vr dt ∂ρr,V dhg,r ∂ρr,V dpr + 1/2 + FV Vr dt ∂hr,V dpr ∂pr ∂ρr,V dhr,in + 1/2 FV Vr dt ∂hr,V

m ˙ r,in − m ˙ r,V −T P =


B.2.2

dFV dFT P [ρr,V − ρg,r ] Vr + [(¯ γ − 1)(ρg,r − ρf,r )] Vr dt dt dρf,r ∂ρr,V dhg,r ∂ρr,V dρg,r dpr ∂¯ γ + FV 1/2 + + FT P γ¯ + (1 − γ¯ ) + (ρg,r − ρf,r ) Vr dt ∂hr,V dpr ∂pr dpr dpr ∂pr ∂ρr,V dhr,in 1/2 FV Vr + dt ∂hr,V dhr,out ∂¯ γ + (ρg,r − ρf,r ) FT P Vr dt ∂hr,out


(B.2.22)

• energiebehoud:

dFV [(ρh)r,V − (ρh)g,r ] Vr dt ∂(ρh)r,V dhg,r ∂(ρh)r,V dpr + 1/2 + − 1 FV V r dt ∂hr,V dpr ∂pr ∂(ρh)r,V dhr,in 1/2 FV Vr + dt ∂hr,V

m ˙ r,in [hr,in − hg,r ] − Qrt,V =

(B.2.23)

− hg,r [m ˙ r,in − m ˙ r,V −T P ]

dFT P [(¯ γ − 1)((ρh)g,r − (ρh)f,r )] Vr dt ∂(ρh)f,r d(ρh)g,r dpr d¯ γ + γ¯ + (1 − γ¯ ) + ((ρh)g,r − (ρh)f,r ) − 1 FT P Vr dt dpr ∂pr dpr dhr,out ∂¯ γ + ((ρh)g,r − (ρh)f,r ) FT P Vr dt ∂hr,out

m ˙ r,out [hg,r − hr,out ] − Qrt,T P =

(B.2.24)


Optellen van deze 2 vergelijkingen resulteert in:

− hg,r [m ˙ r,V −T P − m ˙ r,out ] 91

kF

dFL = − FL dt

dTt,L = Tt,T P − Tt,L dt d(T2,out )L = (T2,out )T P − (T2,out )L k2 dt kt

werd tijdens de simulaties gelijk gesteld aan 0.001. Voor kF , kt en k2 zijn respectievelijk de waarden 1000, 15 en 15 s gekozen.

(B.2.25) (B.2.26) (B.2.27)



92

Bijlage C. Statische zone-specifieke behoudsvergelijkingen verdamper en condensor

93

Bijlage C

Statische zone-specifieke behoudsvergelijkingen verdamper en condensor C.1

Verdamper

Deze sectie geeft een overzicht van de statische zone-specifieke behoudsvergelijkingen van de verdamper. Deze betrekkingen kunnen eenvoudig afgeleid worden uit de 1ste hoofdwet van de thermodynamica of door de tijdsafgeleiden in het dynamisch model uit bijlage B gelijk te stellen aan nul. In elke mode geldt dat:

C.1.1

FL + FT P + FV = 1

(C.1.1)

Q2t,L + Qrt,L = 0

(C.1.2)

Q2t,T P + Qrt,T P = 0

(C.1.3)

Q2t,V + Qrt,V = 0

(C.1.4)

−Q2t,L − FL m ˙ 2 Cp,2 [(T2,out )L − T2,in ] = 0

(C.1.5)

−Q2t,T P − FT P m ˙ 2 Cp,2 [(T2,out )T P − T2,in ] = 0

(C.1.6)

−Q2t,V − FV m ˙ 2 Cp,2 [(T2,out )V − T2,in ] = 0

(C.1.7)

L-TP-V-mode

Buiswand

Secundair flu¨ıdum

Koelmiddel • massabehoud: m ˙r=m ˙ r,out =m ˙ r,in

(C.1.8)


94

• energiebehoud:

C.1.2

m ˙ r (hr,in − hf,r ) − Qrt,L = 0

(C.1.9)

m ˙ r (hf,r − hg,r ) − Qrt,T P = 0

(C.1.10)

m ˙ r (hg,r − hr,out ) − Qrt,V = 0

(C.1.11)

TP-V-mode

Buiswand Q2t,T P + Qrt,T P = 0

(C.1.12)

Q2t,V + Qrt,V = 0

(C.1.13)


(C.1.14)


(C.1.15)



(C.1.16)

• energiebehoud:

C.1.3

m ˙ r (hr,in − hg,r ) − Qrt,T P = 0

(C.1.17)

m ˙ r (hg,r − hr,out ) − Qrt,V = 0

(C.1.18)

L-TP-mode

Buiswand Q2t,L + Qrt,L = 0

(C.1.19)


(C.1.20)


(C.1.21)


(C.1.22)



(C.1.23)

• energiebehoud: m ˙ r (hr,in − hf,r ) − Qrt,L = 0

(C.1.24)

m ˙ r (hf,r − hr,out ) − Qrt,T P = 0

(C.1.25)


C.1.4

95

TP-mode


(C.1.26)


(C.1.27)



(C.1.28)

• energiebehoud: m ˙ r (hr,in − hr,out ) − Qrt,T P = 0

C.1.5

(C.1.29)

L-mode


(C.1.30)


(C.1.31)



(C.1.32)

• energiebehoud: m ˙ r (hr,in − hr,out ) − Qrt,L = 0

C.1.6

(C.1.33)

V-mode

Buiswand Q2t,V + Qrt,V = 0

(C.1.34)


(C.1.35)



96


(C.1.36)

• energiebehoud: m ˙ r (hr,in − hr,out ) − Qrt,V = 0

C.1.7

(C.1.37)


Bij het ontbreken van de V-zone 0 = − FV

(C.1.38)

0 = Tt,T P − Tt,V

(C.1.39)

0 = (T2,out )T P − (T2,out )V

(C.1.40)

werd tijdens de simulaties gelijk gesteld aan 0.001. Bij het ontbreken van de TP-zone 0 = − FT P

(C.1.41)

0 = Tt,x − Tt,T P

(C.1.42)

0 = (T2,out )x − (T2,out )T P

(C.1.43)

x staat voor een andere zone die wel actief is. Dit is de L-zone of de V-zone. werd tijdens de simulaties gelijk gesteld aan 0.001. Bij het ontbreken van de L-zone 0 = − FL

(C.1.44)

0 = Tt,T P − Tt,L

(C.1.45)

0 = (T2,out )T P − (T2,out )L

(C.1.46)

werd tijdens de simulaties gelijk gesteld aan 0.001.

C.2

Condensor

Deze sectie geeft een overzicht van de statische zone-specifieke behoudsvergelijkingen van de condensor. Deze betrekkingen kunnen eenvoudig afgeleid worden uit de 1ste hoofdwet van de thermodynamica of door de tijdsafgeleiden in het dynamisch model uit bijlage B gelijk te stellen aan nul. In elke mode geldt dat:

C.2.1

FV + FT P + FL = 1

(C.2.1)

Q2t,V + Qrt,V = 0

(C.2.2)


(C.2.3)

Q2t,L + Qrt,L = 0

(C.2.4)

V-TP-L-mode

Buiswand


97

Secundair flu¨ıdum −Q2t,V − FV m ˙ 2 Cp,2 [(T2,out )V − T2,in ] = 0

(C.2.5)


(C.2.6)


(C.2.7)


(C.2.8)

• energiebehoud:

C.2.2

m ˙ r (hr,in − hg,r ) − Qrt,V = 0

(C.2.9)

m ˙ r (hg,r − hf,r ) − Qrt,T P = 0

(C.2.10)

m ˙ r (hf,r − hr,out ) − Qrt,L = 0

(C.2.11)

V-TP-mode


(C.2.12)


(C.2.13)


(C.2.14)


(C.2.15)



(C.2.16)

• energiebehoud:

C.2.3

m ˙ r (hr,in − hg,r ) − Qrt,V = 0

(C.2.17)

m ˙ r (hg,r − hr,out ) − Qrt,T P = 0

(C.2.18)

TP-L-mode


(C.2.19)

Q2t,L + Qrt,L = 0

(C.2.20)


98

Secundair flu¨ıdum −Q2t,T P − FT P m ˙ 2 Cp,2 [(T2,out )T P − T2,in ] = 0

(C.2.21)


(C.2.22)


(C.2.23)

• energiebehoud:

C.2.4

m ˙ r (hr,in − hf,r ) − Qrt,T P = 0

(C.2.24)


(C.2.25)

TP-mode


(C.2.26)

Q2t,L + Qrt,L = 0

(C.2.27)


(C.2.28)


(C.2.29)



(C.2.30)

• energiebehoud:

C.2.5

m ˙ r (hr,in − hf,r ) − Qrt,T P = 0

(C.2.31)


(C.2.32)

L-mode


(C.2.33)


(C.2.34)



99


(C.2.35)

• energiebehoud: m ˙ r (hr,in − hr,out ) − Qrt,L = 0

C.2.6

(C.2.36)

V-mode


(C.2.37)


(C.2.38)



(C.2.39)

• energiebehoud: m ˙ r (hr,in − hr,out ) − Qrt,V = 0

C.2.7

(C.2.40)


Bij het ontbreken van de V-zone 0 = − FV

(C.2.41)

0 = Tt,T P − Tt,V

(C.2.42)

0 = (T2,out )T P − (T2,out )V

(C.2.43)

werd tijdens de simulaties gelijk gesteld aan 0.001. Bij het ontbreken van de TP-zone 0 = − FT P

(C.2.44)

0 = Tt,x − Tt,T P

(C.2.45)

0 = (T2,out )x − (T2,out )T P

(C.2.46)

x staat voor een andere zone die wel actief is. Dit is de L-zone of de V-zone. werd tijdens de simulaties gelijk gesteld aan 0.001. Bij het ontbreken van de L-zone 0 = − FL

(C.2.47)

0 = Tt,T P − Tt,L

(C.2.48)

0 = (T2,out )T P − (T2,out )L

(C.2.49)

werd tijdens de simulaties gelijk gesteld aan 0.001.

Bijlage D. Wiskundige technieken

100

Bijlage D

Wiskundige technieken In deze bijlage wordt een overzicht gegeven van de wiskundige technieken die in dit huidig werk gebruikt zijn.

D.1

De regel van Leibniz

Om tijdsafgeleiden te integreren, waarbij de integratiegrenzen tijdsfuncties zijn, wordt de regel van Leibniz gebruikt: βZ2 (t)

β1 (t)

D.2

∂ d f (t, z) d z = ∂t dt

βZ2 (t)

f (t, z) d z + f (t, β1 )

∂ ∂ β1 − f (t, β2 ) β2 ∂t ∂t

β1 (t)

Het algoritme van Newton-Raphson [27]

Gegeven is een algebra¨ısch stelsel bestaande uit N vergelijkingen Fi , waaruit de N onbekenden xi moeten opgelost worden: F(x) = 0 i = 1, 2, ..., N (D.2.1) met: 0 F = F1 ... F2 0 x = x1 ... x2

(D.2.2) (D.2.3)

Dit stelsel wordt iteratief opgelost, zodat: xnew = xold + δx

(D.2.4)

δx = −J−1 F

(D.2.5)

met: J de Jacobiaanse matrix, die gedefinieerd wordt als: Jij =

∂Fi ∂xj

i, j = 1, 2, ..., N

(D.2.6)

Bijlage E. Werking simulatieprogramma

101

Bijlage E

Werking simulatieprogramma E.1

Werking steady-state simulatieprogramma

Het steady-state simulatieprogramma wordt opgeroepen via het script main.m. Dit script leest eerst de nodig data in, daarna wordt het sequentiëel algoritme uit hoofdstuk 5 afgelopen. De programma-structuur is weergegeven in figuur E.1. De ingelezen data omvat: • de relaxatieparameters (SOR.dat); • de randvoorwaarden (input.dat); • de initiële waarden (init.dat); • de correcties voor de warmte-overdracht (heat cor.dat); • de grootte van de inactieve zone-lengtes (len cor.dat); • de constructieve details van elke component en de regeling van compressor en klep (read data.m). Deze data kan door de gebruiker aangepast worden. Eenmaal alle nodige gegevens ingelezen zijn, lost het programma 1 voor 1 elke component op. Dit gebeurt door achtereenvolgens de componentfuncties evap.m, comp.m, cond.m en tev.m op te roepen. Hierin wordt het Newton-Raphson algoritme doorlopen. De matrix J en de vector F (zie bijlage D)voor een gegeven iteratiestap worden berekend in de componentspecifieke functie g.m. Deze spreekt daarvoor een andere component-specifieke functie f.m aan. f.m bevat het model van de component. De iteraties worden ook visueel weergegeven. Als een component is opgelost wordt dit gesignaleerd in het command window. Eenmaal elke component is opgelost, wordt de totale lading in het systeem berekend met charge.m. Hieruit volgt de correctie voor de verdamperdruk. Ten slotte worden de nieuwe waarden geplot en vergeleken met de waarden uit de vorige iteratiestap. Indien het verschil voldoende klein is stopt het algoritme, anders begint het een nieuwe berekening.

E.2

Werking dynamisch simulatieprogramma

Het dynamisch simulatieprogramma wordt opgeroepen via de functie main.m met als argumenten het begin- en eindtijdstip van de simulatie. Deze functie leest eerst de nodige data in, daarna wordt het tijdsverloop berekend met de ingebouwde MATLAB-solver ode15s. De programma-structuur is weergegeven in figuur E.2.


Figuur E.1: Programma-structuur steady-state simulatieprogramma

Figuur E.2: Programma-structuur dynamisch simulatieprogramma

102


103

De ingelezen data omvat: • het tijdsverloop van de randvoorwaarden (evap BC.dat en cond BC.dat); • het tijdsverloop van de instelwaarde van de PLR (PLR.dat); • de initiële waarden (init.dat); • de correcties voor de warmte-overdracht (heat cor.dat); • de grootte van de inactieve zone-lengtes (zone.dat); • de constructieve details van elke component en de regeling van compressor en klep (read data.m); • de volgconstanten voor de pseudo-vergelijkingen (tracking.dat). Deze data kan door de gebruiker aangepast worden. Het tijdsverloop wordt gegeven aan de hand van een lijst waarden. Het eerste getal van de lijst waarden duidt de tijdstap aan tussen de waarden. Eenmaal alle nodige gegevens ingelezen zijn, wordt ode15s opgeroepen. Deze berekent M en F aan de hand van de functies give F en give M, met als argumenten het tijdstip t en de toestandsvector y (zie vergelijking (6.2.4)). Om de rekentijd te verkleinen zijn ook de opties ’MvPattern’ en ’JPattern’ gespecifieerd, respectievelijk in de functies JPattern en MvPattern. Na elke geslaagde iteratie wordt de functie handle output opgeroepen. Deze slaat de bekomen data op en geeft ze ook weer op een grafiek. De output van het programma omvat: • prog t: de tijd; • prog x E : het tijdsverloop van de toestand en het massadebiet aan inen uitlaat en de interne grootheden van de verdamper; • prog x C : het tijdsverloop van de toestand en het massadebiet aan inen uitlaat en de interne grootheden van de condensor; • prog x B : het tijdsverloop van de temperatuur in het reservoir; • prog mode E : het tijdsverloop van de optredende mode in de verdamper; • prog mode C : het tijdsverloop van de optredende mode in de condensor; • prog sys: het tijdsverloop van de PLR;

Bibliografie

104

Bibliografie [1] N. Abuaf, O. C. Jones en B. J. C. Wu (1983). Critical flashing flows in nozzles with subcooled inlet conditions. Journal of Heat Transfer-Transactions of the Asme, 105(2):379–383. [2] ASHRAE (2011). Centrifugal compressor PresentationNo2(Centri)version2.pdf.

technology.

www.ashrae.gr/

[3] K. Beatty en D. Katz (1948). Condensation of vapors on outside of finned tubes. Chemical engineering progress, 44(5):55–57. [4] S. Bendapudi, J. E. Braun en E. A. Groll (2008). A comparison of moving-boundary and finite-volume formulations for transients in centrifugal chillers. International Journal of Refrigeration, 31(8):1437–1452. [5] J. E. Bourdouxhe, M. Grodent en J. Lebrun (1996). HVAC1KIT - a toolkit for primary HVAC system energy caclulation. Université de Liège. [6] J. E. Braun, J. W. Mitchell, S. A. Klein en W. A. Beckman (1987). Analysis of the Energy Use and Control Characteristics of a Large Variable Speed Drive Chiller System (RP-409). ASHRAE. [7] L. Bromley (1950). Heat transfer in stable film boiling. Chemical engineering progress, 46(5):231. [8] M. W. Browne en P. K. Bansal (1998). Steady-state model of centrifugal liquid chillers. International Journal of Refrigeration, 21(5):343–358. [9] M. W. Browne en P. K. Bansal (2001). An elemental ntu-epsilon model for vapourcompression liquid chillers. International Journal of Refrigeration-Revue Internationale Du Froid, 24(7):612–627. [10] L. Chen, J. P. Chen, J. H. Liu en Z. J. Chen (2009). Experimental investigation on mass flow characteristics of electronic expansion valves with r22, r410a and r407c. Energy Conversion and Management, 50(4):1033–1039. [11] M. C. Comstock en J. E. Braun (1999). Experimental Data From Fault Detection And Diagnostic Studies On A Centrifugal Chiller (RP-1043). ASHRAE. [12] M. De Paepe (2009). Thermische Installaties. Universiteit Gent. [13] J. T. Gravdahl en O. Egeland (1999). Centrifugal compressor surge and speed control. Control Systems Technology, IEEE Transactions on, 7(5):567–579. [14] F. P. Incropera, D. DeWitt, T. L. Bergman en A. Lavine (2006). Fundamentals of Heat and Mass Transfer. John Wiley & Sons, Inc.

Bibliografie

105

[15] X. Jia, C. P. Tso, P. Jolly en Y. W. Wong (1999). Distributed steady and dynamic modelling of dry-expansion evaporators. International Journal of Refrigeration, 22(2):126– 136. [16] W. Jiang, J. Khan en R. A. Dougal (2006). Dynamic centrifugal compressor model for system simulation. Journal of Power Sources, 158(2):1333–1343. [17] T. S. Lee en W. C. Lu (2010). An evaluation of empirically-based models for predicting energy performance of vapor-compression water chillers. Applied Energy, 87(11):3486– 3493. [18] B. Li en A. G. Alleyne (2010). A dynamic model of a vapor compression cycle with shutdown and start-up operations. International Journal of Refrigeration, 33(3):538–552. [19] N. Liang, S. Q. Shao, C. Q. Tian en Y. Y. Yan (2010). Dynamic simulation of variable capacity refrigeration systems under abnormal conditions. Applied Thermal Engineering, 30(10):1205–1214. [20] J. Liu, J. Chen en Z. Chen (2008). Investigation on the choking flow characteristics in electronic expansion valves. International Journal of Thermal Sciences, 47(5):648–658. [21] MathWorks (2011). Differential equations in matlab - how can differential algebraic equations systems be solved in matlab? http://www.mathworks.com/support/tech-notes/ 1500/1510.html. [22] T. L. McKinley en A. G. Alleyne (2008). An advanced nonlinear switched heat exchanger model for vapor compression cycles using the moving-boundary method. International Journal of Refrigeration, 31(7):1253–1264. [23] P. Mithraratne, N. E. Wijeysundera en T. Y. Bong (2000). Dynamic simulation of a thermostatically controlled counter-flow evaporator. International Journal of Refrigeration, 23(3):174–189. [24] M. J. Moran en H. N. Shapiro (2006). Fundamentals of Engineering Thermodynamics. John Wiley & Sons, Inc. [25] M. Morini, M. Pinelli en M. Venturini (2007). Development of a one-dimensional modular dynamic model for the simulation of surge in compression systems. Journal of Turbomachinery, 129(3):437–447. [26] N. B. O. L. Pettit, M. Willatzen en L. Ploug-Sørensen (1998). A general dynamic simulation model for evaporators and condensers in refrigeration. part ii: simulation and control of an evaporator. International Journal of Refrigeration, 21(5):404–414. [27] W. H. Press, S. A. Teukolsky, W. T. Vetterling en B. P. Flannery (2007). Numerical recipes. Cambridge University Press. [28] S. M. Sami en A. Dahmani (1996). Numerical prediction of dynamic performance of vapour-compression heat pump using new hfc alternatives to hcfc-22. Applied Thermal Engineering, 16(8-9):691–705. [29] V. E. Schrock, E. S. Starkman en R. A. Brown (1977). Flashing flow of initially subcooled water in convergent-divergent nozzles. Journal of Heat Transfer-Transactions of the Asme, 99(2):263–268. [30] L. F. Shampine en M. W. Reichelt (1997). The matlab ode suite. SIAM Journal on Scientific Computing, 18(1):1–22.

Bibliografie

106

[31] J. R. Simoes-Moreira (2000). Oblique evaporation waves. Shock Waves, 10(4):229–234. [32] J. R. Simoes-Moreira en C. W. Bullard (2003). Pressure drop and flashing mechanisms in refrigerant expansion devices. International Journal of Refrigeration, 26(7):840–848. [33] J. R. Simoneau en R. C. Hendricks (1979). Two-phase choked flow of cryogenic fluids in converging-diverging nozzles. Technical report. [34] R. F. Stearns, R. R. Johnson, R. M. Jackson en C. A. Larson (1951). Flow measurement with orifice meters. D. Van Nostrand Company, Toronto. [35] M. A. Tahat, G. A. Ibrahim en S. D. Probert (2001). Performance instability of a refrigerator with its evaporator controlled by a thermostatic expansion-valve. Applied Energy, 70(3):233–249. [36] D. Thirumalai en R. E. White (2000). Steady-state operation of a compressor for a proton exchange membrane fuel cell system. Journal of Applied Electrochemistry, 30(5):551–559. [37] J. R. Thome (2010). Engineering Data Book III. Wolverine Tube, Inc. [38] R. Tirnovan, S. Giurgea, A. Miraoui en M. Cirrincione (2008). Surrogate modelling of compressor characteristics for fuel-cell applications. Applied Energy, 85(5):394–403. [39] J. Vierendeels (2009). Stromingsmechanica. Universiteit Gent. [40] F. Q. Wang, G. G. Maidment, J. F. Missenden en R. M. Tozer (2007). A novel special distributed method for dynamic refrigeration system simulation. International Journal of Refrigeration, 30(5):887–903. [41] G. L. Wedekind, B. L. Bhatt en B. T. Beck (1978). System mean void fraction model for predicting various transient phenomena associated with 2-phase evaporating and condensing flows. International Journal of Multiphase Flow, 4(1):97–114. [42] M. Willatzen, N. B. O. L. Pettit en L. Ploug-Sørensen (1998). A general dynamic simulation model for evaporators and condensers in refrigeration. part i: moving-boundary formulation of two-phase flows with heat exchange. International Journal of Refrigeration, 21(5):398–403. [43] H. Yasuda, S. Touber en C. Machielsen (1983). Simulation-model of a vapor compression refrigeration system. Ashrae Journal-American Society of Heating Refrigerating and AirConditioning Engineers, 25(5):61–61. [44] F. W. Yu en K. T. Chan (2007). Modelling of a condenser-fan control for an air-cooled centrifugal chiller. Applied Energy, 84(11):1117–1135. [45] W. J. Zhang en C. L. Zhang (2006). A generalized moving-boundary model for transient simulation of dry-expansion evaporators under larger disturbances. International Journal of Refrigeration, 29(7):1119–1127. [46] W. J. Zhang, C. L. Zhang en G. L. Ding (2009). On three forms of momentum equation in transient modeling of residential refrigeration systems. International Journal of Refrigeration, 32(5):938–944. [47] W. J. Zhang, C. L. Zhang en G. L. Ding (2009). Transient modeling of an air-cooled chiller with economized compressor. part i: Model development and validation. Applied Thermal Engineering, 29(11-12):2396–2402. [48] L. Zhao en M. Zaheeruddin (2005). Dynamic simulation and analysis of a water chiller refrigeration system. Applied Thermal Engineering, 25(14-15):2258–2271.

Lijst van figuren

107

Lijst van figuren 1.1

Koelmachine . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1

2.1 2.2 2.3 2.4

Basis compressie-koelcyclus [12] . . . . . . Schematische weergave van de proefstand De klep uit de proefstand (naar [11]) . . . Testsequentie . . . . . . . . . . . . . . . .

. . . .

4 5 6 7

3.1 3.2

10

3.3 3.4 3.5 3.6 3.7 3.8

Standaard thermostatische expansieklep . . . . . . . . . . . . . . . . . . . . . . Analogie stroming doorheen een klep en stroming doorheen een convergentdivergent gevolgd door dumpdiffusie [20] . . . . . . . . . . . . . . . . . . . . . . Het stromingspatroon en het drukverloop in een klep [20] . . . . . . . . . . . . Centrifugaalcompressor met inlaatschoepen [2] . . . . . . . . . . . . . . . . . . Compressie-koelcyclus van een koelmachine met inlaatschoepen [44] . . . . . . . Finite-volume methode [23] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Optredende modes in verdamper en condensor volgens [18] . . . . . . . . . . . . Schakelcriterium op basis van enthalpie [26] . . . . . . . . . . . . . . . . . . . .

4.1 4.2 4.3 4.4 4.5 4.6 4.7 4.8 4.9 4.10 4.11 4.12 4.13

Krachtenevenwicht van de klepafsluiter . . . . . . . Schematische weergave van de rotor [5] . . . . . . . Datafit compressor . . . . . . . . . . . . . . . . . . Compressorkarakteristiek . . . . . . . . . . . . . . Infinitesimaal volume dVt in een pijp . . . . . . . . Infinitesimaal volume dVr in het koelmiddel . . . . Infinitesimaal volume dV2 in het secundair flu¨ıdum De verdamper in de TP-V-mode . . . . . . . . . . Lineair verloop enthalpie in de verdamper . . . . . Het schakelcriterium in de verdamper . . . . . . . De condensor in de V-TP-L-mode . . . . . . . . . . Lineair verloop enthalpie in de condensor . . . . . Het schakelcriterium in de condensor . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

24 27 29 30 31 33 33 35 39 40 41 44 45

5.1 5.2 5.3 5.4 5.5 5.6 5.7 5.8 5.9 5.10 5.11

Schematische weergave van het oplossingsalgoritme van het systeemmodel . . Schematische weergave van het oplossingsalgoritme van een componentmodel Schakelcriterium voor verdamper in steady-state model . . . . . . . . . . . . Simulatieresultaten inen uitlaattoestand verdamper met ± 5 % foutgrenzen Simulatieresultaten inen uitlaattoestand condensor met ± 5 % foutgrenzen . Simulatieresultaten systeemparameters met ± 10 % foutgrenzen . . . . . . . . Simulatieresultaten systeemparameters uitgezet ten opzichte van de PLR . . Het verloop van het elektromechanisch rendement in functie van de PLR . . . Simulatieresultaten verdamperen condensordruk in functie van de PLR . . . Invloed PLR op toestandsgrootheden en massadebiet . . . . . . . . . . . . . . Entropietoename tijdens compressie . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . .

50 51 52 54 54 55 55 57 57 59 59

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . . . . . . . . . . .

. . . .

. . . . . . . . . . . . .

. . . .

. . . . . . . . . . . . .

. . . .

. . . . . . . . . . . . .

. . . .

. . . . . . . . . . . . .

. . . .

. . . . . . . . . . . . .

. . . .

. . . . . . . . . . . . .

. . . .

. . . . . . . . . . . . .

. . . .

. . . . . . . . . . . . .

. . . .

. . . . . . . . . . . . .

. . . .

. . . . . . . . . . . . .

. . . .

. . . . . . . . . . . . .

. . . .

. . . . . . . . . . . . .

. . . .

. . . . . . . . . . . . .

. . . .

11 12 13 14 18 19 20

Lijst van figuren 5.12 5.13 5.14 5.15

6.1 6.2 6.3 6.4 6.5 6.6 6.7 6.8 6.9 6.10 6.11

108

Invloed temperatuur koelwater en ijswater op COP . . . . . . . . . . . . . . . . Invloed lading koelmiddel op testgeval 1 . . . . . . . . . . . . . . . . . . . . . . Invloed correctiefactoren op de convectiecoëfficiënt in testgeval 1 . . . . . . . . Verloop zone-groottes in condensor bij dalende druk en constante inlaatenthalpie en massadebiet . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

60 61 62

Simulatieresultaten verdamperen condensordruk . . . . . . . . . . . . . . . . . Simulatieresultaten inen uitlaatenthalpie van verdamper en condensor . . . . Simulatieresultaten uitlaattemperatuur koelen ijswater . . . . . . . . . . . . . Simulatieresultaten massadebiet compressor en expansieklep en tijdsverloop PLR Overeenkomst tussen opgemeten en gesimuleerd transiënt gedrag condensordruk Afwijking tussen opgemeten en gesimuleerd transiënt gedrag massadebiet . . . Invloed te kleine tijdsconstante van het reservoir . . . . . . . . . . . . . . . . . Invloed te grote tijdsconstante van het reservoir . . . . . . . . . . . . . . . . . . Tijdsverloop PLR om schakelen te testen . . . . . . . . . . . . . . . . . . . . . Tijdsverloop zone-groottes en temperaturen in de verdamper tijdens schakelen Tijdsverloop toestandsgrootheden en massadebiet verdamper tijdens schakelen

67 67 68 68 69 70 71 72 72 73 73

63

E.1 Programma-structuur steady-state simulatieprogramma . . . . . . . . . . . . . 102 E.2 Programma-structuur dynamisch simulatieprogramma . . . . . . . . . . . . . . 102

Lijst van tabellen

109

Lijst van tabellen 2.1 2.2 2.3

Eigenschappen verdamper en condensor . . . . . . . . . . . . . . . . . . . . . . Eigenschappen compressor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Testsequentie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5 6 8

4.1 4.2

Model van de klep . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Model van de compressor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

25 28

5.1 5.2

Overeenkomst met de testomstandigheden uit [37] . . . . . . . . . . . . . . . . Correctiefactoren op de convectiecoëfficiënt in iedere zone . . . . . . . . . . . .

58 61

6.1

Inlaatvoorwaarden om schakelen te testen . . . . . . . . . . . . . . . . . . . . .

72

A.1 A.2 A.3 A.4

Nomenclatuur . . Vervolg tabel A.1 Absolute fouten . Vervolg tabel A.3

78 79 80 81

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

. . . .

Dynamische thermodynamische modellering van koelmachines

Recommend Documents