Dynamic Simulation of an Organic Rankine Cycle—Detailed Model of a Kettle Boiler

: Organic Rankine Cycles (ORCs) are nowadays a valuable technology to produce electricity from low and medium temperature heat sources, e.g., in geothermal, biomass and waste heat recovery applications. Dynamic simulations can help improve the ﬂexibility and operation of such plants, and guarantee a better economic performance. In this work, a dynamic model for a multi-pass kettle evaporator of a geothermal ORC power plant has been developed and its dynamics have been validated against measured data. The model combines the ﬁnite volume approach on the tube side and a two-volume cavity on the shell side. To validate the dynamic model, a positive and a negative step function in heat source ﬂow rate is applied. The simulation model performed well in both cases. The liquid level appeared the most challenging quantity to simulate. A better agreement in temperature was achieved by increasing the volume ﬂow rate of the geothermal brine by 2% over the entire simulation. Measurement errors, discrepancies in working ﬂuid and thermal brine properties and uncertainties in heat transfer correlations can account for this. In the future, the entire geothermal power plant will be simulated, and suggestions to improve its dynamics and control by means of simulations will be provided.


Introduction
Organic Rankine Cycles (ORCs) have been increasingly applied to produce electricity from low and medium temperature sources in the last decade. The reported global installed capacity has currently surpassed 2700 MW e , and an additional 523. 6 MW e have been planned [1]. Typical application fields for ORC are geothermal, biomass, solar thermal and low grade waste heat recovery plants [2][3][4]. The size of ORC plants can vary from few kW e to some MW e [4,5]. Several studies report the advantages of organic fluids in comparison to water for low and medium temperature small-size power plants [5][6][7]. As an example, in [8,9] the organic fluid resulted to be better performing than water and air for a heat source temperature below 300 • C.
Given the modest temperature of the source, the electrical efficiency of ORCs is generally less than 20% [2,4]. The remaining energy is either dissipated through piping, compression and expansion machines, heat losses, auxiliary consumption and heat transfer to the condenser. The heat rejected to the condenser can be recovered to provide thermal energy for district heating, commercial buildings and industrial processes. In this case, the ORC is designed to work in combined heat and power (CHP) mode, providing electricity and heat with a relatively high energy utilization factor [10][11][12]. A comparison of the different configurations for ORC in CHP mode is provided in [13]. A simple ORC system consists of four main elements: a preheater/evaporator, where the working fluid is preheated and vaporized at saturated state or slightly superheated, an expansion machine, where the thermal energy of the fluid is converted into rotational energy of the expander shaft, a condenser, where the fluid is cooled down and condensed back to liquid state, and a pump, which brings the condensed liquid to the evaporation pressure to start the process cycle again. A recuperator might be applied, being advantageous mainly for dry fluids. Sub-, trans-and supercritical ORCs with pure organic fluids and their mixtures have also been investigated [14][15][16][17][18][19][20][21]. Supercritical ORCs might offer higher system efficiencies in dependency of the temperature source, at the expenses of higher investment and maintenance costs. The reduced temperature difference at the evaporator causes higher investment costs, because of the larger heat transfer surface and the higher pressure that the evaporator has to withstand [15,22]. For dry fluids, superheating has no positive effect and should be avoided, since it increases the complexity of the evaporator and the thermal stresses on the expansion machine [23].
The evaporator is the crucial component that links the ORC to the heat source/heat transfer medium, and its design is therefore of major relevance for the thermodynamic, hydraulic and economic performance of the ORC. For mobile applications, the heat exchanger has to be additionally optimized in terms of volume and weight [24][25][26].
In the following, the main configurations of evaporators for stationary ORC are considered. Fin-and-tube heat exchangers are mainly applied for waste heat recovery from gas turbine and internal combustion engines [27][28][29]. The extended surface of the shell side improves the heat transfer for the exhaust gas, whereas the working fluid shows a heat transfer coefficient of at least two orders of magnitude higher and flows inside the tubes [28]. The pressure drop of the exhaust gas on the shell side has a negative impact on the performance of the topping cycle and should be minimized [30][31][32]. Shell-and-tube heat exchangers are applied or assumed for ORC optimization in different studies [33][34][35][36]. This type of heat exchanger can be deployed up to very high pressures and temperatures, has a relatively simple geometry and its design and manufacturing procedures are very well-known and established [37]. Plate heat exchangers might be advantageous for low pressure, low temperature applications thanks to their compactness, effectiveness, easiness to be cleaned, maintained or extended [38]. As a major drawback, plate heat exchangers increase the pressure losses, especially for exhaust gas waste heat recovery, because of the narrow channels. A comparison between shell-and-tube and plate heat exchangers is carried out in [39].
In large size subcritical ORCs, kettle boilers are often used as evaporators [6]. The preheating of the organic fluid generally occurs in once-through shell-and-tube heat exchangers, whereas the evaporation takes place in the kettle boiler. The heating fluid, e.g., thermal brine for a geothermal power plant, passes through the tubes in a single-or multi-pass scheme. If no superheating is desired, the tube bundle is completely submerged in the liquid pool. The boiler disposes of a freeboard at the top of the kettle for separation of the vapor and liquid phase [37]. The vapor is then extracted from the top of the boiler, and fed to the expander. The effective separation between liquid and vapor in a kettle boiler gives more flexibility in power plant operation. In fact, in case of sudden variation of the input source, the boiler responds with a variation in pressure and/or level, avoiding that droplets get entrained to the turbine inlet [40]. In off-design conditions, some tubes may lay above the liquid level, leading to some degree of superheating and hence ensuring that no droplets are carried out to the expander inlet. In case of tube exposed to the vapor, tube material limits should not be overcome as a result of the higher wall temperature. The dynamic flexibility gives a significant advantage to kettle boilers with respect to tube once-through boilers. In the latter some degree of superheating is always necessary for safety reason, and a more demanding control is required to ensure pure vapor conditions at the expander inlet [40]. The pool boiler configuration is also recently being mentioned for mixtures of organic fluids [40]. An accurate evaporator design has also to consider the response of the heat exchanger in transient conditions. Dynamic simulations of ORC power plants are, in fact, gaining significant interest for a number of reasons: Energies 2017, 10, 548 3 of 28

