Modelleren en Simuleren van. een Zonnecentrale

FACULTEIT INGENIEURSWETENSCHAPPEN

Vakgroep Elektrische Energie, Systemen en Automatisering Voorzitter: Prof. dr. ir. J. Melkebeek

Modelleren en Simuleren van een Zonnecentrale door Mathieu Vandenbulcke

Promotor: Prof. dr. ir. R. De Keyser Scriptiebegeleiders: Ir. Clara Mihaela Ionescu en Ir. Manuel Galvéz-Carrillo Scriptie ingediend tot het behalen van de academische graad van burgerlijk natuurkundig ingenieur optie toegepaste natuurkunde Academiejaar 2005-2006

Voorwoord In deze thesis tracht ik uit te leggen hoe een simulator van een zonnecentrale geprogrammeerd kan worden. De bedoeling is dat deze simulator gebruikt wordt bij het testen van geavanceerde regelaars die de centrale moeten besturen. Hiermee hoop ik een, zij het bescheiden, bijdrage te leveren tot de ontwikkeling van nieuwe technologieën inzake duurzame energie. Ik zou hierbij ook iedereen willen bedanken die geholpen heeft bij het verwezenlijken van dit eindwerk. In de eerste plaats gaat mijn dank vooral uit naar mijn thesisbeleidster Clara Ionescu, die altijd klaar stond om te antwoorden op al mijn vragen. Verder wil ik haar ook bedanken voor haar hulp bij mijn tussentijdse presentaties. Ook een woord van dank aan mijn promotor Robin De Keyser. Door zijn feedback tijdens de tussentijdse presentaties wist ik mijn simulator telkens te verbeteren. Verder zou ik nog mijn tweede thesisbegeleider Manuel Galvéz-Carrillo willen bedanken die mij, wegens zijn doctoraat over dezelfde zonnecentrale in het Plataforma Solar de Almeria, bijkomende inzichten kon verschaffen. Daarnaast wil ik ook mijn vriendin Inge bedanken, van wie ik veel steun heb gekregen in deze laatste stresserende weken. Ik hoop dat ik voor haar volgend jaar hetzelfde kan doen, tijdens haar thesisjaar. Ook collega Matthias Marescaux zou ik willen bedanken voor het kritisch nalezen van dit eindwerk en het helpen de spreekwoordelijke puntjes op de i te zetten. Dit eindwerk zou er grafisch ook veel minder aantrekkelijk uitgezien hebben mocht het niet i

ii de professionele tekeningen van Jan Moerman bevatten, waarvoor ik hem dan ook wens te bedanken.

Toelating tot bruikleen De auteur geeft de toelating deze scriptie voor consultatie beschikbaar te stellen en delen van de scriptie 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 scriptie.

Mathieu Vandenbulcke, 27 mei 2006.

Overzicht Er valt een toenemende interesse op merken voor milieuvriendelijke energieproductie. Een mogelijke manier om ’groene energie’ te produceren is via een thermische zonnecentrale. Deze verwarmt een thermische olie d.m.v. geconcentreerd zonlicht. De hete olie kan daarna gebruikt worden in verschillende industriële processen waar warmtebronnen nodig zijn, zoals een ontziltingsinstallatie of een elektriciteitscentrale. Deze thesis is als volgt opgebouwd: hoofdstuk 1 bevat een korte inleiding. In hoofdstuk 2 wordt een beschrijving gegeven van de zonnecentrale die gemodelleerd wordt. In hoofdstuk 3 wordt de probleemstelling geformuleerd, er treden immers enkele fysische verschijnselen op die de modellering bemoeilijken. Vanaf hoofdstuk 4 begint de eigenlijke modellering. In dit hoofdstuk wordt de centrale gemodelleerd via niet-uniforme tijdsmonsters. In deze methode wordt gebruik gemaakt van een gedistribueerd parameter model. Daarna worden de resultaten besproken. In hoofdstuk 5 wordt de tweede methode beschreven die gebruikt wordt om de centrale te modelleren. Deze methode werkt via Proper Orthogonal Decomposition (POD). Eerst wordt een kleine theoretische uiteenzetting gegeven over de POD-methode en vervolgens wordt de methode toegepast op de centrale. Daarna worden opnieuw de resultaten besproken. In hoofdstuk 6 wordt tot slot een korte besluitvorming geformuleerd.

iii

Extended Abstract

iv

v

Modelling and simulating of a solar power plant Mathieu Vandenbulcke Supervisor(s): Prof. dr. ir. Robin De Keyser, ir. Clara Ionescu, ir. Manuel Galvéz-Carrillo Abstract— This article explains how a simulator is constructed for the solar power plant located at Plataforma Solar de Alemeria. Two techniques are explained. The first one using non-uniform sampling time and the second one using Proper Orthogonal Decomposition (POD). Keywords—Solar power plant, non-uniform sampling time, POD.

I. I NTRODUCTION

I

T is obvious that fossile fuels won’t be available forever and that there is an increasing need for renewable energy sources. A source that has been and will be around for quite some time is the sun. Solar irradiation can be converted to useful energy in different ways. The method described in this thesis uses parabolic mirrors to concentrate sunlight on a system of parallel pipes. The thermal oil that flows through these pipes gets heated exstensively and is fed to the top of a storage tank. The heat that is present in the oil can be used in different industrial processes where there’s a need for heat sources. E.g. in an electrical power plant where steam has to be generated or in a desalination plant. II. P LANT DESCRIPTION The plant is located at Plataforma Solar de Almeria, in the South-East of Spain and can basically be devided into three large sections. The first section feeds oil from the bottom of the storage tank to the system of parallel loops and is called the ’feeding section’. The second part of the plant consists of ten parallel loops that collect the sunlight and heat the oil. This section will be called the ’heating section’. The last section returns the hot oil from the ’heating section’ to the top of the storage tank and is called the ’return section’.

from the sunlight as possible and to maintain the temperature of the oil that flows into the storage tank at a desired level. In order to do this, there is only one control variable: the flow of the oil. The influence of the flow on the temperature of the oil is simple: a higher flow will lead to a lower temperature, because the oil has got less time to get heated by the sunlight. III. N ON - UNIFORM TIME SAMPLES The first methode uses non-uniform time samples [1] to simulate the plant. This is inspired by the fact that the flow of the oil can vary substantially during operation. Equation 1 describes the heat received and stored into the metal tubes, while equation 2 describes the amount of heat received by the oil:

ρm Cm Am

∂Tm ∂t

= DIrr(t) − Hl G [Tm − Ta ] −LHt [Tm − Tf ]

∂Tf ∂Tf ρf Cf Af + ρf Cf q(t) ∂t ∂x

= LHt [Tm − Tf ]

These two equations only hold for the heating section. The feeding and return section can be regarded as a pure time delay: these sections are thoroughly insulated and heat losses can be neglected.

Fig. 2. Results for the temperature at the end of the returnsection for the 15th of July. Fig. 1. Schematical drawing of the plant.

The oil flow is delivered by a pump situated in the feeding section. The objective of the plant is to extract as much energy M. Vandenbulcke is a student of Ghent University (UGent), Gent, Belgium. E-mail: mathieu [email protected].

The results obtained with this first simulator are quite good. The error normally stays below 5%, except during the first moments and when the flow is rather low. The first kind of errors are due to the fact that the initial temperature distribution along the system of parallel pipes is unknown. The second error is probably caused by an inaccurate heat transfer coefficient Ht .

(1) (2)

vi IV. POD- METHOD The second simulator is constructed using the POD-method [2], [3]. This method produces a number of orthonormal base functions φi (x) which represent the spatial distribution of an important parameter in the system, in this case the temperature. These base functions can be derived from measured or calculated data. To describe the variations in time, each base function is multiplied by its own time-varying coefficient ai (t). In this way the time- and place dependant temperature of the system can be n P ai (t)φi (x). The advantadescribed by the approximation : i=1

ge of using the POD-method is the fact that only a low number of base functions φi (x) needs to be considered to describe the system. To find the base functions φi (x), data concerning the spatial distribution of temperature along the system, is needed. Since this data is not available, it needs to be calculated. In order to do is, a state-space model of the system is constructed [4]. The state-space model is derived from a first order partial differential equation which describes the conservation of energy in the plant: ρf cf

DT (x, t) = S(x, t) Dt

Fig. 3. Results for the temperature at the end of the returnsection for the 15th of July using POD.

(3)

D ∂ ∂ where Dt = ∂t + vx (t) ∂x is a substantial derivative and S(x, t) denotes a time and place depandent heat source, in this case the solar irradiation. This can be cast into a state space model [5]:

T(t) = A−1 (t)T(t − 1) + A−1 (t)u(t)

(4)

T(t) refers to the temperature distribution of the system at time t and u(t) is a vector with input signals. Substituting T by n P its decomposition ai (t)φi (x) = Φn a(t) leads to:

Fig. 4. Results for the temperature at the end of the returnsection for the 16th of July using POD.

i=1

V. C ONCLUSIONS Φn a(t) = A−1 (t)Φn a(t − 1) + A−1 (t)u(t)

(5)

Projection of equation 5 on the orthonormal set of base functions Φn results in:

a(t) = Φn T A−1 (t)Φn a(t − 1) + Φn T A−1 (t)u(t)

(6)

Equation 6 gives a simple expression for the time varying coefficients ai (t). The result of this method can be seen in figures 3 and 4. The results for the 15th of July look very good, even for low values of the flow. The results for the 16th of July on the other hand are less accurate. Only in the second half of the experiment there is a match between the real and simulated data. A possible explanation for this could be the lack of a flow dependant heat transfer coefficient Ht in equation 3, which was present in the first simulator.

This article has presented two techniques for obtaining a model of the solar power plant located at Plataforma Solar de Alermia. For the first technique some improvements could be made regarding the heat transfer coefficient Ht . The second technique, using the POD-method, is substantially faster than the first one. For the second technique, the use of a flow dependant heat transfer coefficient should lead to better results. R EFERENCES [1] R. N. Silva, L. M. Rato, and J.M. Lemos. Observer based time warped control of distributed collector solar fields. Improving human potential programma, proceedings of the 2nd users workshop, pages 103–110, February 2002. [2] Anindya Chatterjee. An introduction to proper orthogonal decomposition, December 1999. [3] Patricia Astrid. Reduction of Process Simulation Models: a proper orthogonal decomposition approach. PhD thesis, Technische Universiteit Eindhoven, November 2004. [4] Oscar Mauricio Agudelo Manozca. Control of the temperature profile of an one-dimensional bar using a pod approach. presentation, January 2006. [5] Yang, Cao, Chung, and Morris. Applied Numerical Methods Using Matlab. John Wiley & Sons Inc, 2005.

Inhoudsopgave . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

i

Overzicht . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

iii

Extended Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

iv

Voorwoord

1 Inleiding

1

2 Beschrijving van de centrale

7

3 Probleemstelling

14

4 Modelleren via niet-uniforme tijdsmonsters

16

4.1

4.2

4.3

4.4

Vergelijking . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

16

4.1.1

Eigenschappen van de olie . . . . . . . . . . . . . . . . . . . . . . . . .

18

4.1.2

Bepalen van Ht . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

18

4.1.3

Onderstellingen . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

19

Verdelen in secties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

20

4.2.1

Toevoersectie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

20

4.2.2

Verhittingssectie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

22

4.2.3

Afvoersectie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

24

Niet-uniforme tijdsmonsters . . . . . . . . . . . . . . . . . . . . . . . . . . . .

26

4.3.1

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

26

Implementatie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

30

4.4.1

Toevoersectie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

33

4.4.2

Verhittingssectie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

33

4.4.3

Afvoersectie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

36

Berekening van tk

vii

Inhoudsopgave 4.4.4 4.5

4.6

viii Beginwaarden . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

37

Resultaten . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

38

4.5.1

15 juli: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

39

4.5.2

16 juli: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

42

4.5.3

28 juli: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

43

4.5.4

Conclusies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

45

Experimentele bepaling van de lengte van de afvoersectie. . . . . . . . . . . .

46

5 Methode II: Proper Orthogonal Decomposition 5.1

5.2

5.3

52

Theorie . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

53

5.1.1

Bepalen van de basisfuncties . . . . . . . . . . . . . . . . . . . . . . .

54

5.1.2

Bepalen van de tijdsafhankelijke coëfficiënten . . . . . . . . . . . . . .

57

Praktijk . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

59

5.2.1

Vergelijking voor het temperatuursprofiel . . . . . . . . . . . . . . . .

60

5.2.2

Opstellen van het toestandsmodel . . . . . . . . . . . . . . . . . . . .

62

5.2.3

Berekenen van de basisfuncties . . . . . . . . . . . . . . . . . . . . . .

64

5.2.4

Bepalen van de tijdsafhankelijk coëfficiënten . . . . . . . . . . . . . . .

66

5.2.5

Uitbreiding naar de volledige centrale . . . . . . . . . . . . . . . . . .

67

Resultaten . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

68

5.3.1

15 juli . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

68

5.3.2

16 juli . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

70

5.3.3

28 juli . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

71

5.3.4

Conclusies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

72

6 Besluit

73

A Dode tijden

75

B Wiskundige bijlage

78

B.1 Definities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

78

B.2 Oplossing extremumprobleem . . . . . . . . . . . . . . . . . . . . . . . . . . .

79

Inhoudsopgave C Figuren

ix 81

C.1 Figuren bij hoofdstuk 4 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

81

C.2 Figuren bij hoofdstuk 5 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

92

D Matlabcode

100

Bibliografie

121

Hoofdstuk 1

Inleiding Het tijdperk waarin onze energievoorzieningen tegemoet gekomen worden door hoofdzakelijk fossiele brandstoffen lijkt genoodzaakt te eindigen. Hiervoor vallen twee redenen te onderscheiden: vooreerst is de voorraad fossiele brandstoffen eindig. De tijden waarin per jaar meer olie werd gevonden dan geconsumeerd zijn over. Als het huidige verbruikerspeil aanhoudt, zal er voldoende olie zijn tot in 2050. Voor aardgas zijn de vooruitzichten iets beter, geschat wordt dat er nog 75 jaar aardgas beschikbaar zal zijn. Steenkool is nog het meest voorhanden. Daar wordt gesproken over voorraden die minstens 200 jaar aanspreekbaar zijn. Een tweede oorzaak is de milieu-impact. Begrippen als ’het broeikaseffect’, ’C02 -uitstoot’ en ’zure regen’ worden allemaal verbonden met het verbruik van fossiele brandstoffen en zijn helemaal niet wenselijk. Naast de conventionele energiecentrales bestaan er ook nog kerncentrales, die anno 2006 goed zijn voor 55% van stroomvoorziening in België. Het controversiële karakter van kernenergie leidt ertoe dat het bouwen van nieuwe kerncentrales op veel weerstand stuit. Nieuwe, alternatieve en vooral duurzame energiebronnen lijken dus een noodzaak. Veel wordt verwacht van kernfusie. Kernfusie heeft enkele ontegensprekelijke voordelen tegenover fossiele brandstoffen en kernfissie 1 . Zo worden geen milieu-onvriendelijke gassen geproduceerd, zoals CO2 . Een kernfusiereactor brengt geen grote gevaren met zich mee, wat niet kan gezegd worden van kernfissiecentrales. Zelfs bij zware ongevallen is er weinig kans 1

kernfissie of kernsplitsing, het fysisch proces in de huidige kerncentrales, waarbij energie wordt vrijgemaakt

door het splitsen van zware atomen

1

Hoofdstuk 1. Inleiding

2

op ernstige schade. Daarenboven geldt: de grondstoffen voor een fusiereactor zijn praktisch overal ter wereld en in overvloed beschikbaar. Zo kan deuterium (waterstof, maar met één neutron:

2 H)

gevonden worden in water en lithium (6 Li) in gesteenten. Dit zorgt ervoor

dat geografische conflicten, ontstaan door het feit dat sommige landen bijna alle grondstoffen bezitten om energie te produceren en andere niet, in de toekomst bijna niet meer zullen voorkomen. Maar dit is echter toekomstmuziek, gezien commerciële kernfusiereactoren ten vroegste over 40 jaar beschikbaar zullen zijn. En dit cijfer is nog vrij optimistisch. Om de periode tussen vandaag en het kernfusietijdperk te overbruggen is er zeker nog nood aan alternatieve energiebronnen. Hiernaar werd in het verleden reeds veel onderzoek verricht. Er kunnen drie grote groepen onderscheiden worden: zonne-energie, windenergie (eigenlijk ook een gevolg van het zonlicht) en hydro-energie. De laatste groep doelt vooral op stuwdammen. In dit geval is de energie eigenlijk ook afkomstig van het zonlicht, denk maar aan de watercyclus. Deze stuwdammen hebben als nadeel dat ze een grote oppervlakte in beslag nemen. Vele valleien werden reeds getransformeerd tot meren om de turbines in de stuwdammen te laten draaien. Het principe van een hydro-centrale is vrij duidelijk: de potentiële energie van water op een hoger punt wordt omgezet naar kinetische energie door het naar een lager punt te laten lopen. De kinetische energie die hierbij ontstaat wordt gebruikt om een turbine te laten draaien en zodus elektriciteit te produceren. Ook de wind kan omgezet worden in elektrisch energie. Hiervoor worden hoge windturbines gebruikt. Deze manier van energie-opwekking krijgt de laatste tijd veel belangstelling. Denk maar aan de vele windmolenparken die gebouwd worden2 . Het ge¨ınstalleerd vermogen aan windenergie bedroeg eind 2005 167 MW in België (bron: the European Wind Energy Association, www.ewea.org). De laatste groep betreft de energie geproduceerd uit zonlicht. Zonlicht kan enerzijds meteen omgezet worden in elektriciteit. Dit is het geval bij fotovolta¨ısche zonnecellen. Deze cellen worden hoofdzakelijk gemaakt uit halfgeleidermaterialen zoals silicium en germanium. Voorbeelden hiervan zijn cellen met kristallijn silicium. Het rendement van deze cellen blijkt te stagneren rond de 16 à 17 %. Dit betekent dat slechts 16 à 17% van het invallende zonlicht kan omgezet worden in elektriciteit. Daarnaast heb je ook de dunne film cellen. Deze bestaan 2

Kluizendok in haven van Gent, Villers-le-Bouillet en in de toekomst op de fameuze Thorntonbank


3

uit amorf silicium dat gedeponeerd wordt op een substraat. Hun voordeel is dat de productie veel goedkoper is, gezien het bestaan van goedkope depositie technieken. Het rendenment van deze cellen ligt echter lager dan die van de kristallijne cellen en bedraagt 6-10%. Een recentere vorm van fotovolta¨ısche cellen maakt geen gebruik meer van halfgeleiders maar van polymeren. Deze polymere cellen bestaan uit kunststoffen die zich gedragen als halfgeleiders. Aangezien deze technologie nog in zijn kinderschoenen staat, liggen de rendementen ook eerder laag: 2-3%. Het voordeel van deze technologie is dat de energie gegenereerd in deze fotovolta¨ısche cellen rechtstreeks gebruikt kan worden of gestockeerd kan worden in batterijen. Anderzijds kan het zonlicht in een eerste stadium aangewend worden als warmtebron. In dit geval kan de geproduceerde warmte gebruikt worden om, via de creatie van stoom, een generator aan te drijven en zo elektriciteit te produceren. De warmte kan ook opgeslagen worden in een materiaal dat het een tijdlang kan vasthouden. Zo kan de energie ook aangewend worden als er even geen zon is. Dit is het geval voor de zonnecentrale die besproken wordt in dit eindwerk. In dergelijke zonnecentrales wordt olie opgewarmd via zonlicht en gestockeerd in een grote opslagtank. De energie die vervat zit in de warmte van de olie kan later bijvoorbeeld aangewend worden voor het creëren van stoom. Dit eindwerk past in het kader van het onderzoek over dergelijke zonnecentrales. Er wordt getracht de energie-extractie uit het zonlicht zo optimaal mogelijk te laten verlopen. Verder wordt er ook geëist dat de temperatuur van de olie gecontroleerd kan worden, om o.a. oververhitting te voorkomen. Het controleren van de olietemperatuur is echter niet evident: de olie warmt immers sneller op onder een wolkeloze hemel dan onder een bewolkte hemel. Om toch een vooropgestelde temperatuur te krijgen, moet de olie rondgepompt worden met een variërend debiet: sneller als er veel zon is, trager als er weinig zon is. Dit vereist een geavanceerd regelsysteem en om dit regelsysteem te testen is er nood aan een computermodel dat de echte centrale simuleert. Het genereren van een dergelijk computermodel is het eigenlijke onderwerp van deze thesis. Om dit computermodel te schrijven, is er nood aan een wiskundig model dat het fysisch proces beschrijft. Dit wiskundig model kan op drie verschillende manieren bekomen worden. Een eerste manier is via een theoretische analyse. Deze methode vertrekt van de drie bekende behoudswetten: behoud van massa, behoud van impuls en behoud van energie. Deze drie


4

wetten worden toegepast op het proces wat leidt tot een aantal differentiaalvergelijkingen. Het voordeel van deze methode is dat het wiskundige model kan gevonden worden, nog voor het fysisch systeem in werkelijkheid bestaat. Zo kunnen reeds simulaties uitgevoerd worden van een centrale, nog voor ze gebouwd is. Een tweede voordeel van een theoretische analyse is dat het zogenaamde ’white models’ betreft. De term ’white’ slaat op het feit dat alle parameters in de vergelijking een duidelijke fysische betekenis hebben. De bekomen vergelijking schept een duidelijk fysisch inzicht in het proces. Een nadeel van deze methode is dan ook dat er heel wat fysisch inzicht nodig is om het wiskundig model op te stellen. Zo moet door de ingenieur besloten worden welke de belangrijkste fysische verschijnselen zijn. Er kan geen rekening gehouden worden met elk fysisch proces dat optreedt. Dit zou immers tot een te complex model leiden. Er moet dus ergens een compromis gevonden worden tussen nauwkeurigheid en efficiëntie. Een tweede manier is via experimentele analyse. In deze analyse wordt gebruik gemaakt van de informatie die in het in- en uitgangssignaal van het proces zit. Deze informatie, meestal onder de vorm van een Bodeplot, wordt opgemeten en onderworpen aan een identificatieprocedure, wat leidt tot een wiskundig model van het systeem. Dit identificatieproces kan op twee verschillende manieren gebeuren. Heeft men reeds een idee over de vorm van de vergelijking, dan kan een parametrische identificatie toegepast worden. Enkele populaire voorbeelden van parametrische identificatie algoritmen zijn: de least squares methode, instrumental variables methode en de prediction error methode. Hierbij wordt het model vooropgesteld en worden de waarden van de parameters gevonden via identificatie. Is de structuur van het wiskundig model daarentegen niet op voorhand beschikbaar, dan dient men over te stappen op niet-parametrische identificatie. Enkele voorbeelden van niet-parametrische identificatie algoritmen zijn: transiënt analyse, frequentie analyse, correlatie analyse en spectraalanalyse. Het grote voordeel van een experimentele analyse is dat er eigenlijk niets hoeft gekend te zijn over het fysisch proces. Leg een ingang aan aan het systeem, registreer wat eruit komt en er kan een model opgesteld worden. Wat er zich in het systeem afspeelt, hoeft niet perse gekend te zijn. Vandaar dat de verzamelnaam ’black models’ hiervoor frequent gebruikt wordt. Identificatie procedures kunnen niet enkel gebruikt worden voor het opstellen van een volledig model, maar ook voor het schatten van ongekende parameters in een model. Het nadeel is dat er tijdrovende tests moeten uitgevoerd worden.


5

