Diese Arbeit wurde vorgelegt amLehrstuhl fur Mathematik (MathCCES)

Simulation von Direktverdampfung in SolarthermischenKraftwerken mittels Finite Volumen Verfahren

Simulation of Direct-Steam Generation in SolarThermal Power Plants using Finite Volume Methods

ProjektarbeitComputational Engineering Science

September 2017

Vorgelegt von Benjamin Terschanski, 343974Presented by benjamin.terschanski@rwth-aachen.de

Aron Zingler, 344126Aron.Zingler@rwth-aachen.de

Fynn Kepp, 344700fynn.kepp@rwth-aachen.de

Erstprufer Prof. Dr. Martin FrankFirst examiner Lehrstuhl fur Mathematik (MathCCES)

RWTH Aachen University

Koreferent Dr. rer. nat. Pascal RichterCo-supervisor Lehrstuhl fur Mathematik (MathCCES)

RWTH Aachen University

Eigenstandigkeitserklarung

Hiermit versichern wir, dass wir diese Projektarbeit selbstandig verfasst und keine an-deren als die angegebenen Quellen und Hilfsmittel benutzt haben. Die Stellen unsererArbeit, die dem Wortlaut oder dem Sinn nach anderen Werken entnommen sind, habenwir in jedem Fall unter Angabe der Quelle als Entlehnung kenntlich gemacht. Dasselbegilt sinngemaß fur Tabellen und Abbildungen. Diese Arbeit haben wir in dieser odereiner ahnlichen Form noch nicht im Rahmen einer anderen Prufung vorgelegt.

Aachen, im September 2017

II

Contents

List of Figures V

List of Tables VII

1. Introduction 1

2. Optical model 22.1. Problem description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22.2. Deterministic Raytracer Modelling . . . . . . . . . . . . . . . . . . . . 22.3. Plant components and geometry . . . . . . . . . . . . . . . . . . . . . . 22.4. Mirror alignment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32.5. Mirror curvature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32.6. Mirror surface discretization . . . . . . . . . . . . . . . . . . . . . . . . 32.7. Optical effects and losses . . . . . . . . . . . . . . . . . . . . . . . . . . 5

2.7.1. Shadowing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52.7.2. Blocking . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

2.8. Secondary reflector . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62.9. Deterministic error modeling . . . . . . . . . . . . . . . . . . . . . . . . 72.10. Probability integration along the circumference . . . . . . . . . . . . . 82.11. 3D Extension . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

2.11.1. Sun position . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102.11.2. Offset in z direction . . . . . . . . . . . . . . . . . . . . . . . . . 102.11.3. 3D-Errorcone . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

2.12. Computation of irradiative power . . . . . . . . . . . . . . . . . . . . . 122.13. Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132.14. Case study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

2.14.1. Test setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142.14.2. Consistency verification with SolTrace . . . . . . . . . . . . . . 15

2.15. Accuracy against Ray count . . . . . . . . . . . . . . . . . . . . . . . . 162.15.1. Circumferential power distribution . . . . . . . . . . . . . . . . 172.15.2. Changes in the mirror position . . . . . . . . . . . . . . . . . . . 172.15.3. Changes in the azimutal component . . . . . . . . . . . . . . . . 182.15.4. Solution on a 3D discretization . . . . . . . . . . . . . . . . . . 19

3. Thermodynamical Model 203.1. Prior Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 203.2. Finite Volume Method . . . . . . . . . . . . . . . . . . . . . . . . . . . 213.3. Solvers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

3.3.1. Lax-Friedrichs . . . . . . . . . . . . . . . . . . . . . . . . . . . . 223.3.2. Rusanov . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 223.3.3. HLL . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 233.3.4. HLLX . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

III

3.3.5. HLLC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 243.3.6. Godunov . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

3.4. Analysis of computation time and accuracy . . . . . . . . . . . . . . . . 273.4.1. Sod’s shock tube test . . . . . . . . . . . . . . . . . . . . . . . . 283.4.2. Accuracy versus number of cells . . . . . . . . . . . . . . . . . . 293.4.3. Accuracy versus CPU runtime . . . . . . . . . . . . . . . . . . . 31

3.5. Matlab implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . 323.5.1. Lax-Friedrichs . . . . . . . . . . . . . . . . . . . . . . . . . . . . 323.5.2. Rusanov . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 323.5.3. HLL . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 323.5.4. HLLX . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 333.5.5. HLLC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 333.5.6. Godunov-Suliciu . . . . . . . . . . . . . . . . . . . . . . . . . . 33

3.6. Parallelization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 333.7. Source Terms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34

4. Closure of the two-phase flow model 364.1. Mixture Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 364.2. Prior Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 364.3. Noble-Able Stiffened-Gas Equation of State . . . . . . . . . . . . . . . 38

4.3.1. Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 384.3.2. Two-phase region . . . . . . . . . . . . . . . . . . . . . . . . . . 394.3.3. Thompson Coefficient in NASG . . . . . . . . . . . . . . . . . . 40

4.4. IAWPS calculation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 414.4.1. IAWPS ζ calculation . . . . . . . . . . . . . . . . . . . . . . . . 414.4.2. Two-phase region . . . . . . . . . . . . . . . . . . . . . . . . . . 414.4.3. Calculation of (T, ρ′, ρ′′ and α) . . . . . . . . . . . . . . . . . . . 424.4.4. Detection of two-phase region . . . . . . . . . . . . . . . . . . . 42

4.5. Methods under consideration . . . . . . . . . . . . . . . . . . . . . . . . 434.5.1. Look-up table . . . . . . . . . . . . . . . . . . . . . . . . . . . . 444.5.2. Locally fit parameters . . . . . . . . . . . . . . . . . . . . . . . 454.5.3. Comparison of methods . . . . . . . . . . . . . . . . . . . . . . 46

4.6. Threadment of the liquid and vapour region . . . . . . . . . . . . . . . 48

5. Test Case 495.1. Boundaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 515.2. Discussion Test Results . . . . . . . . . . . . . . . . . . . . . . . . . . . 51

6. Conclusion 536.1. Outlook . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

A. Appendix 54

References 58

IV

List of Figures

1. PE1 power plant, Calasparra, Spain . . . . . . . . . . . . . . . . . . . . 12. Definition of the the coordinate system and the geometry of the 2D

projection. In this work we will describe the orientation / deviation ofthe plant from the local north axis with δ . . . . . . . . . . . . . . . . . 3

3. Visualization of the normal projection | ~PQ| for an incident vector ~i . . 44. Sample distribution of the thermal power associated with a single mirror

in the x-y plane among N=5 discrete rays . . . . . . . . . . . . . . . . 45. Gaussian surface discretization for N rays . . . . . . . . . . . . . . . . 56. Shadowing of a surface S (marked in red) . . . . . . . . . . . . . . . . . 67. Blocking of rays reflected on surface B (marked in red) . . . . . . . . . 68. Visualisation of the standard deviation in 2D and 3D . . . . . . . . . . 89. Parameters in the computation of the hit probability for an error cone

onto the segment i. The probability distribution is integrated over theangular deviation with respect to the ideal ray from angle ρ to angle ω 9

10. Introduction of concepts for the 3D extension. Figure 10a defines theazimuth angle τ and the altitude (also: elevation) angle ε . . . . . . . . 10

11. Steps in computing the energy distribution along the z-axis. The cu-mulative distribution function converges against the Heaviside step-function for lim σ → 0 (green) . . . . . . . . . . . . . . . . . . . . . . 12

12. Sample parameters for the computation of the heat flux on the axialsegment i with bounds a and b under the hit probability P(z). Thesegment position relative to the hit point of an ideal ray is characterizedusing the bound-angles γ1 and γ2 . . . . . . . . . . . . . . . . . . . . . 13

13. Simplified Class diagram (only basic functionalities are shown to presentthe modeling framework) . . . . . . . . . . . . . . . . . . . . . . . . . . 14

14. Reference geometry following the geometry of the PE1 power plant . . 1515. Visualization of the 2D ray tracing algorithm to verify the geometry of

reflection for different incident angles . . . . . . . . . . . . . . . . . . . 1517. Absolute deviation from a simulated ideal solution using 10000 rays per

mirror for different raycounts. The mirror configuration corresponds tothe layout shown in Section 2.14.1 . . . . . . . . . . . . . . . . . . . . . 16

18. Qualitative plot of the circumferential power distribution from 0◦ to360◦ in counter-clockwise direction with the zero angle being defined in18a . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

19. The plots show the variation of the offset in z direction as a single mirroris moved away from the position directly under the tube (X Pos=0) toan offset of 9m under a fixed sun-position. . . . . . . . . . . . . . . . . 18

20. Variation of the transition zone over the course of a sunrise on june 21st2017 at PE1, Calasparra, Spain. The simulation was run using threesample mirrors at x=[0, 3.9375, 7.875] [m] of equal width w = 0.75m . 19

V

21. Irradiation condition Qopt at the absorber tube simulated using the ge-ometry 14. Discretization is done along the circumference (angular) andthe z-Axis (axial). The shift along the axial direction is evident for moredominant azimutal components (Fig. 21b) . . . . . . . . . . . . . . . . 20

22. Propagation of information at a cell face for the HLL solver . . . . . . 2323. Propagation of information at a cell face for the HLLC solver . . . . . . 2524. Propagation of information at a cell face in the exact system . . . . . . 2625. Plot of the density vs. the length for 50 cells . . . . . . . . . . . . . . . 2926. Accuracy vs. number of cells for the different solvers . . . . . . . . . . 3027. Accuracy vs. number of cells for the different solvers . . . . . . . . . . 3128. Pipe wall acts as a buffer for heat . . . . . . . . . . . . . . . . . . . . 3529. Flowchart of the phase detection algorithm . . . . . . . . . . . . . . . 4330. Different scalings of the ρ, e plane . . . . . . . . . . . . . . . . . . . . . 4431. Methods that iterate on the IAWPS equations more accurate over the

whole interval. Methods based on stiffened gas equations have high errorin regions of higher temperature . . . . . . . . . . . . . . . . . . . . . . 46

32. Accuracy comparison in the α, T -plane . . . . . . . . . . . . . . . . . . 4734. Execution speed. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4733. Relative error of look-up method . . . . . . . . . . . . . . . . . . . . . 4835. Pressure varies heavily for nearly incompressible region . . . . . . . . . 4936. The simulation results show a clear similarity with results from ColSimQT 5237. Relative error of NASG method with fitted parameters . . . . . . . . . 5438. Relative error of NASG method with global parameters . . . . . . . . . 5439. Relative error of look-up method . . . . . . . . . . . . . . . . . . . . . 5540. Relative error of IAWPS95 Method . . . . . . . . . . . . . . . . . . . . 5541. Absolute error of NASG method with fitted parameters . . . . . . . . . 5642. Absolute error of NASG method with global parameters . . . . . . . . 5643. Absolute error of look-up method with iterations . . . . . . . . . . . . . 5744. Absolute error of look-up method . . . . . . . . . . . . . . . . . . . . . 57

VI

List of Tables

1. Deviation of the numerical solutions from the analytical solution forN=50 cells . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

2. Parameters and statistics of surface fit in the liquid region . . . . . . . 493. Parameters of the simulated pipe . . . . . . . . . . . . . . . . . . . . . 504. Values chosen for boundary values . . . . . . . . . . . . . . . . . . . . . 50

VII

1. Introduction

The fact that fossil resources are limited and most traditional forms of power generationare running on borrowed time motivates the rapid growth of alternative energy sourcessuch as wind, geothermal and solar energy over the last decades.One of multiple approaches to the construction of solar thermal power plants is theFresnel power plant. It consists of rows of mirrors concentrating solar radiation ontoan absorber tube to heat a working fluid (see Fig. 1).

Figure 1: PE1 power plant, Calasparra, Spain