•
To increase the energy conversion efficiency, power plants must be able to respond very fast and in an optimal way to variations in heat source or ambient conditions (e.g., waste heat recovery systems); • Control strategies can be improved, e.g., shifting from classical PID controllers [41,42] to advanced model-based or model-predictive strategies, leading to higher energy recovery [43][44][45]; • The presence of hot spots in heat exchangers during transients can be limited, avoiding thermal degradation of the working fluid [46]; • The increasing complexity in the dynamics of the electricity grid, associated with the increasing penetration of variable-source renewables (i.e., wind and photovoltaics), requires a higher contribution from controllable power plants (biomass, geothermal) for security margin, load reserve and grid stability; • Renewables-based ORCs could be used as stand-alone systems in remote areas and operated under strong dynamic conditions [47]; • The design of systems operated most of the time in off-design conditions (e.g., waste heat recovery) could be improved; stresses and plant lifetime can be increased [48,49]; • Start-ups, shutdowns and emergency situations (turbine trips, safety valve lifting) can also take help from dynamic simulations [50].
Dynamic models of ORC have been mainly developed in the Modelica ® language (Modelica Association, Linköping, Sweden) [51]. The equation-based, object-oriented and a-causal nature of the modelling language makes it very suitable for simulation of nonlinear thermo-fluid systems, such as ORC power plants [47,52]. The different plant components are developed as "individual objects" which can be connected to each other. The objects are generally compiled in libraries. Different libraries for thermo-fluid systems are available, either open-source or commercial [53]. Several software codes based on this modelling language are available, such as Dymola, SimulationX, JModelica or OpenModelica [54]. Thermodynamic and transport properties of the fluid are generally accessed by means of external fluid property libraries, e.g., NIST Reference Fluid and Transport Properties Database (REFPROP), FluidProp, CoolProp or TILMedia [55][56][57][58]. The focus of these dynamic tools is generally on the overall plant performance and the interaction among components rather than on the detailed description of the single elements, which is generally carried out with more demanding tools, such as Computational Fluid Dynamics (CFD) or Finite Element Methods (FEM) [47].
The broad application of kettle boilers for middle and large size ORCs, together with the increasing interest in ORC dynamics and simulation, require a valid model for the dynamic simulation of this component. The complex geometrical arrangement makes also current models unsuitable for this scope. In this paper, a dynamic model of a kettle evaporator is developed and validated with measured data. The model has the advantage of high simulation speed, good accuracy and a high flexibility for different heat exchanger geometries.
The dynamic modelling of heat exchangers and in particular evaporators for ORC systems is discussed in the next section. A closer insight into the proposed model for the kettle boiler is gained in Section 2. Validation and results are reported in Section 3, followed by a discussion in Section 4. A brief summary and conclusions are found at the end of this paper.

Dynamic Modelling of Evaporators for Organic Rankine Cycles
Dynamic models of heat exchangers (and not only) are based on the principle laws of mass, energy and momentum conservation, and integrated by empirical correlation to account for specific properties or characteristics of the system. The models are typically discretized over the length of the heat exchangers (1-D). The choice of the state variables for fluid property computation is discussed in [59,60]. The state variables have a strong impact on how the differential-algebraic equations (DAE) are solved. Most of the tools nowadays use pressure and specific enthalpy as state variables. The state of the solid wall of heat exchangers is defined by the wall temperature. The way the equations are The FV method is based on a static discretization of the heat exchanger in a given number of cells of equal volume. According to the degree of subcooling at the inlet and/or superheating at the outlet of the heat exchanger, the fluid in a cell can be in liquid, vapor or two-phase. On one hand, a higher number of cells is generally desired for higher accuracy of the solution. On the other hand, the computational time increases with the number of cells. A trade-off has to be found between the accuracy of the solution and the computational time. A major problem occurs with the FV method applied to a phase-change heat exchanger. This can be explained looking at the continuity equation: (1) where is the fluid density, the cell volume, the incoming and the outgoing fluid mass flow rate. Since the saturated liquid line shows a large discontinuity in density when the fluid passes from the liquid to the two-phase region, the time derivative of the density can become very high.
If the inlet mass flow rate is given, a fast density change can cause peaks in , resulting in flow reversal, chattering or non-solvable solutions. The simulation can become extremely slow or even fail [61]. MB refers to fast low-order dynamic models that subdivide the heat exchanger in a liquid, vapor and two-phase control volume. If one of the regions is not present, the corresponding volume is neglected. The volumes dynamically change their length according to the thermodynamic conditions in the heat exchangers. Techniques to account dynamically for the presence of each control volume are also available [62]. In the two-phase region, the average density ̅ is computed in MB methods by means of the average void fraction, according to [63]: The FV method is based on a static discretization of the heat exchanger in a given number of cells of equal volume. According to the degree of subcooling at the inlet and/or superheating at the outlet of the heat exchanger, the fluid in a cell can be in liquid, vapor or two-phase. On one hand, a higher number of cells is generally desired for higher accuracy of the solution. On the other hand, the computational time increases with the number of cells. A trade-off has to be found between the accuracy of the solution and the computational time. A major problem occurs with the FV method applied to a phase-change heat exchanger. This can be explained looking at the continuity equation: where ρ is the fluid density, V the cell volume, . m in the incoming and . m out the outgoing fluid mass flow rate. Since the saturated liquid line shows a large discontinuity in density when the fluid passes from the liquid to the two-phase region, the time derivative of the density dρ dt can become very high. If the inlet mass flow rate . m in is given, a fast density change can cause peaks in . m out , resulting in flow reversal, chattering or non-solvable solutions. The simulation can become extremely slow or even fail [61].
MB refers to fast low-order dynamic models that subdivide the heat exchanger in a liquid, vapor and two-phase control volume. If one of the regions is not present, the corresponding volume is neglected. The volumes dynamically change their length according to the thermodynamic conditions in the heat exchangers. Techniques to account dynamically for the presence of each control volume are also available [62]. In the two-phase region, the average density ρ is computed in MB methods by means of the average void fraction, according to [63]: where γ is the average void fraction, and ρ sat,l and ρ sat,v are the densities of saturated liquid and vapor. Assumptions on the average void fraction have to be made. The validity of the assumption has an impact on the accuracy of the dynamic simulation [63]. The MB and FV models have been compared in [59,63,64]. The MB model requires less computation time than FV, even though the accuracy is typically lower. The MB model could preferably be applied for online control systems, rather than system simulation [59,60]. The computation of the amount of working fluid present in the cycle (also called charge) is critical for both models because of the absence of high-accuracy correlations for the void fraction in the evaporator and condenser. The MB model has however shown worse results between the two [64].
Another modeling approach can be used when a fluid changes phase in a large shell or pool, and the heating/cooling medium flows in a tube bundle. The evaporating/condensing fluid can be represented by means of two volumes (TV) in thermal non-equilibrium. The TV are modelled with a lumped approach, as in the MB model. Unlike the MB model, however, the mass transfer between the TV is a function of the vapor quality in each of the TV, and not a result of the mass and momentum equation at the cell boundaries. As an example, in case of the evaporator, when the vapor quality in the liquid phase is bigger than zero, the liquid is vaporizing, so mass is transferred from the liquid to the vapor region. These mass transfers are independent from the outer direction of flow, which is pressure driven. A TV model for a shell-and-tube condenser is included in the ThermoSysPro library in Modelica ® and discussed in [65].
An accurate model of the evaporator is essential to reproduce the dynamic behavior of ORC. As an example, most of the control strategies are based on the degree of superheating at the evaporator outlet, or on the evaporation pressure [41]. The thermal inertia and time response of the evaporator have a crucial impact on how source fluctuations affect the dynamics and control of the power plant. While for small-scale ORC once-through heat exchangers of simple geometry are applied, kettle boilers are very often used in larger scale ORC. Because of the more complex geometry of the kettle, FV and MB cannot represent properly the dynamics of this heat exchanger. In the present work, a dynamic model of a kettle boiler is developed by means of a combination between TV (for the shell-side) and FV method (for the tube/wall side) and validated against measurements from an analogous heat exchanger located in a geothermal power plant in the Munich area. The present work is based on the commercial TIL-Library (Version 3.2.2, TLK-Thermo GmbH, Braunschweig, Germany [66]). The basic models have been taken from this library, and new components have been developed based on the existing library. The resulting combined TV/FV model can be used for design, off-design and dynamic simulations, and applied for testing and development of basic or advanced control strategies for middle and large-scale ORC power plants.