Een derde manier om een wiskundige vergelijking voor een fysisch proces op te stellen maakt gebruik van een combinatie van de vorige twee methodes, theoretische en experimentele analyse. Hier kan bijvoorbeeld een wiskundig model vooropgesteld worden op basis van een theoretische analyse, maar waarvan geen uitdrukking gevonden kan worden voor alle fysische parameters. Een experimentele analyse kan dan aangewend worden om de missende parameters te schatten. Gezien de modellen gevonden via deze methode gebruik maken van zowel de ’white model’ als van de ’black model’ filosofie, worden ze om begrijpelijke redenen ook wel ’grey models’ genoemd. Zoals reeds gezegd wordt in dit eindwerk een model opgesteld van een zonnecentrale. Gezien de fysische processen die zich hier in afspelen goed gekend zijn, wordt begonnen met een theoretische analyse. Er zijn echter enkele parameters waaraan niet onmiddellijk een waarde kan toegekend worden. Om deze te schatten wordt gebruik gemaakt van een experimentele analyse. Vandaar dat het model opgesteld in deze thesis een ’grey model’ is. In de literatuur werden reeds heel wat modellen voor de zonnecentrale voorgesteld. Over het algemeen kunnen al deze modellen bestempeld worden als een ’black model’ of als een ’grey model’. De ’black models’ kunnen verder nog opgedeeld worden in zij die bekomen werden via parametrische en zij die bekomen werden via niet-parametrische identificatie. Bij degene bekomen via parametrische identificatie is reeds één en ander gekend over de structuur van het wiskundig model. Sommige van deze voorgestelde modellen beschrijven de fysische werking van de centrale met slechts één enkele vergelijking. Deze worden ook geconcentreerde parameter modellen genoemd. Voorbeelden hiervan zijn te vinden in [1]. Andere modellen delen de centrale op in verschillende subsecties en definiëren voor elke subsectie een aparte vergelijking. Deze modellen worden gedistribueerde parameter modellen genoemd. Voorbeelden hiervan zijn te vinden in [1],[2],[3] en [4]. In beide gevallen komt het erop neer de waarde van de parameters te vinden via identificatie methodes. Bij de modellen bekomen via niet-parametrische identificatie gaat men ervan uit dat niets gekend is over de werking van het systeem. Hier worden enkel de in- en uitgang van het systeem bekeken. Een niet parametrische identificatietheorie, die men voor de zonnecentrale placht te gebruiken, is die van de artificiële neurale netwerken (ANN). Hier wordt een verband gelegd tussen in- en uitgang, ge¨ınspireerd op de werking van biologische neurale netwerken.


6

Net zoals mensen en dieren bepaalde gedragspatronen moeten aanleren, moet ook het artificiële neurale netwerk dit doen. Het programma moet met andere woorden een leerproces ondergaan. Hiervoor wordt beroep gedaan op de data gegenereerd door het werkelijke proces. Voorbeelden zijn te vinden in [5] en [6]. Een andere methode van niet-parametrische identificatie betreft de modellen opgesteld via zogenaamde ’fuzzy logic’. Zoals de naam laat vermoeden wordt er in deze algoritmen geen gebruik gemaakt van strikte binaire waarden zoals ’goed-fout’ en ’1-0’, maar van meer ’troebele’ of ’fuzzy’ waarheden zoals ’relatief goed’, ’normaal’, ’eerder slecht’, enz. Voorbeelden hiervan zijn te vinden in [7]. Een derde methode van niet-parametrische identificatie heet Proper Orthogonal Decomposition of kortweg POD. Deze methode werkt aan de hand van orthogonale basisfuncties, die worden afgeleid uit gemeten of berekende data. Deze basisfuncties geven de ruimtelijke distributie van een bepaalde fysische grootheid, zoals bijvoorbeeld de temperatuur, doorheen het proces, weer. Om de tijdsafhankelijkheid van de fysische grootheid te berekenen, wordt het oorspronkelijke fysische proces, over het algemeen beschreven door een differentiaalvergelijking, geprojecteerd op deze basisfuncties volgens de Galerkin projectie. Uit deze projectie volgen dan een aantal tijdsafhankelijke coëfficiënten die de variatie van de basisfuncties in de tijd weergeven. Deze methode is zeer gelijkend op Fourieranalyse, met dat verschil dat de hierbesproken basisfuncties over het algemeen geen sinussen of cosinussen zijn. Het grote voordeel van POD is dat het aantal basisfuncties meestal heel beperkt is (orde 10), waardoor de simulatie weinig berekeningen vergt [8],[9],[10],[11]. Verder kan er ook nog een onderscheid gemaakt worden tussen algoritmen die gebruik maken van uniforme tijdsmonsters en algoritmen die gebruik maken van niet-uniforme tijdsmonsters. Het feit dat het debiet varieert in de tijd heeft sommigen ertoe aangezet niet-uniforme tijdsmonsters te gebruiken. Het toepassen ervan leidt tot een zeer eenvoudige implementatie van het variërend debiet. Een voorbeeld van het gebruik van niet-uniforme tijdsmonsters is te vinden in [12].

Hoofdstuk 2

Beschrijving van de centrale De zonnecentrale bevindt zich op het Plataforma Solar de Almer´ıa (PSA1 ), het Europees centrum voor zonne-energie applicaties in het zuidoosten van Spanje. Op deze site bevinden zich meerdere centrales (zie figuur 2.1), waarvan het ACUREX-veld onderwerp vormt van dit eindwerk. Dit veld is een zogenaamd gedistribueerd collector systeem (distributed collector system of DCS). Dit betekent dat de zonne-energie op verschillende plaatsen wordt omgezet in thermische energie. Daarnaast bestaan ook systemen waarbij de omzetting op slechts één plaats gebeurt. Dergelijke systemen noemt men centrale ontvanger systemen (central receiver system of CRS).

Figuur 2.1: Overzicht van PSA.

1

Meer info over PSA op www.psa.es

7

Hoofdstuk 2. Beschrijving van de centrale

8

De bedoeling van de zonnecentrale is olie op te warmen door middel van zonlicht. De verwarmde olie wordt gestockeerd in een grote opslagtank en kan gebruikt worden in verschillende industriële processen, waar veel gebruik gemaakt wordt van warmtebronnen. Twee veel geciteerde voorbeelden zijn: het creëren van stoom om elektriciteit op te wekken en energietoevoer in een ontziltingsinstallatie. Om het zonlicht te collecteren wordt op het ACUREX-veld gebruik gemaakt van parabolische spiegels (figuur 2.2). Hun parabolisch profiel zorgt ervoor dat het invallende licht geconcentreerd wordt in het brandpunt. Daar de zon overdag van oost naar west beweegt, hebben de spiegels één rotationele vrijheidsgraad om tijdens de operationele fase steeds naar de zon gericht te zijn (zie figuur 2.3). Een eenvoudig systeem van twee fotodiodes met daartussen een klein plaatje zorgt hiervoor, zoals te zien is in figuur 2.4. Is de spiegel niet volledig naar de zon gericht dan werpt het plaatje zijn schaduw op een van de fotodiodes, die daardoor minder spanning zal produceren. Hierdoor genereren de diodes een verschillende spanning en met dit spanningsverschil kan de positie gecorrigeerd worden. Deze werkingsfase wordt de ’Track toestand’ genoemd. Om oververhitting te vermijden is het ook mogelijk de spiegels enkele graden van de zon weg te richten (Desteer toestand). Op het eind van elke dag of bij serieuze defecten wordt de ’Stow toestand’ ingenomen. Hierbij worden de spiegels volledig naar de grond gericht.

Figuur 2.2: Parabolische spiegels zorgen voor concentratie van het zonlicht.


9

Figuur 2.3: De spiegels volgen de zon bij haar beweging van oost naar west.

Figuur 2.4: Fotodiodes met daartussen een plaatje.

In het brandpunt van de parabolische spiegels loopt, evenwijdig met de as van de spiegels, een metalen pijpleiding omhuld door een doorzichtige, glazen buis. Deze glazen buis isoleert de metalen leiding van de wind. Daar het glas doorzichtig is wordt het zonlicht nauwelijks belemmerd. Er zijn in totaal 20 rijen spiegels, waarvan er in figuur 2.2 één te zien is. Deze 20 rijen vormen samen 10 parallelle lussen of ’loops’ met elk een lengte van 172 m. Elke van deze loops bevat 48 parabolische spiegels en elke spiegel heeft een nuttige oppervlakte of apertuur van 1, 83 × 3, 05 m2 . Het totale collectoroppervlak bedraagt bijgevolg 10 × 48 × 3, 05 m × 1, 83 m ≈ 2680 m2 . De buitendiameter van de metalen pijpleiding is 3,18 cm, wat leidt tot een concentratiefactor van het zonlicht van

1,83m 0,0318m

≈ 57, 5. Deze concentratiefactor zorgt


10

ervoor dat de metalen pijpleiding in het brandpunt van de spiegels sterk verhit wordt. Door de leiding stroomt een olie die wordt opgewarmd door de verhitte pijpleiding. Aan het begin van elke loop bevindt zich een klep die de loop kan afsluiten in geval van een noodsituatie zoals een lek. De olie wordt voortgestuwd door één pomp (zie figuur 2.5), waarvan het debiet geregeld kan worden. Om veiligheidsredenen is dit debiet beperkt tussen 2 l/s en 10 l/s. In figuur 2.5 is er ook een buffertank en een driewegenklep te zien. Deze laatste kan aangesproken worden om de olie enkele malen door het collectorsysteem te laten lopen, tot ze een voldoend hoge temperatuur bereikt heeft om de opslagtank te betreden.

Figuur 2.5: Schematisch overzicht van de zonnecentrale.

Het volledige buizensysteem kan onderverdeeld worden in drie secties (figuur 2.6). In de eerste sectie wordt olie vanuit de bodem van de opslagtank naar het collectorveld gepompt. Deze leiding heet de toevoersectie en is ge¨ısoleerd om het warmteverlies te beperken. In de tweede sectie wordt de leiding gesplitst in tien parallelle leidingen of loops. Elke loop is, zoals gezegd, voorzien van collectoren die de olie opwarmen. Dit deel wordt de verhittingssectie genoemd. Om de warme olie terug naar de top van de opslagtank te brengen wordt gebruik gemaakt van een derde leiding, de afvoersectie, die ook ge¨ısoleerd is. De centrale heeft een thermisch piekvermogen van zo’n 1, 3 M Wt bij invallend zonlicht van 900 W/m2 . Dagelijks produceert de centrale gemiddeld een 6, 5 M W ht aan thermische energie. Het volume van de opslagtank (zie figuur 2.7) bedraagt 140 m3 waardoor er bij een olietemperatuur van zo’n


11

290°C ongeveer 2, 3 M W h kan in opgeslagen worden. Door de stratificatie-eigenschappen2 van de opslagtank en het thermoclinisch effect3 van de olie kan de warme temperatuur van de bovenste olielagen verschillende dagen behouden worden. De onderverdeling in drie secties wordt meer gedetailleerd besproken in hoofdstuk 4.

Figuur 2.6: Onderverdeling in toevoer-, verhittings- en afvoersectie.

2 3

Stratificatie eigenschappen: olielagen op verschillende temperatuur gaan niet met elkaar mixen. Thermoclinisch effect: de temperatuur van de olie in de opslagtank stijgt sterk van onder naar boven toe.


12

Figuur 2.7: Foto van de opslagtank.

Het data acquisitie systeem bestaat uit een centrale computereenheid die alle data vergaart. In totaal wordt rekening gehouden met meer dan 150 variabelen. De belangrijkste hiervan zijn de in- en uitgangstemperatuur van het buizensysteem, het debiet, de intensiteit van het zonlicht, de temperatuur op het eind van elke loop4 , de omgevingstemperatuur en de tijd. In het model kan niet met al deze variabelen rekening gehouden worden, dit zou immers teveel rekenwerk vragen. De variabelen die gebruikt worden door het model zijn onderverdeeld in: • Input: het debiet van de olie in m3 /s • Output: de temperaturen op het eind van elke loop in de verhittingssectie en de temperatuur op het eind van de afvoersectie in °C • Storingen: warmteflux van de zon in W /m2 , de temperatuur van de olie die uit de opslagtank vloeit (= olietemperatuur aan het begin van de toevoersectie) in °C en de omgevingstemperatuur in °C. 4

Op het eind van elke loop is een temperatuursensor aangebracht.


Figuur 2.8: Model

13

Hoofdstuk 3

Probleemstelling Het grote nadeel van een zonnecentrale als deze, is de onvoorspelbaarheid van het weer. Het is eenvoudig in te zien dat onder een wolkeloze hemel de olie veel sneller zal opwarmen dan onder een zwaar bewolkte hemel. Denk maar aan de plotselinge warmte die je voelt, als de zon opeens vanachter een wolk verschijnt. Daarnaast ondergaat het zonlicht ook meer regelmatige veranderingen, zoals de dag-nacht en zomer-winter cyclus. Deze kunnen, in tegenstelling tot de vorige, exact berekend worden [3]. Niet alleen het zonlicht bevat een onvoorspelbare factor, ook de omgevingstemperatuur kan een invloed hebben. De warmteverliezen worden beperkt door de aangebrachte isolatie, maar volledig isoleren kan nooit, waardoor de warmte van de olie gedeeltelijk weglekt naar de omgeving. Een derde belangrijke factor is de ingangstemperatuur van het systeem. Hoe lager die ingangstemperatuur, hoe meer er moet opgewarmd worden om eenzelfde eindtemperatuur te krijgen. Er is immers een veel grotere hoeveelheid energie nodig om olie van bijvoorbeeld 60°C op te warmen tot 300°C dan olie van 120°C. De sterke invloed van de ingangstemperatuur kan men zien in figuur C.2 in bijlage C. Om de invloed van al deze fluctuaties te controleren en zo de uitgangstemperatuur van de olie op een gewenste waarde te houden, is er slechts één input: het debiet. De invloed van het debiet kan als volgt begrepen worden: onderstel dat het niet bewolkt is en dat de olie wordt rondgepompt met een welbepaald debiet. Op deze manier wordt de temperatuur van de olie tussen in- en uitgang van het systeem met een bepaald aantal graden verhoogd. Als er plots een wolk voor de zon verschijnt, zal de opwarming van de olie bij eenzelfde debiet fors 14

Hoofdstuk 3. Probleemstelling

15

dalen. Om dezelfde temperatuurstijging te krijgen, moet de olie trager rondgepompt worden. Zo heeft de olie meer tijd om op te warmen. Het nadeel van deze debietvariaties is dat ze een tijdsafhankelijke dode tijd (het tijdsverschil tussen het ogenblik dat de olie het systeem betreedt en het ogenblik dat dezelfde olie het systeem terug verlaat) induceren:

V = τ (t) q(t)

(3.1)

• V : volume van de olie in het systeem in m3 • q(t): debiet van de olie in m3 /s • τ (t): tijdsafhankelijke dode tijd in s Het is duidelijk dat er een geavanceerd regelsysteem nodig is om met al deze factoren rekening te houden. Dit is het onderwerp van de doctoraatsstudie van ir. Manuel Galvéz-Carrillo. In het eerste deel van zijn proefschrift wordt gezocht naar een degelijk model voor predictive control. In het tweede deel wordt dit model gebruikt in een predictive controller ontworpen via de EPSAC1 of NEPSAC (alnaargelang het gevonden model lineair dan wel niet-lineair is) om een regelaar te ontwerpen. Om de werking van deze regelaar te kunnen testen is er nood aan een simulatiemodel van de zonnecentrale. De constructie van deze simulator is het eigenlijke onderwerp van dit eindwerk.

1

Extended Prediction Self-Adaptive Control

Hoofdstuk 4

Modelleren via niet-uniforme tijdsmonsters In de eerste methode wordt gebruik gemaakt van niet-uniforme tijdsmonsters om op een eenvoudige manier het variërend debiet te modelleren. Zoals in de inleiding reeds vermeld werd, wordt er in dit eindwerk uitgegaan van een gedistribueerd parameter model. Dit betekent dat het volledige buizennetwerk wordt opgedeeld in een groot aantal subsystemen (cfr discretisatie). Voor elk subsysteem wordt vervolgens de bijhorende vergelijking opgelost.

4.1

Vergelijking

Als men het systeem in meer detail bekijkt, kunnen er twee verschillende fysische processen onderscheiden worden. In een eerste fase wordt de metalen leiding opgewarmd door het geconcentreerde zonlicht en in een tweede fase geeft het metaal haar warmte af aan de olie. In de literatuur wordt het systeem daarom vaak beschreven door twee vergelijkingen [1],[3]:

ρm Cm Am

∂Tm (x, t) = DIrr(t) − Hl G [Tm (x, t) − Ta (x, t)] − LHt [Tm (x, t) − Tf (x, t)] (4.1) ∂t ρf Cf Af

1

∂Tf (x, t) ∂Tf (x, t) + ρf Cf q(t) = LHt [Tm (x, t) − Tf (x, t)] 1 ∂t ∂x

(4.2)

de temperatuursafhankelijkheid van ρf , Cf , Ht en Hl wordt niet expliciet vermeld, om de notatie te

verlichten.

16

Hoofdstuk 4. Modelleren via niet-uniforme tijdsmonsters •

m:

17

metaal

• f : olie • ρ: dichtheid in kg/m3 • C: specifieke warmtecapaciteit in J/kg K • A: transversale oppervlakte in m2 • T (x, t): tijds- en plaatsafhankelijke temperatuur in K • : optische efficiëntie (dimensieloos) • D: breedte van de spiegels in m • Irr(t): warmteflux van de zon in W /m2 • Hl : warmteverliescoëfficiënt tussen metaal en omgeving in W /m2 • G: buitendiameter van de pijpleiding in m • L: binnendiameter van de metalen leiding in m • Ht : warmtetransfercoëfficiënt tussen metaal en olie in W /m2 K • q(t): debiet van de olie in m3 /s

Hierbij berekent de eerste vergelijking de temperatuur van de metalen leiding. In het rechterlid van vergelijking 4.1 slaat de eerste term op de warmte die het metaal van de zon ontvangt, de tweede term staat voor de warmte die het metaal verliest aan de omgeving en de derde term staat voor de warmte die het metaal doorgeeft aan de olie. De tweede vergelijking legt de temperatuur van de olie vast. In het rechterlid is enkel een term aanwezig die het warmtetransport tussen metaal en olie bepaalt. Bovenstaande vergelijkingen zijn enkel geldig in dat deel van het systeem waar de collectoren aanwezig zijn. In de andere delen (toevoer- en afvoersectie) dient in vergelijking 4.1 de term DIrr(t) weggelaten te worden, aangezien hier geen verhitting voorkomt. Verder moet ook de warmteverliescoëfficiënt Hl aangepast worden, aangezien deze stukken sterk ge¨ısoleerd zijn.

Hoofdstuk 4. Modelleren via niet-uniforme tijdsmonsters

4.1.1

18

Eigenschappen van de olie

De olie die gebruikt wordt om de warmte te transporteren is Santotherm 55. Twee kenmerkende eigenschappen van deze olie zijn haar lage warmtegeleidbaarheid en haar dichtheid die sterk temperatuursafhankelijk is. Door de laatste eigenschap weegt warme olie veel minder dan koude olie, waardoor er in de opslagtank lagen ontstaan met verschillende temperatuur (cfr. thermoclinisch effect). Hierna volgen twee temperatuursafhankelijke eigenschappen van de olie die gebruikt zullen worden: • dichtheid: ρf (T ) = 903 + 0, 672 T (kg/m3 ) • specifieke warmtecapaciteit: Cf = 1820 + 3, 478 T (J/kg K)

4.1.2

Bepalen van Ht

Ht geeft het tempo aan waarmee de olie wordt opgewarmd door de verhitte leiding bij een gegeven temperatuursverschil tussen olie en leiding. Het is als het ware een convectiecoëfficiënt (zelfde dimensie). Men kan aanvoelen dat dit tempo afhankelijk zal zijn van het debiet: hoe groter het debiet, hoe sneller het warmtetransport zal gebeuren. Vergelijk dit met de gedwongen convectie die gebruikt wordt om een processor sneller af te koelen. Daar zorgt een ventilator ervoor dat de luchtstroom met een grotere snelheid langs de processor passeert. Op deze manier wordt de afkoeling versneld. Daarnaast is Ht ook een functie van de temperatuur van de olie. Ht kan daardoor geschreven worden als:

• Ht = Hv (T ) q(t)0.8 (W /m2 K) • Hv (T ) = 2.17 106 − 5.01 104 T + 4.53 102 T 2 − 1.64T 3 + 2.1 10−3 T −4 (W /m5 K) Met q(t) opnieuw het debiet in m3 /s en Hv een polynoom die ge¨ıdentificeerd werd via een least squares methode [3]. Voor de waarden van de overige parameters in vergelijkingen 4.1 en 4.2 wordt verwezen naar sectie 4.4.


4.1.3

19

Onderstellingen

Omdat het niet nodig is elk fysisch proces tot in het kleinste detail te beschrijven, dit leidt immers tot een veel te complex model, worden er een paar vereenvoudigende veronderstellingen gemaakt: 1. ten eerste wordt het warmtetransport in de olie zelf verwaarloosd. Twee naburig gelegen punten in de olie op een verschillende temperatuur zullen geen warmte uitwisselen, noch door convectie, noch door conductie. Aan deze eis is grotendeels voldaan, gezien de lage warmtegeleidingscoëffiënt van de olie: 0.1923 − 1.3 10−4 T (W /m◦ C), vergelijk dit met 0,6 W /m°C voor water en 50 W /m°C voor staal. 2. ten tweede wordt de olie onsamendrukbaar ondersteld. Hierdoor zal bij een debietswijziging, ge¨ınduceerd door de pomp, het debiet in elk punt van het systeem op hetzelfde moment en op eenzelfde manier veranderen. Er treedt met andere woorden geen vertraging op in de debietsverandering tussen een punt net na de pomp en een punt ver van de pomp verwijderd. 3. het debiet in elk van de tien loops is hetzelfde. Stel dat de pomp een debiet levert van 10 l/s dan zal, aangezien er 10 parallelle loops zijn, door elk van de loops olie vloeien aan een tempo van 1 l/s. In werkelijkheid zal het debiet in loop 1 hoger zijn dan in loop 2, omdat deze eerste dichter bij de pomp ligt dan de tweede. Het debiet in loop 2 zal hoger zijn dan dat in loop 3, enz. 4. de eigenschappen van elke loop worden identiek beschouwd. In werkelijkheid zullen deze licht verschillen van loop tot loop, de spiegels zullen bijvoorbeeld niet allemaal dezelfde reflectiviteit hebben. Daarnaast zullen ook de warmteverliezen in elke loop licht verschillen. Voor de eenvoud worden deze verschillen genegeerd. 5. het warmteverlies van de toevoer- en afvoersectie wordt verwaarloosd; deze leidingen zijn sterk ge¨ısoleerd.


4.2

20

Verdelen in secties

In hoofdstuk 2 werd het systeem reeds opgedeeld in drie secties: de toevoersectie, de verhittingssectie en de afvoersectie. In de volgende paragrafen wordt elke sectie apart behandeld. Er wordt soms verwezen naar punten A, Bi , Ci en D waarbij A verwijst naar het begin van de toevoersectie, Bi naar het begin van loop i in de verhittingssectie, Ci naar het einde van loop i van de verhittingssectie en D naar het einde van de afvoersectie, zoals ook kan gezien worden in figuur 4.1:

Figuur 4.1: Definiëren van locaties A, Bi , Ci en D.

4.2.1

Toevoersectie