To maximize the plant’s efficiency an accurate prediction of the flow conditions insideof the absorber tube is essential. This is a rather complex problem and the optimiza-tion of such a comprehensive model is still subject to current research.

This report documents a project work simulating a two phase flow of water in anabsorber tube of a Fresnel solar thermal power plant. The written simulation softwareconsists of three components, which are going to be put together in the end to formthe complete model of one absorber tube:

1. An optical model to compute the optical heat flux onto the absorber tube

2. A finite volume flow solver for the internal flow conditions

3. A thermodynamical model for the thermodynamic closure of the finite volumemodel

An efficient deterministic 3D ray tracing algorithm will be described and the irradiationconditions on a 3D discretization of the tube surface will be determined in Section2. Several finite volume algorithms (Lax-Friedrichs, Rusanov, HLL, HLLX, HLLC,Godunov) will be implemented and tested in terms of accuracy and efficiency in Section3. Different approaches for the computation of the thermodynamic state in the two-phase regime are implemented and compared in Section 4.

1

2. Optical model

2.1. Problem description

The flow conditions inside the tube and therefore all structural considerations as wellas the overall plant efficiency are governed by the distribution of irradiation on theabsorber tube. The source terms of the Euler equations to be solved for the internalflow include the heat source term, thus an adequate description of the heat transferonto the absorber tube is essential.For a given plant the geometry of the plant as well as input parameters such as the sunposition as a function of daytime are given. To determine the energy fraction reflectedonto the tube surface a ray tracing algorithm needs to be designed to solve for the heatflux on a discretization of the tube surface.

Related work Several codes have been developed to efficiently simulate the irradiativeheat flux for different use cases. Examples include codes like the SPRAY and STRAL[2]software for heliostat receivers developed by the DLR and the general-purpose opticalMonte-Carlo ray tracing software SolTrace[31]. A comparison of existing software hasbeen done by Osorio et al.[17]. Since the SolTrace software allows for user definedplant geometries and sun shapes it will be used as a validation tool in the developmentof our custom code.

2.2. Deterministic Raytracer Modelling

For the sake of consistent optimization it is desirable to rely on a deterministic approx-imation of the heat flux onto the tube surface which will yield averaged results directlyas opposed to running a Monte Carlo simulation multiple times.A deterministic optical model will be described in this section using a 2D simplifica-tion that will be extended for arbitrary 3D sun positions. We will first explain plantcomposition and physical effects taken into account before presenting a fast approachto simulate 3D offsets.

2.3. Plant components and geometry

The Fresnel power plant consist of n parallel rows of Fresnel mirrors whose inclinationcan be aligned independently along their longitudinal axis. All rows reflect sunlightonto the absorber tube parallel to the mirrors to heat a working fluid. A secondaryreflector covers the tube to re-reflect rays missing the tube. In the following model thecoordinate system shown in Figure 2a will be used where the z-axis is aligned with thelongitudinal component of the tube. The x-axis is defined as the South-axis for theazimutal component in the description of the sun position. The geometric descriptionfollows the approach described by Pino et al.[20].

2

(a) Global coordinate system (b) Mirror alignment geometry

Figure 2: Definition of the the coordinate system and the geometry of the 2D projec-tion. In this work we will describe the orientation / deviation of the plantfrom the local north axis with δ

2.4. Mirror alignment

For a given point in time the sun’s position is parametrized in the x-y plane using thesun angle α. As shown in Figure 2b the mirror alignment can be characterized by theangle between the x-axis and the mirror normal ~n. Our goal is to set this angle suchthat the reflected ray is directed at the tube center. Using basic physics and geometrywe find from the law of specular reflection that

γ = α +arctan2(PosTube.Y − PosMirror.Y, PosTube.X − PosMirror.X)

2(1)

2.5. Mirror curvature

To maximize efficiency mirrors are not flat surfaces but slightly curved. We accountfor curvature effects by modeling the shape of the mirror surface with a highly dilatedparabola. If we define a focal length f we can find the shape function accordingly using

y(x) =1

4 · f· x (2)

For real plants the focal length is typically large and curvature small.

2.6. Mirror surface discretization

It is customary to express the irradiative heat flux in [Wm

] emitted by the sun in termsof the Direct Normal Irradiance (DNI). For a single mirror surface it is thenpossible to compute the total irradiative power associated with the mirror surfaceQDNI by computing the product of the DNI with a normal projection | ~PQ| of themirror surface with respect to the incident ray as shown in Figure 3.

3

Figure 3: Visualization of the normal projection | ~PQ| for an incident vector ~i

For a single mirror we thus obtain

QDNI,mirror = DNI · widthmirror ·cos(γ − α)

2(3)

where γ is the mirror normal angle and α the projected 2D sun incident angle asintroduced in Figure 2b.A more precise prediction of the power distribution is achieved by distributing this netenergy contribution of the mirror among discrete rays which are traced independently.Figure 4 shows how the total heat flux associated with a mirror projection in 2D issplit among discrete rays.

(a) Total reflected heat flux (b) Discrete heat flux

Figure 4: Sample distribution of the thermal power associated with a single mirror inthe x-y plane among N=5 discrete rays

The discretization of evaluation points for the computation of reflections on the mirrorsurface is done using a Gaussian integration where the ray hit points on the mirror aregiven by the abscissae of the Legendre polynomial Pn(x).

4

Following the Gauss-Legendre Quadrature quadrature rule we associate a weight ωiwith all rays such that

ωi = − 2

(1− x2i )[P

′n(xi)]2

(4)

Qmirror,ray,i = QDNI · ωi (5)

Qmirror,xy ≈N∑i=1

Qmirror,ray,i (6)

Figure 5 shows the discretization scheme for different surface resolutions. The choiceof Legendre abscissae for the position of discretization points has the advantage thatsurface edges have a higher discretization resolution (see Figure 5b) which facilitatesthe accurate simulation of optical losses.

(a) N=11 (b) N=101

Figure 5: Gaussian surface discretization for N rays

2.7. Optical effects and losses

For a realistic modeling of the irradiated heat flux several optical losses have to beconsidered. Since the mirrors are aligned parallel to the tube it is sufficient to considershadowing and blocking in the x-y plane:

2.7.1. Shadowing

Shadowing occurs when incident rays are obstructed from reaching the surface of amirror by another geometric entity. Figure 6 shows how optical losses occur at theedge of a mirror surface.

5

Figure 6: Shadowing of a surface S (marked in red)

2.7.2. Blocking

Blocking occurs when mirrors are located in the path of reflected rays (see Figure 7).

Figure 7: Blocking of rays reflected on surface B (marked in red)

2.8. Secondary reflector

The shape of the secondary reflector shape follows, for positive abscissas, a parabolicequation y = −a · x2 + b · x + c with a > 0. It is extended to negative abscissasby symmetry along the vertical axis. The mirror is set in a box which protects itagainst adverse weather conditions. It is also a way here to improve the plant efficiencyby optimizing isolation. The user has to provide several parameters to define thesecondary’s geometry:

1. Apex point coordinates

2. Isolation width

3. Secondary width

4. Distance above the absorber-tube

6

2.9. Deterministic error modeling

A key feature of the deterministic approach is the way that uncertainties in the modelare treated. Several parameters affect the ray tracing directly but can not be deter-mined analytically:

Mirror surface roughnessBecause of micro-scale roughness of the mirror surface incident rays will not be reflectedin a perfectly specular manner. The deviation from an ideal reflection is expressed us-ing the deviation angle σsurface

Tracking errorsErrors are committed by rotating Fresnel collectors because of mechanical limitationsof the tracking motors. This uncertainty is expressed in the parameter σtracking.

Sun shapeUp to now, ideal incident rays have been esteemed coming from the center of the sunand reaching the earth in parallel. Real incident rays are emitted at a random locationon the sun surface and reach the Earth with non-uniform directions. Thus we needto quantify the probability that one incident ray comes from a certain location on thesun. Each point on a Fresnel collector can be lit by a ray coming from the center of thesun (ideal ray), its extremity, or any other area as well. The deviation angle betweenreal and ideal rays is written as σsun which is assumed to follow a Gaussian distribution.

In our model all uncertainties are taken into account and the probability for a ray todeviate from the ideal specular reflection is assumed to behave according to a Gaussiandistribution (see Figure 8a) whose standard deviation is calculated from the σerrorConementioned above as

σerrorCone =√σ2

tracking + σ2surface + σ2

sun (7)

For 3D calculations we can simplify the expression since σtracking → 0 if we considerthe error in the axial direction of the mirror row. This is because mirrors only rotatearound one axis / in the x-y plane.A physical representation of the probability distribution can be found in the conceptof the error cone which can also be extended to the 3D case: The reflection of realrays will deviate from that of the ideal ray in all directions. σerrorCone can be regardedas the angle of a cone around the the ideal reflection in which 66.6% of the rays lieaccording to the Gaussian distribution we assumed. This is illustrated in Figure 8b.

7

(a) The directional probability of a real rayscenters around the ideal reflection

(b) The error cone in the 3D plane

Figure 8: Visualisation of the standard deviation in 2D and 3D

2.10. Probability integration along the circumference

Based on the stochastic model introduced in Section 2.9 we can discretize the tubearound its circumference. A reflection scenario is depicted in Fig 9 where one cir-cumferential section i has boundaries with deviation angles ρ and ω. If the ray hasa total energy Qtot,i we find the fraction sensed at circumferential segment i from thecumulative distribution function of the normal distribution:

P(Segmenti) =

∫ ω

ρ

e− x2

2σerrorCone

σerrorCone√

2πdx (8)

which can be approximated using the fundamental theorem of calculus as

P(Segmenti) ≈1

2(erf(

b√2σerrorCone

)− erf(a√

2σerrorCone)) (9)

with the error function

erf(x) =2√π

∫ ∞0

e−t2

dt (10)

8

Figure 9: Parameters in the computation of the hit probability for an error cone ontothe segment i. The probability distribution is integrated over the angulardeviation with respect to the ideal ray from angle ρ to angle ω

2.11. 3D Extension

In previous considerations we considered the sun to be aligned in the x-y plane of theplant and thus we could safely neglect the offset of the reflected ray into the z direction.In real scenarios the azimutal component of the sun’s position is not negligible and weexpect the reflection vector to have a z-component. This results in some fraction ofthe tube not sensing any irradiation as highlighted in green in Figure 10b.As suggested by Pino et al.[20] one may reuse the circumferential distribution and thetotal 2D irradiative power Qtotal,xy =

∑Mirrors Qmirror,xy obtained from the 2D-model to

solve for the 3D distribution of irradiation onto the absorber tube.

9

(a) Definition of azimuth and elevation[28] (b) Tube-offset (green) for nonzero azimuth

Figure 10: Introduction of concepts for the 3D extension. Figure 10a defines the az-imuth angle τ and the altitude (also: elevation) angle ε

2.11.1. Sun position

It is customary to describe the hemispherical position of celestial objects in terms ofazimuth and elevation angle. By using the coordinates as defined in 1.2.1 we can thenfind the solar incident angle to be

~i =

cos(ε) · cos(−τ + δ)

sin(ε)

cos(ε) · sin(−τ + δ)

(11)

where τ and ε describe the azimuth and elevation of the sun’s 3D position as definedin Figure 10a and δ is the orientation of the plant (see Figure 2a).

2.11.2. Offset in z direction

Since the mirrors are aligned in parallel with the absorber tube it is sufficient tocalculate the projection of the sun’s 3D position onto the x-y plane of the power plantand proceed the 2D calculations with the angle obtained. From the definition of theincident angle we find the projected sun angle α to be

α =