TV/FV Model for the Kettle Boiler
In this section, the theory and equations necessary for the modeling of the kettle boiler are described. It is very important to limit the complexity of the model, in order to allow for its integration in a closed-loop cycle simulation, without incurring in extremely large simulation time or simulation failures.

Two-Volume Cavity
The shell-side of the kettle boiler is modelled by means of a TV cavity (TVC), shown in Figure 2. The lower volume (Region 1, grey) represents the vaporizing liquid, whereas the upper volume (Region 2, white) describes the vapor region. The laws of conservation for mass, momentum and energy are developed for each (lumped) volume, which are in general in thermal non-equilibrium. The TVC can exchange mass according to the vapor quality of the volumes. If the vapor quality in Energies 2017, 10, 548 6 of 28 the Region 1 is larger than zero, a mass flow will be transferred from Region 1 to 2, since vapor is produced. If the vapor quality in the Region 2 drops below one, condensation occurs and mass is transferred to Region 1. The non-equilibrium among the TVC is particularly suitable if, for example, the vapor produced by the evaporator is superheated, or if the liquid at the inlet is subcooled. The two regions can transfer heat through their interface, whose surface changes dynamically according to the liquid level. The liquid level depends on the region occupied by the TV. The sum of the two volumes is equal to the shell-side volume of the heat exchanger and remains constant.

Two-Volume Cavity
The shell-side of the kettle boiler is modelled by means of a TV cavity (TVC), shown in Figure 2. The lower volume (Region 1, grey) represents the vaporizing liquid, whereas the upper volume (Region 2, white) describes the vapor region. The laws of conservation for mass, momentum and energy are developed for each (lumped) volume, which are in general in thermal non-equilibrium. The TVC can exchange mass according to the vapor quality of the volumes. If the vapor quality in the Region 1 is larger than zero, a mass flow will be transferred from Region 1 to 2, since vapor is produced. If the vapor quality in the Region 2 drops below one, condensation occurs and mass is transferred to Region 1. The non-equilibrium among the TVC is particularly suitable if, for example, the vapor produced by the evaporator is superheated, or if the liquid at the inlet is subcooled. The two regions can transfer heat through their interface, whose surface changes dynamically according to the liquid level. The liquid level depends on the region occupied by the TV. The sum of the two volumes is equal to the shell-side volume of the heat exchanger and remains constant.  Four differential equations are used to describe the TVC, and the corresponding four differential states are assumed to be the specific enthalpies in the liquid and vapor, the pressure at the interface and the volume of liquid. The input port is located at the bottom (below the Region 1), and the outlet port at the top (above Region 2), so that only vapor is sent to the expander, as in the real heat exchanger. In this case, the mass balances in each of the TV are: where m is the volume mass, the subscript 'l' refers to the liquid and 'v' to the vapour volume, 'in' to the input port at the liquid volume, and 'out' to the output port connected to the vapor volume. The densities of the liquid and vapor volume (resp. ρ l and ρ v ) are computed taking into account the vapor quality of the volumes, and are in general different from their values at saturation. Since only the specific enthalpy h and the pressure p are the state variables, the time derivative of the density is extended as: The evaporation and condensation mass flow rates are defined as: .
. where x is the steam quality and x re f is a reference value. The factors C cond and C evap are two tuning parameters that have to be chosen in order to achieve a proper dynamics with respect to the experimental data. They depend on the fluid, and on the evaporator geometry. A higher value leads to a faster response of the system, and lower value to a slower.
x v,re f and x l,re f should be theoretically 1 and 0, but they might differ slightly for a better numerical performance. The momentum balance can be considered as steady; liquid head and pressure drops at the nozzles (∆p drop,in and ∆p drop,out ) can be included as well: where g is the acceleration of gravity and l l the liquid level, computed from the bottom point of the evaporator. The energy balance for the TV is: where U l and U v are the internal energies, h sat, v and h sat, l are the saturation specific enthalpies of the condensing and evaporating flows at the pressures p l and p v , and ∑ . Q l and ∑ .
Q v account for the heat transferred by respectively the liquid and vapor volume with the tube bundle, with the wall of the heat exchanger and at the interface with the other volume itself. Combining Equations (12) and (13), with the definition of internal energy U = H − pV and including Equations (3)-(6), the energy balances become:

Heat Transfer Network and Correlations
The heating fluid flows on the tube-side and vaporizes the organic fluid on the shell side. In this work, the tube-side pathway is discretized in a number of cells along its main flow direction, according to the FV method. Parallel tubes of the same pass are lumped in the same cells. Each cell is connected to a wall cell representing the solid material of the pipe, and transfers heat either to Region 1 or 2 (or a fraction of both) of the TVC, according to the liquid level and the position of the cell. The heat transfer scheme is shown in Figure 3. In general, the heat transferred by each cell i can be written as: where α cell,i is the cell heat transfer coefficient, A i the cell heat transfer area, T cell,i the cell temperature and T w,cell,i the temperature of the tube wall at the cell side. The heat transfer between the wall cell and its interface with the cell i is: . where T w,i is the average temperature of the wall and R w,i the thermal resistance of the wall cell. For the solid wall, the following dynamic equation is valid in each cell i: where c w is the specific heat capacity (for stainless steel assumed at 500 J/kgK), ρ w the density (for stainless steel assumed as 7800 kg/m 3 ) and V w,i the volume of the wall cell.
. Q l,wall,i and .
Q v,wall,i refer to the heat transfer rate for the wall cell with the liquid and vapor volume. If the cell is completely submerged . Q v,wall,i = 0, if it is completely exposed to the vapor . Q l,wall,i = 0. The heat transfer between the tube and the interface with Region 1 and/or 2 is: where T w,lv,i is the temperature of the wall at the interface with Region 1 and/or 2. The heat transferred by Region 1 with each single wall cell i is described as: .
where α l,i is the heat transfer coefficient of the evaporating liquid and A l,i is the external heat transfer area of the cell in contact with the liquid. The heat transferred by Region 2 with each wall cell is analogously: where α v,i is the heat transfer coefficient of the vapor and A v,i is the external heat transfer area of the cell in contact with the vapor. The heating fluid considered in this paper is geothermal water.
To compute the heat transfer coefficient α cell,i , the equation proposed in McEagle [37] is used: where u i is the liquid velocity inside the pipe and d i the internal diameter of the pipe. The McEagle correlation has been developed specifically for water and suggested by [37]. In absence of further information on the characteristics of the geothermal water, this correlation will be used.
where is the liquid velocity inside the pipe and the internal diameter of the pipe. The McEagle correlation has been developed specifically for water and suggested by [37]. In absence of further information on the characteristics of the geothermal water, this correlation will be used. In a kettle boiler, the evaporation occurs in a quiescent pool of saturated liquid. The convective effects are in this situation almost negligible. The evaporation regime is called "nucleate boiling", since the process takes place with the formation of bubbles, which detach when close to the vapor-liquid interface. The bubbles are generally formed in vapor filled cavities of rough tube surfaces [67]. A condition of stable nucleate boiling is given in [37], where a maximum heat flux at the tube bundle is  In a kettle boiler, the evaporation occurs in a quiescent pool of saturated liquid. The convective effects are in this situation almost negligible. The evaporation regime is called "nucleate boiling", since the process takes place with the formation of bubbles, which detach when close to the vapor-liquid interface. The bubbles are generally formed in vapor filled cavities of rough tube surfaces [67]. A condition of stable nucleate boiling is given in [37], where a maximum heat flux at the tube bundle is defined. Since the heat transfer coefficient is dependent on the nature and condition of the heat exchange surface, a correlation with high accuracy for every system cannot be provided [37]. The heat transfer coefficient for pool or nucleate boiling can be described by means of several correlations. The Mostinski equation [68], which is a function of the cell heat flux . q l.i = . Q l,wall,i /A l,i and the pressure p l of the liquid, has a relative simple form: where p c is the critical pressure of the fluid. The pressures have to be expressed in bar. Another proposed correlation, developed from statistical multiple regression techniques for organic fluids is the Stephan and Abdelsalam [69]: where λ sat,l is the thermal conductivity and a sat,l is the thermal diffusivity of the saturated liquid, T sat is the saturation temperature and ∆h Lv is the latent heat of vaporization at the pressure p l . The bubble departure diameter d bub is defined as: (25) where β is equal to 35 • for every fluid and σ is the surface tension.
The Cooper equation is also function of the heat flux . q l.i and liquid pressure p l [69]: where ε is the pipe roughness in µm (for stainless steel ε = 1.5 µm) and M is the molecular weight of the organic fluid in [g/mol]. The discussed correlations (23)- (26) for nucleate boiling are compared in Figure 4 for a common organic working fluid, in the range of pressures 5 < p l < 20 bar and heat flux equal to . q l.i = 13 kW m 2 . Stainless steel is considered as pipe material. It can be seen that the Mostinski equation appears more conservative than the other two. The Cooper and Stephan equations show a very similar pattern, but the Stephan is slightly more conservative and has a more complex structure. Because of its simplicity, only the Cooper correlation is used in the following.
The film-boiling equation can be used for the heat transfer at cells that are not submerged under the liquid level and exchange heat with the vapor that surrounds them. The heat transfer is limited in this case by the conduction through the film of vapor and can be described by the Bromley equation [70]: where λ v is the thermal conductivity of the vapor, c pv is the specific heat at constant pressure of the vapor, d o the external diameter of the pipe and η v the dynamic viscosity of the vapor. The wall cell thermal resistance R w,i is defined as: where d i is the tube internal diameter, λ w is the tube thermal conductivity (for stainless steel assumed at 15 W/mK), L i is the length of the wall cell and N p i is the number of parallel tubes per pass. where is the pipe roughness in μm (for stainless steel = 1.5 μm) and is the molecular weight of the organic fluid in [g/mol].
The discussed correlations (23)- (26) for nucleate boiling are compared in Figure 4 for a common organic working fluid, in the range of pressures 5 < < 20 and heat flux equal to . = 13 .
Stainless steel is considered as pipe material. It can be seen that the Mostinski equation appears more conservative than the other two. The Cooper and Stephan equations show a very similar pattern, but the Stephan is slightly more conservative and has a more complex structure. Because of its simplicity, only the Cooper correlation is used in the following.
where is the thermal conductivity of the vapor, is the specific heat at constant pressure of the vapor, the external diameter of the pipe and the dynamic viscosity of the vapor. The wall cell thermal resistance , is defined as: where is the tube internal diameter, is the tube thermal conductivity (for stainless steel assumed at 15 W/mK), is the length of the wall cell and is the number of parallel tubes per pass.

Geometry of the Kettle Boiler
The kettle boiler used for the model validation is, according to the standards of the Tubular Exchanger Manufacturers Association, Inc. (TEMA) [71], of the type NKN with no baffles. The tube bundle extends horizontally for the entire heat exchanger and is welded to the shell in the front and rear side to form a box (for this reason, it is also called "box-type" [72]). As mentioned, the boiler has a freeboard used to vaporize the fluid on the shell side in subcritical conditions, and separate the liquid droplets that might be carried by the outgoing vapor stream itself. That is why such evaporator has a kettle diameter significantly bigger than the bundle one (usually one third larger than the inner bundle diameter). The bundle diameter is also called "shell diameter". To minimize the amount of liquid extracted by the vapor flow, an impingement-type baffle is usually also present at the vapor outlet nozzle.