Dit deel brengt olie vanuit de bodem van de opslagtank naar het begin van elke loop van het collectorsysteem. De toevoerbuis bestaat uit één leiding met een diameter van 10,16 cm. De afstand tussen het begin van de toevoersectie en het begin van een loop van de verhittingssectie, is voor elke loop verschillend (zie figuur 4.2). Daardoor wordt in het model de leiding verdeeld in tien parallelle leidingen, met elk een verschillende lengte:


Figuur 4.2: Geometrie in werkelijkheid.

Figuur 4.3: Geometrie gebruikt in de simulator.

21


22

Deze parallellisering van de toevoersectie blijkt weinig invloed te hebben op het dynamisch gedrag van het systeem. Uit figuren C.1 tot 4.13 en 4.16 tot 4.18 blijkt er een goede overeenkomst te zijn tussen het dynamisch gedrag in werkelijkheid en in de simulator. In de toevoersectie worden warmteverliezen verwaarloosd en is er geen verhitting aanwezig aangezien er zich geen collectoren bevinden. De toevoersectie kan dus beschouwd worden als louter een dode tijd. Aangezien deze sectie opgedeeld werd in tien parallelle leidingen met verschillende lengte, zal elke leiding zijn eigen dode tijd τi1 (t) (i=1...10 en 1 : toevoersectie) hebben, die verder nog debiet- en dus tijdsafhankelijk is. Een uitdrukking voor deze tijdsafhankelijke dode tijd τi1 (t) valt eenvoudig op te stellen:

τi1 (t)

:

Vi1

Zt =

0

0

q(t ))dt

(4.3)

t−τi1 (t)

• τi1 (t): tijdsafhankelijke dode tijd voor leiding i, in de toevoersectie, in s • Vi1 : volume van leiding i in de toevoersectie in m3 • q(t): debiet in m3 /s

τi1 (t) geeft dus weer hoe lang het duurt om het volume Vi1 op te vullen met olie rekeninghoudend met het tijdsafhankelijk debiet. De temperatuur in het punt B van leiding i, TBi , is dus eenvoudigweg de temperatuur in punt A, een tijd τi1 (t) vroeger: TBi (t) = TAi (t − τi1 (t))

(4.4)

De praktische berekening van τi1 (t) in uitdrukkingen 4.4 en ??.

4.2.2

Verhittingssectie

Vervolgens is de verhittingssectie aan de beurt. In dit deel van het systeem zijn de collectoren aanwezig en wordt de olie opgewarmd. Naast de opwarming is er ook een warmteverlies aanwezig waarmee rekening gehouden moet worden. De berekening van de olietemperatuur gebeurt in twee stappen: eerst moet de temperatuur van de metalen leiding bepaald worden, vervolgens kan overgegaan worden tot de berekening van olietemperatuur.


23

Voor het beschrijven van de opwarming van de metalen leiding wordt gebruik gemaakt van de vooropgestelde vergelijking 4.1. Er wordt eerst een kleine wijziging doorgevoerd in deze vergelijking; uitdrukking 4.1 is geldig per lengte-eenheid van de leiding. Daarom wordt ze vermenigvuldigd met de lengte van de verhittingssectie, zodat de energiebalans per loop verkregen wordt. Aangezien in de verhittingssectie elke loop dezelfde lengte heeft, wordt telkens met eenzelfde waarde vermenigvuldigd.

ρm Cm Vm

∂Tm (x, t) Ht = Amirr Irr(t) − Hl Gs [Tm (x, t) − Ta (x, t)] − Amf [Tm (x, t) − Tf (x, t)] ∂t π (4.5)

• Vm : volume van één loop in m3 • Amirr : apertuur van de spiegel voor elke loop in m2 • Gs : buitendiameter pijpleiding maal lengte verhittingssectie in m2 • Amf : contactoppervlak tussen metaal en olie per loop in m2

Deze continue partiële differentiaalvergelijking dient gediscretizeerd te worden, zoals reeds werd aangehaald in het begin van dit hoofdstuk. De partiële afgeleide naar de tijd wordt vervangen door: ∂Tm (x, t) Tm (x, t) − Tm (x, t − ∆t) → ∂t ∆t

(4.6)

De discrete vorm van vergelijking 4.5 ziet er als volgt uit:

  Amirr Irr(t − ∆t) − Hl Gs [Tm (x, t − ∆t) − Ta (x, t − ∆t)] ∆t   Tm (x, t) = Tm (x, t−∆t)+ ρm Cm Vm −A Ht [T (x, t − ∆t) − T (x, t − ∆t)] m mf π f (4.7) De temperatuur van een stuk pijpleiding in de verhittingssectie is duidelijk samengesteld uit vier bijdragen (in respectievelijke volgorde): 1. de temperatuur van hetzelfde stuk leiding als een tijd ∆t ervoor


24

2. de temperatuurstijging gecreëerd door het geconcentreerde zonlicht 3. de temperatuurdaling, toe te schrijven aan warmteverlies naar de omgeving 4. de temperatuurverandering, te wijten aan de warmte-uitwisseling met de olie

Voor de berekening van de olietemperatuur wordt gebruik gemaakt van een andere vergelijking dan uitdrukking 4.2: dTf (x, t) Amf Ht = [Tm (x, t) − Tf (x, t)] dt ρf Cf Vf π

(4.8)

• Vf : volume van de olie in één loop van de verhittingssectie

De tijdsafgeleide in vergelijking 4.8 moet ook gediscretizeerd worden, wat leidt tot volgende iteratieve formule:

Tf (x, t) = Tf (x, t − ∆t) +

Amf Ht [Tm (x, t − ∆t) − Tf (x, t − ∆t)] ρf Cf Vf π

(4.9)

Het enige wat voor opwarming van de olie zorgt, is de warmte van de metalen leiding. Vergelijkingen 4.7 en 4.9 lijken eenvoudig te implementeren, ware het niet dat ze geen rekening houden met het tijdsvariërend debiet. De invoering van het tijdsafhankelijk debiet is gebaseerd op het gebruik van niet-uniforme tijdsmonsters, terug te vinden in [12]. Dit wordt uitvoerig besproken in deel 4.3.

4.2.3

Afvoersectie

De laatste sectie brengt de warme olie terug naar de top van opslagtank. De ge¨ısoleerde leiding die hiervoor zorgt, heeft een binnendiameter van 10,16 cm. De lengte van de leiding wordt verderop experimenteel bepaald, aangezien hierover geen nauwkeurige informatie beschikbaar is. Evenals bij de toevoersectie worden ook hier de warmteverliezen verwaarloosd waardoor ook de afvoersectie beschouwd kan worden als louter een dode tijd. Ook hier wordt in het model de leiding opgesplitst in 10 parallelle leidingen. Dat dit geen grote veranderingen in de dynamica van systeem teweegbrengt, werd reeds aangehaald in het stuk over de toevoersectie.


25

Door deze splitsing ontstaan er opnieuw tien verschillende dode tijden τi3 (t) (i=1...10 en 3 : afvoersectie ): voor elke leiding één. De dode tijden kunnen op exact dezelfde wijze berekend worden als vorige keer:

τi3 (t)

:

Vi3

Zt =

0

q(t )dt

0

(4.10)

t−τi3 (t)

• τi3 (t): tijdsafhankelijke dode tijd voor leiding i in de afvoersectie in s • Vi3 : volume van leiding i in de afvoersectie in m3 • q(t): debiet in m3 /s

De temperatuur in het punt D van leiding i, TDi , valt opnieuw eenvoudig te berekenen als de temperatuur in punt C, een tijd τi3 (t) vroeger: TDi (t) = TCi (t − τi3 (t))

(4.11)

Op deze manier berekent het model tien temperaturen voor het eind van de afvoersectie. In werkelijkheid is er natuurlijk maar één temperatuur, aangezien er maar één leiding is. Daarom wordt het gemiddelde genomen van de tien temperaturen die het model genereert. Zo wordt slechts één waarde bekomen en kan het model beter gevalideerd worden. 10

1 X TD (t) = TDi (t) 10

(4.12)

i=1

De praktische berekening van τi3 (t) in uitdrukkingen 4.9 en 4.10 wordt hierna behandeld.


4.3

26

Niet-uniforme tijdsmonsters

Figuur 4.4: buis

Beschouw een rechtlijnige, lange buis waardoor olie gepompt wordt. Het debiet van de olie, q(t), wordt geregeld door een pomp en is tijdsafhankelijk. Onderstel dat er in de olie een rood balletje aanwezig is, dat een verwaarloosbare massa heeft en dus even snel als de olie beweegt. De tijdsafhankelijke plaats van het balletje is weergegeven in figuur 4.5: De buis is opgedeeld in 6 segmenten, elk met eenzelfde lengte l. Het is zeer belangrijk dat elk segment dezelfde grootte heeft. Telkens het balletje een nieuw segment betreedt, wordt het tijdstip waarop dit gebeurt, tk (k=0...5), onthouden. Aangezien het debiet varieert in de tijd, zal de afstand tussen twee opeenvolgende tijdstippen, tk −tk−1 , eveneens variëren. Ze zijn met andere woorden niet uniform. Uit de grafiek volgt duidelijk dat het balletje op tijdstip t0 segment 1 betreedt, op tijdstip t1 segment 2, enzovoort tot het op tijdstip t5 segment 6 betreedt. De olie in segment 6 op een tijdstip tk is met andere woorden de olie in segment 1 op tijdstip tk−5 . Dit voorbeeldje illustreert dat de tijdsafhankelijke dode tijden τi1 (t) en τi3 (t) uit vergelijkingen 4.3, 4.4, 4.9 en 4.10 eenvoudig berekend kunnen worden, als de tijdstippen tk gekend zijn.

4.3.1

Berekening van tk

Om de tk ’s van het buizensysteem van de zonnecentrale te berekenen, moet dit buizensysteem ook opgedeeld worden in even grote segmenten. Voor de eenvoud worden enkel de tk ’s van de toevoersectie berekend. De werkwijze voor de verhittings- en de afvoersectie is exact dezelfde. Herinner u eraan dat in het model alledrie de secties opgedeeld zijn in tien parallelle leidingen. Elk van deze tien parallelle leidingen moet dus opgedeeld worden in


27

Figuur 4.5: Verplaatsing van het balletje.

segmenten. Aangezien elke leiding in het model van de toevoersectie een verschillende lengte heeft (zie figuur 4.3), zal elke leiding een verschillend aantal segmenten bevatten. Dit aantal wordt bepaald door het volume van elke leiding te delen door het volume van het segment. Dit segmentvolume Vsegm wordt vooropgesteld en moet voldoen aan een aantal voorwaarden: • zo mag het niet te klein zijn, anders zijn er teveel segmenten en duren de berekeningen te lang. • anderzijds mag het ook niet te groot zijn, want dan zijn er te weinig segmenten om nauwkeurige resultaten te krijgen. • verder moet het segmentvolume Vsegm zo gekozen worden, dat het liefst een geheel aantal keer in elke leiding past.

Een goede waarde voor het Vsegm kan bekomen worden door te kijken naar de geometrie van de toevoersectie. In het model wordt de toevoersectie verdeeld in tien parallelle leidingen, elk met een verschillende lengte en dus verschillend volume. Het volumeverschil ∆V tussen twee opeenvolgende leidingen is constant en gelijk aan 7,7 liter (zie figuur 4.3). Intu¨ıtief kan aangevoeld worden dat het segmentvolume best niet groter is dan dit volumeverschil ∆V . Stel dat dit wel het geval is en dat Vsegm bijvoorbeeld tweemaal ∆V bedraagt, dan zullen er


28

in het model, doordat elk volume afgerond wordt tot een geheel aantal keer Vsegm , zeker twee opeenvolgende leidingen zijn met eenzelfde aantal segmenten. Dit leidt tot twee dezelfde dode 1 , wat fouten in de uitkomst kan induceren. Dit fenomeen wordt verderop nog tijden τi1 en τi+1

verduidelijkt. Nu de toevoersectie in segmenten opgedeeld is, kan overgegaan worden op de berekening van de tk ’s. Om dit te kunnen doen, wordt een denkbeeldig balletje, met verwaarloosbare massa, in het begin van de toevoersectie geplaatst. De positieverandering van dit balletje en dus die van de olie kan gevolgd worden door een klok in te voeren. De periode van de klok is constant en gelijk aan Tklok . Tijdens elke periode van de klok worden de volgende bewerkingen uitgevoerd: • Tklok wordt vermenigvuldigd met het debiet op dat moment. Noem dit product Vtemp . • Vtemp wordt vervolgens vergeleken met het Vsegm .

Is Vtemp groter of gelijk aan Vsegm , dan beweegt het balletje tijdens Tklok zeker door gans segment 1 en is t1 = Tklok . t1 wordt onthouden en de klok wordt op nul gesteld. Is daarentegen Vtemp kleiner dan Vsegm , dan tikt de klok verder. Opnieuw wordt Tklok vermenigvuldigd met het debiet op dat moment. Deze nieuwe Vtemp wordt bij de vorige opgeteld. Vervolgens wordt deze som vergeleken met het segmentvolume. Is de som groter dan Vsegm , dan heeft het balletje twee klokperiodes Tklok nodig om één segment te doorlopen en is t1 = 2Tklok . Daarna wordt de klok op nul gesteld. Is de som daarentegen kleiner dan Vsegm , dan tikt de klok opnieuw verder enzovoort. Eenmaal een waarde voor t1 gevonden is, kan op dezelfde manier een waarde voor t2 gevonden worden. Het komt er dus op neer een goede waarde voor Tklok te vinden. Deze waarde mag niet te klein zijn, anders zal de berekening te lang duren. Noch mag ze te groot zijn, anders treden er onnauwkeurigheden op. Voor een goede eerste schatting wordt opnieuw gekeken naar de geomtrie van de leidingen. Het volume van één segment Vsegm werd gekozen op 7,7 liter. Het maximumdebiet in het systeem bedraagt 9 l/s(2 ). Aangezien er tien parallelle leidingen aanwezig zijn, zal het debiet in elke leiding 0,9 l/s bedragen. Eén segment wordt dus bij 2

Het debiet overschreed tijdens de experimenten nooit 9 l/s.

Hoofdstuk 4. Modelleren via niet-uniforme tijdsmonsters dit debiet opgevuld in

7,7 0,9

29

= 8, 55s. Aangezien 9 l/s het maximumdebiet was, zal 8,55 s de

minimum tijd zijn die de olie nodig heeft om één segment te vullen. Een eerste schatting voor Tklok is 10 s . Er treedt echter een probleem op bij een Tklok van 10 s. Er werd reeds berekend dat bij een debiet van 9 l/s, het ongeveer 8,55 s duurt om Vsegm te vullen: 7, 7 = 8, 55 ≈ Tklok 0, 9

(4.13)

Stel nu dat het debiet verlaagt van 9 l/s naar 8 l/s. Bij dit debiet doet de olie er eveneens 1 Tklok over om door een segment te vloeien: 7, 7 = 9, 625 ≈ Tklok 0, 8

(4.14)

Met deze Tklok zal het model dus geen onderscheid maken tussen een debiet van 9 l/s en een debiet van 8 l/s. In tabel 4.3.1 wordt voor verschillende waarden van Tklok weergegeven hoeveel klokperiodes de olie nodig heeft om één segment te vullen bij verschillende waarden van het debiet:

debiet per leiding

tijd om segment te vullen

aantal klokperiodes

(l/s)

(s)

Tklok = 10s

Tklok = 2s

Tklok = 1s

0,2

38,5

4

20

39

0,3

25,66

3

13

26

0,4

19,25

2

10

20

0,5

15,4

2

8

16

0,6

12,833

2

7

13

0,7

11

2

6

11

0,8

9,625

1

5

10

0,9

8,55

1

5

9

Tabel 1

Op deze manier vindt het model met een Tklok van 10 s slechts vier waarden voor de dode tijd van één segment, nl 4Tklok , 3Tklok , 2Tklok en Tklok , terwijl er in werkelijkheid toch acht


30

verschillende waarden voor het debiet zijn. Dit is natuurlijk niet erg nauwkeurig. Een mogelijkheid om dit probleem op te lossen is Tklok verlagen. Pas vanaf een waarde voor Tklok van 1 s, worden acht verschillende waarden voor de dode tijd gevonden. Dat een goede keuze voor Tklok belangrijk is om goede resultaten te verkrijgen kan ook gezien worden in figuur C.2 van bijlage C.

4.4

Implementatie

Nu beide vergelijkingen, die het fysisch proces beschrijven, omgezet zijn naar hun discrete vorm en er een methode gevonden werd om het variërend debiet in rekening te brengen, kan overgegaan worden tot de implementatie ervan in matlab. In figuur 2.8 van hoofdstuk 2, werd getoond dat de ingang van het systeem het debiet van de olie is, de uitgang de temperatuur op het eind van de centrale is en verder dat er drie storingen aanwezig zijn: de warmteflux van de zon, de ingangstemperatuur van de olie en de omgevingstemperatuur. Om het model te kunnen testen zijn er waarden nodig voor deze variabelen. Hiervoor kan beroep gedaan worden op data die verzameld werd gedurende drie dagen: 15, 16 en 28 juli 2003. Voor elk van deze drie dagen is er data beschikbaar over: • het tijdstip • het debiet • de temperatuur aan het begin van de toevoersectie • de temperatuur op het eind van de afvoersectie • de omgevingstemperatuur • de warmteflux van het invallend zonlicht • de temperatuur van de olie op het eind van elke leiding3 in de verhittingssectie.

Een eerste probleem dat opduikt is de bemonsteringstijd van deze data. Die bedraagt 30 s. Dus pas om de 30 seconden wordt een nieuwe waarde van de variabelen gegeven. Aangezien 3

leiding 1 was afgesloten wegens problemen, er waren dus slechts 9 leidingen operationeel.


31

een goede periode voor de interne klok 1 s bedraagt, is er elke seconde een waarde voor de variabelen nodig. Dit probleem kan eenvoudig opgelost worden door een zero order hold te gebruiken: elk monster van het data-acquisitie systeem wordt gedurende 30 klokperiodes vastgehouden. Vervolgens wordt berekend in hoeveel segmenten elke leiding in elke sectie opgedeeld moet worden. Hiervoor moeten eerst de volumes van elke leiding in elke sectie gekend zijn. De lengtes en diameters van elke sectie zijn weergegeven in de tabel 4.4:

toevoersectie

verhittingssectie

afvoersectie

lengte (m)

35.12

172

37.694

binnendiameter (cm)

10.16

2.969

10.16

Tabel 2

De toevoer- en afvoersectie worden in het model gesplitst in tien parallelle leidingen, met elk een verschillende lengte. Dit leidt tot de volgende waarden (leiding 1 wordt afgesloten):

m3

toevoersectie

verhittingssectie

afvoersectie

leiding 2

0.041

0.103

0.043

leiding 3

0.049

0.103

0.052

leiding 4

0.058

0.103

0.061

leiding 5

0.067

0.103

0.070

leiding 6

0.076

0.103

0.079

leiding 7

0.085

0.103

0.087

leiding 8

0.094

0.103

0.096

leiding 9

0.103

0.103

0.105

leiding 10

0.112

0.103

0.114

Tabel 3: volumes van de leidingen in elke sectie.

4

Deze waarde blijkt niet volledig correct te zijn, daarom wordt ze in een volgende hoofdstuk experimenteel

bepaald


32

toevoersectie

verhittingssectie

afvoersectie

leiding 2

5

13

6

leiding 3

6

13

7

leiding 4

8

13

8

leiding 5

9

13

9

leiding 6

10

13

10

leiding 7

11

13

11

leiding 8

12

13

13

leiding 9

13

13

14

leiding 10

15

13

15

Tabel 4: aantal segmenten in elke leiding.

De waarden voor het aantal segmenten worden gestockeerd in de matrix N. Nu wordt de tijdsiteratie met tijdsindex j gestart. Aangezien Tklok 1 seconde bedraagt, zal de tijdsiteratie gebeuren in stappen van 1 seconde. t0 wordt gelijkgesteld aan 0. Vervolgens wordt berekend hoelang het duurt vooraleer de olie één segment vult, gegeven het huidige debiet. Dit gebeurt door Tklok te vermenigvuldigen met het ogenblikkelijke debiet, q(j). De waarde van deze vermenigvuldiging, Vtemp , wordt vergeleken met het volume van één segment, Vsegm . Is Vtemp kleiner dan Vsegm , dan loopt de tijdsiteratie j verder naar j + 1. Opnieuw wordt Tklok vermenigvuldigd met het debiet op dat moment, q(j+1). Dit product wordt opgeteld bij Vtemp . Deze nieuwe Vtemp wordt opnieuw vergeleken met Vsegm . Stel nu dat de nieuwe Vtemp groter dan of gelijk is aan dan Vsegm , dan wordt t1 gelijkgesteld aan j+1 en Vtemp wordt teruggebracht op nul. t1 wordt bewaard in de niet-uniforme tijdsmonster vector TIJD= [t0 t1 t2 ...]. Opnieuw wordt een onderscheid gemaakt in de berekeningen voor de drie secties. Maar eerst worden nog enkele parameters gedefinieerd.

Temperatuursafhankelijke parameters Aangezien enkele van de parameters in de vergelijkingen 4.7 en 4.9 temperatuursafhanklijk zijn moeten deze eerst geëvalueerd worden. Het gaat hier in het bijzonder om de parameters Ht , ρf en Cf . Deze worden voor elk oliesegment apart berekend.


33

Bepalen van Hl en en de overige parameters Er moet nog een waarde toegekend worden aan de parameters Hl en die voorkomen in vergelijkingen 4.1, 4.5, 4.7 en 4.16. Deze parameters worden op een experimentele manier bepaald. Door ’trial and error’ blijken de beste waarden: • Hl = 28 W/m2 • = 0, 7 Verder moet er ook nog een waarde toegekend worden aan de parameters ρm en Cm : • ρm = 7850 Kg/m3 • Cm = 419 J/Kg K

4.4.1

Toevoersectie

Aangezien de toevoersectie louter als dode tijd beschouwd kan worden en de waarden voor tk tot het huidige tijdstip bekend zijn, is het eenvoudig om te temperatuur op het eind van de toevoersectie te berekenen. Formule 4.4 wordt als volgt ingevoerd:

TBi (t) = TAi (t − τi1 (t)) → TBi (tk ) = TAi (tk−l )

(4.15)

waarbij l aangeeft in hoeveel segmenten leiding i in de toevoersectie verdeeld is. Deze waarde is terug te vinden in de matrix N.

4.4.2

Verhittingssectie

Hierna worden de temperatuursafhankelijke coëfficiënten ρf , Cf en Ht geëvalueerd. Voor de evaluatie op tk wordt gekeken naar de waarde van de olietemperatuur op het vroegere tijdstip tk−1 , zoals ook kan gezien worden in uitdrukkingen 4.7 en 4.9. ∆t is gelijk aan het verschil tk − tk−1 . Eens deze waarden gekend zijn, hebben we alle waarden van de rechterleden in uitdrukkingen 4.7 en 4.9 (in de veronderstelling dat de iteratie al enige tijd bezig is en er reeds waarden bestaan voor Tm (x, t − ∆t) en Tf (x, t − ∆t)). Voor de implementatie moet nog onderscheid gemaakt worden tussen de verschillende loops (index i; i=1..9) en er moet


34

nog ge¨ıtereerd worden over alle segmenten van elke loop (index n). De implementatie van uitdrukking 4.7 is eenvoudig, de vergelijking kan zó overgenomen worden:

 Tm (n, tk , i) = Tm (n, tk−1 , i)+

Amirr Irr(tk−1 ) − Hl Gs [Tm (n, tk−1 , i) − Ta (n, tk−1 , i)]