{arctan

~iy~ix, for |τ | ≤ 90

180 + arctan~iy~ix, else

(12)

From that the correct mirror alignment and all optical losses can be computed as ex-plained earlier.

Since the geometry of the plant (position of the mirrors, height of the absorber tube) isknown we can find an approximate value for the z-component of the 3D unit reflection

10

vector ~r using a Householder reflection

~r = 2 · ~n|~n|· 〈 ~n|~n|,~i

|~i|〉 −

~i

|~i|(13)

where ~n is the normal vector of the reflection plane / mirror surface and ~i is the 3Dincident vector.If we assume the effects of the tube radius r and the distance of the reflection pointon the mirror surface from the rotational axis of the mirror d to be negligible withrespect to the total length of a reflected ray |~r| we may approximate the z-offset for areflection on mirror i as

zoffset = ~rz ·ytubeCenter − ymirror,i

~ry(14)

By doing so we significantly reduce the number of computational steps for large num-bers of rays by avoiding to track offsets individually on a per-ray basis. This comes ofcourse at the expense of a slight decrease in accuracy.

2.11.3. 3D-Errorcone

In order to achieve deterministic behaviour in the presence of uncertainties we extendthe concept of the error cone introduced for the 2D model into the third dimension.For a certain z-offset we directly obtain the approximate z-coordinate of an ideallyreflected ray along the tube. Since the mirrors have same length Ltube as the absorbertube we may simply state

zonTube =

{zoffset, for zoffset ≥ 0

Ltube + zoffset, else(15)

Figure 11a shows the error cone projection along the tube. Since we assume the errorto be normally distributed with standard deviation σ we find the probability for anaxial segment of the tube to experience irradiation from the mirror-row starting atorigin O to behave according to the cumulative distribution function of the normaldistribution as introduced in Section 2.10. The more accurately we can predict thephysical behaviour and the smaller our error angle σ the more narrow our band ofpossibly hit segment gets. As one would expect the distribution converges against astep function (as shown in 11b where all irradiated potential is sensed at the tubeinstantaneously for an ideally reflected ray.

11

(a) Error cone in z-direction (b) CDF of the normal distribution

Figure 11: Steps in computing the energy distribution along the z-axis. The cumulativedistribution function converges against the Heaviside step-function for limσ → 0 (green)

2.12. Computation of irradiative power

The concepts explained in the sections above may now be combined to obtain an ex-pression for the irradiative heat flux Q in [W ] and the irradiative heat flux density q in[W/m2] sensed at the absorber tube. Our ray tracing method allows for an accurateprediction of irradiation conditions along the circumference as computed in the 2Dplane and we already described the power associated with mirrors in the 2D plane inSection 2.6.

Now the question remains how to translate this instantaneous 2D distribution into a3D distribution. Since we already described a hit probability along the third dimensionof the absorber tube in Section 2.11.3 we may now directly obtain Q along the thirddimension by integrating the probability along the z-axis of the tube.

Using the probability distribution for the mirror-row j Pj(z) (see also: Figure 12) wefind the 3D heat flux for a segment i as described in Figure 12a to be

Qtot,ij = Qxy,Mirror ·∫ b

a

Pj(z)dz ≈ Qxy,Mirror · Pij · (b− a) (16)

In this equation the integral over the continuous hit probability for the rays reflectedby mirror j, Pj(z) = Pj(γ) is approximate using the discrete averaged hit probability

Pij(z) =P (γ1) + P (γ2)

2(17)

where the z-position is expressed using the segment boundary angles γ1 and γ2 asshown in Figure 12.

12

(a) Continuous hit probability P (b) Averaged hit probability P

Figure 12: Sample parameters for the computation of the heat flux on the axial segmenti with bounds a and b under the hit probability P(z). The segment positionrelative to the hit point of an ideal ray is characterized using the bound-angles γ1 and γ2

Since the surface area for any cylindrical axial segment of the tube is readily computedwe may finally find the surface heat flux in [W/m2] from

qtot,ij =Qtot,ij

πr2 · (b− a)(18)

where the indices i and j refer to an axial segment and a mirror row respectively.Parameters a and b describe the bounds of the axial segment i along the axial directionsuch that (b − a) is the length of the segment along this direction or the height of acylindrical tube segment with radius r.

2.13. Implementation

The goal of the implementation was to achieve functional segregation based on the realplant components which could naturally be achieved using an object oriented approach.A simplified class diagram in Figure 13 reveals key classes and functionalities thatmodel the power plant behaviour in our FREDRay (Fresnel Deterministic Raytracer)software. A helper class FPPConfiguration is used to either conveniently set modelparameters or import them from a settings file in the json format. A time step can besimulated using only few intuitive function calls:

FPPConfiguration p lantConf ig ;FresnelPowerPlant fpp ( p lantConf ig ) ;fpp . traceRaysFPP ( ) ;fpp . computeZDistr ibut ion ( ) ;cout << fpp . pr intEnergy ( ) << endl ;

A container class HeatDataContainer is provided for the export of simulation data.Crucial data such as the irradiation on the discretized tube grid may be exported to atext file for a convenient usage as a boundary condition in the solver for the internalflow.

13

Figure 13: Simplified Class diagram (only basic functionalities are shown to presentthe modeling framework)

2.14. Case study

Please note that all following simulations have been done using aDirectNormalIrradianceof 1. The net power distribution shown in the plots below is therefore to be regardedas a qualitative validation.

2.14.1. Test setup

Having implemented all components of the 3D ray tracing software we may finallypresent the results of several testcases. First of all we validated the geometry of thereflection using a geometric configuration similar to that of the PE1 power plant[16].Figure 14 shows a sample image of the projection of ray traces in the x-y plane. Blockedor shaded rays are shown in red.

14

Figure 14: Reference geometry following the geometry of the PE1 power plant

Figure 15: Visualization of the 2D ray tracing algorithm to verify the geometry ofreflection for different incident angles

2.14.2. Consistency verification with SolTrace

The accuracy and physical consistency of the model has been verified using the existingMonte-Carlo ray tracing software SolTrace.The validation testcase uses a mirror configuration of 16 equidistant mirrors of equalwidth w = 0.75m and an aborbertube diameter d = 0.1m. The ”ideal” solution istaken from a Monte-Carlo simulation with 160.000 rays. As depicted in Figure 16aour code shows excellent agreement with SolTrace even for lower ray counts. Thecomputation time is shown in Figure 16b with simulations being run on the sameprocessing environment. Note that the computation in SolTrace is done in paralleland in multiple runs (n=20 for the following plots) to achieve an average value withlow variance.

15

(a) Accuracy for increasing ray counts(b) Computation time for increasing ray

counts

2.15. Accuracy against Ray count

The number of rays per mirror is the parameter with the highest influence on computa-tion time and we are thus interested and the degree of accuracy achieved by a certainnumber of rays with respect to an ”ideal solution”. The comparison is done usingthe total irraditive power on a segment. Simulations for several elevation angles haveshown that a deviation below 1 percent is generally achieved for less than 100 rays permirror as shown in Figure 17. It is notable that a high resolution of the mirror surfaceespecially increases accuracy when the tube diameter is significantly smaller than theprojected mirror surface. This is true for the real implementation of the power plant.

Figure 17: Absolute deviation from a simulated ideal solution using 10000 rays permirror for different raycounts. The mirror configuration corresponds to thelayout shown in Section 2.14.1

16

2.15.1. Circumferential power distribution

For a given sun position (ε = 60◦, τ = 0◦ in Figure 18) we may now obtain a circum-ferential distribution of irradiative power for every axial segment. If the azimuth angleis zero the circumferential distribution will be equal for all axial segments.

(a) Visualisation of the angular direction(b) Qualitative power distribution around the

circumfence

Figure 18: Qualitative plot of the circumferential power distribution from 0◦ to 360◦

in counter-clockwise direction with the zero angle being defined in 18a

For a very large number of rays the distribution converges against a continuous curveas shown in Figure 18b.

2.15.2. Changes in the mirror position

For a fixed sun position we may vary the mirror’s x-Position relative to the tubeto observe the development of the z-offset as we move the mirror further away fromtube. Figure 19 shows how changes in the z direction become more visible for largerdeviations from the origin straight below the tube and for smaller elevation angles.

17

Figure 19: The plots show the variation of the offset in z direction as a single mirroris moved away from the position directly under the tube (X Pos=0) to anoffset of 9m under a fixed sun-position.

2.15.3. Changes in the azimutal component

Having obtained a solution for the offset along the tube we can now quantify the irra-diation behaviour with emphasis on the transition region. Since azimuth and elevationare not varying independently over the course of a real sun trajectory we chose tosimulate a sunrise with real sun data (changing elevation and azimuth). Figure 20clearly shows how independent rows of mirrors cause distinct offsets on the tube whichresults in a stairway-shaped step function. The effect becomes less pronounced forhigher elevation angles ε, that is during solar high noon.Quantitatively it is striking that the overall offset is in the range of a few meters atmost. For industrial plants where the absorber tube can easily be several hundreds ofmeters long this might barely corresponds to a single axial cell in a 1D flow discretiza-tion. Thus the effect of the z-offset on the flow solution can be expected to be rathernegligible.

18

Figure 20: Variation of the transition zone over the course of a sunrise on june 21st2017 at PE1, Calasparra, Spain. The simulation was run using three samplemirrors at x=[0, 3.9375, 7.875] [m] of equal width w = 0.75m

2.15.4. Solution on a 3D discretization

The integration of the 3D solver allows for a computation of heat fluxes on a fulldiscretization (axial and circumferential) of the tube. Figure 21 shows the 3D distri-bution on a discretization along the circumference and along the tube’s axial direction.The distribution along the circumference corresponds to the one simulated in Section2.15.1. It is clearly visible how increasing the azimuth angle shifts the distribution tothe right in the direction of increasing z-offset.

19

(a) Azimuth 5 ◦, Elevation 60 ◦ (b) Azimuth 10 ◦, Elevation 60 ◦

Figure 21: Irradiation condition Qopt at the absorber tube simulated using the geometry14. Discretization is done along the circumference (angular) and the z-Axis(axial). The shift along the axial direction is evident for more dominantazimutal components (Fig. 21b)

3. Thermodynamical Model

We need to describe the flow condition of a two phase flow inside of the absorbertube. It is important to have a fast simulation code, as it may also be used in futureapplications e.g. the influence of cloud movements on the state of the flow. This meansthat the simulation should be able to be run in real time.It is also important to have conservation of mass, momentum and energy. This is whywe are using a finite volume method. The problem is reduced to just one dimension, asthis is enough to describe the problem in adequate accuracy while obtaining a suitableruntime.

3.1. Prior Work

Up until this point code for simulating the thermodynamics inside an absorber tube ofa solar power plant mostly consists of two or three dimensional finite volume code [14][21]. We break the problem down to one dimension and trade accuracy for runtime.Another approach also solves the problem in one dimension and couples the finitevolume code directly with the water equations and the solar irradiation [23]. We usethis approach and use even more exact water equations and a deterministic ray tracermodel combined with finite volume code optimized for this specific problem to obtaina faster runtime. Since it is uncertain which of the known solvers is going to be the

20

most accurate given a certain computation time for this problem, we need to comparedifferent solvers.

3.2. Finite Volume Method

For modeling the thermodynamics of a two-phase flow inside of the absorber tubewe solved the Euler equations, a set of hyperbolic differential equations, using a onedimensional finite volume method. In the two-phase region, we employ a HomogeneousEquilibrium Model where the two phases are assumed to be in thermodynamical andmechanical equilibrium and they are spatially homogeneous. These are the assumptionsthat enable us to use a three equation model for the two-phase region. The threeequation system of the Euler equations is given by

∂u

∂t+∂f(u)

∂x= s(u) (19)

with the state vector

u =

ρρvρE

,

the analytical flux function

f(u) =

ρvρv2 + p

(ρE + p)v

,

the source terms s(u), which are going to be presented in section 3.7 and the Jacobianmatrix

A(u) =∂f(u)

∂u=

0 1 0

−v2 + (p)ρ − (p)e · E−v2

ρ2v + −(p)e·v

ρ(p)eρ

v((p)ρ − (p)e · E−v2

ρ− p

ρ− E) E + p−(p)e·v2

ρv + v+(p)e

ρ

(20)

with the three eigenvalues v − c, v and v + c.ρ represents the density, v the velocity, e the inner energy, E the total energy, c thespeed of sound and p the pressure.

By discretizing in time and space we get the update formulation

un+1i = uni −

∆t

∆x· (Fn

i+ 12− Fn

i− 12) (21)

where the superscript with a lowercase n represents the time step and the subscripti = 1, ..., N represents the cell index. The cells with the index 0 and N + 1 are theghost cells, which do not get updated by the formula above, but by the specification

21

of the boundary conditions.This is the conservative form for a finite volume method, with F being the numericalflux function, representing the flow over the left (subscript i − 1

2) or right (subscript

i+ 12) cell boundary. Every solver used in this project work can be written in conser-

vative form. The pseudo-code of the algorithm looks like this:

while (t < tend)compute ∆t via CFL condition;

for i = 1, ..., Nupdate ui;

apply source terms;

update boundary conditions;

t+=∆t;end

The time step ∆t is computed via the CFL condition, meaning

∆t ≤ ∆x

λmax(22)

where ∆x is the length of one cell and λmax is the largest absolute eigenvalue of theJacobian matrix A(u).

3.3. Solvers

3.3.1. Lax-Friedrichs

The classical Lax-Friedrichs method is forward in time and centered in space [23]. Itsnumerical flux function is given by

FLF(uL,uR) =1

2(f(uL) + f(uR)− ∆x

∆t· (uR − uL)). (23)

When writing both of the numerical flux functions down, we notice that terms eliminateeach other, so we can write the update formulation as

un+1i = uni −

1

2

∆t

∆x· (f(uR)− f(uL)− ∆t

∆x· (uni−1 − 2uni + uni+1)). (24)

3.3.2. Rusanov

The Rusanov or local Lax-Friedrichs method uses the maximum wave speed at everycell interface [23]. Its numerical flux function is given by

FRusanov(uL,uR) =1

2(f(uL) + f(uR)− λmax

local · (uR − uL)) (25)

withλmax

local = max{λmax(|A(uL)|), λmax(|A(uR)|)}. (26)

22

x

tSRSL

0

u∗

uRuL

Figure 22: Propagation of information at a cell face for the HLL solver

3.3.3. HLL

The HLL method uses the fastest wave speeds in each direction of flow λminlocal and λmax

local

of the system and solves a Riemann problem [23] [29]. The HLL flux function is givenby

FHLL(uL,uR) =

f(uL) if 0 ≤ SL

F∗ if SL ≤ 0 ≤ SR

f(uR) if 0 ≥ SR

(27)

withSL = min{0, λmin

local(A(uL)), λminlocal(A(uR))} (28)

andSR = max{0, λmax

local(A(uL)), λmaxlocal(A(uR))} (29)

representing the propagation speeds of the information in our system.

Considering we are simulating a subsonic flow, we notice that the smallest eigen-value of our Jacobian matrix, which is given by λmin = v − c is always going to besmaller than 0, as v is always going to be smaller than c.Thus SL < 0 and SR > 0 holds, meaning that the first and the third case are notneeded, which leaves us at

FHLL(uL,uR) = F∗ (30)

With the Rankine-Hugoniot conditions we then obtain the numerical flux function

FHLL(uL,uR) =SRf(uL)− SLf(uR) + SLSR(uR − uL)

SR − SL

. (31)

3.3.4. HLLX

The HLLX method [26] uses the same amount of information about the system asHLL does, meaning it only needs the slowest and the fastest wave speeds λmin

local andλmax

local. Its flux function is a weighted combination of the known flux functions fromthe Lax-Friedrich and the HLL method and additionally uses the Lax-Wendroff flux

23

function. Lax-Wendroff is a non-monotone second order method in space and time.Its numerical flux function is given by

FLW(uL,uR) = f(1

2(uL + uR)− 1

2

∆t

∆x(f(uR)− f(uL)). (32)

Additionally to the mentioned three numerical flux functions we need to compute f ,which is given by

f =1

2(f(uL) + f(uR))

as well as νmin and νmax, which are the global minimum and maximum of the Courantnumber, as known given by νi = λi

∆t∆x

and the weighting coefficient α, given by

α =νmax − νmin − ||νmax| − |νmin||

(νmax − νmin)2. (33)

We now can compute the flux function of the HLLX method (using equations (23),(31) and (32)) by

FHLLX(uL,uR) = FHLL + α{|νminνmax|(FLF(uL,uR)− f)− (|νmax|+ |νmin|)(FHLL(uL,uR)− f + (FLW(uL,uR)− f)}.

(34)

3.3.5. HLLC

The HLLC method is a modified version of the HLL method, which also takes thecontact (this is what the ”C” stands for) on surfaces and shear waves into account. [29]Its numerical flux function is given by

FHLLC(uL,uR) =

f(uL) if 0 ≤ SL

F*L if SL ≤ 0 ≤ S∗

F*R if S∗ ≤ 0 ≤ SR

f(uR) if 0 ≥ SR

(35)

with SL and SR being the velocities defined in the equations (28) and (29) and S∗being a velocity we are looking to compute. Because we only are considering subsonicflow we can neglect the first and the fourth case, as described in the HLL section. Forcomputing F*L and F*R we once more use the Rankine-Hugoniot conditions and get

F*L = f(uL) + SL(u*L − uL), (36)

with u*L being our state vector in the intermediate state,

u*L = ρL(SL − vL

SL − S∗)

1ρ*Lv*L

EL + (S∗ − vL)(S∗ + pLρL(SL−vL)

)

(37)

24

x

tSRSL

S∗

0

u*Ru*L

uRuL

Figure 23: Propagation of information at a cell face for the HLLC solver

in which we have a flux function given by

f(uL) =

ρ*Lv*L

ρ*Lv2*L + p*L

(ρ*LE*L + p*L)v*L

. (38)

We can compute f(uR) analogously.

We notice that we have four unknowns in our analytical flux function, while we onlyhave three equations, which means that we need one equation more. We use theassumption that

p*L = p*R = p∗. (39)

Additionally we know that v*L = v*R = v∗ with v∗ = S∗ holds.When we rearrange equation (36) we get

SLu*L − F*L = SLuL − f(uL), (40)

where the left side only contains unknowns, while the right side only contains knownvariables. We split up equation (40) into its vector components, while we only use thefirst and the second component and get

SLρ*L − ρ*Lv*L = SLρL − ρLvL (41)

SLρ*Lv*L − (ρ*Lv2*L + p*L) = SLρLvL − (ρLv

2L + pL) (42)

Out of this two equations we get

ρ*L = pL + ρL(SL − vL)(S∗ − vL). (43)

Analogously we can do this for the right side.We then use equation (39) and receive our needed velocity S∗ after rearranging:

S∗ =pR − pL + ρLvL(SL − vL)− ρRvR(SR − vR)

ρL(SL − vL)− ρR(SR − vR), (44)

which only requires known variables. We now know the if-decider from equation (35),as we know SL (28), SR (29) and S∗ (44), u∗ (37) and F*L as well as F*R (36). Fromthere we can compute the numerical HLLC flux function using (35).

25

x

tv + cv − c v

0

uL

u1 u2

uR

Figure 24: Propagation of information at a cell face in the exact system

3.3.6. Godunov

The Godunov method follows another approach compared to the the solvers mentionedabove, as we had approximate solvers of an exact system and now are using an ap-proximated model of the system, which we can solve exactly [23].The numerical flux function for Godunov is given by

FGodunov(uL,uR) =1

∆t

∫ tn+1

tnf(u(x, t)) dt (45)

with x being the value in space between the left and the right cell and tn+1 = tn+∆t. Asseen in figure 24 the exact solution of the Riemann problem is constant with progressingtime, meaning the integral can be solved exactly, which also results in an exact solutionof the approximated model.The numerical flux function thus is given as analytical flux function of the Riemannsolution uRiemann,

F(uL,uR) = f(uRiemann(0;uL,uR)) (46)

with uRiemann being either u1 or u2, depending on v being greater or smaller than 0.

To find u1 and u2 we extend our two-phase flow model using Suliciu-relaxation toobtain a linearly degenerate model as done in [23]. We now have a state vector u andan analytical flux function f(u) containing four components – the three componentsof the Euler equations and one component representing the pressure of the relaxatedmodel.

u =

ρρvρEρπ

f(u) =

ρv

ρv2 + p(ρE + p)vρvπ + a2v

When computing the Jacobian matrix A(u) := ∂f(u)

∂uof this relaxated homogeneous

equilibrium system we still want the eigenvalues to be the same as in the homogeneousequilibrium system. This is given if for the parameter a ≥ ρc holds, as the eigenvalues

26

are given by

λ1 = v − a

ρ≤ λ2 = v ≤ λ3 = v +

a

ρ.

Thus, we choose a = max{ρLcL, ρRcR}. It is important to use the eigenvalue v + aρ

forcomputing the time step.

Using the Rankine-Hugoniot conditions at the discontinuities we get eight equationscoupling our known states uL and uR with the unknown states u1 and u2, which con-tain the unknowns ρ1, v1, u1, and π1 and analogously for state u2. From there we cansolve for our unknowns and receive the solutions

u1 =

ρ1

ρ1v1

ρ1E1

ρ1π1

and u2 =

ρ2

ρ2v2

ρ2E2

ρ2π2

with

v∗ = v1 = v2 =vR + vL

2− pR − pL

2a

π∗ =πr + πl

2− avR − vL

2

ρ1 =aρL

a+ ρL(v∗ − vL)

ρ2 =aρR

a+ ρR(vR − v∗)

u1 = uL +π2∗ − π2

l

2a2

u2 = uR +π2∗ − π2

r

2a2

Now we can compute u1 and u2, set uRiemann by using the known if-decider and computethe analytical flux function of uRiemann, which results in the numerical flux function.

3.4. Analysis of computation time and accuracy

As introduced in section 3.2 every solver presented in section 3.3 can be written inconservative form. This means that the code of every solver – in both Matlab andC++ – looks quite similar to the pseudocode mentioned in section 3.2, containing awhile-loop over the time and a for-loop iterating over the cells. When parallelizingthe different solvers, the biggest part of the speed up comes from parallelizing thefor-loop, which also is the point where the solvers differ most from each other.It is important to note though that the tendencies of the computation times of thesolvers relative to each other remain the same independently from the level of paral-lelization, programming language and system they are run on, as long as those three

27

parameters are constant for every solver while measuring the runtime. Because of thisit is possible to do the test with the unparallelized code in Matlab, whose result weare going to present to you in the following.

3.4.1. Sod’s shock tube test

Sod’s shock tube problem is a commonly used test for numerical fluid solvers. Thelength of the tube is one meter and the simulation time is 0.2 seconds. For setting theinitial conditions, the tube gets split in two equal halves. The primitive variables onthe left and right hand side are:ρL

pL

vL

=

1.01.00.0

and

ρR

pR

vR

=

0.1250.10.0

,

which leads us to the initial conditions for our state vector u

uL =

ρL

ρLvL

ρLEL

=

1.00.02.5

and uR =

ρR

ρRvR

ρRER

=

0.1250.00.25

.

For boundary conditions we set the values of the ghost cells as the values of theirneighbor cells:

u0 = u1and uN+1 = uN .

The equations for the pressure p and the speed of sound c are given by the stiffenedgas equations of state.

After executing Sod’s shock tube test we obtain results for the density, pressure, ve-locity and specific internal energy. In figure 25 you can see the density over the lengthof the tube for N = 50 cells.

28

Figure 25: Plot of the density vs. the length for 50 cells

3.4.2. Accuracy versus number of cells

In figure 25 it is noticeable that the Lax-Friedrichs and Rusanov methods are heav-ily smearing out, as they do not use much information about the system itself. Thismakes them fast but not as exact solvers. Rusanov though is more accurate than Lax-Friedrichs, as it is using the local wave speeds instead of one single global one per timestep.The HLL solver is better, because it uses two wave speeds per interface; even better isthe HLLC solver, which was expectable as it is taking the contact into account. Also ituses an additional wave speed, which means it also uses one state more, approximatingthe solution more precisely. HLLX differs a lot from the analytical solution as it is amixture between the Lax-Friedrichs, HLL and the unstable Lax-Wendroff method. Alittle more exact than the HLL method is the Godunov-Suliciu solver.

For comparing the different methods we computed the 2-norm as well as the infinitynorm of the error and the root-mean-square deviation (RMSD). The result is shown intable 1. The error εi for the value of one cell i is given by

εi = |ρnumericali − ρanalytical

i |.

The 2-norm of the error is defined as

||ε||2 =

√√√√ N∑i=1

ε2i ,

the infinity norm as||ε||∞ = max{ε1, ..., εN}

and the RMSD as

||ε||2 =

√√√√ 1

N

N∑i=1

ε2i .

29

Figure 26: Accuracy vs. number of cells for the different solvers

The accuracy is defined as 1− εrelative with

error measure HLLX Lax-Friedrichs Rusanov HLL HLLC Godunov-Suliciu2-norm 0.3210 0.3441 0.3101 0.2454 0.2273 0.2389infinity-norm 0.1149 0.1167 0.1155 0.0903 0.0863 0.0899RMSD 0.0454 0.0487 0.0439 0.0347 0.0321 0.0338

Table 1: Deviation of the numerical solutions from the analytical solution for N=50cells

εrelative =||ε||2||ρ||2

being the relative error and ||ρ||2 being the 2-norm of the vector, which contains theanalytical value of the density for every cell. When we plot the accuracy versus thenumber of cells as done in figure 26, we can see that the most precise solver for anynumber of cells is the HLLC solver. Almost as precise are the Godunov-Suliciu andthe HLL method.

30

Figure 27: Accuracy vs. number of cells for the different solvers

3.4.3. Accuracy versus CPU runtime

To choose the solver to work with later on, we have to find the solver that is most accu-rate for a given computation time. This is not necessarily the most accurate solver fora given number of cells. The reason for this is that the solvers need different executiontimes for one time step hence for the whole computation. Thus a faster but less precisesolver can perform a computation with a higher number of cells in the same time aslower but more accurate solver needs for a computation with a lower number of cells.This can result in the ”less accurate” solver obtaining a higher accuracy for a fixedCPU runtime.Because of this we plotted the accuracy versus the computation time. The result isobservable in figure 27. You can see that the Godunov solver is the best choice regard-ing the computation time, which is why we wrote and parallelized it in C++.

It is important to mention that the implementations of the different solvers in Mat-lab are structured similarly to make sure we obtain a fair comparison of the differentmethods.

31

3.5. Matlab implementation

For the implementation of the solvers in Matlab we tried out different approaches andcompared them to make sure we optimize the solvers as much as possible before com-paring the runtimes. In this subsection we present you the most important approaches.

3.5.1. Lax-Friedrichs

When writing the update formulation of the Lax-Friedrichs method as done in equa-tion (24) instead of using the formulation of each numerical flux function as done inequation (23), we reduce the number of necessary floating point operations. This ap-proach is especially faster than at first computing all numerical flux functions, becausememory accesses need significantly more time than floating point operations.

Additionally we computed the left and middle analytical flux function before the for-loop and the right analytical flux function inside of the for-loop, handing over thevalues of the middle analytical flux function to the left one in the next step in spaceand analogously the right to the middle one.This also speeds the code up, as we neither have to compute the left, middle andright analytical flux function in every step in space nor have we to precompute everyanalytical flux function before the for-loop, which would not be as slow as the firstapproach, but which would be more costly in memory requirements. By doing this, wesave up to 15% computation time, depending on the number of cells. This method isgoing to be used in the following solvers as well.

3.5.2. Rusanov

As the right analytical flux function in step in space i becomes the left analyticalflux function in step i+ 1, one could think that handing over the right numerical fluxfunction to the left one in the next step in space and thus only computing each differentnumerical flux function once could be an idea to save computation time. This is notthe case, because the memory accesses need more time than the saved floating pointoperations would need. Thus the numerical flux functions get computed separately inevery step.

3.5.3. HLL

From a implementation point of view there is nothing special to mention for the HLLapproach besides keeping in mind that we do not have to use the if-decider as writtenin equation (27), but can directly use equation (31) to compute the analytical fluxfunction.

32

3.5.4. HLLX

From a implementation point of view there is nothing special to mention for the HLLXapproach.

3.5.5. HLLC

The reminder for the HLL solver also holds for the HLLC solver. Additionally wecompute F*L and F*R not beforehand but inside the if-case to save a little bit ofcomputation time. We also know that SL < 0 and SR > 0 always holds, meaningthat our only deciding factor in the if-case is S∗, thus we do not have to do as muchcomparisons.

3.5.6. Godunov-Suliciu

For the standard Godunov-Suliciu solver there is nothing special to mention. However,since we found that this solver is the most efficient one of our solvers we tried a newimplementation approach, in which we precompute the numerical flux function at everycell face first and then we update every cell. The floating point operations we save bydoing this by far compensates the time we need for the memory accesses. We did notuse this approach in the accuracy plot pictures 26 and 27, but used it from now on asit is even faster.

3.6. Parallelization

In order to parallelize the Godunov-Suliciu solver in C++ we

• vectorized the computation of the numerical flux functions and the updatingprocess itself as described in section 3.5.6,

• used the Eigen library to optimize vector operations and

• used OpenMP to parallelize the for-loops for computing and applying the bound-ary conditions, computing a and λmax, updating the cells and applying the sourceterms.

33

3.7. Source Terms

A fractional-step method is employed for the consideration of source terms. Twosub-steps are solved sequentially. First the evolution step solves the homogeneoushyperbolic PDE from the current value uni [23]

∂tu + ∂xf(u) = 0 .

The resulting value un+1,−i is used in solving the ODE with the same time step as the

PDE above.∂tu = s(u)

which is done here with a simple Euler step.

un+1i = un+1,−

i + ∆t · s(un+1,−i ) .

The source terms are calculated as

s(u) =

0−Fr

−Fru+QtoWater

πr2in∆x

Two effects are captured in the source terms: friction and heat flux.

Friction: The friction density in our simulation is calculated as detailed in [3] as

Fr =fρu|u|

2dn

where f is the Darcy friction factor, dh is the hydraulic diameter with dh =4Apipe

Pwetwith

Pwet as the wetted perimeter. For a homogeneous flow within a cylindrical pipe thissimplifies to dh = 4πr2

2πr= Din and can thus be identified with the inner diameter of the

pipe. The friction factor is computed in dependency of the Renould number (Re):

f(Re) =

1, if 0 ≤ Re < 64

flam(Re) = 64Re, if 64 ≤ Re < 2200 laminar

flam(2200)w(Re) + (1− w(Re))fturb(3000), if 2200 ≤ Re < 3000 transient

fturb(Re) else

The weighting factor w(Re) is calculated as

w(Re) =3000(Re−2200)

(3000− 2200) Re.

Lastly fturb(Re) is given by:

1

fturb(Re)= −2log10

(ε

3.7Din

+2.51

Re

[1.14− 2log10

(ε

Din

+21.25

Re0.9

)])were ε is the surface roughness and Din is the pipe inner diameter.

34

Figure 28: Pipe wall acts as a buffer for heat

Heat Transfer: The pipe wall acts as a heat buffer. As illustrated in Figure 28the temperature of the pipe is increased by the sun rays focused on the pipe. Atthe same time heat is lost by radiation and transferred to the water. Following theprocedure in [23] the heat loss by thermal radiation is modeled using a polynomial ofthe temperature in Celsius:

Qloss = 0.16155W

K m(T − 273.15K) + 6.4407 · 10−9 W

K4m(T − 273.15K)4∆x .

The heat added from the sun rays Qext is calculated using the simulation from theoptical model described earlier.

The amount of heat transferred from the pipe to the water is calculated accordingto Fourier’s First Law [22]:

QtoWater = κ(T − TWater)(Ain)

with Ain the contact area between water and pipe:

Ain = 2πrin∆x .

The thermal transmittance κ is calculated as:

κ =λwallα

λwall + ln(routrin

)rinα

.

The entity α describes the heat transfer coefficient

α = 0.0235λ

2rin

Re0.8 Pr0.4

with Pr the Prandtl number, λ the thermal conductivity of the fluid and λwall thethermal conductivity of the wall.

35

We assume a constant temperature in radial direction and that heat capacity anddensity of the pipe are constant over the pipe. If we model the heat transfer only inradial direction [22] the temperature of the pipe wall T can be updated with

T n+1 = T n + ∆t (Qext −Qloss −QtoWater)1

ρWallVWallcpWall

.

4. Closure of the two-phase flow model

A thermodynamic closure for the finite volume solver is needed. We require highaccuracy to avoid error cumulation and fast computation for the large number ofqueries generated within.

4.1. Mixture Model

In the two-phase region, we employ a Homogeneous Equilibrium Model (HEM) wherethe two phases are assumed to be in thermodynamic equilibrium and they are spatiallyhomogeneous. This is the assumption that enables the usage of the three equation Eulerequations for the two-phase region.

Following [3] the thermodynamic properties of the pseudo fluid are modeled as anvapour fraction weighted average of the two phases. Properties like thermal conduc-tivity and viscosity collapse to a single value as a mass-weighted average:

µ = αµvap + (1− α)µliq (47)

Speed of Sound is averaged based on Wallis formula:

1

ρc2=

α

ρvapc2vap

+1− αρliqcliq

(48)

Density and inner energy are the primitive variables of the Euler equations, butare rarely used as primary variables for equations of state (EOS). Since we enter thetwo-phase region in our simulation and the vapour fraction of our system is unknown,we have to reliably find for a given density and inner energy of the whole mixture,the corresponding vapour fraction, temperature and pressure for an equilibrium state.This is complicated by the fact that most equations of state give separate equationsfor different phases. The treatment of this problem will be the topic of the remainderof this section.

4.2. Prior Work

Prior Work in this field has focused on developing accurate but fast approximationsto the IAWPS95 equations [7] (which are formulated as a function of temperatureT and density ρ rather than density ρ and inner energy e), or efficiently computingthe vapour fraction of the system. In [18] the authors derived an adaptive version of

36

the stiffened gas equation for single-phase liquid water in multiple temperature andpressure intervals. First the inner energy of the system was fitted to a bilinear functionof pressure and specific volume to determine the coefficients of the model,

e = Apv +Bv + C (49)

A =1

γ − 1(50)

B =γp∞γ − 1

(51)

C = q. (52)

Then the model parameter cv was fit to the temperature data

T = D(e− q − p∞v) with D =1

cv. (53)

This method works well for liquid water at high pressures, but has “unacceptably highrelative error levels“ near the saturation line. We conducted a fitting approach similarto the presented liquid region fit, but using an additional parameter (the co-volume)with the pseudo linear regression idea presented in the aforementioned paper and alsowith the least squares fitting method for stiffened gas parameters presented in [12].While good data fits were possible for low temperature ranges, large errors occurredfor high temperatures as the vapour phase became difficult to fit to the model. Oneexplanation might be that the model will predict a monotonic relationship betweentemperature and inner energy, but since the inner energy of saturated steam decreasemonotonically with temperature above 240 Celsius the model displayed unphysicalbehavior. An negative value for the coefficient γ would for example lead to an imaginaryspeed of sound, so the parameters are also constrained to predict realistic behavior.Meanwhile, even with the fitted parameters, high errors occurred at temperaturesabove 500K for the approach.

In [11] and the corresponding IAPWS release, a spline based interpolation methodwas presented, to allow fast computations on water properties with inverse functionswith high accuracy. At the present time the software to construct those splines Fluid-Spline is not available to the public and the manual construction of the splines needscareful placement of the sampling points with different transformations and samplingdensities chosen at different parts of the corresponding diagram.

Chiapolino et al [4] presented a framework in the HEM context (extended to multi-component flows in [5]) for calculating the changes in vapour fraction t (mass-transport)not immediately in each time step, but by integrating it into the Euler equations, con-verging it over multiple time steps. For this the system is extended by

∂ραl∂t

+ div(ραlu) = 0. (54)

The last term is later relaxed to

∂ραl∂t

+∂(ραlu)

∂x= ρλ(gg − gl). (55)

37

Additionally an algorithm is determined to solve this equation robustly and how togenerate robust starting values. One such idea is to compute the mass fractions of theliquid phase according to current estimate of volume and inner energy and use thoseestimates as starting values for the liquid fraction (Yl = 1− α).

Defining

Y ml (p) =

v − vg(p)vl(p)− vg(p)

(56)

Y el (p) =

e− eg(p)el(p)− eg(p)

(57)

r =Y ml (p)− Y prior

l

Y el (p)− Y prior

l

(58)

(59)

an estimate is given by Yl for r < 0, Y el for r > 1 and Y m

l else. Equilibrium pressurecan also be approximated.

4.3. Noble-Able Stiffened-Gas Equation of State

Le Metayer et al. [12] presented an approach for a convex equation of state for water.Based on experimental data for the saturation curve and a liquid phase reference point,several model parameters are fitted with a least square method to the thermodynamicmodel.

4.3.1. Equations

In the following we use τ as the specific volume and e as the specific inner energy. TheNobel-Abel Stiffened-Gas equation of state (NASG-EOS) for one phase is:

p(τ, e) = (γ − 1) · (e− q)(τ − b)

− (γ · p∞)

Equations for other thermodynamic properties are derived in [12]:

e(p, τ) =p+ γp∞γ − 1

(τ − b) + q (60)

τ(p, T ) =(γ − 1)cvT

p+ p∞+ b (61)

h(p, T ) = γcvT + bp+ q (62)

g(p, T ) = (γcv − q′)T − cvT ln(

T γ

(p+ p∞)γ−1

)+ bp+ q (63)

c(p, τ)2 =γτ 2(p+ p∞)

(τ − b)(64)

38

Wherein γ,p∞,cv,b,q and q′ are parameters to be determined.From these equations we can also derive an equation for T:

T (τ, e) =e− p∞ · (τ − b)− q

cv. (65)

With this formulas we are able to compute all thermodynamic properties starting fromτ and e. It is important, however, to note that these equations are valid for a singlephase and the two phases have to be coupled to calculate two-phase behavior.

4.3.2. Two-phase region

With the assumption of a homogeneous system(p′ = p′′ ,g′ = g′′ ,T ′ = T ′′), the followingcorrelation between vapour pressure and temperature can be derived:

ln (p+ p∞,g) = A+B + E · p

T+ Cln(T ) +D ln(p+ p∞,l) (66)

With:

A =cp,l − cp,g + q′g − q′l

cp,g − cv,g

B =ql − qgcp,g − cv,g

C =cp,g − cp,lcp,g − cv,g

D =cp,l − cv,lcp,g − cv,g

E =bl − bgcp,g − cv,g

Since we assumed thermodynamic equilibrium, we have the following system of equa-tions that we can solve for (p, T and the vapour fraction a):

ln (p+ p∞,g) = A+B + E · p

T+ Cln(T ) +Dln(p+ p∞,l)

τ = τl(T, p) + a · (τg(T, p)− τl(T, p))e = el(T, p) + a · (eg(T, p)− el(T, p))

This system is solved with a quasi-Newton-method described in [13]. The viscosity,surface tension and thermal conductivity are implemented as specified in the corre-sponding IAPWS releases [19][10][6]. As recommended, the viscosity calculation ne-glects the critical enhancement, as it is only measurable in the extreme proximity ofthe critical point [6]. The thermal conductivity calculation in contrasts, has a criticalenhancement term, that can in general not be neglected. However, as we will show

39

below, the form of the NASG-EOS guarantees that this term vanishes: The IAPWSrelease [10] defines:

∆X (T , ρ) = ρ

[ζ(T , ρ)− ζ(TR, ρ)

TR

T

](67)

ζ(T , ρ) =

(∂ρ

∂p

)T

(68)

But we can show that ζ is anti-proportional in T using the fact that we assumed ahomogeneous mixture

ρ =1

τ=

1

τg(p, T )a+ τl(p, T )(1− a)=

1

τ(p, T )

leads to

∂ρ

∂p=∂ρ

∂τ· ∂τ∂p

= − 1

τ(p, T )2·(∂τg∂p

a+∂τl∂p

(1− a)

)Since in equation 67, the same density for both calculations of ζ is used, it follows thatτ is also constant and thus we only have to show that δτ

δpis anti-proportional in T .

δτ(p, T )

δp= −(γ − 1)cvT

(p+ p∞)2= − (τ − b)2

(γ − 1)cvT

where we substituted p + p∞ = (γ−1)cvT(τ−b) from equation 61. Since ζ(T , ρ) is anti-

proportional, equation 67 will always yield ∆X = 0 and correspondingly:

λ2(T , ρ) = ΛρcpT

µZ(y) = 0,

because y = kζ(T , ρ) < 1.2 · 10−7 and thus Z(y) is set to zero as per specification.

4.3.3. Thompson Coefficient in NASG

The Joule-Thompson coefficient is defined as(dTdp

)H

. From the NASG-EOS we can

write: h(p, T ) = γcvT + bp+ q which leads to T = h−bp−qγcv

and finally µ = −bγcv

.

40

4.4. IAWPS calculation

For exact calculations, both for validation and for model regression, the IAPWS for-mulation 1995 (Rev. 2016) was implemented [7]. The formulation is based on anaccurate equation for the specific Helmholtz free energy f . This property is expressedas a function of (normalized) temperature and density (τ = T

Tcrit, δ = ρ

ρcrit).

f(ρ, T ) = RTΦ(δ, τ) = RT (Φo(δ, τ) + Φr(δ, τ)) . (69)

All other properties can be calculated from appropriate derivatives of this function,where explicit but algorithmic formulas exist for most of the relevant thermodynamicalproperties. For example for the isobaric heat capacity:

cp = (∂h/∂T )p = R

(−τ 2 (Φo

ττ + Φrττ ) +

(1 + δΦrδ − δτΦr

δτ )2

1 + 2δΦrδ + δ2Φr

δδ

), (70)

where each of the partial derivatives is given as an expression in terms of τ and δ.For example Φr

δ is given by:

Φrδ =

7∑i=1

nidiδdi−1τ ti +

51∑i=8

nie−δci [δdi−1τ ti (di − ciδci)

]+

54∑i=52

niδdiτ tie−αi(δ−εi)

2−βi(τ−γi)2[diδ− 2αi(δ − εi)

]

+56∑i=55

ni

[∆bi

(Ψ + δ

∂Ψ

∂δ

)+∂∆bi

∂δδΨ

](71)

where Ψ and ∆ are also functions of δ and τ . The full formulation and the coefficientsare available in [7].

4.4.1. IAWPS ζ calculation

In the IAWPS calculation, the calculation of ζ is necessary. In [10] ζ is defined as:

ζ =(∂ρ∂pT

)= − 1

v2· 1

( ∂p∂v )T. Overlined quantities are scaled with a reference value

(e.g. p = p/p∗). The term(∂p∂v

)T

= −p · βp = −RTρ(1 + δ · Φrδ) · ρ ·

[1 +

δ·Φrδ+δ2·Φrδδ

1+δ·Φrδ

]is equal to −ρ2RT · [1 + 2δ · Φr

δ + δ2 · Φrδδ]. We can use this to calculate ζ utilizing(

∂ρ∂p

)T

= p∗

ρ∗

(∂ρ∂p

)T

4.4.2. Two-phase region

In the two-phase region, the phase-equilibrium condition (Maxwell criterion) has to besatisfied. For the IAWPS95 equations, this means that the following equations have tobe satisfied:

Jl =pσRTρ′

= 1 + δ′Φrδ(δ′, τ) (72)

41

Jv =pσ

RTρ′′= 1 + δ′′Φr

δ(δ′′, τ) (73)

pσRT

(1

ρ′′− 1

ρ′

)− ln

(ρ′

ρ′′

)= Φr(δ′, τ)− Φr(δ′, τ). (74)

4.4.3. Calculation of (T, ρ′, ρ′′ and α)

In the finite volume solver, the primitive variables give direct access to the densityof the mixture ρ and the inner energy e. Again, thermodynamic properties of thepseudo-fluid are modeled as an vapour fraction weighted average of the two phasesand temperatures in both phases are assumed to be equal (T ′ = T ′′, τ ′ = τ ′′′). Giventhe density of the mixture ρ and inner energy e, we want to determine the primaryvariables for the IAWPS calculations, namely the temperature T and ρi for both phasesas well as the composition as the (mass) vapour fraction α. With knowledge of allthese quantities, we could evaluate g(T, ρi) for any thermodynamical property g forboth phases i ∈ {liquid, gas} and combine the results according to the vapour fraction:

g (e, ρ) = αg (T (e, ρ), ρv(e, ρ)) + (1− α) g (T (e, ρ), ρl(e, ρ)) (75)

To compute the quantities that are used in this procedure, the following system ofhighly nonlinear equations has to be solved (see 72,73,74):

Jl(T, ρ′)ρ′ − Jv(T, ρ′′)ρ′′ = 0

Jl(T, ρ′)− Jv(T, ρ′′)− ln(

ρ′

ρ′′)− Φr(δ′, τ)− Φr(δ′, τ) = 0

e− e′(1− α) + e′′α = 0

(76)

Solving this system can be stated as solving f(T, ρ′, ρ′′) = 0 because α is given asα = τ−τ ′

τ ′′−τ ′ with τ the specific volume (ρ = 1τ). Since the Jacobian of f is infeasi-

ble to obtain analytically, the above mentioned Broyden algorithm seemed promisingas it avoids direct estimation of the Jacobian. Unfortunately for this problem, theimplemented Broyden and Newton algorithms did not converge because of nearly sin-gular Jacobian approximations and too large step sizes leading the algorithm to tryto evaluate functions with unphysical inputs e.g. negative temperature. In order toaccomplish a robust convergence, two steps where taken: First a Eigens implementa-tion of the Powell Hybrid Dogleg algorithm, a port of the MINPACK-1 [15] was usedfor its employment of a trust region and overall more reliable convergence in for thissystem of equations. Second an initial guess was calculated using the NASG-EOS asdescribed above with the coefficients from the original paper [12] and an initial pointof T = 423.15K, α = 0.5 and p = 476.160kPa.

4.4.4. Detection of two-phase region

Since the water in our model will start in the pure liquid region, cross the two-phaseregion and become overheated steam, we had to check whether or not our assumption

42

Input ρ, e

Assume 1 Phase. Solve e− eIAPWS(ρ, T0) = for T0

ρ ∈ [ρ′′(T0), ρ′(T0)]? Vapour PhaseLiquid Phase

Find initial guess from NASG

Solve Two-Phase Problem for (α, T, ρ′, ρ′′)

ρ < ρ′′ρ > ρ′

true

Figure 29: Flowchart of the phase detection algorithm

to be in the two-phase region is valid. If we are not in the two-phase region, thecalculation is simplified, since we know the vapour fraction to be zero for the liquidregion and one for the overheated steam. Additionally, we do not have to converge theequilibrium condition present in both the NASG and IAWPS95 formulations. Sincewe know that for water the specific volume and energy of the liquid phase will alwaysbe lower than that of the vapour phase, we used that fact to check for the region weare in:

For a temperature T which fulfills h(T, ρ) − e − p(T,ρ)ρ

= 0 i.e. works for the as-

sumption that only one phase is present (α ∈ {0; 1}), it is determined if this state isconsistent with not being on the saturation line. Saturation properties are quickly cal-culated from [27] which gives saturation properties as long polynomials of saturationtemperature. If ρ′(T ) ≤ ρ ≤ ρ′′ and e′(T ) ≤ e ≤ e′′(T ) than we can conclude that theassumption was wrong and we are indeed on the saturation line (but not necessarily atthat temperature), otherwise the actual phase can easily be determined and the calcu-lated temperature is the right one, so that (T, ρ) are known for further computation.The algorithm is visualized in Figure 29.

4.5. Methods under consideration

In this section we want to compare the methods under considerations and assess theiraccuracy and computational speed. The methods are:

• IAWPs Calculations (I)

• IAWPS Calculations with starting value from a look-up table (II)

• Interpolation between results from IAWPS from a look-up table (III)

• Stiffened Gas with parameters for the whole region (IV)

43

• Stiffened Gas with locally fitted parameters (V)

0 1 2 3

·106

10−2

100

102

Inner Energy in J/kg

Den

sity

inkg/

m3

(a) Two-Phase Region in log(ρ) vs. e scaling

[

104 105 106 107

10−2

100

102 ← 275K

← 305K← 335K← 365K← 395K← 425K← 485K← 545K

← 635K

Inner Energy in J/kg

Vol

um

ein

m3/k

g

(b) Two-Phase Region in log(τ) vs. log(e)scaling

0 1 2 3

·106

0

200

400

600

800

1,000← 393K← 453K← 513K← 573K

← 633K

Inner Energy in J/kg

Den

sity

inkg/

m3

(c) Two-Phase Region in ρ vs. e

103 104 105 106 1070

100

200

Inner Energy in J/kg

Vol

um

ein

m3/k

g

(d) Two-Phase Region in τ vs. log(e) scaling

Figure 30: Different scalings of the ρ, e plane

We recognize that method (I) and (IV) will have a disadvantage, since in normaloperation they could be given good initial values from the last time step of the finitevolume solver. These are used more as an guideline on how big the improvement on theother methods is to determine the best method before coupling with the finite volumesolver.

4.5.1. Look-up table

The look-up table was created based on sample points (100x100 in this comparison)taken from a regular grid on the log(ρ) vs. e plane (see Figure 30a). Different scal-ings are compared in Figure 30a and following. The scaling in Figure 30d was clearly

44

unusable for regular sampling. The scaling in Figure 30c makes it hard to differenti-ate between different temperatures at low densities. The plot in Figure 30b has theisotherms very close together at high temperatures and increased spacing at lowertemperatures. The sampled points were then inverted with the implementation of theIAWPS95 equations and saved this yielded (ρ′, ρ′′, e′, e′′, T, p, h′′, h′ and x). More pointscan then be generated efficiently by varying x and mixing the vapour and liquid prop-erties accordingly. All points are saved inside a kd-tree structure which enables fastlook-up and nearest neighbour search.

Interpolation Since pressure is the most important property we will be calculating forwater in the finite volume solver and the most varying over the simulation, interpolationis not just left at a nearest neighbour basis, but threaded more accurately with thin-plate spline interpolation largely taken from [8]. The function has the form

p(e, ρ) = α1 + α2x+ α3y +n∑i=1

wiU (‖Pi − (e, ρ)‖) . (77)

The base function can be written as U(r) = r2log(r) if r > 0 and Pi describes aninterpolation/control point. The linear system to solve for w and α is described as

Ki,j = U(Pi − Pj) (78)

P =(~1n ~e ~ρ

)(79)(

~p~03

)=

(K PPT 03x3

)(~w~α

)(80)

where ~1n is a vector with n ones, 03x3 is a square matrix with three rows and onlycontains zeros. Other entities written in special vector notation are column vectorsconsisting of scalars consistent with the symbol of the vector, ~w for example consistsof all wi.

4.5.2. Locally fit parameters

The methods of [18], [25] and [12] were tested for generating locally fitted parametersof the NASG equations. As mentioned in [18], the first method failed to generatesatisfactory results on or near the saturation line and was therefore discarded. Thesecond method uses just two points, but allows fast determination of the parameters,but due to taking only two points into account and the usage of derivatives, is verysusceptible to noise in the data. Additionally, the time of precomputation was a minorconcern, as long as it stayed in a reasonable realm.

As an result the method of [12] was adapted, but without the assumption that theco-volume of the gaseous phase is negligible (bg 6= 0):

1. p∞l is guessed, p∞g is assumed to be zero

2. τ(T, p) is fitted with least squares to find (cp − cv) and b for both phases

45

3. h(T, p) is fitted with least squares to find (cp and q) for both phases

4. p∞l is then iterated by Ridders’ method to correctly predict the speed of sound

4.5.3. Comparison of methods

The accuracy in the two-phase region and computational cost of the methods intro-duced in the beginning of this section were evaluated. Data from [24] was used (5significant digits). Figure 32 shows the ability of the methods to correctly backcom-pute the thermodynamical state (T, α, psat).

We see that the stiffened gas equations tend to overestimate the vapour fraction.Furthermore, they tend to overestimate the temperature with increasing correct tem-perature, which is amplified at higher vapour fractions. Up to 623K, the fitted NASGmethod generally outperforms the version with global parameters, but becomes unsta-ble at temperatures beyond this point. All methods based on the IAWPS95 formulationrecover a good approximation of the initial state for the whole temperature range.

In Figures 31a and 31b the previous assessment is confirmed. The NASG equationsonly return valid results up to about 500K, but in this range, they are still muchconsiderably less accurate than the IAWPS95 equations.

We also see that the methods based on the IAWPS95 equations that involve iterationyield very accurate result in terms of relative error (less than 0.3%) and only hassignificant absolute error above 550K with a maximal error less than 50kPa.

The purely look-up based method performs very well with less than 3% error in theworst case. Figure 33 shows that this worst case only happens at low temperaturesand when the vapour fraction is below 20%. Corresponding figures are available for the

Temperature [K]

300 400 500 600ma

x.

Pre

ssu

re E

rro

r (r

ela

tive

)

0

0.05

0.1

0.15

0.2IAWPS-95

95-Look-Up

95-Look-Up (refined)

NASG

NASG (fitted)

(a) Comparison of relative pressure error forthe different methods in the Two-Phaseregion over temperature.

Temperature [K]300 400 500 600

ma

x P

ressu

re E

rro

r [P

a] ×10

4

0

2

4

6

8

10IAWPS-95

95-Look-Up

95-Look-Up (refined)

NASG

NASG (fitted)

(b) Comparison of absolute pressure error forthe different methods in the Two-Phaseregion over Temperature

Figure 31: Methods that iterate on the IAWPS equations more accurate over the wholeinterval. Methods based on stiffened gas equations have high error in regionsof higher temperature

46

300 400 500 600 700

Temperature [K]

0

0.2

0.4

0.6

0.8

Massfr

action o

f V

apour

Accuracy Comparision of Water EOSs

Target

IAWPS-95

NASG

NASG (fitted)

95-Look-Up

95-Look-Up (refined)

Figure 32: Accuracy comparison in the α, T -plane

other methods in the appendix, where for the NASG methods for clarity only valuesbelow 0.5 relative and 1.0 · 107Pa are visualized in the contour colors .

In Figure 34 we compare the computational cost of the methods.

Method

Execution tim

e [s

]

0

1

2

3

43.5655

0.0140

0.6843 0.6057

1.1491

Execution time for 180 points

IAWPS95

95Look-Up

95-Look-Up (refined)

NASG

NASG (fitted)

Figure 34: Execution speed.

We observe that using the look up dataas starting values, while generating equiv-alent results, is about five times faster incomputation making it an attractive op-tion. We also see that the fitted NASGwhile improving on the global variant inaccuracy is much slower in computation.This might be in part due to the instabil-ity discovered above leading to excessiveiterations, but part of the performancepenalty comes from having to search andset the parameters from the database. Weconclude, that the fitted NASG approachdid no yield enough benefit over the usingglobal parameters, but increased accuracyfor temperatures between 400K and 550K especially for higher vapour fractions.

The simple look-up method is more than 50 times faster than using iterations forrefinement. If small relative errors are acceptable this method is the most promisingand accuracy can be further improved by using more sampling points. This is themethod we implemented.

47

Figure 33: Relative error of look-up method

4.6. Threadment of the liquid and vapour region

Liquid Water has to be separately sampled and in a small density region since aspictured in Figure 35

it becomes incompressible and τ does not vary much with changes in p. This makesthe computation of pressure from density and inner energy very hard with just pre-computed and locally interpolated points, since p depends very strongly on the specificvolume. One approach is to use the lookup table to calculate the temperature andthen use the accurate function p(T, ρ) since it is still possible to compute temperatureaccurately with the look up method. Still we found that fitting a surrogate functionthat takes into consideration many points, but is restricted to the problematic region isthree to four times faster than using the IAWPS formulation for pressure, even if tem-perature is known. For this the MATLAB curve fitting tool was used to fit a analyticfunction to the data obtained from [30] and [1]. The resulting equation is summarizedbelow.

p [MPa](ρ

[kg

m3

]= x, e

[kJ

kg

]= y)

=p00 + p10x+ p01y + p20x2 + p11xy + p02y

2 + p21x2y + p12xy

2 + p03y3

48

0

2

2

4

p [

Pa]

Pressure sampled as a function of and e

e [J/kg]

1

1000

[kg/m 3]

800600400200

Figure 35: Pressure varies heavily for nearly incompressible region

The degree of the multidimensional polynomial was chosen because the adjusted R-square value did not increase after this point, indicating that the model is not describingthe data better for what would be expected with an additional parameter. Coefficientsand statistics of the polynomial fit can be found in Table 2.

Parameter 00 10 01 20 11 02 21 12 03Value 6923 -15.79 -9.003 0.008863 0.01443 0.003607 -5.454e-6 -2.707e-6 -4.911e-7

R 1 Adj-R 1 RMSE 0.05491 RMSE(validation) 0.15473

Table 2: Parameters and statistics of surface fit in the liquid region

5. Test Case

We validated our solver against ColsimQT [9] from Fraunhofer ISE. The parametersof the pipe are summarized in Table 3. ColsimQT assumes constant pressure over the

49

pipe, ignores the momentum equation and uses the IAWPS84 equations for the waterproperties. The equations used are shown next.

Mass balance:∂ρ

∂t+∂(ρv)

∂x= 0

Energy balance:

ρ∂h

∂t+ ρv

∂h

∂x=Qint

Vint

Heat transfer model:

cW · ρW ·∂TW∂t

=Qext − Qint

VW.

Parameter Value

Pipe Length ` 1000 m

Pipe inner diameter 2rin 0.125 m

Pipe outer diameter 2rout 0.140 m

Heat capacity wall cW 540 J/(kg K)

Thermal conductivity wall λW 38 W/(m K)

Density pipe ρW 7500 kg/m3

Table 3: Parameters of the simulated pipe

For boundary conditions, we considered a prescribed mass flow with a known en-thalpy at the entry and a given pressure for the exit. The values for each of theboundary conditions can be found in Table 4. Since enthalpy is not a preliminaryvariable in our formulation for the flow simulation, the method to incorporate thisboundary condition is described next.

Parameter Value

Mass Flow at Entry m 3.5 kg/s

Pressure at Exit p 70.0 bar

Enthalpy at Entry h 944.96 kJ/kg

Heat from irridation qW 8000 W/m

Table 4: Values chosen for boundary values

50

5.1. Boundaries

We use ghost cells for the implementation of the boundary conditions. For the inflowwe are given the mass flow m and enthalpy h and we extrapolate pressure from thefirst cell. Then we find density ρ and inner energy e from

htarget!

= e+p(ρ, e)

ρ(81)

ptarget!

= p(ρ, e) . (82)

The velocity e can be recovered from m = ρuA where the cross-section area of the pipeis constant and known.

For the outflow we are given pressure and extrapolate mass flow and enthalpy fromthe last cell. Otherwise we proceed analogous to the inflow case.

5.2. Discussion Test Results

The test results of the simulation with 100 cells after 1600s are shown in Figure 36 indirect comparison with the output from ColSimQT for 1000 cells. While 100 cells isnot a fine space discretization, the CFL condition causes finer space discretization toreduce time step length leading to quadratic dependence of runtime on number of cells.Additionally the high speed of sound of water in the order of 1000 m/s also reducesthe upper bound of the time step. In this test increasing the number of time steps isespecially expensive since in each time step the equations 81 and 82 must be solvediteratively twice which is why we did not use as many cells in our simulation as in thereference ColSimQT.

The recovered results suggest that the constructed model is able to realisticallyreconstruct the behaviour of the flow inside the tube. We recover parts of the wellknown Rankine cycle in the temperature profile (see Figure 36a). It can be clearly seenthat water in the two-phase region is present in the middle of the pipe between thepoints where the temperature profile has discontinuous slope despite the inner energysteadily rising. Figure 36b shows good parity of the simulation results in the primaryvariable density. The fast reduction in density after entering the two-phase regionillustrates the phenomena that the vapour density enters the pseudo-fluids density isa linearly, but has orders of magnitude lower density.

The difference in the velocity plot Figure 36c can be explained by the friction thatis not considered in ColSimQT and mainly effects the velocity of the fluid. It alsohas to be considered that there exist small differences in the IAWPS95 and IAWPS84equations.

In summary we conclude that our model produces realistic results that are consistentwith ColSimQT for a test case where two-phase behaviour is exhibited.

51

(a) Temperature profile after 1600s (b) Density profile after 1600s

(c) Velocity profile after 1600s (d) Inner energy profile after 1600s

Figure 36: The simulation results show a clear similarity with results from ColSimQT

52

6. Conclusion

In the presented report we have:

• Implemented a deterministic 3D ray tracer for simulation of solar heat added toa pipe

• Implemented a fast explicit finite volume solver with an accurate water model

• Investigated multiple possibilities for a fast and accurate water model

• Validated the model against ColSimQT

6.1. Outlook

It may be fruitful to conduct further work building on the results of this project. Longtime step methods or implicit methods may be investigated for a faster thus closerto real time computation, since the currently small time steps must be used due tothe CFL condition. Combining this computation speed with the deterministic natureof the ray tracer, model based control or optimization might be a possible endeavor.For the water model, an accurate equation of state in density and inner energy, orbackwards equations for the IAWPS equations would simplify the use of the Eulerequations for water considerably. Lastly embedding the single pipe simulation in thewhole process cycle or simulating a network of pipes would pose an interesting field fora future project.

The ray tracing algorithm itself may be easily parallelized since in principle tracinga ray can be done independently for every mirror. However, using the object orientedmodeling described in section 1.4, we require every process to hold its on copy of theplant geometry data in form of objects of the class Mirror, Secondary and Absorber-Tube. Since the computational time for a single time step does not justify the overheadof a parallel initialization one may find it more suitable to apply parallelism to opti-mization tasks where simulations are done on a ”per-moment” basis. In this scenarioa nearly perfect speedup can be expected.

Since the optical model is deterministic in nature one may find it suitable for appli-cation in a broad range of nonlinear-optimization scenarios. Example problems include(but are not limited to) optimization of the mirror position on arbitrary ground to-pographies, optimization of mirror width and other parameters of the plant geometry.

53

A. Appendix

Figure 37: Relative error of NASG method with fitted parameters

Figure 38: Relative error of NASG method with global parameters

54

Figure 39: Relative error of look-up method

Figure 40: Relative error of IAWPS95 Method

55

Figure 41: Absolute error of NASG method with fitted parameters

Figure 42: Absolute error of NASG method with global parameters

56

Figure 43: Absolute error of look-up method with iterations

Figure 44: Absolute error of look-up method

57

References

[1] C. F. Beaton. Nbs/nrc steam tables, 2012. URL http://dx.doi.org/10.1615/

AtoZ.s.steam_tables.

[2] Boris Belhomme, Robert Pitz-Paal, Peter Schwarzbozl, and Steffen Ulmer. A newfast ray tracing tool for high-precision simulation of heliostat fields. Volume 131,08 2009.

[3] Ray Alden Berry, John William Peterson, Hongbin Zhang, Richard Charles Mar-tineau, Haihua Zhao, Ling Zou, and David Andrs. Relap-7 theory manual. Tech-nical report, Idaho National Laboratory (INL), 2015.

[4] Alexandre Chiapolino, Pierre Boivin, and Richard Saurel. A simple phase transi-tion relaxation solver for liquid vapor flows. International Journal for NumericalMethods in Fluids, 83(7):583–605, 2017. ISSN 1097-0363. doi: 10.1002/fld.4282.URL http://dx.doi.org/10.1002/fld.4282. fld.4282.

[5] Alexandre Chiapolino, Pierre Boivin, and Richard Saurel. A simple and fast phasetransition relaxation solver for compressible multicomponent two-phase flows.Computers & Fluids, 150:31 – 45, 2017. doi: 10.1016/j.compfluid.2017.03.022.URL https://hal.archives-ouvertes.fr/hal-01502389.

[6] R Cooper and B Dooley. Release on the iaps formulation 2008 for the viscosityof ordinary water substance. The International Association for the Properties ofWater and Steam, IAPWS, Berlin, Germany, 2008.

[7] B Dooley. Revised release on the iapws formulation 1995 for the thermodynamicproperties of ordinary water substance for general and scientific use. The Inter-national Association for the Properties of Water and Steam, IAPWS, Dresden,Germany, 2016.

[8] Jarno Elonen. Thin plate spline interpolation, 2005. URL https://elonen.iki.

fi/code/tpsdemo.

[9] T. Hirsch, M. Berger, and Gabriel Morin. Softwarevergleich zur warmetrager-regelung, projekt fresdemo technischer abschlussbericht zu ap 4.1. Technical re-port, Fraunhofer ISE, 2007.

[10] ML Huber, RA Perkins, DG Friend, JV Sengers, MJ Assael, IN Metaxa, K Miya-gawa, R Hellmann, and E Vogel. New international formulation for the thermalconductivity of h2o. Journal of Physical and Chemical Reference Data, 41(3):033102, 2012.

[11] Matthias Kunick, Hans-Joachim Kretzschmar, Francesca di Mare, and UweGampe. Fast calculation of steam and water properties with the spline-basedtable look-up method (sbtl). J. Eng. Gas Turbines & Power, in preperation.

58

[12] Richard Le Metayer, OlivierSaurel. The noble-abel stiffened-gas equation of state.Physics of Fluids, 28(4):046102, 2016. doi: 10.1063/1.4945981.

[13] Dong-Hui Li and Masao Fukushima. A derivative-free line search and global con-vergence of broyden-like method for nonlinear equations. Optimization Methodsand Software, 13(3):181–201, 2000. doi: 10.1080/10556780008805782.

[14] MA Moghimi, KJ Craig, and Josua P Meyer. A novel computational approachto combine the optical and thermal modelling of linear fresnel collectors using thefinite volume method. Solar Energy, 116:407–427, 2015.

[15] Jorge J More, B. S Garbow, and Kenneth E Hillstrom. User guide for minpack-1.Argonne National Laboratory, 1980.

[16] National Renewable Energy Laboratory (NREL). Reference data for the puerto er-rado 1 thermosolar power plant. URL https://www.nrel.gov/csp/solarpaces/

project_detail.cfm/projectID=46.

[17] Tiago Osorio, Pedro Horta, Marco Larcher, Ramon Pujol, Julian Hertel, De Wetvan Rooyen, Anna Heimsath, Simon Schneider, Daniel Benitez, Antoine Frein, andAlice Denarie. Ray-tracing software comparison for linear focusing solar collectors.1734:020017, 05 2016.

[18] John W. Peterson. Accurate curve fits of iapws data for high-pressure, high-temperature single-phase liquid water based on the stiffened gas equation of state.CoRR, abs/1311.0534, 2013. URL http://dblp.uni-trier.de/db/journals/

corr/corr1311.html#Peterson13.

[19] T Petrova and RB Dooley. Revised release on surface tension of ordinary watersubstance. The International Association for the Properties of Water and Steam,IAPWS, Moscow, Russia, 2014.

[20] F.J. Pino, R. Caro, F. Rosa, and J. Guerra. Experimental validation of an opticaland thermal model of a linear fresnel collector system. Applied Thermal Engineer-ing, 50(2):1463 – 1471, 2013. ISSN 1359-4311. doi: http://dx.doi.org/10.1016/j.applthermaleng.2011.12.020. URL http://www.sciencedirect.com/science/

article/pii/S1359431111007174. Combined Special Issues: ECP 2011 and IM-PRES 2010.

[21] Yu Qiu, Ya-Ling He, Ze-Dong Cheng, and Kun Wang. Study on optical andthermal performance of a linear fresnel solar reflector using molten salt as htfwith mcrt and fvm methods. Applied Energy, 146:162–173, 2015.

[22] Pascal Richter. Modelling and simulation of direct steam generation in the ab-sorber tubes of solar thermal power plants. Master’s thesis, 2011.

[23] Pascal Richter. Simulation and Optimization of Solar Thermal Power Plants.PhD thesis, PhD thesis, RWTH Aachen University, 2017.

59

[24] S.I. Sandler. Chemical, Biochemical, and Engineering Thermodynamics. NumberBd. 1 in Chemical, Biochemical, and Engineering Thermodynamics. John Wiley& Sons, 2006. ISBN 9780471661740. URL https://books.google.de/books?

id=4MXDAgAAQBAJ.

[25] Richard Saurel, Fabien Petitpas, and Remi Abgrall. Modelling phase transition inmetastable liquids: application to cavitating and flashing flows. Journal of FluidMechanics, 607:313–350, 2008. doi: 10.1017/S0022112008002061.

[26] Birte Schmidtmann and Manuel Torrilhon. A hybrid riemann solver for largehyperbolic systems of conservation laws. arXiv preprint arXiv:1607.05721, 2016.

[27] Levelt Sengers and B. Dooley. Revised supplementary release on saturation prop-erties of ordinary water substance. The International Association for the Proper-ties of Water and Steam, IAPWS, Russia, St. Petersburg, 1992.

[28] M.H.M. Sidek, W WAN HASAN, Mohd Zainal Abidin Ab Kadir, Suhaidi Shafie,Mohd Amran Mohd Radzi, Siti Ahmad, and Mohammad Hamiruce Marhaban.Gps based portable dual-axis solar tracking system using astronomical equation.pages 245–249, 03 2015.

[29] Eleuterio F Toro. Riemann solvers and numerical methods for fluid dynamics: apractical introduction. Springer Science & Business Media, 2013.

[30] Ohio Univesity. Properties of compressed liquid water. URL https://www.ohio.

edu/mechanical/thermo/property_tables/H2O/H2O_Compressed.html.

[31] Tim Wendelin, Aron Dobos, and Allan Lewandowski. Soltrace: A ray-tracingcode for complex solar optical systems. Contract, 303:275–3000, 2013.

60