Geometry of the Kettle Boiler
The kettle boiler used for the model validation is, according to the standards of the Tubular Exchanger Manufacturers Association, Inc. (TEMA) [71], of the type NKN with no baffles. The tube bundle extends horizontally for the entire heat exchanger and is welded to the shell in the front and rear side to form a box (for this reason, it is also called "box-type" [72]). As mentioned, the boiler has a freeboard used to vaporize the fluid on the shell side in subcritical conditions, and separate the liquid droplets that might be carried by the outgoing vapor stream itself. That is why such evaporator has a kettle diameter significantly bigger than the bundle one (usually one third larger than the inner bundle diameter). The bundle diameter is also called "shell diameter". To minimize the amount of liquid extracted by the vapor flow, an impingement-type baffle is usually also present at the vapor outlet nozzle.
The geometry of the heat exchanger can be seen in Figure 5. The heating medium (for instance, geothermal water) flows inside the tubes, in a single or multi-pass arrangement. The tube path is discretized in cells to represent the model discretization discussed in Section 2.2. Parallel tubes of the same pass are lumped in the same cells. Figure 5 also shows the liquid and vapor volumes of the working fluid (resp. green and white regions). The flow directions are highlighted: the heating fluid (geothermal brine) enters and leaves the evaporator from the left, whereas the organic fluid enters the evaporator as liquid from the bottom and leaves it as vapor from the top.
Region 1 and 2 occupy a larger or smaller volume according to the liquid level, and according to their volume, the tubes transfer more or less heat with each of the regions. To correlate the occupied volume by Region 1 and 2 with the liquid level, the geometry of the kettle boiler has to be defined. The cross-section of kettle boiler at the medium kettle length is shown in Figure 6. Different geometric parameters are defined in the following. The kettle is characterized by a kettle radius R s , corresponding to the largest internal diameter of the evaporator. The kettle has an internal shell (defined by the dash line in Figure 6a), characterized by a radius R. The tube bundle is found within the shell. The bundle consists of a given number of tubes N tubes , having a certain outer diameter d o , a wall thickness t w and length L. As an example, the cross-section of a single tube is shown in grey in Figure 6a. The shell is redrawn in Figure 6b together with the tube bundle (in grey) and subdivided in four segments, representing four tube passes. The cross-sectional area occupied by the tube bundle over the shell is (grey area of Figure 6b): Energies 2017, 10, 548 10 of 27 The geometry of the heat exchanger can be seen in Figure 5. The heating medium (for instance, geothermal water) flows inside the tubes, in a single or multi-pass arrangement. The tube path is discretized in cells to represent the model discretization discussed in Section 2.2. Parallel tubes of the same pass are lumped in the same cells. Figure 5 also shows the liquid and vapor volumes of the working fluid (resp. green and white regions). The flow directions are highlighted: the heating fluid (geothermal brine) enters and leaves the evaporator from the left, whereas the organic fluid enters the evaporator as liquid from the bottom and leaves it as vapor from the top. Region 1 and 2 occupy a larger or smaller volume according to the liquid level, and according to their volume, the tubes transfer more or less heat with each of the regions. To correlate the occupied volume by Region 1 and 2 with the liquid level, the geometry of the kettle boiler has to be defined. The cross-section of kettle boiler at the medium kettle length is shown in Figure 6. Different geometric parameters are defined in the following. The kettle is characterized by a kettle radius , corresponding to the largest internal diameter of the evaporator. The kettle has an internal shell (defined by the dash line in Figure 6a), characterized by a radius . The tube bundle is found within the shell. The bundle consists of a given number of tubes , having a certain outer diameter , a wall thickness and length . As an example, the cross-section of a single tube is shown in grey in Figure 6a. The shell is redrawn in Figure 6b together with the tube bundle (in grey) and subdivided in four segments, representing four tube passes. The cross-sectional area occupied by the tube bundle over the shell is (grey area of Figure 6b    Region 1 and 2 occupy a larger or smaller volume according to the liquid level, and according to their volume, the tubes transfer more or less heat with each of the regions. To correlate the occupied volume by Region 1 and 2 with the liquid level, the geometry of the kettle boiler has to be defined. The cross-section of kettle boiler at the medium kettle length is shown in Figure 6. Different geometric parameters are defined in the following. The kettle is characterized by a kettle radius , corresponding to the largest internal diameter of the evaporator. The kettle has an internal shell (defined by the dash line in Figure 6a), characterized by a radius . The tube bundle is found within the shell. The bundle consists of a given number of tubes , having a certain outer diameter , a wall thickness and length . As an example, the cross-section of a single tube is shown in grey in Figure 6a. The shell is redrawn in Figure 6b together with the tube bundle (in grey) and subdivided in four segments, representing four tube passes. The cross-sectional area occupied by the tube bundle over the shell is (grey area of Figure 6b): (a) (b)   The remaining free cross-sectional area of the shell is (white area of Figure 6b): An important quantity is the ratio of the free cross-sectional area of the shell (white area of Figure 6b) to total cross-sectional area of the shell (white plus grey area of Figure 6b): The larger σ, the larger the volume available for the working fluid on the shell side. In a first approximation, it can be assumed that the local σ, i.e., the area ratio considering only a portion of the shell, is uniform and equal to the shell σ of Equation (31). As mentioned, the tube-side flow can be arranged in multiple passes N passes . The inner diameter of each tube is d i = d o − 2t w . The heat transfer area for each tube cell from the inner side (used in Equation (16)) is: For each pass, the cross-sectional area of the tube bundle is: First pass : Subsequent passes 2 ≤ k ≤ N passes : where dR = 2 R N passes . All the tube cells that are part of the same pass k have the same cross-sectional area A sk . The process to determine the contact areas A l,i and A v,i for each cell i with the TV as a function of the liquid level l l is shown in Appendix A.
The equations that relate the liquid volume V l in the evaporator and the liquid level l l are defined in Appendix B and the results from Equation (A23) are plotted in Figure 7 for two kettle evaporators. The evaporators are located in an ORC geothermal power plant in Sauerlach, Germany (see Section 3.1). For confidentiality issues, the quantities have been normalized to the manufacturer data. The discrepancy between the proposed volume method and the values given by the manufacturer is shown in Table 1. The errors are 2.29% and −0.37% for the two evaporators.
An important quantity is the ratio of the free cross-sectional area of the shell (white area of Figure 6b) to total cross-sectional area of the shell (white plus grey area of Figure 6b): The larger , the larger the volume available for the working fluid on the shell side. In a first approximation, it can be assumed that the local , i.e., the area ratio considering only a portion of the shell, is uniform and equal to the shell of Equation (31). As mentioned, the tube-side flow can be arranged in multiple passes . The inner diameter of each tube is = − 2 . The heat transfer area for each tube cell from the inner side (used in Equation (16)) is: For each pass, the cross-sectional area of the tube bundle is: 2 ≤ ≤ : All the tube cells that are part of the same pass have the same cross-sectional area . The process to determine the contact areas , and , for each cell with the TV as a function of the liquid level is shown in Appendix A. The equations that relate the liquid volume in the evaporator and the liquid level are defined in Appendix B and the results from Equation (A23) are plotted in Figure 7 for two kettle evaporators. The evaporators are located in an ORC geothermal power plant in Sauerlach, Germany (see Section 3.1). For confidentiality issues, the quantities have been normalized to the manufacturer data. The discrepancy between the proposed volume method and the values given by the manufacturer is shown in Table 1. The errors are 2.29% and −0.37% for the two evaporators.

Validation and Results
In this section, the model of the evaporator is validated and compared with available measurements from the geothermal power plant in Sauerlach. Section 3.1 describes the power plant configuration and which measurements are available. The simulation set-up is compared in Section 3.2. In Section 3.3 the case studies are described and the results are shown in Section 3.4.

Geothermal CHP Plant in Sauerlach
The geothermal power plant in Sauerlach is one of the largest geothermal power plants in Germany, with an electrical output of approximately 5 MW e and a thermal output of around 4 MW th . The plant has 10.99% of net electrical efficiency. At the design point, the geothermal brine has a temperature of 140 • C, a nominal volume flow rate of 110 L/s. The plant configuration is depicted in Figure 8.

Validation and Results
In this section, the model of the evaporator is validated and compared with available measurements from the geothermal power plant in Sauerlach. Section 3.1 describes the power plant configuration and which measurements are available. The simulation set-up is compared in Section 3.2. In Section 3.3 the case studies are described and the results are shown in Section 3.4.

Geothermal CHP Plant in Sauerlach
The geothermal power plant in Sauerlach is one of the largest geothermal power plants in Germany, with an electrical output of approximately 5 MWe and a thermal output of around 4 MWth. The plant has 10.99% of net electrical efficiency. At the design point, the geothermal brine has a temperature of 140 °C, a nominal volume flow rate of 110 L/s. The plant configuration is depicted in Figure 8. The plant consists of a dual-pressure ORC that makes use of the split concept [73]. The high temperature (HT, or high-pressure) cycle has one evaporator and two preheaters. The low temperature (LT, or low-pressure) loop has one evaporator and one preheater. The thermal water is pumped from the extraction probe (1) to the high temperature evaporator (HTE), and then flows to the high temperature preheater of the high temperature cycle (2  3). At this point the still hot thermal brine is fed to the low temperature evaporator (3  4). Then the flow is split into lines (5a) and (5b) that go respectively to the low temperature preheater of the high temperature cycle and the only preheater of the low temperature loop. After being cooled down to around 45 °C in (6) and (7) the thermal water is fed back to the injection probes at point 8. The HT cycle and the LT cycle are completely separated in terms of organic fluid process. The high and low temperature condensers The plant consists of a dual-pressure ORC that makes use of the split concept [73]. The high temperature (HT, or high-pressure) cycle has one evaporator and two preheaters. The low temperature (LT, or low-pressure) loop has one evaporator and one preheater. The thermal water is pumped from the extraction probe (1) to the high temperature evaporator (HTE), and then flows to the high temperature preheater of the high temperature cycle (2 → 3). At this point the still hot thermal brine is fed to the low temperature evaporator (3 → 4). Then the flow is split into lines (5a) and (5b) that go respectively to the low temperature preheater of the high temperature cycle and the only preheater of the low temperature loop. After being cooled down to around 45 • C in (6) and (7) the thermal water is fed back to the injection probes at point 8. The HT cycle and the LT cycle are completely separated in terms of organic fluid process. The high and low temperature condensers are air-cooled, by means of induced-draught fans. The condensate from the air cooled condensers is then collected in the condensate tanks, from which the pump drives the liquid back to the preheaters (14 → 9 and 19 → 15). The turbines are connected via a gearbox to a single synchronous generator, since they have different rotational speed (around 2683 and 1457 rpm for the HT and LT turbines respectively). They are single-stage axial turbines. In the start-up phase, a bypass valve controls the flow passing through each turbine. The plant is in fact a CHP plant. Different bypass valves are available through the thermal water pathway, to allow a proper control of the heat and electrical load. In particular, an extraction for the district heating system is placed at point 3, but also a bypass directly at point 1 can be used. In this way a combination of parallel and series (at least in a small portion) configuration of the power and the heat production systems is achieved, which guarantees more flexibility.
The present work focuses on the HTE (in red box in Figure 8). The measurement available for the thermal brine are the volume flow rate and the pressure at point 1 and the temperatures at points 1 and 2. For the organic fluid, the temperatures at points 11 and 12 and the pressure in the evaporator (ca. 12). The volume flow rate of organic fluid is measured at the outlet of the pump (point 9). By means of the temperature at point 14 and pressure at point 9, the mass flow rate is computed. Since between points 9-11 the organic fluid is at the liquid phase, the mass flow rate at the HTE inlet is assumed equal to the one at the pump outlet. A summary of the available measurements is provided in Table 2.

Simulation Set-Up
The model of the kettle boiler is simulated in Dymola. The standard solver DASSL is used, with 10 −4 tolerance. The model has been developed combining present models of the TIL library (heating fluid and wall cell) and creating a new TVC model (see Section 2). The fluid properties are computed by means of REFPROP. The model has been assembled in a single component called "KettleBoiler".
The main geometric features of the boiler, the initial states of both fluids and some other settings that can be modified are summarized in Table 3. The reference vapor qualities for evaporation and condensation are set at resp. 7.5 × 10 −3 and 0.9975 for numeric reasons. Both mass transfer coefficients C evap and C cond are set at 1 s −1 , as the impact of these parameters is limited and the set values could guarantee a best match for the simulation. The pressure drops at the nozzles ∆p drop,in and ∆p drop,out in Equations (9)-(11) are neglected. Since the pressure difference between the top and the bottom of the liquid level is lower than 0.5% the absolute pressure, it is assumed that p l = p v . The simulation set-up is shown in Figure 9. The blue line refers to the thermal brine and the green to the organic fluid. The inlet thermal flow rate and temperature of the thermal brine are set as measured. For simplicity, its pressure is fixed at 10 bar. For the organic fluid the inlet temperature is set as measured. The mass flow rate at the inlet is set as the calculated mass flow rate for the pump outlet. The evaporator is connected to a turbine, modelled through the common Stodola's equation: where π = p v p out,t is the turbine pressure ratio and k = 0.0107 m −2 is the turbine characteristics. The outlet pressure of the turbine is also set as measured. The pressure and liquid level in the evaporator, together with the HTE outlet temperatures for both fluids are therefore a result of the simulation. x l, re f -Reference vapor quality for condensation is the turbine pressure ratio and = 0.0107 m is the turbine characteristics.
The outlet pressure of the turbine is also set as measured. The pressure and liquid level in the evaporator, together with the HTE outlet temperatures for both fluids are therefore a result of the simulation.  Figure 9. Kettle boiler: simulation layout.

Case Studies
Two case studies are evaluated for the validation of the evaporator model, based on measurement data from Sauerlach. Case 1 refers to measurements taken on 23 May 2015, whereas Case 2 refers to measurements taken on 14 May 2015. A step function in inlet volume flow rate of thermal brine is

Case Studies
Two case studies are evaluated for the validation of the evaporator model, based on measurement data from Sauerlach. Case 1 refers to measurements taken on 23 May 2015, whereas Case 2 refers to measurements taken on 14 May 2015. A step function in inlet volume flow rate of thermal brine is applied in both cases: in Case 1, the step is negative and has an amplitude of 2.98%, whereas in Case 2 the volume flow rate is increased by 3.26%. The two measured profiles are shown in Figure 10. The values have been normalized for confidentiality issues. It can be seen that the inlet temperature of the geothermal brine does not vary significantly, as it commonly is in deep geothermal power plants.
drop in the measured pressure of approximately 1.5%, as shown in Figure 11a ("Measured" line). The same drop is represented in the simulation, which has a starting offset of 0.9% ("Base" line) that is however kept over the entire simulation, so that the error remains within 1% for the entire simulation. The outlet temperature of the working fluid (WF) and of the thermal brine (HF) are lower than measured, even though within 1% ("Base" line in Figure 11b). In Case 2, the change in pressure is opposite to Case 1, and increases of 1.7% ("Measured" line in Figure 12a).

Simulation Results
The simulation is run separately for the two cases since the boundary conditions are different. All the values except the time scale have been normalized for confidentiality issues. The boiler operates in regime of stable nucleate boiling (mentioned in Section 2.2), given the limited temperature difference and the heat flux at around 50% of the critical heat flux.
In Case 1, the drop in thermal brine volume flow rate at t = 1400 s (see Figure 10a) causes a drop in the measured pressure of approximately 1.5%, as shown in Figure 11a ("Measured" line). The same drop is represented in the simulation, which has a starting offset of 0.9% ("Base" line) that is however kept over the entire simulation, so that the error remains within 1% for the entire simulation. The outlet temperature of the working fluid (WF) and of the thermal brine (HF) are lower than measured, even though within 1% ("Base" line in Figure 11b). In Case 2, the change in pressure is opposite to Case 1, and increases of 1.7% ("Measured" line in Figure 12a). applied in both cases: in Case 1, the step is negative and has an amplitude of 2.98%, whereas in Case 2 the volume flow rate is increased by 3.26%. The two measured profiles are shown in Figure 10. The values have been normalized for confidentiality issues. It can be seen that the inlet temperature of the geothermal brine does not vary significantly, as it commonly is in deep geothermal power plants.

Simulation Results
The simulation is run separately for the two cases since the boundary conditions are different. All the values except the time scale have been normalized for confidentiality issues. The boiler operates in regime of stable nucleate boiling (mentioned in Section 2.2), given the limited temperature difference and the heat flux at around 50% of the critical heat flux.
In Case 1, the drop in thermal brine volume flow rate at = 1400 s (see Figure 10a) causes a drop in the measured pressure of approximately 1.5%, as shown in Figure 11a ("Measured" line). The same drop is represented in the simulation, which has a starting offset of 0.9% ("Base" line) that is however kept over the entire simulation, so that the error remains within 1% for the entire simulation. The outlet temperature of the working fluid (WF) and of the thermal brine (HF) are lower than measured, even though within 1% ("Base" line in Figure 11b). In Case 2, the change in pressure is opposite to Case 1, and increases of 1.7% ("Measured" line in Figure 12a).  The simulation follows the same pattern, even though an offset lower than 1% is observable ("Base" line). The outlet temperature of both fluids follow the measured value within a 1% difference (cf. "Measured" and "Base" line in Figure 12b). A sensitivity analysis on the effect of the increase in volume flow rate of thermal water is carried out and shown in Figure 13. To account for differences between measurement and simulation over the entire simulation time, the relative root mean square error RRMSE is considered: where is the number of time intervals, is the measured value and is the respective simulated value. Figure 13 shows that an increase of 2% of volume flow rate of thermal brine leads to the minimum RRMSE in temperature. In fact, by increasing the thermal brine volume flow rate of 2% over the entire simulation, the RRMSE for the outlet temperatures goes below 0.2% in both Cases 1 and 2. The time behavior in pressure and temperature with the 2% increase in thermal brine volume flow rate can be observed in the "+2%" line in respectively Figures 11 and 12.  The simulation follows the same pattern, even though an offset lower than 1% is observable ("Base" line). The outlet temperature of both fluids follow the measured value within a 1% difference (cf. "Measured" and "Base" line in Figure 12b). A sensitivity analysis on the effect of the increase in volume flow rate of thermal water is carried out and shown in Figure 13. To account for differences between measurement and simulation over the entire simulation time, the relative root mean square error RRMSE is considered: where N is the number of time intervals, x is the measured value and x i is the respective simulated value. Figure 13 shows that an increase of 2% of volume flow rate of thermal brine leads to the minimum RRMSE in temperature. In fact, by increasing the thermal brine volume flow rate of 2% over the entire simulation, the RRMSE for the outlet temperatures goes below 0.2% in both Cases 1 and 2. The time behavior in pressure and temperature with the 2% increase in thermal brine volume flow rate can be observed in the "+2%" line in respectively Figures 11 and 12. Figure 14 shows the measured and simulated liquid level in Cases 1 and 2. It can be seen that the agreement is not as accurate as for pressure and temperatures. In Case 1, the simulated level shows an offset below 5% relative to the measured data (cf. "Measured" and "Base" line in Figure 14a). When the flow rate of thermal brine is increased by 2% over the entire simulation ("+2%" line), the offset becomes negative. Referring to the maximum level (i.e., evaporator full of liquid), the difference from the measured value remains below 3% in the "Base" line and within −7% in the "+2%" line. In Case 2 (Figure 14b), the offset in the "Base" line is larger than in Case 1, but it is reduced to less than 10% (relative to the measured value) when the thermal water flow rate is increased by 2% ("+2%" line). Referring to the maximum evaporator level, the differences are lower than 12% in the "Base" line and 2% in the "+2%" line. The effect of the increase in volume flow rate of thermal brine on the RRMSE for the evaporator liquid level is depicted in Figure 15a. In Case 1, the offset would be minimum for no increase of thermal brine volume flow rate, whereas for Case 2 the increase should be around 2%.
From Figure 14a, a growing tendency in simulated evaporator level over time can be observed ("Base" line), which does not appear in the measured value ("Measured" line). The growing tendency is made noticeable by the significant simulation time (2000 s), and analyzed in Figure 15b, where the level increase in the simulation over the average value of the measured level is shown. The level Energies 2017, 10, 548 18 of 28 increase drops as the volume flow rate of thermal brine is increased. The effect can also be seen by comparing the "Base" and "+2%" line in Figure 14a. The difference becomes even smaller for higher volume flow rate of thermal brine. The growing tendency of the liquid level in Case 1 seems to be generated by a non-sufficiently accurate inlet mass flow rate of organic fluid in the evaporator. between measurement and simulation over the entire simulation time, the relative root mean square error RRMSE is considered: where is the number of time intervals, is the measured value and is the respective simulated value. Figure 13 shows that an increase of 2% of volume flow rate of thermal brine leads to the minimum RRMSE in temperature. In fact, by increasing the thermal brine volume flow rate of 2% over the entire simulation, the RRMSE for the outlet temperatures goes below 0.2% in both Cases 1 and 2. The time behavior in pressure and temperature with the 2% increase in thermal brine volume flow rate can be observed in the "+2%" line in respectively Figures 11 and 12. (a) (b)   Figure 14 shows the measured and simulated liquid level in Cases 1 and 2. It can be seen that the agreement is not as accurate as for pressure and temperatures. In Case 1, the simulated level shows an offset below 5% relative to the measured data (cf. "Measured" and "Base" line in Figure  14a). When the flow rate of thermal brine is increased by 2% over the entire simulation ("+2%" line), the offset becomes negative. Referring to the maximum level (i.e., evaporator full of liquid), the difference from the measured value remains below 3% in the "Base" line and within −7% in the "+2%" line. In Case 2 (Figure 14b), the offset in the "Base" line is larger than in Case 1, but it is reduced to less than 10% (relative to the measured value) when the thermal water flow rate is increased by 2% ("+2%" line). Referring to the maximum evaporator level, the differences are lower than 12% in the "Base" line and 2% in the "+2%" line. The effect of the increase in volume flow rate of thermal brine on the RRMSE for the evaporator liquid level is depicted in Figure 15a. In Case 1, the offset would be minimum for no increase of thermal brine volume flow rate, whereas for Case 2 the increase should be around 2%. From Figure 14a, a growing tendency in simulated evaporator level over time can be observed ("Base" line), which does not appear in the measured value ("Measured" line). The growing tendency is made noticeable by the significant simulation time (2000 s), and analyzed in Figure 15b, where the level increase in the simulation over the average value of the measured level is shown. The level increase drops as the volume flow rate of thermal brine is increased. The effect can also be seen by comparing the "Base" and "+2%" line in Figure 14a. The difference becomes even smaller for higher volume flow rate of thermal brine. The growing tendency of the liquid level in Case 1 seems to be generated by a non-sufficiently accurate inlet mass flow rate of organic fluid in the evaporator.
This quantity was in fact approximated as the product of the volume flow rate measured at the pump outlet and the density computed from the measured pressure at the pump outlet and the temperature before the pump. In addition to this, the measured data are available only at a minute sampling time between two consecutive measurements. As an example, Figure 14a shows the liquid  This quantity was in fact approximated as the product of the volume flow rate measured at the pump outlet and the density computed from the measured pressure at the pump outlet and the temperature before the pump. In addition to this, the measured data are available only at a minute sampling time between two consecutive measurements. As an example, Figure 14a shows the liquid level in the evaporator ("Sine" line), when this is subjected to a sinusoidal inlet mass flow rate of organic fluid oscillating around the same average value as the measured one for 0 < t < 1400 s, with a period of 180 s and an amplitude of 2 kg/s. The liquid level also oscillates around an average value, with no growing tendency. This confirms the hypothesis of the influence of inaccuracies on the inlet flow rate of the organic fluid on the growing trend for the liquid level. For t > 1400 s, the liquid level suddenly increases because of the negative step in volume flow rate of thermal brine (Figure 10a), which can no longer lead to evaporation of the same amount of organic fluid. If the inlet mass flow rate of organic fluid continues to sinusoidally oscillate around the same average value and the evaporation rate drops, the liquid level grows as shown in Figure 14a.

Discussion
In the present analysis, a dynamic model of a kettle boiler was developed in Modelica ® language, by partially making use of models provided by the TIL-library and partially creating new components or features. The boiler model can be effectively useful for larger scale ORC plants, where kettle boilers are very often applied. The model was tested against available measurements from a geothermal CHP plant in Sauerlach, Germany. The pressure of the organic fluid and the outlet temperatures of both organic and heating fluid could be reproduced within a 1% error in case of a negative (Case 1) or positive (Case 2) step change in inlet volume flow rate of the heating fluid.
In Section 3, the pressure difference between top and bottom of the liquid volume in the evaporator was neglected. The assumption of negligible pressure difference in the liquid level is justified by the fact that the pressure difference is typically below 0.5% of the absolute pressure at the top of the volume. Nevertheless, if such pressure difference is considered, no significant pressure and temperature difference of the vapor is found. The average pressure in the liquid volume would be slightly higher because of the liquid column, and this can lead to a variation of liquid level to 2% with respect to the simplified case. The heat transfer coefficient for pool nucleate boiling (with the Cooper correlation) increases by less than 1%. The accuracy on the liquid level in the evaporator could not in general be as high as for the evaporator pressure and outlet temperatures. In the attempt to explain the reasons for the liquid level offset (and the improved temperature matching) when increasing the thermal water flow rate by 2% over the entire simulation time, the fluid properties of the thermal water and the working fluid were analyzed, but no clear explanation through these parameters could be found. It might be assumed that the improvement achieved by increasing the volume flow rate of thermal brine could be caused by a deviation in latent heat of vaporization for the working fluid. If the volume flow rate of thermal brine is increased and the amount of working fluid that evaporates does not change, the latent heat of vaporization of the working fluid should be higher in the simulation than in the real case. Using the Peng-Robinson equation instead of REFPROP, the calculated heat of vaporization becomes even smaller, in contrast with the hypothesis. Future work should focus on analyzing the working fluid used in Sauerlach and comparing to the REFPROP database to gain major information on the differences in fluid properties.
Because of its simplicity and better agreement with the results, the Cooper correlation was used in this work. The Stephan and Mostinski correlations (Equations (23) and (24)) would lead to a higher liquid level because the lower heat transfer coefficient requires a higher heat transfer area to vaporize the same amount of fluid.

Summary and Conclusions
Dynamic simulations of ORC systems have been recently gaining significant interest. A dynamic model of a kettle boiler was presented in this paper. The model was developed as a combination of a discretized finite volume method on the tube side, and two volumes in non-equilibrium for the evaporating fluid on the shell side. The simulation results were compared to measurements available from a geothermal CHP plant in Sauerlach, Germany. Two cases were investigated: one case where the evaporator undergoes a negative step variation in inlet volume flow rate of thermal brine, and one case where a positive step variation is instead applied. The boundary conditions were set as measured. The simulation results showed a good agreement with the measurements, for what regards the evaporator pressure and outlet temperatures of the working and heating fluids. An increase in 2% of the inlet volume flow rate of the heating fluid over the entire simulation could lead to a better agreement with the temperature measurements. The liquid level was the most challenging quantity to compare, because of measurement uncertainties and possible discrepancies in fluid properties. The model can reproduce with high calculation speed the dynamic behavior of the kettle boiler, keeping the error in pressure and outlet temperatures below 1%. The model offers a high flexibility in geometry of the heat exchanger, being able to handle a multi-pass scheme. It can be implemented for future cycle optimization, in both design and off-design conditions. In future work, the entire cycle in Sauerlach will be simulated, and validated with additional measurements. The fluid properties should also be analyzed more thoroughly. As a result of the simulations, solution for an improved dynamics and control of the power plant will be proposed. The results will be valid not only for the power plant in Sauerlach, but also for other similar power plants.
First pass : Subsequent passes 2 ≤ k ≤ N passes : When w f k = 1, the pass k is completely submerged, when w f k = 0 it is completely dry. If only partially submerged, the pass has 0 < w f k < 1. It is clear, that all cells of the same pass have the same area fraction w f . In addition, when a pass k is partially or completely submerged, all previous passes 1, 2, ..., k − 1 must be completely submerged, i.e., w f p = 1 for 1 ≤ p < k − 1. The concept is illustrated in Figure A1. Attention must be paid to the fact that in this figure the cross-sectional area of each pass A sk and the submerged area A l refer only to the tube cross-sectional area. In other words, the free cross-sectional area of the shell (white area of Figure 6b) is excluded.
Given the pipe submerged fraction w f , the transfer area for a tube cell with the liquid can be computed. For a discretization of each pass into N cells , given the length L of the single pipe and given N p pipes per pass, the ratio of the heat transfer surface A k with respect to the fluid cross-sectional area A cs,k for the i-th cell in the pass k is:

Appendix B
To compute the liquid volume as a function of the liquid level , the evaporator is divided axially in 5 regions, as shown in Figure A2. The internal and external regions are cylindrical, and the two regions in between (2 and 4) are conic, with inclined axis, since the bottom line follows the

Appendix B
To compute the liquid volume V l as a function of the liquid level l l , the evaporator is divided axially in 5 regions, as shown in Figure A2. The internal and external regions are cylindrical, and the two regions in between (2 and 4) are conic, with inclined axis, since the bottom line follows the horizontal plane of the ground where the evaporator is mounted.
For the cylindrical regions, the following areas are defined: And, with R s the kettle inner radius: The liquid volume of regions 1 and 5, each of the length L edge , is: For the liquid volume of region 3, a demister at the top of the kettle has to be considered ( Figure A2): where V d is the demister volume. The liquid volume of region 3, given the length L int , becomes: For the conic regions 2 and 4, each cone of length L cone is subdivided axially in N cylinders of equal length (see Figure A3): Fixing an x-axis on the horizontal bottom line of the cone starting from the end of bigger radiusR s , the position of each sub-cylinder j = 1, 2, .., N is: For the conic regions 2 and 4, each cone of length is subdivided axially in N cylinders of equal length (see Figure A3): Fixing an x-axis on the horizontal bottom line of the cone starting from the end of bigger radius , the position of each sub-cylinder = 1, 2, . . , is: Given the slope of the upper side of the cone as: The radius of each cylinder is: In this way, the volume of each small cylinder is computed as for regions 1, 3 and 5, substituting the radius R or R s with R cylind,j : If l l < 2R cylind,j A cylind,j = R cylind,j 2 acos 1 − l l R cylind,j + l l − R cylind,j R cylind,j 2 + l l − R cylind,j 2 (A20) else A cylindr,j = πR cylind,j 2 (A21) The volume of the conic region is therefore: The liquid volume as a function of the liquid level is therefore: V l = V l,int + 2V l,cone + 2V l,edge (A23) Energies 2017, 10, 548 24 of 27 Figure A3. Discretization of the conic region.