∆t  ρm Cm Vm −A Ht [T (n, t , i) − T (n, t , i)] m mf π k−1 f k−1

 (4.16)

Deze vergelijking drukt uit dat de temperatuur van het metaal in segment n van leiding i op een tijdstip tk gelijk is aan de temperatuur van datzelfde stuk metaal op het vorige tijdstip tk−1 (eerste term in het rechterlid), vermeerderd met de warmte ontvangen van de zon (eerste term binnen de rechte haken), verminderd met de warmte die weggelekt is naar de omgeving(tweede term binnen de rechte haken) en gecorrigeerd voor de warmte-uitwisseling met de olie (derde term binnen de rechte haken). Voor de implementatie van vergelijking 4.9 komt het grote voordeel van de niet-uniforme tijdsmonsters naar voren. Vergelijkingen 4.8 en 4.9 zijn eigenlijk enkel geldig als de olie stilstaat; ze houden geen rekening met de verplaatsing van de olie. De niet-uniforme tijdsmonsters werden in de eerste plaats ingevoerd, om dit te verhelpen. Stel dat het rode oliesegment in figuur 4.6 wordt gevolgd. Onderstel dat dit oliesegment zich op tijdstip tk−1 in het segment n − 1 bevindt. Het gebruik van niet-uniforme tijdsmonsters zorgt ervoor dat ditzelfde segment zich op tijdstip tk in segment n zal bevinden. Intu¨ıtief kan aangevoeld worden dat de temperatuur van het rode olievolume op tijdstip tk in segment n, de temperatuur van datzelfde olievolume zal zijn op het vroegere tijdstip tk−1 in het segment n − 1, vermeerderd met de warmte die de olie ontvangen heeft van het metaal in segment n − 1 (zie figuren 4.64.7. Of in computertaal:

Tf (n, tk , i) = Tf (n − 1, tk−1 , i) +



Amf Ht [Tm (n − 1, tk−1 , i) − Tf (n − 1, tk−1 , i)] (4.17) ρf Cf Vf π


Figuur 4.6: Een deel uit loop i van de verhittingssectie op tijdstip tk−1 .

Figuur 4.7: Een deel uit loop i van de verhittingssectie op tijdstip tk .

35


36

Nu kunnen de temperaturen TCi (i=2...10) in het punt C berekend worden. Deze temperaturen zijn voor elke loop i gelijk aan de temperatuur van de olie in het laatste segment van loop i in de verhittingssectie:

TCi (tk ) = Tf (ni , tk , i) (i = 2...10)

met ni het aantal segementen in loop i van de verhittingssectie. Deze waarde is terug te vinden in de matrix N.

4.4.3

Afvoersectie

Het laatste deel van het traject bestaat uit de afvoersectie. Net als de toevoersectie kan ook deze sectie beschouwd worden als een sectie die enkel een dode tijd invoert. Formule 4.15 kan dus rechstreeks overgenomen worden, mits aanpassing van de indices:

TDi (t) = TCi (t − τi3 (t)) → TDi (tk ) = TCi (tk−l0 )

(4.18)

0

Waar l aangeeft in hoeveel segmenten loop i van de afvoersectie verdeeld is. Ook deze waarde werd gestockeerd in matrix N. Op deze manier worden 9 (normaal 10 maar loop 1 is afgesloten) temperaturen verkregen op het eind van de afvoersectie. Om deze te kunnen vergelijken met de eenwaardige, werkelijk eindtemperatuur wordt het gemiddelde van de berekende uitgangstemperaturen genomen: 9

TD (t) =

1X Di (t) 9

(4.19)

i=1

Op deze manier worden de gesimuleerde waarde voor de uitgangstemperatuur bekomen. Het probleem dat hier optreedt, is dat er in elk van de vier vergelijkingen 4.15, 4.16, 4.17 en 4.18 verwezen wordt naar een temperatuurswaarde op een vorig tijdstip tk−1 ,tk−l of tk−l0 . Dus voor de berekening van de temperatuur op de eerste tijdstippen t0 , t1 , t2 ,... kan het zijn dat wordt verwezen naar tijdstippen met een negatieve index, wat natuurlijk niet kan. In de volgende paragraaf wordt uitgelegd hoe dit opgelost wordt.


4.4.4

37

Beginwaarden

Wanneer het programma gestart wordt, is er nood aan enkele beginwaarden. Bij de berekening van vergelijking 4.15 wordt verwezen naar de temperatuur op tijdstip tk−l . Indien het programma zich op tijdstip t0 bevindt, bestaat er nog geen waarde voor de temperatuur op dit tijdstip. TBi kan dus niet berekend worden. Daarom worden de eerste l waarden van TBi ge¨ınitialiseerd op TAi (0), de ingangstemperatuur5 op tijdstip 0. Dus:

TBi (t0 ) = TBi (t1 ) = TBi (t2 ) = ... = TBi (tl−1 ) = TAi (t0 )

De volgende stap in de berekening is het uitrekenen van vergelijking 4.16. In deze vergelijking wordt telkens verwezen naar de temperatuur van het metaal en de olie en naar de warmteflux van de zon op een vorig tijdstip tk−1 . Hier treedt enkel een probleem op bij de berekening van Tm op het moment t0 . Dan kan niet verwezen worden naar t−1 . De temperatuur van Tm moet dus ook ge¨ınitialiseerd worden. Er wordt hier ook gekozen om als beginwaarde de ingangstemperatuur te nemen:

Tm (n, t0 , i) = TAi (t0 ) (i = 2...10; n = 1...ni )

met ni het aantal segmenten in loop i in de verhittingssectie. Hierna wordt de temperatuur van de olie in de verhittingssectie berekend. Formule 4.17 wordt hiervoor uitgerekend. Hier treedt een dubbel probleem op: er wordt enerzijds verwezen naar de temperatuur op een vorig tijdstip (tk−1 ) en anderzijds naar de temperatuur in een vorig segment (n − 1). Als beginwaarde voor de temperatuur van de olie wordt ook hier de ingangstemperatuur gekozen:

Tf (n, t0 , i) = TAi (t0 ) (i = 2...10; n = 1...ni )

Het tweede probleem, de verwijzing naar het segment n − 1 als n = 1 levert niet eenmalig, maar voor alle tijdstippen tk een probleem. Aangezien de verhittingssectie voorafgegaan 5

Deze waarde wordt gemeten en hoeft dus niet berekend te worden


38

wordt door de toeversectie (zie figuur 2.6), kan verwezen worden naar het laatste segment van de toevoersectie voor de berekening van de olietemperatuur in het eerste segment van de verhittingssectie (n = 1). De temperatuur van de olie op het einde van de toevoersectie speelt hier dus de rol van het segment n − 1 als n = 1. Dit leidt tot de volgende uitdrukking voor de randvoorwaarde:

Tf (0, tk , i) = TBi (tk ) (i = 2...10)

In de laatste sectie, de afvoersectie, treedt net hetzelfde probleem op als in de toevoersectie. 0

Bij het uitrekenen van vergelijking 4.18 wordt er de eerste l tijdstippen verwezen naar een 0

temperatuur die niet gekend is. Hier dienen de eerste l waarden voor TDi ge¨ınitialiseerd te worden. Als beginwaarde wordt opnieuw gekozen voor de ingangstemperatuur:

TDi (t0 ) = TDi (t1 ) = TDi (t2 ) = ... = TDi (tl0 −1 ) = TAi (t0 )

Op deze manier worden temperatuurswaarden voor alle olie- en metaalsegmenten bekomen op alle tijdstippen. Dit impliceert dat het berekenen van de temperatuursafhankelijke coëfficiënten Ht , ρf en Cf geen problemen oplevert. Nu alles ge¨ımplementeerd is en de begin- en randvoorwaarden gedefinieerd zijn, kan overgegaan worden tot het simuleren.

4.5

Resultaten

Nu het volledige programma af is kan overgegaan worden tot simuleren. Hiervoor kan beroep gedaan worden op data gegenereerd op drie verschillende testdagen: 15, 16 en 28 juli. Aangezien de temperatuur in de centrale gemeten wordt op zowel het einde van de verhittingssectie (9 waarden, voor elke loop 1) als op het einde van de afvoersectie (één waarde), kan er op deze twee punten een vergelijking gemaakt worden tussen de werkelijke temperatuur en de temperatuur bekomen via de simulatie. Omdat het wat omslachtig is om de negen gemeten temperatuurswaarden op het einde van de verhittingssectie te vergelijking met de negen gesimuleerde waarden, wordt telkens het gemiddelde genomen. Met andere woorden, de negen


39

gemeten temperatuurswaarden worden vervangen door het gemiddelde van die negen waarden. Hetzelfde wordt gedaan met de gesimuleerde waarden. Deze gemiddelde waarden wordt in de volgende grafieken aangeduid met Tgem . De fout in de grafieken is gedefinieerd als het verschil tussen de werkelijke en de gesimuleerde waarde, genormeerd op het temperatuursbereik van de werkelijke waarden. De eenheid van de tijdsas in onderstaande grafieken is één tijdsmonster. Een tijdsmonster slaat op de bemonsteringsperiode van klok en bedraagt Tklok of dus 1 seconde.

4.5.1

15 juli:

Eerst worden de resultaten bekeken als de gegevens van 15 juli gesimuleerd worden. In figuur C.16 zijn de resultaten te zien van de temperatuur op het einde van de verhittingssectie en in de tweede grafiek die op het einde van de afvoersectie.

zonlicht (W/m²)

1000

10 debiet (l/s)

800 600

5

400

2000

4000

6000 8000 tijdsmonster

0

10000 12000 14000

0

2000

4000


10000 12000 14000

160 250 150 140 200 130 120 0

2000

4000


10000 12000 14000

10

Tgem (°C)

ingangstemperatuur (°C)

200 0

150

fout in %

100 5 50 0

0

2000

4000


10000 12000 14000

reële temperatuur gesimuleerde temperatuur 0

2000

4000


10000 12000 14000

Figuur 4.8: Resultaten voor de temperatuur in punt C op 15 juli. 6

Voor een groter formaat van de figuren wordt verwezen naar bijlage C


40

De berekende waarden aan het begin van de simulatie blijken niet goed overeen te stemmen met de werkelijkheid. Dit komt omdat de eerste temperatuurswaarden die worden gemeten in punt C, afkomstig zijn van de olie die nog in de toevoer- en verhittingssectie aanwezig was. Deze olie heeft een volledige nacht in de leidingen doorgebracht en is daarom vrij koud. De eerste temperatuurswaarden die berekend worden door de simulator, starten met de eerste waarde van de ingangstemperatuur. Daarom begint de gesimuleerde curve op een hogere temperatuur en begint ze ook onmiddellijk te stijgen. Verder kan rond het 2900ste tijdsmonster een klein dipje (omcirkeld) opgemerkt worden in de reële temperatuur. Na nauwkeuriger onderzoek van de data, bleek dit afkomstig te zijn van loop 9 (zie figuur C.1 in bijlage C). In het temperatuursverloop op het einde van loop 9 is rond dit tijdstip een sterke daling van de temperatuur merkbaar. Wellicht werden de spiegels van deze loop even van de zon weggedraaid. Helemaal op het einde van het experiment treedt er ook een grote fout op. Dit komt omdat op het eind van de dag de spiegels volledig naar beneden gericht worden om ze beschermen tegen stof en regen. Daardoor neemt de temperatuur van de olie op het eind zo sterk af. Dit verdraaien van de spiegels wordt niet gesimuleerd. Buiten deze fenomenen komt de gesimuleerde temperatuur redelijk goed overeen met de werkelijke. De fout blijft bijna overal onder de 5%. In de foutplot werd bij 5% een horizontale lijn getekend om te zien waar de simulator een fout genereert groter dan 5%. Na de 10 000 tijdsmonsters treedt er echter een grote fout op, groter dan de 5% limiet. Opvallend is dat het debiet op dit tijdstip laag is. De tendens dat er grotere fouten optreden bij een lager debiet zal op volgende testdagen ook merkbaar zijn, wat leidt tot het vermoeden dat de enige coëfficiënt die van het debiet afhankelijk is, Ht , niet bepaald nauwkeurig is. In de figuur 4.9 staan de resultaten voor de uitgangstemperatuur van het systeem, met andere woorden de temperatuur in punt D:


zonlicht (W/m²)

1000

41

10 debiet (l/s)

800 600

5

400

2000

4000


0

160

260

150

240

140

220

130 120 0

2000

4000


10000 12000 14000

10

fout in %

0

10000 12000 14000

5

uitgangstemperatuur (°C)


200 0

2000

4000

10000 12000 14000

200 180 160 140 120 100 reële temperatuur gesimuleerde temperatuur

80 0


0

2000

4000


10000 12000 14000

60

0

2000

4000


10000 12000 14000

Figuur 4.9: Resultaten voor de temperatuur in punt D op 15 juli.

De werkelijke temperatuur in punt D start hoger dan die in punt C. Ook dit is vrij logisch: de afvoersectie is veel sterker ge¨ısoleerd dan de verhittingssectie. Hierdoor is de temperatuur van de olie die ’s nachts in de afvoersectie aanwezig was, minder sterk afgenomen dan de temperatuur van de olie in de verhittingssectie. Vervolgens is er een sterke daling waarneembaar in de werkelijke temperatuur in punt D. Op dit ogenblik vloeit de olie die ’s nachts in de verhittingssectie aanwezig was en dus sterker is afgekoeld, langs de temperatuursensor in punt D. Daarna stijgt de olietemperatuur sterk tot de werkelijke en gesimuleerde temperatuur elkaar goed benaderen. Dit gebeurt rond het 2000ste tijdsmonster. Als er gekeken wordt naar de grafiek die de fout weergeeft, kan gezien worden dat die licht toeneemt in vergelijking met figuur C.1. De fout komt nu enkele keren boven de 5%. Dit is te wijten aan een lichte tijdsshift die optreedt tussen de werkelijke en de gesimuleerde waarden. De gesimuleerde waarden komen eerder in tijd dan de werkelijke. Deze shift komt nog beter tot


42

uiting in de resultaten van 28 juli, zoals zal blijken. Aangezien deze shift niet aanwezig is in de berekening van de temperaturen in punt C, mag besloten worden dat ze volledig te wijten is aan de afvoersectie. Later komt aan bod hoe deze shift weggewerkt kan worden. Daarnaast treden er opnieuw fouten op voor de beginwaarden, die nooit goed benaderd kunnen worden, en voor de waarden na de 10 000 tijdsmonsters. Maar er werd reeds gesuggereerd dat deze laatste fout waarschijnlijk te wijten is aan een verkeerde uitdrukking voor Ht .

4.5.2

16 juli:

Hierna volgen de resultaten voor 16 juli. Opnieuw eerst de waarden in punt C en erna die in punt D.

debiet (l/s)

10

800 600


400 0

0.5

1 tijdsmonster

1.5

5

0

2

0

0.5

4

x 10

1 tijdsmonster

1.5

2 4

x 10

200 250 150 200

100

50

0

0.5

1 tijdsmonster

1.5

2 4

x 10

10

Tgem (°C)

zonlicht (W/m²)

1000

150

fout in %

100 5 50 0

0

0.5

1 tijdsmonster

1.5

2 4

x 10


0.5

1 tijdsmonster

1.5

2 4

x 10

Figuur 4.10: Resultaten voor de temperatuur in punt C op 16 juli.

Hier gelden dezelfde opmerkingen in verband met de slechte resultaten voor de beginwaarden als voor 15 juli. Uit de grafiek van de fout blijkt dat de resultaten vrij goed zijn; de fout


43

overstijgt nergens de 5% lijn, behalve voor de beginwaarden en helemaal op het einde. Er komen opnieuw mindere resultaten voor bij lage waarden van het debiet, zo rond het 14 000ste en 16 000ste tijdsmonster. Een mogelijke oorzaak is reeds aangehaald: de debietafhankelijke parameter Ht .

debiet (l/s)

10

800 600


400 0

0.5

1 tijdsmonster

1.5

0

2

0

0.5

4

200

1 tijdsmonster

1.5

2 4

x 10

260 240

150

220 100

50

0

0.5

1 tijdsmonster

1.5

2 4

x 10

10 fout in %

5

x 10

5


zonlicht (W/m²)

1000


80 0

0

0.5

1 tijdsmonster

1.5

2 4

x 10

60

0

0.5

1 tijdsmonster

1.5

2 4

x 10


Uit de foutgrafiek van figuur 4.11 blijkt duidelijk dat die dramatisch is toegenomen. Deze toename kan opnieuw verklaard worden door het optreden van een tijdsshift. Hoe dit verholpen kan worden, komt later aan bod.

4.5.3

28 juli:

De derde meetwaarden werden bekomen op 28 juli. Hieronder zijn de resultaten van deze data zichtbaar. Er wordt opnieuw eerst gekeken naar de temperatuur in het punt C en vervolgens die in punt D.


zonlicht (W/m²)

1000

44

10 debiet (l/s)

800 600

5

400 0.5

1 1.5 tijdsmonster

2

0

2.5

0

0.5

4

x 10

1 1.5 tijdsmonster

2

2.5 4

x 10

200 250 150 100 200 50 0

0

0.5

1 1.5 tijdsmonster

2

2.5 4

x 10

10

Tgem (°C)


200 0

150

fout in %

100 5 50 0

0

0.5

1 1.5 tijdsmonster

2

2.5


4

x 10

0.5

1 1.5 tijdsmonster

2

2.5 4

x 10

Figuur 4.12: Resultaten voor de temperatuur in punt C op 28 juli.

Voor de resultaten van 28 juli zien we opnieuw de fout in de eerste waarden voor de temperatuur. Ongeveer halverwege het experiment, zo rond de 10 000 tijdsmonsters, is er een dip in de ingangstemperatuur (omcirkeld). Deze komt duidelijk naar voor in zowel de werkelijke als de gesimuleerde eindtemperatuur. Dit onderstreept de sterke invloed van de ingangstemperatuur op het systeem. Het feit dat de dip ook voorkomt in de gesimuleerde waarden, is een teken dat het model goed reageert op deze storing.


zonlicht (W/m²)

1000

45

10 debiet (l/s)

800 600

5

400 0.5

1 1.5 tijdsmonster

2

0

200

260

150

240

100

220

50 0

0

0.5

1 1.5 tijdsmonster

0.5

4

2

2.5 4

x 10

10 fout in %

0

2.5 x 10



200 0

5

2

2.5 4

x 10


80 0

1 1.5 tijdsmonster

0

0.5

1 1.5 tijdsmonster

2

2.5 4

x 10

60

0

0.5

1 1.5 tijdsmonster

2

2.5 4

x 10


In bovestaande figuur komt de tijdsshift pas echt goed naar voor. Door de vele pieken in de eindtemperatuur vanaf 15 000 tijdsmonsters kan een schatting gemaakt worden van de grootte van de shift: ongeveer 120 s of 2 minuten.

4.5.4

Conclusies

Voor alle drie de testdagen zijn de resulaten in het begin, ongeveer de eerste 2500 tijdsmonsters, niet erg nauwkeurig. Dit komt omdat de temperatuursdistributie van de olie in het systeem, voor de centrale wordt opgestart, niet gekend is. Dat de tijd die de olie nodig heeft om gans het systeem te doorlopen ongeveer 1800 tijdsmonsters bedraagt bij een debiet van 2 l/s, lijkt deze hypothese deels te bevestigen. Het lage debiet van 2 l/s wordt in het begin bijna altijd gebruikt om de pomp niet te overbelasten aangezien de olie in het begin koud en dus erg visceus is. Voorts werden er steeds grote fouten ontdekt als het debiet laag is.


46

Dit leidt tot het vermoeden dat de enige coëfficiënt die van het debiet afhangt, Ht , niet erg nauwkeurig is. Een derde opmerking is dat er bij de temperatuursberekeningen in punt D een tijdsshift optreedt van ongeveer 2 minuten. In wat volgt zal de lengte van de afvoersectie op een experimentele manier berekend worden om deze shift weg te werken. De voordelen van het model zijn dat de simulator slechts een kleine steady sate fout heeft. Dit is zichtbaar in de resultaten van 15 juli, tussen de 2000 en de 4000 seconden, aangezien gedurende deze periode de storingen constant zijn. Verder reageert het model ook goed op veranderingen in de intensiteit van het zonlicht. Dit is zichtbaar in de resultaten van 16 juli. Hier neemt de intensiteit van het zonlicht vrij sterk toe tussen de 5000 en de 10 000 seconden. Dit genereert geen grote fouten in de berekende resultaten. En als laatste reageert de simulator blijkbaar ook goed op veranderingen van de ingangstemperatuur. Dit is zichtbaar in de resultaten van 28 juli.

4.6

Experimentele bepaling van de lengte van de afvoersectie.

Er wordt nu getracht om op een experimentele manier de lengte van de afvoersectie te berekenen. De gegevens over de lengte ervan, lijken niet zeer nauwkeurig. Deze experimentele bepaling gebeurt als volgt: 1. Er wordt begonnen met de data van 28 juli te analyseren. Hieruit kan gevonden worden hoeveel tijd de olie nodig heeft om van het einde van loop 2 naar het einde van de afvoersectie te reizen. 2. Vervolgens wordt gekeken bij welk debiet dit gebeurt. 3. De lengte van het afvoerstuk kan dan eenvoudig berekend worden met de volgende formules (de precieze binnendiameter van de afvoersectie is wel gekend): ∆t q(t) = Vf Vf = lA ⇒l=

• ∆t: tijdsverschil in seconden

∆t q(t) A


47

• q(t): het debiet in m3 /s • Vf : volume van de olie in m3 • A: doorsnede van de afvoersectie in m2 • l: lengte van de afvoersectie in m

Figuur 4.14: detail van de afvoersectie.

Aangezien het einde van loop 2 zich het dichtst bij het einde van de afvoersectie bevindt (loop 1 is afgesloten), zal een temperatuursverandering op het einde van de afvoersectie eerst en vooral te wijten zijn aan een temperatuursverandering in T2 7 . Daarom wordt gekeken naar het tijdsverschil tussen een temperatuursverandering in T2 en Tout (=TD ), dat ge¨ınduceerd wordt door een debietverandering. Dit tijdsverschil wordt ∆t genoemd en kan gehaald worden uit de figuur 4.15 (Hier is de eenheid van de tijdas de bemonsteringstijd van de data en die bedraagt 30 s): 7

Dit is de temperatuur geregistreerd door de temperatuursensor op het einde van loop 2


48

255

250

245

temperatuur (°C)

240

235

230

225

220

215

210

einde afvoersectie einde loop 2 500

520

540

560


620

640

660

680

Figuur 4.15: Detail van T2 en Tout op 28 juli.

Als gekeken wordt naar het debiet dat op de overeenkomstige tijdstippen geleverd werd, kan de volgende tabel opgesteld worden:

wijziging

wijziging

∆t

q(t)

l

T2 (K)

Tout (K)

(s)

(m3 /s)

(m)

258

263

150

5 10−3

92,51

270

281

330

5 10−3

203,52

395

399

120

8 10−3

118,41

405

412

210

8 10−3

207,22

519

525

180

7 10−3

155.42

567

573

180

9 10−3

199,82

597

603

180

8 10−3

177,62

652

658

180

7,5 10−3

166,52

Tabel 5: Experimentele bepaling van de lengte van de afvoersectie.

De eerste kolom geeft het tijdstip weer van een temperatuurswijziging op het eind van loop 2, met als tijdseenheid de bemonsteringstijd van de data, nl. 30 s. De tweede kolom geeft het tijdstip aan van de overeenkomstige temperatuurswijziging op het einde van de afvoersectie.


49

Het gemiddelde van alle berekende lengtes bedraagt 165,13 m en deze waarde wordt nu gebruikt in de simulator. Dit levert de volgende resultaten voor de temperatuur op het eind van de afvoersectie (punt D) voor 15, 16 en 28 juli:

zonlicht (W/m²)

1000

10 debiet (l/s)

800 600

5

400

2000

4000


0

160

260

150

240

140

220

130 120 0

2000

4000


10000 12000 14000

10

fout in %

0

10000 12000 14000

5



200 0

2000

4000

10000 12000 14000


80 0


0

2000

4000


10000 12000 14000

60

0

2000

4000


10000 12000 14000

Figuur 4.16: Resultaten voor de temperatuur in punt D op 15 juli, met de nieuwe lengte voor de afvoersectie.


debiet (l/s)

10

800 600


400 0

0.5

1 tijdsmonster

1.5

0

2

0

0.5

4

200

1 tijdsmonster

1.5

2 4

x 10

260 240

150

220 100

50

0

0.5

1 tijdsmonster

1.5

2 4

x 10

10 fout in %

5

x 10

5


zonlicht (W/m²)

1000

50


80 0

0

0.5

1 tijdsmonster

1.5

2 4

x 10

60

0

0.5

1 tijdsmonster

1.5

2 4

x 10



zonlicht (W/m²)

1000

51

10 debiet (l/s)

800 600

5

400 0.5

1 1.5 tijdsmonster

2

0

200

260

150

240

100

220

50 0

0

0.5

1 1.5 tijdsmonster

0.5

4

2

2.5 4

x 10

10 fout in %

0

2.5 x 10

5



200 0

2

2.5 4

x 10


80 0

1 1.5 tijdsmonster

0

0.5

1 1.5 tijdsmonster

2

2.5 4

x 10

60

0

0.5

1 1.5 tijdsmonster

2

2.5 4

x 10


De tijdsshift is praktisch verdwenen. Het positieve effect van de nieuwe lengte van de afvoersectie wordt pas duidelijk als gekeken wordt naar de procentuele fout (linksonder in de figuren 4.16, 4.17 en 4.18). Als de beginwaarden genegeerd worden, is deze enkel nog groter dan 5% voor 15 juli, maar dat is hoogstwaarschijnlijk te wijten aan een slechte waarde voor Ht .

Hoofdstuk 5

Methode II: Proper Orthogonal Decomposition Een veelgebruikte methode om systemen waarin bewegende vloeistoffen of gassen voorkomen te modelleren is via CFD-technieken (Computational Fluid Dynamics). Deze methode verdeelt het systeem in een groot aantal1 segmenten. Op elk van deze segmenten worden de wetten van behoud van energie, momentum en massa toegepast. Dit leidt voor elk segment tot een aantal differentiaalvergelijkingen. Het is duidelijk dat via deze methode heel wat berekeningen moeten worden uitgevoerd. Tegenwoordig is dat meestal geen probleem, gezien de grote rekenkracht van moderne computers. Doch, in sommige gevallen is het gewenst de rekenintensiteit te verlichten, bijvoorbeeld in het geval van ’real time prediction’, ’model based predictive control’ (MBPC), ... Voor deze applicaties moet de simulator veel sneller werken dan het werkelijke proces omdat de simulator de toekomstige uitgang moet voorspellen. Vandaar de nood aan technieken die op een snellere manier de uitgang van een proces kunnen reproduceren. Ondermeer hiervoor werd de POD-methode ontwikkeld. Het aantal berekeningen kan via deze methode soms met een factor 100 gereduceerd worden [11]. De tweede simulator steunt op de POD-methode. Zoals reeds aangehaald in de inleiding gebruikt ze een aantal orthonormale basisfuncties die uit gemeten of berekende gegevens kunnen gedistilleerd worden. Deze gegevens worden best gemeten (of berekend) in een situatie waarbij het systeem sterk geëxciteerd wordt teneinde zoveel mogelijk informatie te verzamelen 1

Dit aantal kan oplopen tot 108 .

52

Hoofdstuk 5. Methode II: Proper Orthogonal Decomposition

53

in de basisfuncties. De orthonormale basisfuncties geven de ruimtelijke distributie van een bepaalde fysische grootheid, zoals bijvoorbeeld de temperatuur, weer. Om verder de tijdsafhankelijkheid van het systeem te beschrijven wordt elke basisfunctie vermenigvuldigd met zijn tijdsafhankelijke coëfficiënt. Deze coëfficiënten kunnen gevonden worden door de differentiaalvergelijking, die het systeem beschrijft, op de basisfuncties te gaan projecteren. Uit deze projectie ontstaat een set van vergelijkingen die samen met een aantal rand- en beginvoorwaarden uitdrukkingen genereren voor de tijdsafhankelijke coëfficiënten [8],[9],[10],[11].

5.1

Theorie

Omdat het misschien allemaal wat magisch overkomt is een stukje theorie hier op zijn plaats. Omdat dit werk echter niet theoretisch is opgevat, blijft dit inleidend stukje theorie beperkt tot de belangrijkste begrippen. Voor een wiskundige, maar zeer begrijpelijke uitleg wordt verwezen naar [11]. Een belangrijk theorema hierin luidt als volgt: Theorema 5.1 Als (χ, (., .)) een separabele Hilbert ruimte2 is met orthonormale basis {φi }i∈I , dan kan elk element f ∈ χ geschreven worden als:

f=

X

(f, φi )φi

(5.1)

i∈I

met χ een lineaire vectorruimte en I een aftelbare verzameling.

Hierin stelt (f, φi ) de projectie van de functie f op φi voor. We kunnen deze expansie nu toepassen op de tijds- en plaatsafhankelijke temperatuur T (x, t) van het systeem. Om dit te kunnen doen moet de ruimte χ waartoe de functie T (x, t) behoort, voldoen aan de volgende voorwaarden: • de ruimte χ moet een lineaire ruimte zijn, voorzien van optelling en vermenigvuldiging. • de ruimte χ moet voorzien zijn van het inwendig product: χ × χ → R, zodat het begrip orthonormaliteit geformuleerd kan worden. 2

Zie bijlage B voor een definitie van separabele Hilbert ruimte.


54

• de ruimte χ moet een separabele Hilbertruimte zijn.

Aan deze drie voorwaarden is voldaan [11] en theorema 5.1 is van toepassing. De functie T (x, t) kan als volgt voorgesteld worden:

T (x, t) =

X

ai (t)φi (x)

(5.2)

i∈I

met voor alle t ∈ T: ai (t) = (T (x, t), φi (x)) de tijdsafhankelijk coëfficiënten. Beschouw nu Tsnap als Tsnap = {T(t) ∈ χ | t ∈ T} Met T(t) wordt de temperatuursdistributie over het ganse systeem op tijdstip t bedoeld en T ⊆ R stelt het tijdsdomein voor. Tsnap is dus een verzameling van tijden plaatsafhankelijke temperaturen, die worden bekomen uit metingen of berekeningen. Uit deze Tsnap moet de set van orthonormale basisfuncties {φi }i∈I gevonden worden.

5.1.1

Bepalen van de basisfuncties

Alnaargelang de Hilbertruimte χ eindig of oneindig dimensionaal is, zal de sommatie in 5.2 gebeuren over een eindig of oneindig aantal basisfuncties. Gezien een computer slechts kan werken met een eindig aantal variabelen, worden de functies φ(x)i voorgesteld door vectoren met een eindig aantal elementen φ(k)i , k = 1...N . Dit impliceert dat de ruimte die kan opgespannen worden door deze vectoren φi (k) maximaal N dimensionaal kan zijn. Bijgevolg gebeuren de sommaties in wat volgt over een eindig aantal basisfuncties. Verder moeten de signalen ai (t) in de computer ook voorgesteld worden door vectoren met een eindig aantal elementen ai (t) waarbij t nu een discrete waarde voorstelt (t = 1...M ). Dit leidt tot de volgende nieuwe uitdrukking voor 5.2:

T (k, t) =

N X i

ai (t)φi (k), t = 1...M, k = 1...N

(5.3)


55

Bedoeling van de POD methode is het reduceren van het aantal basisfuncties φi , maar op een zodanige manier dat de expansie 5.3 nog nauwkeurig blijft. Op een meer wiskundige manier wordt dit als volgt geformuleerd: Probleem 5.1 Gegeven de expansie

N P

ai (t)φi (k) van de functie T (k, t) in de Hilbertruimte

i 0

χ waartoe T (k, t) behoort, vind T (k, t) = 0 d T (k, t) , T (k, t) minimaal is.

n P

ai (t)φi (k) waarvoor geldt dat n < N , zo dat

i

met de volgende definitie: 0 0 Definitie 5.1 d T (k, t) , T (k, t) = kT (k, t) − T (k, t) k 0

Er wordt dus een T (k, t) gezocht waarbij de sommatie in expansie 5.3 over n loopt met 1 P n < N . In het extreme geval dat n = 1 wordt expansie 5.3 dus ai (t)φi (k) = a1 (t)φ1 (k). i=1

De ruimtelijke component van de tijden plaatsafhankelijke temperatuur T (k, t) wordt in deze expansie dus enkel gegeven door φ1 . Deze eerste basisfunctie moet dus de ruimtelijke distributie van de temperatuur zo goed mogelijk weergeven en dit voor alle tijdstippen t = 1...M . Om dit te bekomen wordt de volgende extremumvergelijking opgesteld: (T (k, t), φ1 )2 φ1 = max φ1 ∈χ kφ1 k2

(5.4)

Waarbij h i staat voor het tijdsgemiddelde. De projectie van φ1 op de functie T (k, t) moet dus gemaximaliseerd worden en dit voor alle tijdstippen t. De oplossing van dit extremumprobleem wordt afgeleid in bijlage B. Het resultaat van de berekening is eenvoudig; φ1 moet voldoen aan volgende eigenwaardevergelijking:

h(T (k, t), φ1 )χ T (k, t)i = λ1 φ1

(5.5)

Zoals reeds eerder vermeld worden de basisfuncties φi bekomen uit gemeten of berekende data. Deze data werd verzameld in de snapshot matrix Tsnap die dus niets anders is dan een matrixvoorstelling van de functie T (k, t) geëvalueerd in een aantal punten k = 1...N op een aantal tijdstippen t = 1...M :




T (1, 1)

T (1, 2)

···

  Tsnap =  T (2, 1) T (2, 2) · · ·  .. ..  ..  . . .  T (N, 1) T (N, 2) · · ·

T (1, M ) T (2, M ) .. .

56

       

T (N, M )

Vergelijking 5.5 kan dus geschreven worden als: D

E (Tsnap , φ1 )χ Tsnap = λ1 φ1

(5.6)

Dit kan korter geschreven worden in operatorvorm, waarbij de operator C gedefinieerd wordt als C(φi ) = h(Tsnap , φi )χTsnap i:

C(φ1 ) = λ1 φ1

(5.7)

λ1 is dus een eigenwaarde van de operator C met bijhorende eigenvector φ1 die de eerste basisvector vormt van de expansie 5.3. De tweede basisfunctie φ2 wordt uit dezelfde extremumvergelijking 5.4 gevonden, maar met de restrictie dat φ2 orthonormaal moet zijn met φ1 . Ditzelfde geldt nu voor alle basisfuncties; ze voldoen allemaal aan de extremumvergelijking 5.4 en ze moeten allemaal orthonormaal zijn. Als nu de operator C zelftoegevoegd is, dan kan aangetoond worden dat al zijn eigenvectoren orthogonaal (en na normering orthonormaal) zijn en dat zijn eigenwaarden voldoen aan de ongelijkheid:

λ1 ≥ λ2 ≥ ... ≥ λN

(5.8)

De gezochte basisfuncties zijn dus eenvoudigweg de eigenvectoren operator C. Uit de manier waarop de basisfuncties φi geconstrueerd worden is het duidelijk dat de basisfunctie met de meeste informatie over de ruimtelijke distributie van de temperatuur, basisfunctie φ1 is. Basisfunctie φ2 is die met het op één na meeste informatie over de ruimtelijke distributie van de temperatuur enzovoort. Het is echter niet de bedoeling om alle gevonden basisfuncties φi te gaan gebruiken, maar slechts een deel ervan. Het feit dat niet alle φi ’s nodig zijn voor


57

een nauwkeurige voorstelling van de functie T (x, t) is juist het grote voordeel van de PODmethode. Er is echter een criterium nodig om te weten welke functies φi nog als belangrijk mogen beschouwd worden en welke niet. Hiervoor wordt gekeken naar de eigenwaarden λi die volgen uit eigenwaardevergelijking 5.5. Er werd reeds gezegd dat de relevantie van de basisfuncties φi afneemt bij toenemende index i. Samen met ongelijkheid 5.8 levert dit volgend criterium: n P

Pn =

i=1 N P

λi (5.9) λi

i=1

Als alle basisfuncties in gebruik worden genomen is n = N en wordt Pn gelijk aan 1. Hoe minder basisfuncties beschouwd worden, hoe minder nauwkeurig decompositie 5.3 en hoe lager Pn zal worden. Pn kan beschouwd worden als een maat voor de nauwkeurigheid van decompositie 5.3.

5.1.2

Bepalen van de tijdsafhankelijke co¨ effici¨ enten

Er wordt nu getracht een uitdrukking te vinden voor tijdsafhankelijke coëfficiënten ai (t). Definieer daartoe de volgende partiële differentiaalvergelijking voor de functie T (x, t) : X × T → R:

M (T ) = D(T )

(5.10)

Met D een differentiaaloperator die enkel afgeleiden naar de plaats bevat en M een differentiaaloperator die enkel afgeleiden naar de tijd bevat: P

∂ T M (T ) = M0 + M1 ∂T ∂t + ... + MP ∂tP ∂ 0 D(T ) = D0 + D1 ∂T ∂x + ... + DP

P

0

T

∂xP

(5.11)

0

Definieer nu de restoperator R van vergelijking 5.10 als volgt:

R(T ) = M (T ) − D(T )

(5.12)


58

Voor de functie T in vergelijking 5.10 is R(T ) dus duidelijk nul. Stel dat er reeds een PODbasis {φi }N i=1 voor de functie T (x, t) gevonden werd. Als de oplossing T (x, t) van differentiaalvergelijking 5.10 vervangen wordt door een gereduceerde lineaire combinatie van de basisfuncties φi , waarbij de sommatie over n loopt met n ≤ N :

Tn =

n X

ai (t)φi (x) ≈ T (x, t)

(5.13)

i=1

dan zal de restoperator R(Tn ) over het algemeen niet meer nul zijn. Om de originele partiële differentiaalvergelijking 5.10 zo goed mogelijk te benaderen wordt geprobeerd om de tijdsafhankelijke coëfficiënten ai (t) zodanig te kiezen dat de restoperator R(Tn ) nul wordt. Om dit na te streven wordt gebruik gemaakt van een zogenaamde Galerkin Projectie. Definitie 5.2 Gegeven een orthonormale basis {φi }N i=1 van een Hilbertruimte χ en een residuele operator R van het model M (T ) = D(T ). Als T benaderd wordt door een expansie van n P de vorm Tn (k, t) = ai (t)φi (k) dan resulteert een Galerkin projectie in: i=1

(R (Tn ) , φi ) = 0; i = 1...n

(5.14)

Deze definitie zegt dat de projectie van de restoperator R(Tn ) op het gereduceerde aantal basisfuncties {φi }ni=1 ; n ≤ N nul moet zijn. Deze projectie leidt tot een uitdrukking voor de tijdsafhankelijk coëfficiënten ai (t): 0

∂Tn ∂ P Tn ∂Tn ∂ P Tn 0 +...+MP −D −D −...−D R(Tn ) = M (Tn )−D(Tn ) = M0 +M1 (5.15) 0 1 0 P ∂t ∂tP ∂x ∂xP De restoperator R(Tn ) wordt nu geprojecteerd op de gereduceerde basis {φj }nj=1 en Tn wordt voluit geschreven:

(R(Tn ), φj (x)) = n n P P M0 ai (t) (φi (x), φj (x)) + M1

n P ∂ P ai (t) (φi (x), φj (x)) + ... + MP (φi (x), φj (x)) ∂tP i=1 i=1 i=1 0 n n n P P P ∂φi (x) ∂ P φi (x) 0 −D0 ai (t) (φi (x), φj (x)) − D1 ai (t) ai (t) , φj (x) 0 ∂x , φj (x) − ... − DP P i=1

∂ai (t) ∂t

i=1

i=1

∂x

=0 (5.16)


59

Rekeninghoudend met de orthonormaliteit van de basisfuncties: (φi , φj ) = δij leidt dit tot een differentiaalvergelijking voor de coëfficiënten ai (t): ∂ P a (t)

∂a (t)

j MP ∂tPj + ... + M1 ∂t + M0 aj (t) = i h P0 P ∂φi (x) ∂ φi (x) 0 , φ (x) nai (t) D0 + D1 , φ (x) + ... + D 0 j j P ∂x P

i=1

5.2

(5.17)

∂x

Praktijk

De analytische formules en algebra¨ısche basis zijn nu beschikbaar om dit in de praktijk toe te passen. Er bestaat een eenvoudige methode die toelaat de beschreven theorie in praktijk om te zetten [13]. Daartoe wordt een eerste orde partiële differentiaalvergelijking opgesteld die het tijdsafhankelijke temperatuursprofiel van de olie beschrijft. Vervolgens wordt van deze vergelijking een discreet toestandsmodel opgesteld. Dit toestandsmodel genereert het tijdsafhankelijk temperatuursprofiel in het systeem en kan dus gebruikt worden om de waarden voor Tsnap te berekenen. Uit de uitdrukking voor het toestandsmodel kunnen verder de tijdsafhankelijk coëfficiënten ai (t) gevonden worden. Deze methode impliceert wel dat de resultaten via de POD-methode nooit beter kunnen zijn dan de resultaten van het toestandsmodel. De reden dat de POD-methode gebruikt wordt, is dat het aantal berekeningen sterk gereduceerd wordt. Er wordt opnieuw ondersteld dat het volledige buizensysteem is opgedeeld in tien parallelle leidingen zoals te zien in figuur 4.3. Verder wordt opnieuw een onderscheid gemaakt tussen de toevoer-, verhittings- en afvoersectie. Het is belangrijk dat hier ondersteld wordt dat de leidingen dezelfde diameter hebben in elk van deze drie secties. In de methode waarbij nietuniforme tijdsmonsters gebruikt werden, was dit niet het geval (zie tabel 4.4). Daar was het ook niet nodig omdat enkel het volume van elk segment gelijk moest zijn. In de methode waarbij POD-technieken gehanteerd worden, moeten alle segmenten hetzelfde volume en dezelfde lengte ∆x hebben. Dit laatste volgt uit het feit dat er bij de discretizering slechts één ruimtelijke discretisatie geformuleerd wordt. Dit impliceert dat de diameter van elk segment gelijk moet zijn. Daardoor wordt de diameter van de leidingen in de toevoer- en afvoersectie gelijkgesteld aan de diameter van de leidingen in de verhittingssectie. Dit zal de lengte van de leidingen in toevoer- en afvoersectie veranderen, maar dit is niet erg want hun volume blijft


60

ongewijzigd. De nieuwe lengte voor leiding i in de toevoersectie (en een volledige analoge formule voor de afvoersectie) is eenvoudig te berekenen volgens de formule:

Vitoev

Ltoev nieuwi

2 dtoev = π 9 2 verhit 2 d toev = Lnieuwi π 2 ⇓ verhit 2 toev 1 d = Li 9 (dtoev )2 1 Ltoev i

(5.18)

met: • Vitoev : het volume van leiding i in de toevoersectie in m3 • Ltoev : de lengte van leiding i in de toevoersectie in m i • dtoev : de werkelijke diameter van de toevoersectie in m • Ltoev nieuwi : de aangepaste lengte van leiding i in de toevoersectie in m • dverhit : de werkelijke diameter van een leiding in de verhittingssectie in m

Voor de eenvoud wordt de methode eerst enkel toegepast op leiding 2. In een volgende stap wordt overgegaan op de uitbreiding naar de volledige centrale.

5.2.1

Vergelijking voor het temperatuursprofiel

Voor het berekenen van een vergelijking die het temperatuursprofiel langsheen het systeem beschrijft, wordt vertrokken van de wet van behoud van energie [14]:

ρf cf

∂2 DT (x, t) = k 2 T (x, t) + S(x, t) Dt ∂x

• ρf : dichtheid van de olie in kg/m3 • cf : specifieke warmtecapaciteit van de olie in J/kg K

(5.19)


61

• T (x, t): plaatsen tijdsafhankelijk temperatuursprofiel van de olie • k: warmtegeleidingscoëfficiënt van de olie in W /m K • S(x, t): plaatsen tijdsafhankelijke warmtebron in W /m3

Hierin stelt

D Dt

=

∂ ∂t

∂ + vx (t) ∂x een substantiële afgeleide voor. Dit is een tijdsafgeleide die

rekening houdt met de verplaatsing van de olie (vx (t) is de tijdsafhankelijke snelheid van de olie in m/s). De warmtegeleidingscoëfficiënt k wordt opnieuw verwaarloosd en de warmtebron is in dit geval het geconcentreerde zonlicht Irr(t) (in W /m2 ). De differentiaalvergelijking 5.19 is enkel geldig voor de verhittingssectie. In de toevoer- en afvoersectie is immers geen warmtebron aanwezig. In deze twee secties wordt volgende vergelijking gebruikt:

ρf cf

DT (x, t) Dt

= ρf cf

∂T (x, t) ∂T (x, t) + ρf cf vx (t) ∂t ∂x

= 0

(5.20)

Het zonlicht wordt gecollecteerd door parabolische spiegels die een apertuur Amirr en een reflectiecoëfficiënt hebben. Vergelijkingen 5.19 en 5.20 zijn geldig per m3 en worden dus vermenigvuldigd met het volume van de respectievelijke secties: Vtoev ,Vverhit en Vaf v . Dit leidt tot:

∂T (x, t) ∂T (x, t) + ρf cf Vtoev vx (t) ∂t ∂x ∂T (x, t) ∂T (x, t) ρf cf Vverhit + ρf cf Vverhit vx (t) ∂t ∂x ∂T (x, t) ∂T (x, t) ρf cf Vaf v + ρf cf Vaf v vx (t) ∂t ∂x ρf cf Vtoev

= 0

(5.21)

= Irr(t)Amirr

(5.22)

= 0

(5.23)

Waarbij vergelijkingen 5.21, 5.22 en 5.23 geldig zijn voor respectievelijk de toevoer-, de verhittings- en de afvoersectie. Verder kan elk van de drie volumes Vtoev ,Vverhit en Vaf v opgesplitst worden als Vi = A.Li waarbij A (m2 ) de dwarsdoornede en Li (m) de lengte van leiding 2 in sectie i voorstelt (i =toev, verhit of avf). De diameter in elk van de drie secties werd gelijk gekozen dus is de dwarsdoorsnede A in elk van de drie secties ook gelijk. Het product van vx (t) met A vormt een uitdrukking voor het debiet q(t) van de olie in loop 2. Dit leidt tot de volgende vormen voor de partiële differentiaalvergelijking:


∂T (x, t) ∂T (x, t) + ρf cf Ltoev q(t) ∂t ∂x ∂T (x, t) ∂T (x, t) ρf cf Vverhit + ρf cf Lverhit q(t) ∂t ∂x ∂T (x, t) ∂T (x, t) ρf cf Vaf v + ρf cf Laf v q(t) ∂t ∂x ρf cf Vtoev

62

= 0

(5.24)

= Irr(t)Amirr

(5.25)

= 0

(5.26)

Er wordt een randvoorwaarde gedefinieerd voor het begin van het systeem, aan de ingang van de toevoersectie:

T (0, t) = Tin (t)

(5.27)

De ruimtelijke distributie van de temperatuur in het systeem is niet gekend op tijdstip t = 0, daarom wordt ze gelijk gesteld aan een constante, nl. de ingangstemperatuur op tijdstip t = 0:

T (x, 0) = Tin (0)

5.2.2

(5.28)

Opstellen van het toestandsmodel

Vergelijkingen 5.24, 5.25 en 5.26 dienen nu gediscretizeerd te worden om een discreet toestandsmodel te krijgen. Voor de discretizering wordt beroep gedaan op de ’Implicit Backward Euler Method’[15]. Dit leidt tot de volgende uitdrukkingen:

t Tkt − Tk−1 Tkt − Tkt−1 c31 + c21 =0 ∆t ∆x t T t − Tk−1 T t − Tkt−1 c32 k + c22 k = c1 Irrkt ∆t ∆x t T t − Tk−1 T t − Tkt−1 c33 k + c23 k =0 ∆t ∆x

voor k = 1...N1

(5.29)

voor k = N1 + 1...N1 + N2

(5.30)

voor k = N1 + N2 + 1...N1 + N2 + N3 (5.31)

met • Tkt : temperatuur in segment k op tijdstip t in K • c1 = Amirr

Hoofdstuk 5. Methode II: Proper Orthogonal Decomposition • c21 = ρf cf Ltoev

c22 = ρf cf Lverhit

c23 = ρf cf Laf v

• c31 = ρf cf Vtoev

c32 = ρf cf Vverhit

c33 = ρf cf Vaf v

63

• N1 , N2 en N3 het aantal ruimtelijke discretizatiestappen in respectievelijk de toevoer-, verhittings- en afvoersectie

Door de verschillende termen te herschikken in vergelijking 5.29 wordt een uitdrukking bekomen voor Tkt−1 :

c21 ∆t t Tkt − Tk−1 c31 ∆x t t = Tk + at11 Tkt − Tk−1

Tkt−1 = Tkt +

t = (1 + at11 )Ttk − at11 Tk−1

(5.32)

Voor vergelijkingen 5.30 en 5.31 volgen analoge resultaten:

t + at2 Irrt Tkt−1 = (1 + at12 )Ttk − at12 Tk−1

(5.33)

t Tkt−1 = (1 + at13 )Ttk − at13 Tk−1

(5.34)

met • at11 =

c21 ∆t t c31 ∆x q

at12 =

c22 ∆t t c32 ∆x q

at13 =

c23 ∆t t c33 ∆x q

∆t • at2 = − c12 c32

Om T1t te berekenen wordt randvoorwaarde 5.27 gebruikt: T1t−1 = (1 + at11 )T1t − at11 T0t + a2 Irrt

(5.35)

waarbij T0t verwijst naar de ingangstemperatuur Tin van het systeem. Vergelijkingen 5.32, 5.33 en 5.34 kan in matrixvorm gegoten worden:


 t 1 + a11   −1 1 + at11   .. ..  . .    −1 1 + at12   ..  .  

64

  t−1 t t    0  T1 + a1 T0           0    T2t T2t−1          .    . . .. ..    ..           + =   Tt  a Irrt   T t−1     2    N1 +1 N1 +1    .    .. . ..        .. . .    ..          −1 1 + 1at13 TNt 1 +N2 +N3 0 TNt−1 1 +N2 +N3 (5.36) 

T1t

 



De matrices kunnen wat beter gerangschikt worden:  t 1 + a11   −1 1 + at11   .. ..  . .    −1 1 + at12   ..  .  





    t−1 t t   T1 + a1 T0   0             0  T2t T2t−1            .  .. .. .       . .      .   = −   Tt   T t−1  a Irrt       2  N1 +1 N1 +1      .  .. .. ..        . . .      ..        t−1 t t −1 1 + 1a13 TN1 +N2 +N3 TN1 +N2 +N3 0 (5.37) T1t

Dit kan nog korter geschreven worden als:

A(t)T(t) = T(t − 1) + u(t) ⇓ T(t) = A−1 (t)T(t − 1) + A−1 (t)u(t)

5.2.3

(5.38)

Berekenen van de basisfuncties

Er is dus een uitdrukking beschikbaar om het temperatuursprofiel in het systeem op elk tijdstip t = 1...M te berekenen. Hiermee kan de snapshotmatrix Tsnap opgesteld worden. Eens deze gekend is kan vergelijking 5.6 berekend worden. Dit leidt tot de volgende uitdrukking:

C(φi ) = h(Tsnap , φi )χTsnap i =

1 1 Tsnap T φi Tsnap = Tsnap Tsnap T φi = λi φi M M

(5.39)


65

De operator C kan dus zeer eenvoudig berekend worden aan de hand van de matrixvermenigvuldiging:

C=

1 Tsnap Tsnap T M

(5.40)

Het totaal aantal segmenten N1 + N2 + N3 waarin leiding 2 wordt opgedeeld wordt vanaf nu N genoemd. Op de N × N matrix C wordt nu een eigenwaardendecompositie uitgevoerd, wat leidt tot N eigenwaarden λi en N bijhorende eigenvectoren φi . In figuur 5.1 zijn als voorbeeld de eerste drie basisfuncties zichtbaar voor de gegevens van 15 juli:

Figuur 5.1: De eerste drie basisfuncties voor de gegevens van 15 juli.

De essentie van de POD-methode is dat er slechts een klein aantal basisfuncties φi gebruikt worden. Er werd reeds een criterium opgesteld om te bepalen hoeveel basisfuncties nodig zijn om nog nauwkeurige resultaten te krijgen (zie vergelijking 5.9). Als voorbeeld worden in figuur 5.2 de eigenwaarden λi van C geplot voor de gegevens van 15 juli: Als een waarde voor Pn gekozen wordt, kan berekend worden hoeveel basisfuncties er nodig zijn.


66

Figuur 5.2: Logaritmische plot van de eigenwaarden.

5.2.4

Bepalen van de tijdsafhankelijk co¨ effici¨ enten

De tijdsafhankelijke coëfficiënten ai (t) worden nu bepaald via de Galerkin projectie. Daartoe worden de vectoren T(t) en T(t − 1) uit vergelijking 5.38 vervangen door hun decompositie in basisfuncties:

T(t) → T(t − 1) →

n X i=1 n X

ai (t)φi (k) = Φn a(t) = Tn (t) ai (t − 1)φi (k) = Φn a(t − 1) = Tn (t − 1)

(5.41)

i=1

De differentieverelijking 5.38 kan dus geschreven worden als:

Tn (t) = A−1 (t)Tn (t − 1) + A−1 (t)u(t)

(5.42)

Vervolgens wordt de restoperator R(Tn ) van vergelijking 5.42 gedefinieerd: R(Tn ) = Tn (t) − A−1 (t)Tn (t − 1) − A−1 (t)u(t)

(5.43)


67

Uit definitie 5.2 volgt dat de projectie van R(Tn ) van op de set van basisfuncties Φn nul moet zijn. Dit geeft, rekeninghoudend met 5.41:

Φn T Φn a(t) = Φn T A−1 (t)Φn a(t − 1) + Φn T A−1 (t)u(t)

(5.44)

Aangezien de basisfuncties orthonormaal zijn is het matrix product Φn T Φn gelijk aan 1. Φn T A−1 (t)Φn wordt in wat volgt B(t) genoemd en Φn T A−1 (t) wordt vanaf nu korter geschreven als D(t). Dit leidt tot de zeer eenvoudige uitdrukking voor de tijdsafhankelijk coëfficiënten ai (t):

a(t) = B(t)a(t − 1) + D(t)u(t)

5.2.5

(5.45)

Uitbreiding naar de volledige centrale

De voorgaande afleiding is slechts geldig voor één enkele leiding uit het model. Het is de bedoeling dit uit te breiden naar de volledige centrale. Daartoe wordt de centrale gemodelleerd als één enkele leiding, waarvan het volume hetzelfde is als de som van alle leidingen in de centrale. Opnieuw wordt de leiding opgesplitst in een toevoer-, verhittings- en afvoersectie. Om een model van deze nieuwe geometrie op te stellen, kan vertrokken worden van vergelijkingen 5.21, 5.22 en 5.23, mits aanpassing van de parameters Vtoev , Vverhit , Vaf v en Amirr . Er wordt nu immers rekening gehouden met alle spiegels en niet enkel met die van loop 2. Amirr moet dus vermenigvuldigd worden met 9, het aantal loops dat tijdens het experiment in gebruik werd genomen. Omdat alle leidingen in de verhittingssectie hetzelfde volume hebben, kan Vverhit ook gewoon met 9 vermenigvuldigd worden, om dezelfde reden. Vtoev is voor elk van de 9 leidingen verschillend, daarom wordt het gemiddelde van de negen volumes in de toevoersectie genomen. Dit komt overeen met het volume van leiding 6 in de toevoersectie. Vervolgens wordt dit volume ook vermenigvuldigd met 9. Dezelfde redenering geldt voor Vaf v ; als nieuwe waarde voor deze parameter wordt negen keer het volume van leiding 6 in de afvoersectie genomen. De waarde voor de bemonsteringstijd wordt dezelfde gekozen als de bemonsteringstijd van de data, namelijk 30 s. Als ruimtelijke discretisatiestap wordt 25 m genomen. Dit bleek door


68

’gis en mis’ de beste waarde te zijn. Deze waarde voor de discretisatiestap zorgt ervoor dat de leiding wordt opgedeeld in 22 segmenten. De Hilbertruimte van de basisfuncties heeft dus dimensie 22 en er kunnen bijgevolg 22 orthonormale basisfuncties gevonden worden.

5.3

Resultaten

De basisfuncties moeten maar één keer opgesteld worden. Ze worden best berekend door gebruik te maken van een experiment waarbij het systeem goed geëxciteerd wordt (dit is: er worden grote variaties aangebracht in het ingangssignaal en de storingen). Op deze manier komt veel informatie in de basisfuncties terecht. Dit is de reden dat de basisfuncties berekend worden aan de hand van de gegevens van 28 juli, aangezien er die dag grote variaties in het debiet waren. Er moet nog besloten worden hoeveel basisfuncties er in rekening gebracht worden. Daarvoor werd een criterium opgesteld, zie vergelijking 5.9. De nauwkeurigheid van de POD-methode werd er voorgesteld door de parameter Pn ; hoe dichter Pn bij 1 hoe nauwkeuriger. Daarom wordt als waarde voor Pn 1 − 10−5 genomen. Voor de data van 28 juli komt dit overeen met 4 basisfuncties.

5.3.1

15 juli

In figuur 5.33 zijn de resultaten zichtbaar voor de temperatuur op het eind van de afvoersectie berekend via de POD-methode. Ze maken zoals gezegd gebruik van eerste de vier basisfuncties berekend met gegevens van 28 juli. 3

Voor een groter formaat van de figuren wordt verwezen naar bijlage C


zonlicht (W/m²)

1000

69

10 debiet (l/s)

800 600

5

400

50

100

150 200 250 tijdsmonster

300

350

160

0

50

100


300

350

400

250

150 140 130 120 0

50

100


300

350

400

10

fout in %

0

400



200 0

200

150

5 reele temperatuur gesimuleerde temperatuur 0

0

50

100


300

350

400

100 0

50

100


300

350

400

Figuur 5.3: Resultaten met de POD-methode voor 15 juli.

De resultaten blijken vrij goed te zijn. De fout komt nergens boven de 5% behalve tijdens de beginwaarden en helemaal op het einde. De reden voor de afwijkende beginwaarden is dezelfde is als bij de vorige simulator. De temperatuursdistributie van de olie in het systeem is niet gekend op de eerste tijdstippen. De reden voor de afwijkende eindwaarden is eveneens dezelfde als bij de vorige simulator: op het eind van elke dag worden de spiegels naar beneden gericht om ze te beschermen tegen stof en regen. Daardoor neemt de temperatuur van de olie sterk af op het einde. Dit verdraaien van de spiegels wordt niet gesimuleerd. De berekende temperaturen zijn verder ook goed bij lage waarden voor het debiet. Dit in tegenstelling tot de vorige simulator. Verder is de tijdsshift die in de eerste versie van de vorige simulator opdook hier afwezig. Dit komt omdat de experimenteel bepaalde lengte van de afvoersectie (zie sectie 4.6) gebruikt werd. Een laatste opmerking: de simulatie gebeurt zeer snel en neemt


70

ongeveer één seconde in beslag in vergelijking met ongeveer 30 s bij de vorige simulator. Dit werd ook verwacht, aangezien het de bedoeling is van POD om het aantal berekeningen sterk te reduceren. Er wordt nogmaals benadrukt dat de resultaten via de POD-methode niet beter kunnen zijn dan de resultaten van het toestandsmodel. Voor een vergelijking beide resultaten voor alledrie de testdagen wordt verwezen naar bijlage C. De afwijkingen die voorkomen tussen de reële temperatuur en de temperatuur berekend via POD zijn dus vooral te wijten aan de nauwkeurigheid van het toestandsmodel.

5.3.2

16 juli

In figuur 5.4 zijn de resultaten van de POD-methode zichtbaar voor de gegevens van 16 juli.

debiet (l/s)

10

800

600


400 0

100

200


500

200

0

100

200


500

600

240 220

150

200 100

50

0

100

200


500

600

10

fout in %

5

0

600

5


zonlicht (W/m²)

1000

180 160 140 120 100 80 reele temperatuur gesimuleerde temperatuur

60 0

0

100

200


500

600

40

0

100

200


500

600


Zelfde opmerkingen voor de beginwaarden als bij de vorige figuur. Daarnaast is er ook een


71

opmerkelijk grote fout te zien in de eerste helft van het experiment. De reden hiervoor wordt later onderzocht. De tweede helft van de simulatie ziet er redelijk goed uit. De fout blijft steeds een stuk onder de 5%.

5.3.3

28 juli

In figuur 5.5 worden de resultaten getoond van de simulatie van 28 juli.

zonlicht (W/m²)

1000

10 debiet (l/s)

800 600

5

400

100

200


500

600

200

0

100

200


500

600

700

300

150 250 100 50 0

0

100

200


500

600

700

10

fout in %

0

700



200 0

200

150

100

5 50 reele temperatuur gesimuleerde temperatuur 0

0

100

200


500

600

700

0

0

100

200


500

600

700


Opnieuw zijn er sterk afwijkende resultaten voor ruim de eerste helft van de simulatie, maar dit wordt verderop onderzocht. In figuur 5.5 is ook te zien dat de simulator goed reageert op een storing in de ingangstemperatuur. De daling in de werkelijke uitgangstemperatuur, veroorzaakt door een daling in de ingangstemperatuur, is ook aanwezig in de gesimuleerde uitgangstemperatuur (omcirkeld).


5.3.4

72

Conclusies

Vooreerst is de POD-simulator veel sneller dan de vorige. Het duurt ongeveer 1 seconde om alle berekeningen uit te voeren, terwijl de eerste simulator er ongeveer een halve minuut over deed. Verder bleken de resultaten het best voor 15 juli. De fout bleef er telkens onder de 5%, behalve in het begin en op het einde. Er werden ook goede resultaten geboekt voor lage waarden van het debiet, wat minder het geval was met de vorige simulator. Voor 16 en 28 juli waren de resultaten van de simulator sterk afwijkend. Enkel voor de tweede helft van beide experimenten was er een overeenkomst tussen de waarden berekend door de simulator en de werkelijke waarden. De reden hiervoor is waarschijnlijk te wijten aan het feit dat de opwarming van de olie via twee mechanismen afhankelijk is van het debiet. Een eerste mechanisme is vrij logisch; hoe hoger het debiet hoe minder tijd de olie in de verhittingssectie doorbrengt dus hoe minder de olie zal opwarmen. Hiermee wordt rekening gehouden in de POD-methode via het gebruik van de substantiële afgeleide

D Dt .

Deze houdt rekening met het debiet van de

olie. Een tweede mechanisme is minder voor de hand liggend en werd reeds aangehaald in het opstellen van de eerste simulator. Daar werd een warmtetransfercoëfficiënt Ht ingevoerd die op een niet-lineaire manier afhankelijk was van het debiet. Deze kan beschouwd worden als een convectiecoëfficiënt: hoe sneller de olie langs de metalen wand van de pijpleiding vloeit, hoe meer warmte de olie zal opnemen. Dit werd vergeleken met een ventilator die gebruikt wordt om processors sneller af te koelen. Dit mechanisme heeft dus een omgekeerd effect op de opwarming van de olie in vergelijking met het vorige mechanisme. In het toestandsmodel dat werd opgesteld voor de tweede simulator werd de temperatuur van de metalen pijpleiding niet gesimuleerd en ook de warmtetransfercoëfficiënt Ht werd niet gedefinieerd. Aangezien de POD-simulator opgesteld wordt uit het toestandsmodel, zal deze de parameter Ht ook niet bevatten. Dus de gesimuleerde temperatuur zal bij lage debietwaarden waarschijnlijk te hoog zijn. Dit kan een mogelijke verklaring zijn voor de afwijkingen in de eerste helft van de simulaties van 16 en 28 juli.

Hoofdstuk 6

Besluit In dit eindwerk werd op twee manieren een simulator ontworpen voor het ACUREX-veld op het Plataforma Solar de Almeria. De eerste maakte gebruik van een gedistribueerd parameter model met niet-uniforme tijdsmonsters. De reden voor het gebruik van niet-uniforme tijdsmonsters was het variërend debiet. Er werd vertrokken van twee vergelijkingen die enerzijds de temperatuur van de metalen leiding en anderzijds de temperatuur van de olie bepaalden. Er werden enkele vereenvoudigende veronderstellingen gemaakt om het model niet te complex te maken. De resultaten bleken de werkelijkheid vrij goed te benaderen. Er was wel telkens een probleem in het begin en op het eind van de simulatie. Het eerste probleem werd toegeschreven aan het feit dat de initiële temperatuursdistributie in het systeem niet gekend is. Het tweede probleem was een gevolg van het feit dat de spiegels op het einde van het experiment naar beneden werden gericht. Verder leken er ook problemen op te treden bij lage waarden van het debiet. De reden hiervoor is waarschijnlijk een slechte uitdrukking voor de warmtetransfercoëfficiënt Ht . Een voordeel van de eerste methode is dat alle parameters een duidelijke fysische betekenis hebben. Voor sommige diende de waarde op een experimentele manier bepaald te worden, zoals voor Ht , Hl en de lengte van de afvoersectie. Een nadeel van de eerste methode is dat het aantal berekeningen relatief groot is; het simuleren van de gegevens van 15, 16 en 28 juli duurt gemiddeld ongeveer 30 s. De tweede simulator steunde op POD-technieken. Deze maken gebruik van orthonormale

73

Hoofdstuk 6. Besluit

74

basisfuncties die de ruimtelijke distributie van een belangrijke parameter, in dit geval de temperatuur, beschrijven. De basisfuncties worden vermenigvuldigd met tijdsafhankelijke coëfficiënten om de tijdsafhankelijkheid van het systeem weer te geven. Er werd een discreet toestandsmodel afgeleid uit de vergelijking van behoud van energie. Dit toestandsmodel leverde genoeg data om de basisfuncties te kunnen berekenen. Uit dit toestandsmodel konden verder op een eenvoudige manier de tijdsafhanklijke coëfficiënten bepaald worden. De resultaten waren het best voor 15 juli. Er was telkens een goede overeenkomst tussen de werkelijk en de gesimuleerde waarde, behalve in het begin en op het einde. De reden hiervoor is dezelfde als bij de vorige simulator. Voorts waren de resultaten ook goed voor lage waarden van het debiet. Voor 16 en 28 juli daarentegen, bleek er slechts vanaf halverwege het experiment een goede overeenkomst te bestaan tussen de resultaten van de simulator en de werkelijke waarden. Een mogelijke oorzaak hiervan werd toegeschreven aan het ontbreken van een niet-lineaire warmtetransfercoëfficiënt Ht . Een voordeel van deze methode is het lage aantal berekeningen. Daar vooral de eerste basisfuncties het meest representatief zijn voor de temperatuursdistributie in het systeem, moeten ook enkel die in rekening gebracht worden. Het lage aantal berekeningen heeft als gevolg dat de simulaties zeer snel kunnen gebeuren. Zo duurde het simuleren van de gegevens van 15, 16 en 28 juli gemiddeld ongeveer 1 s. Verder moeten de basisfuncties slechts eenmaal berekend worden. Hiervoor wordt best gebruik gemaakt van gegevens verzameld tijdens een experiment waarbij het systeem voldoende geëxciteerd werd, bijvoorbeeld dat van 28 juli. Een nadeel van de POD-methode is dat de simulator minder goed reageert op data van twee van de drie testdagen. De simulator kan waarschijnlijk verbeterd worden door in het toestandsmodel een warmtetransfercoëfficiënt in te voeren die op een niet-lineaire manier van het debiet afhangt.

Bijlage A

Dode tijden In de onderstaande tabel staan de dode tijden die optreden in de verschillende leidingen van de drie secties, voor alle waarden van het debiet. De waarden werden berekend via de formule τ = 1000 Vq met τ de dode tijd in s, V het volume in m3 , q het debiet in l/s en 1000 een omrekenfactor om het volume om te zetten van m3 naar l. Ze zijn dus niet experimenteel maar theoretisch.

(s) debiet

Toevoersectie loop 1

loop 2

loop 3

loop 4

loop 5

loop 6

loop 7

loop 8

loop 9

loop 10

2

138

177

216

255

293

332

371

410

449

488

3

92

118

144

170

196

222

247

273

299

325

4

69

88

108

127

147

166

186

205

224

244

5

55

71

86

102

117

133

148

164

180

195

6

46

59

72

85

98

111

124

137

150

163

7

39

51

62

73

84

95

106

117

128

139

8

38

49

60

71

82

92

103

114

125

136

9

34

44

54

64

73

83

93

103

112

122

(l/s)

Tabel A.1: Dode tijden voor de toevoersectie.

75

Bijlage A. Dode tijden

76

De volumes van alle loops in de verhittingssectie zijn dezelfde. Daardoor is de dode tijd voor elke loop ook dezelfde. (s)

Verhittingssectie

debiet

loop i; i = 1...10

(l/s) 2

514

3

343

4

257

5

206

6

171

7

147

8

128

9

114

Tabel A.2: Dode tijden voor de verhittingssectie.

(s) debiet

Afvoersectie loop 1

loop 2

loop 3

loop 4

loop 5

loop 6

loop 7

loop 8

loop 9

loop 10

2

610

648

687

726

765

804

843

882

921

959

3

406

432

458

484

510

536

562

588

614

640

4

305

324

344

363

383

402

421

441

460

480

5

244

259

275

290

306

322

337

353

368

384

6

203

216

229

242

255

268

281

294

307

320

7

174

185

196

207

219

230

241

252

263

274

8

152

162

172

182

191

201

211

220

230

240

9

135

144

153

161

170

179

187

196

205

213

(l/s)

Tabel A.3: Dode tijden voor de afvoersectie.

Bijlage A. Dode tijden

77

Als laatste worden de dode tijden voor gans het systeem (toevoer-, verhittings- en afvoersectie samen) gegeven per leiding.

(s) debiet

Toevoer-, verhittings- en afvoersectie loop 1

loop 2

loop 3

loop 4

loop 5

loop 6

loop 7

loop 8

loop 9

loop 10

2

1261

1339

1417

1495

1572

1650

1728

1806

1883

1961

3

841

893

945

996

1048

1100

1152

1204

1256

1307

4

631

670

708

747

786

825

864

903

942

981

5

505

536

567

598

629

660

691

722

753

784

6

420

446

472

498

524

550

576

602

628

654

7

360

383

405

427

449

471

494

516

538

560

8

315

335

354

374

393

413

432

451

471

490

9

280

298

315

332

349

367

384

401

419

436

(l/s)

Tabel A.4: Dode tijden voor gans het systeem.

Bijlage B

Wiskundige bijlage B.1

Definities

Definitie B.1 Een deelverzameling V van een genormeerde lineaire vectorruimte is dicht in χ als zijn compactum gelijk is aan χ. Dit betekent dat voor elk element f van χ een v ∈ V gevonden kan worden, zodanig dat voor elke > 0 geldt dat kf − vk < .

Definitie B.2 Een genormeerde lineaire vectorruimte χ is separabel als het een aftelbaar dichte deelverzameling bevat. Een voorbeeld van een aftelbare verzameling is de verzameling van de gehele getallen en een voorbeeld van een onaftelbare verzameling is de verameling van de reële getallen. Definitie B.3 Een rij {fn }∞ n=1 in een genormeerde ruimte χ is een Cauchy rij als

kfn − fm kχ → 0 voor n, m → ∞

Definitie B.4 Een genormeerde, lineaire ruimt χ is compleet als elke Cauchy rij een ophopingspunt heeft in χ.

Definitie B.5 Een Hilbertruimte is een ruimte waarin het lineair product gedefinieerd is en die compleet is. 78

Bijlage B. Wiskundige bijlage

B.2

79

Oplossing extremumprobleem (T (k, t), φ1 )2 φ1 = max φ1 ∈χ kφ1 k2

(B.1)

Definieer om bovenstaand extremumprobleem op te lossen de volgende Lagrangiaan J(φ1 ) met de Lagrangiaanse vermenigvuldiger λ1 : D E J(φ1 ) = |T(t), φ1 |2 − λ1 kφ1 k2 − 1 T

(B.2)

Om het extremumprobleem op te lossen wordt gebruik gemaakt van variatierekening en wordt φ1 vervangen door: φ1 → φ1 + δψ, δ ∈ R; ψ ∈ χ met ψ een arbitraire functie in de Hilbert ruimte χ. Een voorwaarde opdat φ1 een extremum zou zijn van is dat: d J(φ1 + δψ) |δ=0 = 0 dδ

(B.3)

Vergelijking voluit schrijven levert:

d J(φ1 + δψ) |δ=0 = dδ

d

| (T(t), φ1 ) |2 T − λ1 kφ1 k2χ − 1 dδ E i d hD = (T(t), φ1 + δψ)χ (T(t), φ1 + δψ)χ dδ T i d h − λ1 (φ1 + δψ, φ1 + δψ)χ |δ=0 dδ

(B.4)

met (, )χ het inwendig product in de Hilbert ruimte χ. Vergelijking kan verder uitgewerkt worden tot:

Bijlage B. Wiskundige bijlage

80

d d J(φ1 + δψ) |δ=0 = − [λ1 ((φ1 , φ1 ) + (φ1 , δψ) + (δψ, φ1 ) + (δψ, δψ))] |δ=0 dδ dδ d + [h((T(t), φ1 ) + (T(t), δψ)) ((T(t), φ1 ) + (T(t), δψ))i] |δ=0 dδ d h(T(t), φ1 ) (T(t), φ1 ) + (T(t), φ1 ) (T(t), δψ) + (T(t), δψ) (T(t), φ1 )i |δ=0 = dδ d − λ1 ((φ1 , φ1 ) + (φ1 , δψ) + (δψ, φ1 ) + (δψ, δψ)) |δ=0 dδ (B.5)

Door de termen die niet functie zijn van δ zoals (T(t), φ1 ), gelijk te stellen aan nul en door alles te herschikken wordt volgende uitdrukking bekomen: d J (φ1 + δψ) |δ=0 = 2 h(T(t), φ) (T(t), ψ)i − 2λ1 (φ1 , ψ) dψ

(B.6)

Deze laatste variatie gelijkstellen aan nul levert de uitdrukking:

h(T(t), φ1 ) (T(t), ψ)i − λ1 (φ1 , ψ1) = 0

(B.7)

Gebruik makend van de commutativiteit van de operatoren h i en (, ) vinden we:

(h(T(t), φ1 ) (T(t))i − λ1 φ1 ) ψ = 0

(B.8)

En omdat moet gelden voor elke willekeurige functie ψ geldt er:

h(T(t), φ1 ) (T(t))i − λ1 φ1 = 0 h(T(t), φ1 ) (T(t))i = λ1 φ1 (B.9)

Bijlage C

Figuren C.1

Figuren bij hoofdstuk 4

In figuur C.1 (elk tijdsmonster meet 30s) is de oorzaak van het dipje in figuur zichtbaar. Er is een grote daling merkbaar in de temperatuur op het eind van loop 9. Dit kan veroorzaakt worden doordat de parabolische spiegels van deze loop gedurende enkele ogenblikken van de zon werden weggedraaid.

81

Bijlage C. Figuren

82

1000

10

900

debiet (l/s)

zonlicht (W/m²)

8 800 700 600

6 4

500 2

400 300 0

50

100


300

350

0

400

50

100


50

100


300

350

400

250 uitgangstemperatuur (°C)


160

0

150

140

130

200

150

100 reële temperatuur

120 0

50

100


300

350

400

0

300

350

400

Figuur C.1: Gegevens van 15 juli met i.h.b. rechtsonderaan het temperatuursverloop op het eind van loop 9.

In figuur C.2 is de invloed van de ingangstemperatuur op de uitgangstemperatuur merkbaar. Er is een plotse verlaging van de ingangtemperatuur rond het 310ste tijdsmonster (elk tijdsmonster meet 30s). Deze dip komt ook voor in de uitgangstemperatuur en kan niet te wijten zijn aan een debietsverhoging aangezien het debiet gedurende dit tijdstip constant blijft. Het minimum in de dip in de ingangstempetuur komt voor op het 318de tijdsmonster terwijl die in de eindtemperatuur voorkomt op het 338ste tijdsmonster. Dit komt overeen met een dode tijd van 20 30s = 600s. Deze waarde ligt tussen de dode tijd van loop 1 en loop 10 bij een debiet van 5 l/s in tabel A.4.

Bijlage C. Figuren

83

1000

10

900

debiet (l/s)

zonlicht (W/m²)

8 800 700 600

6 4

500 2

400 300 0

100

200


500

600

0

700

150

100

50

0

0

100

200


500

600

700

100

200


500

600

700

250 uitgangstemperatuur (°C)


200

0

100

200


500

600

700

200

150

100

0

Figuur C.2: Gegevens van 28 juli: invloed van de ingangstemperatuur.

In de twee onderstaande figuren is de invloed van Tklok te zien. In de eerste grafiek werd een Tklok (en dus tijdsmonster) van 2 s gebruikt. Ter vergelijking werd in de figuur eronder een Tklok (en dus tijdsmonter) van 1 s gebruikt.

Bijlage C. Figuren

84

debiet (l/s)

10

800

600


400 0

2000


8000

200

0

2000


8000

10000

260 240

150

220 100

50

0

2000


8000

10000

10

fout in %

5

0

10000

5


zonlicht (W/m²)

1000


80 0

0

2000


8000

10000

60

0

2000


8000

Figuur C.3: Gegevens van 16 juli met Tklok = 2s.

10000

Bijlage C. Figuren

85

debiet (l/s)

10

800 600


400 0

0.5

1 tijdsmonster

1.5

0

2

0

0.5

x 10

200

1 tijdsmonster

1.5

2 4

x 10

260 240

150

220 100

50

0

0.5

1 tijdsmonster

1.5

2 4

x 10

10 fout in %

5

4

5


zonlicht (W/m²)

1000


80 0

0

0.5

1 tijdsmonster

1.5

2 4

x 10

60

0

0.5

1 tijdsmonster

1.5

2 4

x 10

Figuur C.4: Gegevens van 16 juli met Tklok = 1s.

In de volgende figuren zijn de resultaten van de simulator bekomen bij het testen van de data van 15, 16 en 28 juli nog eens in het groot afgedrukt (tijdsmonster telkens 1 s).

zonlicht (W/m²)


fout in %

0

5

10

0

120 0

130

140

150

160

200 0

400

600

800

2000

2000

2000

4000

4000

4000




10000 12000 14000

10000 12000 14000

10000 12000 14000

debiet (l/s) 0

50

0

100

150

200

250

0

5

10

Tgem (°C)

1000

2000

2000

4000

4000

10000 12000 14000


10000 12000 14000

reële temperatuur gesimuleerde temperatuur


Bijlage C. Figuren 86

Figuur C.5: Resultaten voor de temperatuur in punt C op 15 juli.

zonlicht (W/m²)


afvoersectie.

fout in %

0

5

10

0

120 0

130

140

150

160

200 0

400

600

800

2000

2000

2000

4000

4000

4000




10000 12000 14000

10000 12000 14000

10000 12000 14000

debiet (l/s) uitgangstemperatuur (°C)

1000

0

60

80 0

100

120

140

160

180

200

220

240

260

0

5

10

2000

2000

4000

4000

10000 12000 14000


10000 12000 14000




Figuur C.6: Resultaten voor de temperatuur in punt D op 15 juli met verbeterde lengte van de

zonlicht (W/m²)


fout in %

0

5

10

50

0

0

100

150

200

400 0

600

800

0.5

0.5

0.5

1 tijdsmonster

1 tijdsmonster

1 tijdsmonster

1.5

1.5

1.5

2

2 4

x 10

4

x 10

4

x 10

2

debiet (l/s) 0

50

0

100

150

200

250

0

5

10

Tgem (°C)

1000

0.5

0.5

1.5

1 tijdsmonster

1.5

2

2 4

x 10

4

x 10


1 tijdsmonster


Figuur C.7: Resultaten voor de temperatuur in punt C op 16 juli.

zonlicht (W/m²)


fout in %

afvoersectie.

0

5

10

50

0

0

100

150

200

400 0

600

800

0.5

0.5

0.5

1 tijdsmonster

1 tijdsmonster

1 tijdsmonster

1.5

1.5

1.5

x 10

4

2

4

x 10

2

4

x 10

2


1000

0

60

80 0

100

120

140

160

180

200

220

240

260

0

5

10

0.5

0.5

1.5

1 tijdsmonster

1.5

2

x 10

4

2

4

x 10


1 tijdsmonster



zonlicht (W/m²)


fout in %

0

5

10

0

50

0

0

100

150

200

200 0

400

600

800

0.5

0.5

0.5

1 1.5 tijdsmonster

1 1.5 tijdsmonster

1 1.5 tijdsmonster

2

2

2

Figuur C.9: Resultaten voor de temperatuur in punt C op 28 juli. x 10

4

2.5

x 10

4

2.5

x 10

4

2.5

debiet (l/s) 0

50

0

100

150

200

250

0

5

10

Tgem (°C)

1000

0.5

0.5

2

1 1.5 tijdsmonster

2

4

x 10

2.5

4

2.5 x 10


1 1.5 tijdsmonster


zonlicht (W/m²)


afvoersectie.

fout in %

0

5

10

0

50

0

0

100

150

200

200 0

400

600

800

0.5

0.5

0.5

1 1.5 tijdsmonster

1 1.5 tijdsmonster

1 1.5 tijdsmonster

2

2

2

x 10

4

2.5

x 10

4

2.5

x 10

4

2.5


1000

0

60

80 0

100

120

140

160

180

200

220

240

260

0

5

10

0.5

0.5

2

1 1.5 tijdsmonster

2

4

x 10

2.5

4

2.5 x 10


1 1.5 tijdsmonster



Bijlage C. Figuren

C.2

92

Figuren bij hoofdstuk 5

In de volgende figuren zijn de resultaten van de simulator bekomen bij het testen van de data van 15, 16 en 28 juli via de POD-methode nog eens in het groot afgedrukt.

zonlicht (W/m²)


fout in %

methode.

0

5

10

0

120 0

130

140

150

160

200 0

400

600

800

50

50

50

100

100

100




300

300

300

350

350

350

400

400

400


1000

0

100 0

150

200

250

0

5

10

50

50

100

100

300

350

400


300

350

400

reele temperatuur gesimuleerde temperatuur



Figuur C.11: Resultaten voor de temperatuur op het eind van de afvoersectie op 15 juli via POD-

zonlicht (W/m²)


fout in %

methode.

0

5

10

50

0

0

100

150

200

400 0

600

800

100

100

100

200

200

200




500

500

500

600

600

600


1000

0

40

60

80

0

100

120

140

160

180

200

220

240

0

5

10

100

100

200

200

500

600


500

600





zonlicht (W/m²)


fout in %

methode.

0

5

10

0

50

0

0

100

150

200

200 0

400

600

800

100

100

100

200

200

200




500

500

500

600

600

600

700

700

700


1000

0

0

50

0

100

150

200

250

300

0

5

10

100

100

200

200

500

600

700


500

600

700





Bijlage C. Figuren

96

In de volgende drie figuren worden de resultaten van het toestandsmodel met die van de POD-methode vergeleken. De basisfuncties worden telkens berekend met de gegevens van de overeenkomstige dag.

zonlicht (W/m²)



300

300

300

350

350

350

400

400

400

0

160

180

200

220

240

260

0

120 0

100



0

50

100

100

5

10

140

0

50

50

debiet (l/s)

1

2

3

120 0

130

140

150

160

200 0

400

600

800

fout in %

15 juli. uitgangstemperatuur (°C)

1000

50

50

100

100



350

300

350

toestandsmodel POD

300

400

400


Figuur C.14: Vergelijking van de resusultaten voor het toestandsmodel en de POD-methode voor

zonlicht (W/m²)


fout in %

16 juli. 300 400 tijdsmonster


500

500

500

600

600

600

0

0

100

120

140

160

180

200

220

240

0

60

200

200


0

100

100

200

80

0

0

100

5

10

0.5

1

1.5

2

50

100

150

200

400 0

600

800


1000

100

100

200

200

500

600


500

600

toestandsmodel POD




zonlicht (W/m²)


fout in %

28 juli.

0

0.5

1

1.5

0

50

0

0

100

150

200

200 0

400

600

800

100

100

100

200

200

200




500

500

500

600

600

600

700

700

700


1000

0

0

50

0

100

150

200

250

300

0

5

10

100

100

200

200



600

700

500

600

700

toestandsmodel POD

500



Bijlage D

Matlabcode

100

Bijlage D. Matlabcode

101

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % Auteur : Mathieu Vandenbulcke % % laatste versie : 06/06/2006 % % % % Dit programma simuleert de werking van het ACUREX-veld % % op het Plataforma Solar de Almeria. De simulator maakt % % gebruik van niet-uniforme tijdsmonsters. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear all; %%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%% % HET VOORAFGAANDE WERK % %%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%% % De gegevens inladen % %%%%%%%%%%%%%%%%%%%%%%% load(’data1507R’); load(’data1607R’); load(’data2807R’); data=data15R; %kies hier data15R, data16R of data28R voor de testdata %van 15,16 of 28 juli %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Initialiseren van vectoren en % % waarden toekennen aan alle parameters % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% lengte=size(data); lengte=lengte(1,1); T_klok=1; V_segm=0.0077;

%in seconden %in m

V_temp=0; klok=0; N=zeros(10,3); aantal=9;

%aantal operationele leidingen 1


r_inner=0.0318/2-0.00211; r_outer=0.0318/2; R_inner=4*2.54/2*0.01; rho_m=7850; c_m=419;

102

%straal binnenkant leiding in verhittingssectie %straal buitenkant leiding in verhittingssectie %straal binnenkant toever- en afvoersectie %dichtheid metaal van de leiding %specifieke warmtecapaciteit metaal van leiding

TIJD=zeros(30/T_klok*lengte,1); Ta=zeros(lengte,1); Irr=zeros(lengte,1); flow1=zeros(lengte,1); V=zeros(9,3);

%vector %vector %vector %vector %matrix

met verschillende t_k’s voor omgevingstemperatuur voor zonlicht voor debiet voor volume van elke sectie

%De temperatuursbereiken van de verschillende loops %worden bewaard voor het berekenen van de fout: maxwaarde=max(data); minwaarde=min(data); bereik2=maxwaarde(7)-minwaarde(7); %temperatuursbereik bereik3=maxwaarde(8)-minwaarde(8); %temperatuursbereik bereik4=maxwaarde(9)-minwaarde(9); %temperatuursbereik bereik5=maxwaarde(10)-minwaarde(10); %temperatuursbereik bereik6=maxwaarde(11)-minwaarde(11); %temperatuursbereik bereik7=maxwaarde(12)-minwaarde(12); %temperatuursbereik bereik8=maxwaarde(13)-minwaarde(13); %temperatuursbereik bereik9=maxwaarde(14)-minwaarde(14); %temperatuursbereik bereik10=maxwaarde(15)-minwaarde(15); %temperatuursbereik bereik11=maxwaarde(4)-minwaarde(4); %temperatuursbereik %vector voor het opslaan van de werkelijke %temperatuur van de olie aan het begin %van de toeversectie: Tin=zeros(30/T_klok*lengte,1); %vector voor het opslaan van de berekende %temperatuur op het eind van de toevoersectie: T1=zeros(30/T_klok*lengte,1); %vectoren voor het opslaan van de werkelijke %temperatuur van de olie in de verhittingssectie: T2echt=zeros(30/T_klok*lengte,1); T3echt=zeros(30/T_klok*lengte,1); T4echt=zeros(30/T_klok*lengte,1); T5echt=zeros(30/T_klok*lengte,1); T6echt=zeros(30/T_klok*lengte,1); T7echt=zeros(30/T_klok*lengte,1); T8echt=zeros(30/T_klok*lengte,1); 2

loop 2 loop 3 loop 4 loop 5 loop 6 loop 7 loop 8 loop 9 loop 10 in punt D


103

T9echt=zeros(30/T_klok*lengte,1); T10echt=zeros(30/T_klok*lengte,1); %vector voor het opslaan van de berekende %temperatuur van de olie in de verhittingssectie: T2=zeros(30/T_klok*lengte,9,50); %vector voor het opslaan van de berekende %temperatuur van het metaal in de verhittingssectie: Tm=zeros(30/T_klok*lengte,9,50); %vector voor het opslaan van de werkelijke %temperatuur op het eind van de avoersectie: Toutecht=zeros(30/T_klok*lengte,1); %vectoren om de fout tussen de werkelijke en berekende %temperatuur op het einde van de verhittings- en op het %einde van de afvoersectie te berekenen: fout22=zeros(30/T_klok*lengte,1); fout23=zeros(30/T_klok*lengte,1); fout24=zeros(30/T_klok*lengte,1); fout25=zeros(30/T_klok*lengte,1); fout26=zeros(30/T_klok*lengte,1); fout27=zeros(30/T_klok*lengte,1); fout28=zeros(30/T_klok*lengte,1); fout29=zeros(30/T_klok*lengte,1); fout210=zeros(30/T_klok*lengte,1); foutTout=zeros(30/T_klok*lengte,1); foutTCavg=zeros(30/T_klok*lengte,1);

%vectoren voor het opslaan van de werkelijke en berekende %gemiddelde olietemperatuur op het einde van de verhittingssectie: TCavg_echt=zeros(30/T_klok*lengte,1); TCavg=zeros(30/T_klok*lengte,1); %vectoren voor het opslaan van de berekende %temperaturen op het eind van de afvoersectie: Touter=zeros(30/T_klok*lengte,1); Tout=zeros(30/T_klok*lengte,9); %varia Irrplot=zeros(30/T_klok*lengte,1); F=zeros(30/T_klok*lengte,1);

3


104

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Berekenen van de volumes van elke % % leiding in elk van de drie secties: % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for k=1:3 for j=1:aantal if k==1 V(j,k)=(R_inner^2)*pi*(35.12+(j)*9.9)/aantal; end if k==2 V(j,k)=(r_inner^2)*pi*172; end if k==3 V(j,k)=(R_inner^2)*pi*(155.23+(j)*9.9)/aantal; % met 155.23 de experimenteel berekende lengte end end end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Berekenen hoeveel V_segm in elke % % leiding van elke sectie gaat: % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for j=1:aantal N(j,1)=round(V(j,1)/V_segm); N(j,2)=round(V(j,2)/V_segm); N(j,3)=round(V(j,3)/V_segm); end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Om de 30 s worden nieuwe waarden van de % % data ingelezen % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for i=1:lengte

4


105

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % De werkelijke gegevens worden gestockeerd % % in de overeenkomstige vectoren. % % Gezien de data slechts om de 30 s werd gemeten % % en de klok van het programma een periode T_klok % % heeft, wordt er een zero-order hold toegpast. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for s=1+30/T_klok*(i-1):30/T_klok*i Tin(s,1)=data(i,3); T2echt(s,1)=data(i,7); T3echt(s,1)=data(i,8); T4echt(s,1)=data(i,9); T5echt(s,1)=data(i,10); T6echt(s,1)=data(i,11); T7echt(s,1)=data(i,12); T8echt(s,1)=data(i,13); T9echt(s,1)=data(i,14); T10echt(s,1)=data(i,15); % Het gemiddelde van de werkelijk temperaturen op het % einde van de verhittingssectie wordt berekend: TCavg_echt(s,1)=(T2echt(s,1)+T3echt(s,1)+T4echt(s,1)+T5echt(s,1) +T6echt(s,1)+T7echt(s,1)+T8echt(s,1)+T9echt(s,1) +T10echt(s,1))/9; %varia Toutecht(s,1)=data(i,4); Irrplot(s,1)=data(i,6); F(s,1)=data(i,2);

end %varia Ta(i,1)=data(i,5); flow1(i,1)=data(i,2)/(aantal*1000); %omzetten naar m/s en delen door %aantal leidingen Irr(i,1)=data(i,6); %%%%%%%%%%%%%%%%%%%%%%%%% % De tijdsiteratie die % % verspringt met T_klok % %%%%%%%%%%%%%%%%%%%%%%%%%

5


106

for s=1+30/T_klok*(i-1):30/T_klok*i %Berekening van de t_k’s: V_temp=V_temp+flow1(i,1)*T_klok; if V_temp>V_segm klok=klok+1; V_temp=0; TIJD(klok,1)=s; end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % De berekening van de temperatuur op % % het eind van de toevoersectie % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for n=1:aantal if klok<=N(n,1) T1(s,n)=Tin(1,1); %beginvoorwaarden invullen else T1(s,n)=Tin(TIJD(klok-N(n,1)),1); end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % De berekening van de temperatuur % % in de verhittingssectie % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for n=1:aantal epsilon=0.7; A_mirr=5.42*6*8/N(n,2); G_s=(2*r_outer*pi*172/N(n,2)); A_mf=2*r_inner*172*pi/N(n,2); V_m=(r_outer^2-r_inner^2)*pi*172/N(n,2); V_f=r_inner^2*pi*172/N(n,2); H_l=29; for k=1:N(n,2) if s==1 H_t=1/pi*(2.17*10^6-5.01*10^4*Tin(1,1) +4.53*10^2*Tin(1,1)^2-1.64*Tin(1,1)^3 +2.1*10^(-3)*Tin(1,1)^4)*(flow1(i,1)^0.8);%beginwaarde else H_t=1/pi*(2.17*10^6-5.01*10^4*T2(s-1,n,k) +4.53*10^2*T2(s-1,n,k)^2-1.64*T2(s-1,n,k)^3 +2.1*10^(-3)*T2(s-1,n,k)^4)*(flow1(i,1)^0.8); end 6


107

if (klok==0)||(klok==1) Tm(s,n,k)=Tin(1,1); %beginwaarden invullen T2(s,n,k)=Tin(1,1); %beginwaarden invullen else delta_t=(TIJD(klok,1)-TIJD(klok-1,1))*T_klok; %randvoorwaarden invullen if k==1 rho_f=903-0.672*(T1(TIJD(klok-1),n)); c_f=1820+3.478*(T1(TIJD(klok-1),n)); Tm(s,n,k)=Tm(TIJD(klok-1),n,k) +delta_t/(rho_m*c_m*V_m)*(epsilon*A_mirr*Irr(i,1) -H_l*G_s*(Tm(TIJD(klok-1),n,k)-Ta(i,1)) -H_t*A_mf*(Tm(TIJD(klok-1),n,k) -T1(TIJD(klok-1),n))); T2(s,n,k)=T1(TIJD(klok-1),n)+ delta_t/(rho_f*c_f*V_f)*H_t*A_mf *(Tm(TIJD(klok-1),n,k) -T1(TIJD(klok-1),n)); %berekening van de temperaturen in het metaal %en de olie in de verhittingssectie: else rho_f=903-0.672*(T2(TIJD(klok-1),n,k-1)); c_f=1820+3.478*(T2(TIJD(klok-1),n,k-1)); Tm(s,n,k)=Tm(TIJD(klok-1),n,k) +delta_t/(rho_m*c_m*V_m)*(A_mirr*epsilon*Irr(i,1) -H_l*G_s*(Tm(TIJD(klok-1),n,k)-Ta(i,1)) -A_mf*H_t*(Tm(TIJD(klok-1),n,k) -T2(TIJD(klok-1),n,k))); T2(s,n,k)=T2(TIJD(klok-1),n,k-1) +delta_t/(rho_f*c_f*V_f)*A_mf*H_t *(Tm(TIJD(klok-1),n,k-1) -T2(TIJD(klok-1),n,k-1)); end end

end end 7


108

%Het gemiddelde van de berekende waarden op %het einde van de verhittingssecties berekenen: for n=1:aantal TCavg(s,1)=TCavg(s,1)+T2(s,n,N(n,2)); end TCavg(s,1)=TCavg(s,1)/aantal;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % De berekening van de temperatuur % % op het eind van de afvoersectie % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for n=1:aantal if klok<=N(n,3)+1 Tout(s,n)=Tin(1,1); %beginwaarden invullen else Tout(s,n)=T2(TIJD(klok-N(n,3)),n,N(n,2)); end end %Het gemiddelde nemen van %alle temperaturen op einde %van de afvoersectie: for n=1:9 Touter(s,1)=Touter(s,1)+Tout(s,n); end Touter(s,1)=Touter(s,1)/aantal; %foutenvectoren voor elke loop afzonderlijk: %fout22(s,1)=100*abs(T2echt(s,1)-T2(s,1,N(1,2)))/bereik2; %fout23(s,1)=100*abs(T3echt(s,1)-T2(s,2,N(2,2)))/bereik3; %fout24(s,1)=100*abs(T4echt(s,1)-T2(s,3,N(3,2)))/bereik4; %fout25(s,1)=100*abs(T5echt(s,1)-T2(s,4,N(4,2)))/bereik5; %fout26(s,1)=100*abs(T6echt(s,1)-T2(s,5,N(5,2)))/bereik6; %fout27(s,1)=100*abs(T7echt(s,1)-T2(s,6,N(6,2)))/bereik7; %fout28(s,1)=100*abs(T8echt(s,1)-T2(s,7,N(7,2)))/bereik8; %fout29(s,1)=100*abs(T9echt(s,1)-T2(s,8,N(8,2)))/bereik9; %fout210(s,1)=100*abs(T10echt(s,1)-T2(s,9,N(9,2)))/bereik10; %foutvector voor de gemiddelde temperatuur op het einde %van de verhittingssectie: foutTCavg(s,1)=100*abs(TCavg_echt(s,1)-TCavg(s,1)); %foutenvector voor de temperatuur op het einde 8


109

%van de afvoersectie: foutTout(s,1)=100*abs(Toutecht(s,1)-Touter(s,1))/bereik11;

end end %%%%%%%%%%%%%%%%%%%%%% % De figuren plotten % %%%%%%%%%%%%%%%%%%%%%% %figure; %hold on; %plot(F,’m’) %plot(tt1,’r’) %plot(Tin,’c’) %plot(Irrplot,’y’) %plot(T2echt,’r’) %plot(T3echt,’r’) %plot(T4echt,’r’) %plot(T5echt,’r’) %plot(T6echt,’r’) %plot(T7echt,’r’) %plot(T8echt,’r’) %plot(T9echt,’r’) %plot(T10echt,’r’) %plot(TCavg_echt,’r’) %plot(Toutecht,’k’) %plot(T2(s,1,N(1,2))) %plot(T2(s,2,N(2,2))) %plot(T2(s,3,N(3,2))) %plot(T2(s,4,N(4,2))) %plot(T2(s,5,N(5,2))) %plot(T2(s,6,N(6,2))) %plot(T2(s,7,N(7,2))) %plot(T2(s,8,N(8,2))) %plot(T2(s,9,N(9,2))) %plot(TCavg) %plot(Tm(s,1,N(1,2)),’k’) %plot(Tm(s,2,N(2,2)),’k’) 9


110

%plot(Tm(s,3,N(3,2)),’k’) %plot(Tm(s,4,N(4,2)),’k’) %plot(Tm(s,5,N(5,2)),’k’) %plot(Tm(s,6,N(6,2)),’k’) %plot(Tm(s,7,N(7,2)),’k’) %plot(Tm(s,8,N(8,2)),’k’) %plot(Tm(s,9,N(9,2)),’k’) %plot(Touter/aantal,’b’);ylabel(’C’) %plot(fout22) %plot(fout23) %plot(fout24) %plot(fout25) %plot(fout26) %plot(fout27) %plot(fout28) %plot(fout29) %plot(fout210) %hold off %figure; %subplot(3,2,1);plot(Irrplot,’b’);ylabel(’zonlicht (W/m)’); xlabel(’tijdsmonster’) %subplot(3,2,2);plot(F,’b’);ylabel(’debiet (l/s)’); xlabel(’tijdsmonster’);ylim([0,10]); %subplot(3,2,3);plot(Tin,’b’);ylabel(’ingangstemperatuur (C)’); xlabel(’tijdsmonster’) %subplot(3,2,[4 6]);hold on;plot(TCavg_echt,’r’); plot(TCavg,’b’);ylabel(’T_g_e_m (C)’);ylim([20,260]); xlabel(’tijdsmonster’);hold off; legend(’reele temperatuur’,’gesimuleerde temperatuur’,’Location’,’SouthEast’); %subplot(3,2,5);hold on;plot(foutTCavg/(max(TCavg_echt)-min(TCavg_echt)),’b’); plot(5*ones(30/T_klok*lengte,1),’k’);xlabel(’tijdsmonster’); ylabel(’fout in %’);ylim([0,10]);hold off; %figure; %subplot(2,1,1);hold on; plot(TCavg,’k’);plot(TCavg_echt,’r’);ylabel(’TC average(C)’) ;ylim([160,260]);xlabel(’timesamples’);hold off

10


111

%subplot(2,1,2);hold on;plot(foutTCavg/(max(TCavg_echt)-min(TCavg_echt)),’k’); plot(5*ones(30/T_klok*lengte,1),’r’);ylabel(’error in %’);ylim([0,10]); hold off; %figure; %subplot(3,2,1);plot(Irrplot,’b’); ylabel(’zonlicht (W/m)’);xlabel(’tijdsmonster’) %subplot(3,2,2);plot(F,’b’);ylabel(’debiet (l/s)’); xlabel(’tijdsmonster’);ylim([0,10]); %subplot(3,2,3);plot(Tin,’b’);ylabel(’ingangstemperatuur (C)’);xlabel(’tijdsmonster’) %subplot(3,2,[4 6]);hold on;plot(Toutecht,’r’);plot(Touter,’b’); ylabel(’uitgangstemperatuur (C)’);ylim([60,260]);xlabel(’tijdsmonster’); hold off;legend(’reele temperatuur’,’gesimuleerde temperatuur’,’Location’,’SouthEast’); %subplot(3,2,5);hold on;plot(foutTout,’b’);plot(5*ones(30/T_klok*lengte,1),’k’); xlabel(’tijdsmonster’);ylabel(’fout in %’);ylim([0,10]);hold off; %figure; %subplot(2,1,1);hold on;plot(Touter,’b’);plot(Toutecht,’r’);ylabel(’Tout(C)’); xlabel(’samples’);ylim([180,250]);hold off; %subplot(2,1,2);hold on;plot(foutTout,’b’);plot(5*ones(30/T_klok*lengte,1),’r’); ylabel(’error in %’);hold off; %figure; %subplot(4,1,[1 2]);hold on;plot(Toutecht,’r’);plot(Touter,’b’); ylabel(’outlet temperature (C)’);ylim([100,250]);xlabel(’time(s)’); hold off;legend(’real temperature’,’simulated temperature’,’Location’,’SouthEast’); %subplot(4,1,3);plot(F,’b’);ylabel(’flow (l/s)’);xlabel(’time (s)’);ylim([0,10]); %subplot(4,1,4);hold on;plot(foutTout,’b’);plot(5*ones(30/T_klok*lengte,1),’k’); xlabel(’time (s)’);ylabel(’error in %’);ylim([0,10]);hold off;

11


112

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % Auteur : Mathieu Vandenbulcke % % laatste versie : 06/06/2006 % % % % Dit programma simuleert de werking van het ACUREX-veld % % op het Plataforma Solar de Almeria. De simulator maakt % % gebruik van POD-technieken. % % % % Nota: Om dit programma te kunnen gebruiken dient de % % file ’berekenbasis2.m’ zich in dezelfde directory als % % deze file te bevinden. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear all; load (’data1507R’); load (’data1607R’); load (’data2807R’);

% data inladen van 15 juli % data inladen van 16 juli % data inladen van 28 juli

%%%%%%%%%%%%%%%%%%% % INPUT: % % vul de gewenste % % waarden in % %%%%%%%%%%%%%%%%%%% par1=1;

% ’fine tuning’ parameter 1

par2=0.51;

% ’fine tuning’ parameter 2

data=data28R;

% selecteer de gewenste data % (data15R, data16R of data28R)

dag=3;

% dag waarop je de basisfuncties % wil berekenen (1=15 juli,2=16 juli % 3=28 juli)

begin=1;

% % % %

einde=length(data);

Nref=1 ;

simulatie wordt gestart op dit tijdsmonster simulatie wordt beeindigd op dit tijdsmonster

% aantal keer dat je het tijdsegment van % de data (30s) wil verfijnen

1


113

dx=25;

% ruimtelijke stap in meter

n=12;

% aantal basisfuncties die % gebruikt moeten worden

%%%%%%%%%%%%% % PROGRAMMA % %%%%%%%%%%%%% nr_of_loop=6;

% er wordt gesteund op de eigenschappen % van loop 6

N=einde-begin;

% aantal tijdstappen

N2=Nref*N;

% aantal tijdstappen na verfijning

L1=35.12+(nr_of_loop-1)*9.9; % lengte van L1_apt=L1*(0.1016/2)^2/9*1/(0.02969/2)^2; % lengte van % aanpassing L2=172; % lengte van

de toevoersectie de toevoersectie, na van de diameter de verhittingssectie

L3=155.23+(nr_of_loop-1)*9.9; % lengte van de afvoersectie L3_apt=L3*(0.1016/2)^2/9*1/(0.02969/2)^2; % lengte van de afvoersectie, na % aanpassing van de diameter dt=(30*N)/N2;

% tijdstap

M1=round(L1_apt/dx); M2=round(L2/dx); M3=round(L3_apt/dx);

% aantal segmenten in toevoersectie % aantal segmenten in verhittingssectie % aantal segmenten in afvoersectie

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % de data aanpassen aan de verfijnde tijdstap % % via een zero-order hold % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for i=1:N for s=(i-1)*Nref+1:i*Nref flow(s,1)=data(begin+i-1,2)/(1000); % debiet irr(s,1)=data(begin+i-1,6); % zonlicht Tin(s,1)=data(begin+i-1,3); % ingangstemperatuur 2


114

Tout_real(s,1)=data(begin+i-1,4); % uitgangstemperatuur Tendloopi(s,1)=data(begin+i-1,nr_of_loop+5); % temperatuur op einde % van gebruikte loop end end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % tijdsonafhankelijk parameters % % definieren % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% rho_f=903-0.672*250; L=0.029;

% dichtheid van de olie % binnendiameter verhitt.sect.

V1=(L/2)^2*pi*9*L1_apt; V2=(L/2)^2*pi*9*L2; V3=(L/2)^2*pi*9*L3_apt;

% totaal volume toevoersectie % totaal volume verhitt. sectie % totaal volume afvoer sectie

c_f=1820+3.478*250; eff=0.69; G=9*5.42*8*6;

% warmtecapaciteit olie % reflectiecoeff. spiegels % volledige apertuur spiegels

c1=eff*G; c21=rho_f*c_f*L1_apt; c22=rho_f*c_f*L2; c23=rho_f*c_f*L3_apt; c31=rho_f*c_f*V1; c32=rho_f*c_f*V2; c33=rho_f*c_f*V3;

%%%%%%%%%%%%%%%%%%%%%%% % POD-basis berekenen % %%%%%%%%%%%%%%%%%%%%%%% phi=berekenbasis2(dag,dx,nr_of_loop,L1_apt,L2,L3_apt,M1,M2,M3,par1,par2); phi_n=phi(:,1:n);

% getrunceerde set basisfuncties

3


115

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % berekening van de tijdsaf% % hankelijk coefficienten a_i(t) % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% a(1:n,1)=phi_n’*ones(M1+M2+M3,1)*data(1,3); % als temperatuursdistributie % op het eerste tijdstip wordt % de ingangstemperatuur op het % eerste tijdstip gekozen for k=2:N2 a11=par1*(c21*dt/(c31*dx)*flow(k-1)); a12=par1*(c22*dt/(c32*dx)*flow(k-1)); a13=par1*(c23*dt/(c33*dx)*flow(k-1)); a2=par2*(-c1*dt/c32);

for i=1:M1+M2+M3 if i<=M1

% opvullen van de matrix A

A(i,i)=1+a11; if (i>1) A(i,i-1)=-a11; end elseif (i>M1)&&(i<=M2) A(i,i)=1+a12; A(i,i-1)=-a12; elseif (i>M2) A(i,i)=1+a13; A(i,i-1)=-a13; end end A_inv=inv(A); B=phi_n’*A_inv*phi_n; D=phi_n’*A_inv;

% inverse van A berekenen % matrix B berekenen % matrix D berekenen

a(1:n,k)=B*a(1:n,k-1)+D*([a11*Tin(k-1,1);zeros(M1+M2+M3-1,1)] -[zeros(M1,1);ones(M2,1)*a2*irr(k-1,1);zeros(M3,1)]); end

4


116

T_POD=phi_n*a;

% Temperatuur berekenen met % POD-methode

%%%%%%%%%%%%%%%%%% % fout berekenen % %%%%%%%%%%%%%%%%%% bovengrens=max(Tout_real(:,1)); ondergrens=min(Tout_real(:,1)); fout=abs(Tout_real(:,1)-T_POD(M1+M2+M3,:)’)/(bovengrens-ondergrens)*100; FOUT=sum(fout)/(N2); %%%%%%%%%%%%%%%%%%% % figuren plotten % %%%%%%%%%%%%%%%%%%% figure; subplot(3,2,1);plot(irr,’b’);ylabel(’zonlicht (W/m)’); xlabel(’tijdsmonster’); subplot(3,2,2);plot(1000*flow,’b’);ylabel(’debiet (l/s)’); xlabel(’tijdsmonster’);ylim([0,10]); subplot(3,2,3);plot(Tin,’b’);ylabel(’ingangstemperatuur (C)’); xlabel(’tijdsmonster’) subplot(3,2,[4 6]);hold on;plot(Tout_real,’r’); plot(T_POD(M1+M2+M3,:),’b’);ylabel(’uitgangstemperatuur (C)’); xlabel(’tijdsmonster’);hold off; legend(’reele temperatuur’,’gesimuleerde temperatuur’,’Location’,’SouthEast’); subplot(3,2,5);hold on;plot(fout,’b’);plot(5*ones(length(fout),1),’k’); xlabel(’tijdsmonster’);ylabel(’fout in %’);ylim([0,10]);hold off;

5


117

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % Auteur : Mathieu Vandenbulcke % % laatste versie : 06/06/2006 % % % % Dit programma simuleert de werking van het ACUREX-veld % % op het Plataforma Solar de Almeria. De simulator maakt % % gebruik van POD-technieken. % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [phi]=berekenbasis2(dag,dx,nr_of_loop,L1_apt,L2,L3_apt,M1,M2,M3,par1,par2) load (’data1507R’); load (’data1607R’); load (’data2807R’);

% data inladen van 15 juli % data inladen van 16 juli % data inladen van 28 juli

if dag==1 data=data15R; end

% POD-basis berekenen uit % gegevens 15 juli





%%%%%%%%%%%%%%%%%%% % INPUT: % % vul de gewenste % % waarden in % %%%%%%%%%%%%%%%%%%% begin=1;

% simulatie wordt gestart op dit % tijdsmonster

einde=length(data);

% simulatie wordt beeindigt op dit % tijdsmonster

Nref=1 ;

% aantal keer dat het het tijdsegment van % de data (30s) wil verfijnen 1


118

%%%%%%%%%%%%% % PROGRAMMA % %%%%%%%%%%%%% N=einde-begin;

% aantal tijdstappen

N2=Nref*N;

% aantal tijdstappen na verfijning

dt=(30*N)/N2;

% tijdstap

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % de data aanpassen aan de verfijnde tijdstap % % via een zero-order hold % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for i=1:N for s=(i-1)*Nref+1:i*Nref flow(s,1)=data(begin+i-1,2)/(1000); % debiet irr(s,1)=data(begin+i-1,6); % zonlicht Tin(s,1)=data(begin+i-1,3); % ingangstemperatuur Tout_real(s,1)=data(begin+i-1,4); % uitgangstemperatuur Tendloopi(s,1)=data(begin+i-1,nr_of_loop+5); % temperatuur op einde % van gebruikte loop end end

for i=1:M1+M2+M3 u(i,1)=data(begin,3);

% % % %

temperatuursdistributie op tijdstip 1 is onbekend en wordt gelijkgesteld aan ingangstemp. op tijdstip 1

end rho_f=903-0.672*250; L=0.029;

% dichtheid van de olie % binnendiameter verhitt.sect.

V1=(L/2)^2*pi*9*L1_apt; V2=(L/2)^2*pi*9*L2; V3=(L/2)^2*pi*9*L3_apt;

% totaal volume toevoersectie % totaal volume verhitt. sectie % totaal volume afvoer sectie

2


119

c_f=1820+3.478*250; eff=0.69; G=9*5.42*8*6;

% warmtecapaciteit olie % reflectiecoeff. spiegels % volledige apertuur spiegels

c1=eff*G; c21=rho_f*c_f*L1_apt; c22=rho_f*c_f*L2; c23=rho_f*c_f*L3_apt; c31=rho_f*c_f*V1; c32=rho_f*c_f*V2; c33=rho_f*c_f*V3;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % berekenen van de temperatuursdistributie % % in het systeem op de verschillende tijdstippen % % door middel van discreet toestandsmodel % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for k=2:N2 a11=par1*(c21*dt/(c31*dx)*flow(k-1)); a12=par1*(c22*dt/(c32*dx)*flow(k-1)); a13=par1*(c23*dt/(c33*dx)*flow(k-1)); a2=par2*(-c1*dt/c32); for i=1:M1+M2+M3 if i<=M1

% opvullen van de matrix A

A(i,i)=1+a11; if (i>1) A(i,i-1)=-a11; end elseif (i>M1)&&(i<=M2) A(i,i)=1+a12; A(i,i-1)=-a12; elseif (i>M2) A(i,i)=1+a13; A(i,i-1)=-a13; end end 3


120

b=[a11*Tin(k-1,1);zeros(M1+M2+M3-2,1);0]+u(1:M1+M2+M3,k-1) -[zeros(M1,1);ones(M2,1)*a2*irr(k-1,1);zeros(M3,1)]; u(1:M1+M2+M3,k)=trid(A,b); end

%%%%%%%%%%%%%%%%%%%%%%%%%%% % berekenen basisfuncties % %%%%%%%%%%%%%%%%%%%%%%%%%%% C=1/N2*u*u’; [V W]=eig(C);

% correlatie matrix % eigenwaardendecompositie

for i=1:M1+M2+M3 lambda(i,1)=real(W(i,i)); end

% % % %

phi=real(V);

door numerieke fouten durven er al eens complexe eigenwaarden te ontstaan, wat niet kan, aangezien C reeel symmetrisch is

% zelfde opmerking voor % eigenvectoren

4

Bibliografie [1] E.F. Camacho, M. Berenguel, and F.R. Rubio. Advanced control of solar power plants. 1997. [2] Tor A. Johansen and Camilla Storaa. An internal energy controller for distributed solar collector fields. Improving human potential programme, proceedings of the 2nd users workshop, pages 111–117, February 2002. [3] E.F. Camacho, M. Berenguel, and F.R. Rubio. Simulation software package for the ACUREX field. Springer. [4] Tor A. Johansen and Camilla Storaa. Energy-based control of a distributed solar collector field. [5] Michaela Sbarciog, Bart Wyns, Clara Ionescu, Robin De Keyser, and Luc Boullart. Neural network models for a solar power plant. [6] Clara Ionescu, Bart Wyns, Mihaela Sbarciog, Luc Boullart, and Robin De Keyser. Comparison between physical modelling and neural network modelling of a solar power plant. [7] C. Pereira and A. Dourade. Application of a neuro-fuzzy network with support vector learning to a solar power plant. Improving human potential programma, proceedings of the 2nd users workshop, pages 71–78, February 2002. [8] Y. C. Liang, H. P. Lee, S. P. Lim, W. Z. Lin, K. H. Lee, and C. G. Wu. Proper orthogonal decomposition and its applications-part i: Theory. Journal of Sound and Vibration, 252(3):527–544, 2002.

121

Bibliografie

122

[9] Y. C. Liang, W. Z. Lin, H. P. Lee, S. P. Lim, K. H. Lee, and H. Sun. Proper orthogonal decomposition and its applications-part ii: Model reduction for mems dynamical analysis. Journal of Sound and Vibration, 256(3):515–532, December 2002. [10] Anindya Chatterjee. An introduction to proper orthogonal decomposition, December 1999. [11] Patricia Astrid. Reduction of Process Simulation Models: a proper orthogonal decomposition approach. PhD thesis, Technische Universiteit Eindhoven, November 2004. [12] R. N. Silva, L. M. Rato, and J.M. Lemos. Observer based time warped control of distributed collector solar fields. Improving human potential programma, proceedings of the 2nd users workshop, pages 103–110, February 2002. [13] Oscar Mauricio Agudelo Man˜ nozca.

Control of the temperature profile of an one-

dimensional bar using a pod approach. presentation, January 2006. [14] James R. Welty, Charles E. Wicks, Robert E. Wilson, and Gregory Rorrer. Fundamentals of Momentum, Heat and Mass Transfer. [15] Yang, Cao, Chung, and Morris. Applied Numerical Methods Using Matlab. John Wiley & Sons Inc, 2005.

Modelleren en Simuleren van. een Zonnecentrale

Recommend Documents