Batch evaporation power cycle: Influence of thermal inertia and residence time

Accepted Manuscript Batch evaporation power cycle: Influence of thermal inertia and residence time Moritz Gleinser, Christoph Wieland, Hartmut Spliethoff PII:

S0360-5442(18)30979-4

DOI:

10.1016/j.energy.2018.05.145

Reference:

EGY 12980

To appear in:

Energy

Received Date: 21 December 2017 Revised Date:

26 March 2018

Accepted Date: 22 May 2018

Please cite this article as: Gleinser M, Wieland C, Spliethoff H, Batch evaporation power cycle: Influence of thermal inertia and residence time, Energy (2018), doi: 10.1016/j.energy.2018.05.145. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Batch Evaporation Power Cycle: Influence of Thermal Inertia and Residence Time Moritz Gleinsera,∗, Christoph Wielanda , Hartmut Spliethoffa,b a Institute

for Energy Systems, Technical University of Munich, Boltzmannstrasse 15, 85748 Garching, Germany Center for Applied Energy Research (ZAE Bayern), Walther-Meissner-Strasse 6, 85748 Garching, Germany

RI PT

b Bavarian

Abstract

M AN U

SC

The transition in the energy market and the growing share of renewable energy sources have been boosting the research in new power cycles. For example, the concept of batch evaporation in the Misselhorn Cycle promises to increase the overall efficiency in low-temperature applications and therefore saves resources. In this paper, a dynamic evaporator model was extended in order to prove the feasibility of the Misselhorn Cycle despite its transient character. In this context, the thermal capacity of the wall material as well as the residence time of the heat source medium were added. The previous, underlying model predicted an improved system efficiency for the Misselhorn Cycle of about 50 % compared to an Organic Rankine Cycle (ORC) at 100 ◦C. Initially, the results of the extended model showed a negative influence of the inertial effects on the possible net power output (advantage over ORC only 10 %). However, an unheated discharge phase and reduced dimensions of the heat exchanger could compensate these drawbacks and achieved results (about 40 % better than ORC) in the same range as the previous, simple model predicted. These findings prove the general practical feasibility of the Misselhorn Cycle.

D

Keywords: Misselhorn Cycle, dynamic simulation, waste heat recovery, trilateral cycle

10

15

EP

In the middle of numerous discussions about renewable energy sources, sustainability and the need to increase efficiency, engineering research questions have become more and more complex recently. In the exemplary field of energy systems, both the fluctuations induced by renewable energy sources [1, 2] as well as the introduction of dynamic operation modes [3, 4] and dynamic processes [5] cause a need for transient analyzing tools. This trend is also reflected in the disproportionately increased listing of “dynamic” and “transient” types of simulation topics in Scopus over the last decade (see Figure 1). One of these dynamic applications is the Misselhorn Cycle, which was first introduced by “Maschinenwerk Misselhorn MWM GmbH” and first published by Gleinser and Wieland [6]. This power cycle uses transient batch evaporation to recover low temperature waste heat. As shown by the authors, the system efficiency can, for instance at 100 ◦C, be increased by about 50 % compared to a continuously operated Organic Rankine Cycle (ORC). While there are numerous well-established techniques to calculate stationary processes, a dynamic application, like the Misselhorn Cycle, often requires complex models and simulations. For example, considering heat transfer, long standing and well-known procedures, as introduced by Bell [7] or Kern [8], provide accurate results for complex shell and tube heat exchangers. However, those methods are only valid for

AC C

5

TE

1. Introduction

∗ Corresponding

author Email addresses: [email protected] (Moritz Gleinser), [email protected] (Christoph Wieland), [email protected] (Hartmut Spliethoff)

Preprint submitted to Energy

May 23, 2018

ACCEPTED MANUSCRIPT

2.5

2500

2

2000

1500

"Simulation" "Dynamic Simulation" "Transient Simulation"

1

RI PT

1.5

1000

0.5

Documents by year

Documents by year

105

500 Status: November 2017

1990

2000

Year

2010

0 2020

SC

0 1980

20

R R steady-state operations. Even modern commercial software packages like Aspen Plus or Ebsilon Professional are not able to map a transient behavior easily. Instead, specialized software packages and libraries for dynamic processes [3, 9] or custom-built models [10, 11] have to be used. The custom-built model by Gleinser and Wieland [6] only focuses on the general thermodynamic context. It is still unclear whether a practical implementation is feasible with regard to the thermal inertia and the oscillating process conditions that occur in this cycle. Therefore, the model used in [6] was extended to cover the thermal inertia of the wall material and include the residence time and heat capacity of the heat source medium. In the present paper, this extended model for the evaporation in a common plate heat exchanger is introduced to analyze the various transient effects. In the Misselhorn Cycle (MWM), the convective heat transfer on the hot side of the heat exchanger is governed by the forced flow of the heat source medium like in most stationary cycles. In contrast, the evaporation of the working fluid can best be described as pool boiling due to the batch-wise character of the Misselhorn Cycle. Moreover, the closed batch evaporation causes a continuously changing pressure as well as a shifting liquid level on the pool boiling side. In particular, the heat capacity of the wall material and the moving liquid level are expected to have a high impact on the process behavior. The evaporator model consists of partial differential equations that are derived from the general equations for mass and energy conservation, equations for heat conduction, and coefficients for heat transfer. As this set of equations is too complex to be solved analytically, the results have to be obtained by numerical simulations. The most straightforward approaches for dynamic flow problems with heat transfer are moving boundary models [12, 13]. The problem is split into several sub-problems, each showing only minor or smooth changes in their dependent variables (e.g. temperature or heat transfer coefficients in the two sections of pre-heating and evaporation in the heat exchanger). The time development of those comparably large sub-problems can then be handled as just one cell per section. For more complex problems, where various changes occur, reasonable split positions cannot always be found. In these cases, fully discretized models based on finite differences [14, 15] or finite volumes [3, 9, 16–18] are recommendable. They can be used on various problems and are more robust but are also computationally expensive for problems with sharp gradients in one or more dependent variables [19]. While a coarse grid would generally suffice for the regions with smooth changes, the region with sharp gradients needs a finer resolution. At last, the uniform grid spacing requires the whole domain to be discretized as finely as the few complex parts. A combination of both moving boundary and finite discretization can then reduce the number of cells. If the problem is split at the sharp step, a medium discretization of the sections in-between yields satisfactory results in most cases [20]. In the case of the Misselhorn Cycle, the surface level of the working fluid clearly separates the two regions of heat transfer to the boiling liquid and heat transfer to the gas phase. However, the effects of the heat

40

45

50

EP

35

AC C

30

TE

D

25

M AN U

Figure 1: Disproportionately increasing topic occurrence for dynamic and transient simulations in Scopus (Title, Abstract, Keywords).

2

ACCEPTED MANUSCRIPT

...

Batch

Expansion Engine

Feed Pump Condenser

RI PT

...

n Evaporators

SC

Figure 2: Process flow diagram of the Misselhorn Cycle [6].

2. System Description

75

80

EP

70

2.1. The Misselhorn Cycle The general setup of the Misselhorn Cycle is quite similar to a common ORC and includes a working fluid cycle with evaporation, an expansion engine, a condenser and a pump (see Figure 2). Instead of the continuous evaporation in the ORC, the MWM cycle uses batch evaporation. In a first phase, the evaporator is filled with working fluid. Then, in the second phase, the inlet valve is closed, the evaporator is heated and an isochoric closed batch evaporation takes place. At the end of this evaporation phase, the working fluid is only partially vaporized. The outlet valve is opened and the built-up pressurized vapor fraction is gradually discharged to an expansion engine in the third and last phase. The liquid fraction remains in the heat exchanger and re-evaporates during the discharge phase due to the pressure drop and further heating. To compensate the batch-wise discharging, at least three evaporators are operated in shifted phases. Due to the isochoric evaporation, the boiling pressure will rise over the second phase and a continuously increasing boiling temperature can be realized. Therefore, the pinch point limitation, known from the ORC, can be avoided. To further improve the temperature match, the hot heat source is fed to the heat exchanger in the last phase. The outlet of the heat source from this stage is then redirected to the heat exchanger operating in the previous phase and further on. The combination of batch evaporation and cascaded heat source flow allows for an almost perfect temperature match of heat source and working fluid. Thus, the Misselhorn Cycle combines all the advantages of a trilateral cycle [22, 23] but avoids the drawbacks of single-phase heat transfer and two-phase expansion. In accordance with the experience from a first experimental model of the cycle, six common, commercially available plate heat exchangers with chevron-type plates were chosen to achieve the necessary heat transfer.

AC C

65

TE

D

60

M AN U

55

conduction and the heat capacity of the plates blur the transition area, and the implementation of a sharp cut at the moving boundary is not viable. Therefore, a fully discretized model was chosen in the present work. The introduced model was used to analyze the influence of the transient effects on the feasibility of the process. In addition, different possible methods of operation were compared to the former basic model with regard to power output and efficiency. A genetic algorithm was deployed to optimize the process parameters with the objective of a maximized system efficiency. The general simulation approach is similar to other dynamic models. However, there seems to be a gap in the targeted applications in all the recent research. On one side, comparable work focuses on the dynamic transitions between different operation points of otherwise steady-state and continuous cycles in a time scale > 10 min [3, 10]. Others cover rapid flash evaporation on very short time scales (< 1 s) [5, 21]. To our knowledge, there is no previous research that covers the effects in a batch evaporation power cycle with a comparable setup which is itself transient and operates in the time scale of about 1 min.

3

ACCEPTED MANUSCRIPT

2.2. Methods of Analysis Several efficiencies and numbers were defined to rate the performance of the Misselhorn Cycle and its different operation points following the definitions in [6]. Depending on the application area, the “performance” of a cycle can, for example, be the power output or one of the following efficiencies. In waste heat recovery applications, the intended target or “performance” is normally a maximized power output from a given heat source. In order to account for the general thermodynamic limitations of low-temperature heat sources, all the efficiencies are based on exergy instead of energy. The specific exergy ex is the maximum amount of work that could theoretically be gained from a given heat source:

RI PT

90

Their low costs, the high surface density and their availability in various types and sizes make them suitable for the Misselhorn Cycle. The cycle sequence with six heat exchangers is filling, evaporation 1-4 and discharge.

SC

85

ex = (h − h0 ) − T0 · (s − s0 ) .

D

100

In this definition, h and s are the specific enthalpy and specific entropy of the heat source. T is the temperature in K. The index 0 denotes the reference condition which is normally set to ambient conditions or the properties of the heat sink (also referred to as “dead state”). The usage of exergetic efficiencies facilitates the illustration of how well a cycle uses the accessible potential of a heat source. The main goal of this cycle is to gain maximum output from a given heat source. Therefore, the cycle should reach high efficiencies in the heat transfer process in the heat exchangers (external second law efficiency ηII,ext ) as well as in the conversion from transferred heat to power (internal second law efficiency ηII,int ). The heat exchanger efficiency, or external second law efficiency, is defined by the ratio of the exergy that ˙ in ) versus the overall exergy that is available from the heat source is actually transferred to the cycle (Ex ˙ av ) (Ex ˙ in Ex m ˙ HS · (exHS,in − exHS,out ) = . (2) ηII,ext = ˙ m ˙ HS · (exHS,in − exHS,0 ) Exav | {z }

M AN U

95

(1)

TE

=0

The index HS hereby marks properties of the heat source. The quality of the cycle itself is expressed by the internal second law efficiency, which is defined by

EP

ηII,int =

110

(3)

The net power of the cycle Pnet is the engine power Poutput reduced by the demand of the feed pump and the auxiliary pumps for heat source and heat sink (Pauxiliary ). This efficiency is used to describe the ability ˙ in ) to power. of the cycle to convert actually transferred exergy (Ex The previous two efficiencies are then combined for an overall rating of the cycle. The resulting second law efficiency (also known as exergetic system efficiency) is defined as the product of Equation (2) and Equation (3) ˙ in Pnet Ex Pnet ηII = ηII,ext · ηII,int = · = (4) ˙ ˙ ˙ av Exav Exin Ex

AC C

105

Poutput − Pauxiliary Pnet = . ˙ in ˙ in Ex Ex

and represents the ratio of the net power output versus the available exergy from the heat source. In accordance with the main goal—maximum output from a given heat source—this exergetic system efficiency is chosen as the objective function for the optimization algorithm. Following the definition in Equation (4), a maximized exergetic system efficiency ηII will also match a maximum net power output Pnet for a given heat source.

4

ACCEPTED MANUSCRIPT

3. Model Description

120

SC

• is suitable to find not only local but global minima/maxima

RI PT

115

R The presented model was entirely built in MATLAB . Although MATLAB may be slower than lowlevel language implementations, it is well known for its fast prototyping ability and the numerous built-in functions and toolboxes which makes it suitable for the engineering development of new models. The thermodynamic states were obtained from the NIST Standard Reference Database REFPROP V9.1 [24] and the therein included equations. The dynamic character of the batch evaporation and the dependency of pressure and temperature trends on the progression of the isochoric evaporation lead to a complex interconnection of the free variables and the timing of the cycle. In order to find a suitable operation point, a simple parameter variation was not promising. A genetic algorithm (GA) was chosen as optimization algorithm because of its broad field of application. Similar to evolutionary development, an initial generation of random sets of variables is optimized by means of combination, mutation and selection. In particular, the genetic algorithm:

• can deal with boundaries (min/max heat exchanger dimensions and mass flows) • is not gradient-based and therefore suitable for discontinuous problems

M AN U

125

• supports the optimization of integer variables (optimization of number of heat exchangers in further research). Like for most optimization algorithms, the disadvantage of the GA are the thousands of function calls that are necessary during the optimization.

145

150

155

AC C

140

EP

TE

135

3.1. System Boundaries and Model Setup In order to maintain a reasonable computation time, the detailed design of the heat exchanger geometry is beyond the scope of this paper. In a first stage of the project, the previous model [6] was focused on the general cycle setup and mode of operation. A 0-D model was used to cover the conservation of mass and energy and stick to basic thermodynamic relations. This current model of the second stage is focused on the effects of the thermally inert mass of the wall material and the residence time of the heat source medium within the heat exchanger. The cycle was implemented as a 1-D model and fixed characteristics like pressure drop and heat transfer correlations were used. The focus lay still on the cycle and the coverage of a wide set of parameters and not on the modeling itself. The results will be used in a planned third stage as initial parameters for a detailed 2-D or 3-D model of the evaporators and also as design parameters for a lab scale test rig. The model will then be evaluated against measured data. With regard to the numerous function calls by the optimization algorithm, the evaporators were only modeled as far as necessary to cover all relevant effects. The exact geometry (plate dimensions and wavy structure) is only needed for the flow velocity and the correlations for the pressure drop. Based on these values, the heat transfer could be calculated. Instead of the computation of the intermediate but otherwise unnecessary properties, overall values for heat transfer and pressure drop were directly set by the optimization algorithm and the reference conditions. In order to comply with realistic values, the boundaries of the optimization were compared to common values for plate heat exchangers and to correlations for heat transfer coefficients (see Section 4.1). The area and mass of the plates were then derived from the overall heat transfer, but the length, height and wavy structure were not defined in detail. Based on these simplifications, the whole plate heat exchanger was represented through the repetition of one symmetric sector with a single, plain plate and the surrounding fluids (see Figure 3). The heat source medium (on the left-hand side) was calculated as one-dimensional forced flow. As there are only minor changes in the temperature and pressure of the heat source, its thermodynamic properties (density ρ, heat capacity cp ) were assumed to be constant. In addition, the assumption of an incompressible fluid leads to a constant mass flow. The momentum balance and detailed mass balance were therefore neglected. The flow rate of the heat source medium also outweighs the neglected axial heat conduction. The equation for energy conservation thus only includes the mass-bound heat transport, the pressure drop (which

D

130

5

ACCEPTED MANUSCRIPT

1D

1D

0D

Heat Transfer

Wall

Working fluid

SC

Heat source

RI PT

Heat Conduction

Heat Transfer

Symmetry

z

Symmetry

M AN U

x

Figure 3: System boundaries for the plate heat exchanger model.

TE

D

160

includes the gravitational pressure drop) and the heat transfer to the wall. The latter one was calculated from an estimated heat transfer coefficient for forced convection in plate heat exchangers (αHS ) and the local temperature difference between the fluid and the wall surface for each cell (∆T ). The wall was also implemented as a one-dimensional part of the model. It was discretized into the same number of cells in z-direction as the heat source and into one layer in x-direction. The temperature was assumed to be uniform in each of these cells. Because of the small plate thickness, the heat transfer from the heat source to the wall and from the wall to the working fluid is limiting over the heat conduction λ in the x-direction. The Biot number α Bi = , (5) λwall /dwall

170

175

EP

AC C

165

expresses this ratio of the heat transfer coefficient to the heat conductivity of the wall material. With a given low number in the order of magnitude of Bi ≤ 0.161 , the temperature distribution across the wall is almost uniform and the error from the 1-D implementation is small. In contrast, along the height of the plate, the small cross-section leads to a fairly limited heat conduction and the heat transfer on the surface is the governing effect. Consequently, the temperature in z-direction was expected to follow the profile of the heat source medium rather than be uniform. However, the primary goal of this model was to account for the global effects of thermal inertia on this transient thermodynamic cycle. In this context, it is most important to cover the conservation of heat, the heat capacity of the wall material, and the heat transfer through the wall. The local effects like the exact temperature profile and the small amount of longitudinal heat conduction were therefore of secondary importance. The heat conduction in z-direction was still included to roughly capture and analyze the effects on the temperature gradient. When looking at the detailed results, one has to keep in mind that the assumption of a uniform cell temperature does not exactly map the real longitudinal temperature profile. The side of the working fluid (the right-hand side) was only calculated as one cell with uniform temperature and pressure distribution and an ideal equilibrium state. Research on flash evaporation in pools of liquid [25, 26] shows that the time scale to reach the equilibrium is small (a few seconds). Furthermore, the rapid mixing by the rising bubbles will prevent any thermal stratification [27, 28]. However, the heat transfer coefficient is dependent on the liquid level, as there may be nucleate pool boiling in the lower liquid phase and only natural convection in the upper gas phase. Therefore, the liquid level was still calculated 1 For

an exemplary heat transfer coefficient αHS = 8000 W/(m2 · K), λ = 50 W/(m · K) and a plate thickness of d = 1 mm.

6

ACCEPTED MANUSCRIPT

185

based on the liquid volume, which, in turn, is derived from the overall mass, the vapor quality and the liquid density of the working fluid. The mass balance of the working fluid had to be calculated in order to cover the transient character of the filling and discharge as well as the changes in the overall density due to the evaporation. One goal of the simulations was to be as universal as possible. Without the definition of a detailed geometry, the following overall relations were used:

RI PT

180

• Instead of geometric dimensions, the size of one heat exchanger was defined by the heat transfer ability k · A. This figure combines the overall boiling heat transfer coefficient k with the overall heat transfer area A of the whole evaporator. It was only used as a combined constant (kA)all = koverall · Aoverall

SC

that was varied by the optimization.

(6)

• The overall area could be calculated on an interim basis from the value of (kA)all with the ratio of the estimated heat transfer coefficients for the side of the heat source αHS , the side of the working fluid αWF and the heat conductivity of the wall material λwall :

M AN U

1 d 1 1 = + + ; (kA)all αHS A λwall A αWF A   1 d 1 A = (kA)all · + + . αHS λwall αWF

(7) (8)

TE

D

• The optimization algorithm varied the plate thickness d and the amplitude a of the chevron structure. The plates are in contact at two opposing peaks of the wavy structure. In between those points of contact, the expanding structure builds flow channels. The average distance of two plates is therefore two times the wave amplitude (2a). The heat source medium in this plate gap is split into two layers, where each of them is associated with one of the neighboring plates. The thickness of one heat source layer is therefore only half of the whole flow channel (see Figure 4). The mass of all plates mwall and the mass of the heat source medium mHS that is enclosed in the heat exchanger could then be calculated from (9) (10)

EP

mwall = A · d · ρwall mHS = A · a · ρHS

with the overall area A, the thickness of the wall d and heat source layer a, and the corresponding densities for the wall material ρwall and heat source medium ρHS .

AC C

• In order to avoid the calculation of the flow velocity (and to avoid the definition of the detailed geometry), an alternative figure to represent the residence time was needed. Independently from the exact number and geometry of the plates, the dwell time of the heat source medium in the heat exchanger is always mHS τ= , (11) m ˙ HS where m ˙ HS is the mass flow rate of the heat source. Based on this definition, the whole heat exchanger could be represented by one single model plate when the dwell time was used instead of flow velocities.

190

• The dimensionless fraction of one cell on the whole system is zfrac =

∆z 1 = L Ncells

(12)

for both the heat source side and the wall. The z-fraction and the cell size ∆z are directly related to the number of cells Ncells . The length L of this model plate only affects the flow velocities in the 7

Qz

m

ACCEPTED MANUSCRIPT

T_wall(z,t)

T_HS(z,t)

T_WF(t)

Q_WF

ρ = const. d

a

Heat source

Wall

Working fluid

z

a

x

RI PT

m

Qz

dz

Q_HS

heat exchanger. The velocity is involved in the calculation of the pressure drop and heat transfer coefficients. As both of them were set by the optimization algorithm, the plate dimension did not have any influence in this model and could be set to L = 1.

M AN U

195

SC

Figure 4: Control volume of the plate heat exchanger.

Following these specification, only a minimal set of geometric parameters ((kA)all , plate thickness and wall gap) was needed in order to obtain all relevant thermal properties of the heat exchanger (thermal inertial mass of wall material and enclosed mass of the heat source medium). 3.2. Energy Balances and Differential Equations Figure 4 shows a control volume of the heat exchanger on which the energy balances are based. The assumptions made in Section 3.1 lead to the following simplified model equations. For the hot side of the system, the general energy balance from the control volume is

EP

TE

D

∂h ∂THS = −m ˙ HS · ∆z − kAHS · zfrac · (THS − Twall ) , (13) ∂t ∂z which includes the energy storage in one cell, the mass-bound energy transport and the heat transfer to the wall. The missing information about the exact dimensions of one cell was compensated by using the zfraction (equation (12)) of the overall values. The overall heat transfer (kA)all was converted to the hot-side heat transfer (kA)HS = αHS · A via the ratio between αHS and αWF from Equation (7):    d 1 (kA)HS = (kA)all · 1 + αHS + . (14) λwall αWF mHS · zfrac · cp ·

Furthermore, the derivative of the enthalpy can be replaced by dh = cp dT + vdp for incompressible fluids [29]. Solving (13) for the time derivative of T and substituting the relevant parts with the relationship from Equation (12) then leads to  

AC C

200

∂THS 1 = − ·L ∂t |τ {z }

constant

 ∂T  HS · −  ∂z

1 ∂p · ρHS · cp ∂z | {z }

 (kA)HS  · (THS − Twall,HS ) , −  mHS · cp | {z }

constant for ∂p/∂z=∆p/L

(15)

constant

which includes the flow-bound heat transport, the influence of the pressure drop and the heat transfer to the wall. The pressure drop gradient could be approximated with the mean, constant pressure drop ∆p/L. For the wall, the basic differential equations is mwall · zfrac · cwall

∂Twall ∂ q˙ = + (kA)HS (THS − Twall ) − (kA)WF (Twall − TWF ) −Across · · ∆z ∂t ∂z | {z } (I)

8

(16)

ACCEPTED MANUSCRIPT

which includes the heat transfer from the heat source to the wall and from the wall to the working fluid. The cold-side heat transfer (kA)WF was derived in analogy to Equation (14). As a third factor, the heat conduction within the wall in z-direction is dependent on the heat flux gradient ∂ q/∂z ˙ and the cross section Across of the wall in z-direction. When solving Equation (16) for ∂Twall /∂t, the right term (I) becomes Across · ∆z 1 ∂2T · ·λ 2. mwall · zfrac cwall ∂z

(17)

RI PT

+

The left factor of Equation (17) can be rewritten with the definition from Equation (12) as Across · L Vwall 1 = = . mwall mwall ρwall

SC

The complete differential equation for the wall then leads to

1 λ ∂ 2 Twall ∂Twall . = · ((kA)HS · ∆THS − (kA)WF · ∆TWF ) + · ∂t mwall · zfrac · cwall ρwall · cwall ∂z 2

(18)

(19)

Q˙ WF =

M AN U

For the working fluid, the energy balance is quite simple due to the assumption of a homogeneous state. The summed up heat flow over all cells of the whole plate L X

(kA)WF · zfrac · ∆TWF (z)

(20)

z=0

205

and the in- and outgoing enthalpy flows during filling and discharge are the only influencing factors on the overall energy balance of the working fluid:

215

3.3. Model Discretization As the partial differential equations cannot be handled analytically, numerical solvers have to be used. The basic structure of the problem had to be analyzed to choose the appropriate method.

AC C

210

EP

TE

D

dUWF =m ˙ WF,in · hWF,in − m ˙ WF,out · hWF,out + Q˙ WF (21) dt The overall density of the working fluid could be obtained with the given volume of the evaporator and the overall mass balance dmWF =m ˙ WF,in − m ˙ WF,out . (22) dt In combination with the internal energy UWF from Equation (21), the thermodynamic state of the working fluid could be calculated by REFPROP [24].

3.3.1. Heat Source The heat source equation (15) describes a flux-conservative initial value problem. The fluid conditions at the inlet (z = 0) were known at all times and propagate in the direction of the mass-flow. The initial conditions for the whole fluid at t = 0 could also be approximated, e.g. by the temperature profile of a stationary process. Based on these initial values, the development of the fluid state could be calculated in time and space step by step. Appropriate algorithms for the hyperbolic differential Equation (15) can be the Lax method, the explicit upwind scheme or fully implicit schemes [30]. Considering the relatively simple flow problem on the heat source side, a basic explicit upwind scheme was used. The heat source domain was discretized in a grid, in which the grid points are positioned on the cell borders (see Figure 5). The values of the heat source flow (temperature and pressure) at the inlet and outlet could thus directly be stored in the grid structure. The average temperature of each cell that was

9

ACCEPTED MANUSCRIPT

j= 1

2

3

4

5

6

THS,avg

Heat Source

Staggerd Grid

heatflux

Wall

n=3 n=2 n=1

Twall

2

3

4

5

RI PT

j= 1

z

Figure 5: Grid over the heat exchanger height z and the time t: a cell-average heat source temperature is computed for heat flux calculations.

SC

used for heat transfer calculations to the wall was computed from THS,avg = 1/2 · (THS,j−1 + THS,j ). The discretization scheme for explicit upwind calculations is defined as n Tjn+1 − Tjn Tjn − Tj−1 = −u · , ∆t ∆z

(23)

220

3.3.2. Wall The circumstances are more complex for the wall, as there are boundary conditions at the bottom and top where the plate was defined as adiabatic, and therefore the heat flow has to be zero. In contrast to flux-conservative initial value problems that can be solved from start to end, in these diffusive initial value problems all conditions must be satisfied at once. General algorithms for this parabolic Equation (19) are the explicit FTCS (forward time centered space), the fully implicit or the Crank-Nicholson-Scheme. However, in the latter two implicit schemes, the whole spatial scope has to be solved at once and the resulting large system of equations cannot be solved step by step. Instead, complex algorithms like the shooting method or relaxation methods are recommended [30]. As discussed earlier (section 3.1), the heat conduction over the length of the plate was expected to be less important than the heat storage in the wall. Explicit schemes tend to be more accurate with regards to energy conservation, whereas implicit schemes are more stable by converging to the equilibrium of the transport phenomenons [30]. With the major scope of this model being the global effects of the heat storage, the explicit FTCS scheme was chosen. As a positive side effect, explicit schemes are also faster regarding computation time. The grid points were positioned in the center of each cell to represent the average cell temperature (see Figure 5). The FTCS scheme is defined by

230

AC C

EP

TE

D

225

M AN U

where j is the spacial index and n is the time index. The velocity of the wave propagation u = L/τ as well as further summands were inserted according to Equation (15).

n n Tjn+1 − Tjn Tj−1 − 2Tjn + Tj+1 =D· , ∆t (∆z)2

(24)

where D is the diffusion coefficient. From Equation (19) the diffusion coefficient follows as D = λ/(ρc). Further summands were added accordingly. 3.3.3. Stability Conditions and Cell Size The stability of the differencing schemes is dependent on the step-size in space and time. From a mathematical point of view, the explicit upwind scheme is stable if the Courant condition is met: |u| ∆t ≤1 ∆z

(25)

From a more practical point of view, this means that the wave propagation |u| · ∆t must not travel further than the cell size ∆z within one time-step. 10

ACCEPTED MANUSCRIPT

The stability of the wall equation is guaranteed for 2D∆t ≤ 1. (∆z)2

RI PT

240

Again, this physically represents that the time-step has to be smaller than the heat conduction time across one cell (∆z)2 /D. With the focus on the heat storage rather than the exact transport profiles, a large cell size with only 10 cells along the heat exchanger was chosen in favor of computing time. The initial dimensions of the heat exchanger (see Section 4.1) led to an estimated flow velocity of about 0.1 m/s. The combination of the cell size ∆x = 0.1 m and a time step of ∆t = 0.1 s therefore met the Courant condition from Equation (25). The cell size of the wall, which is coupled to the discretization of the heat source, also fulfilled the stability condition from Equation (26) with the given time step. Both of the stability conditions were smaller than 1 by the factor of 10−1 and 10−4 , respectively. The numerical calculations were thus guaranteed to be mathematically stable even for some variation of the initial heat exchanger geometry during the optimization.

SC

235

M AN U

4. Results 245

(26)

4.1. Reference Conditions To test the evaporation model, a reference case was defined with the conditions shown in Table 1. The heat source was water with a temperature of 100 ◦C and a pressure of 4 bar, which could be either a Table 1: Reference conditions

Exergy Reference

3.2 kg/s 100 ◦C 4 bar

variable 25 ◦C 4 bar

set by optimization 30 ◦C 8.2 bar

25 ◦C 1 bar

D

Working Fluid from Condenser (R134a)

geothermal source or from industrial waste heat. Cooling water with 25 ◦C and 4 bar was used as a heat sink. The mass flow of the cooling water was dependent on the cooling demand and the pinch point in the condenser (5 K). The heat source mass flow was set to 3.2 kg/s to achieve an available heat flow of 1 MW. The heat exchanger was filled with working fluid coming from the condenser with 30 ◦C and a pressure of 8.2 bar. The amount of heat source medium that was filled into the evaporator was adjusted by the optimization algorithm. The properties of the cycle components are given in Table 2. The expansion engine and the pumps were simplified and modeled by an isentropic expansion and compression based on the mass flow and the enthalpy difference. The isentropic efficiencies were 90 % for the feed and auxiliary pumps and 80 % for the expansion engine. The pressure drop of the hot and cold water was estimated with 0.3 bar for each heat exchanger. Those pressure losses were also included in the auxiliary power demand. The pressure loss of the working fluid in the condenser was set to 0.2 bar.

255

AC C

EP

250

Heat Sink (water)

TE

mass flow, m ˙ temperature, T pressure, p

Heat Source (water)

Table 2: Parameters of the cycle components

Component Property

Value

isentropic pump efficiency isentropic engine efficiency pressure drop water (evaporator & condenser) pressure drop working fluid (condenser)

11

90 % 80 % 0.3 bar 0.2 bar

ACCEPTED MANUSCRIPT

Table 3: Estimated heat transfer coefficients

270

liquid boiling liquid natural convection gas natural convection

10 000 1000 150

RI PT

8000 1000

The heat transfer coefficients were estimated for a common geometry with different correlations from the literature (see Table 3). For the side of the heat source medium, the correlations of Martin [31] resulted in a heat transfer coefficient in the scale of αHS ≈ 8000 W/(m2 · K) for a forced convective, single-phase flow in a plate heat exchanger. Pool boiling and free convection are not very common in plate heat exchangers, and therefore only correlations for generic pool boiling and convection on vertical plates could be found as an alternative. The boiling of the liquid working fluid was estimated by the pool boiling correlations from Rohsenow [32, 33] with αWF ≈ 10 000 W/(m2 · K). For temperatures lower than the limit for the onset of nucleate boiling [34], the equation from Churchill and Chu [35] resulted in αWF ≈ 1000 W/(m2 · K) for natural convection on a vertical plate in the liquid phase. This value was also valid for natural convection on the heat source side. With the fluid properties of the gas phase, the natural convective heat transfer coefficient decreased to αWF ≈ 150 W/(m2 · K). Despite their age, these correlations are still widely recommended by current literature [36–38]. The comparison with results from the Aspen Exchanger Design & Rating tool confirmed the validity of these values. In addition, the overall heat transfer coefficient according to Equation (7) (k ≈ 4000 W/(m2 · K)) lay well within the common literature values. Shah and Sekuli´c [36], for example, give the general range of 3000 W/(m2 · K) to 7000 W/(m2 · K) for plate heat exchangers. 4.2. Test Cases The batch operation of the Misselhorn Cycle causes a large number of independent process parameters compared to those in the stationary operated ORC. Therefore, the genetic algorithm toolbox of MATLAB was used in order to find the optimal set of input parameters for each case. The range of possible values is based on common dimensions of heat exchangers [36] and can be grouped in geometrical and operational parameters (see Table 4).

TE

D

275

liquid forced convection liquid natural convection

SC

265

cold side (αWF in W/(m2 · K))

M AN U

260

hot side (αHS in W/(m2 · K))

Table 4: Process parameters: Default values and boundaries for optimization

EP

Parameter

AC C

plate thickness wave amplitude total volume for working fluid heat transfer ability mass of working fluid per evaporator mass flow of working fluid switch pressure limit discharge volume flow

Unit s (mm) a (mm) V (l) kAall (kW/(m2 · K)) m (kg) m ˙ (kg/s) p (bar) V˙ (m3 /min)

Default 1 3 100 30 40 20 13 5

Boundaries Lower Upper 0.5 1.5 10 10 1 1 9 0.75

1.2 7.0 800 40 40 20 20 40

280

285

The former group includes the dimensions of the heat exchanger: plate thickness d, wave amplitude a and the overall heat transfer (kA)all . In addition, the volume for the built-up working fluid vapor is extended by a pressure storage to a total volume V . If only the free volume of the evaporator was available, the built vapor rapidly increased the pressure and the evaporation stopped after just a few seconds due to the rising boiling temperature. Providing the total volume by a bigger evaporator instead of a simple, adiabatic vessel hardly improved the process, as the volume of the liquid is small and the heat transfer to the gas phase is almost negligible compared to the boiling heat transfer. 12

ACCEPTED MANUSCRIPT

1000

100

800

60

600

40

400

evaporation filling

1

2

3

4

200

0 10

20

30

800

60

600

40

400

20

0 0

40

50

200

0

60

1000

80

disscharge

20

1200

RI PT

80

Temperature Heat Source in Temperature Heat Source out Temperature Working Fluid Average Wall Temperature Heat Flow to Working Fluid Heat Flow to Wall

Heat Flow in kW

120

Temperature in °C

Temperature in °C

100

1200

Heat Flow in kW

Temperature Heat Source in Temperature Heat Hource out Temperature Working Fluid Average Heat Flow

120

0

0

20

60

80

Time in s

SC

Time in s

40

(a) NTU-model

(b) current model (1), standard case

290

The latter group of parameters contains the total mass of working fluid m that is filled in one evaporator during each cycle pass, the mass flow of this filling m ˙ WF,in and the pressure limit of the evaporator p that defines the end of the discharge phase and controls the switching of the cycle phases. The lower boundary of the switching pressure was limited by the engine outlet pressure (8.4 bar). The duration of the discharge phase is also influenced by the discharge volume flow V˙ (and the corresponding mass flow rate m ˙ WF,out ), which results from the engine revolution speed and the engine displacement. For comparison, the optimization was run on the previous model from Gleinser and Wieland [6] (subsequently referred to as the “NTU-model”) with the reference conditions from Section 4.1. The temperature profiles over one cycle pass are shown in Figure 6a, where the continuously rising boiling temperature can be seen. The dotted vertical lines indicate the six phases of the cycle (from left to right: filling, evaporation 1-4, discharge). Furthermore, it is clearly visible how the heat source is fed to the cycle in the discharge phase at constant temperature and gradually cools down during the following phases due to the cascaded flow pattern. The underlying quasi-stationary NTU-approach (number of transfer units with a logarithmic mean temperature difference) did not include heat capacities and the residence time of the heat source. The resulting net power of the cycle was 49.9 kW, which equates to a second law efficiency of 45.8 % (see Table 5). The optimized process parameters of this and the following cases can be found in Table 6. In

EP

300

TE

D

295

M AN U

Figure 6: Temperature profiles over one cycle pass (see the web version for colored figures)

Table 5: Cycle performance for the different simulation approaches

AC C

Model

net power, Pnet (kW) second law efficiency, ηII (%) net energy per pass, Enet (MJ) cycle duration, t (s)

NTU

ORC (1)

49.9 45.8 2.9 58.1

33.7 30.9 -

37.2 34.1 3.1 83.4

Current model (2) (3) 46.8 42.9 3.4 73.1

46.1 42.3 3.3 72.5

(4) 52.6 48.3 3.4 64.7

addition, the figures of an ORC under the reference conditions are given in Table 5 as well. 305

4.2.1. Standard Case (1) The current extended model was first tested with the original cycle setup to provide a standard case and reference for further adjustments. In Figure 6b, the temperature profiles over one cycle pass are shown. The blue lines (see also the web version for colored figures) mark the in- and outlet temperature of the 13

ACCEPTED MANUSCRIPT

Table 6: Optimized process parameters

Parameter

Unit

NTU (1)

Cycletime: 0.8s, HX: 1

0.8

0.6

0.4

0.6 1.5 260 13.5 40 20 11.7 4.35

1 1.5 104 19.7 40 20 11.7 4.88

0.9 1.5 86 21.7 40 20 11.9 4.95

0.9 0.4 135 39.6 40 20 11.4 5.5

M AN U

Dimensionless z-Direction (-)

1

236 39.5 36 19 10.8 5.55

RI PT

s (mm) a (mm) V (l) kAall (kW/(m2 · K)) m (kg) m ˙ (kg/s) p (bar) V˙ (m3 /min)

SC

plate thickness wave amplitude total volume for working fluid heat transfer ability mass of working fluid per cycle mass flow of working fluid switch pressure limit discharge volume flow

Current model (2) (3) (4)

Temperature Heat Source in Temperature Wall Temperature Working Fluid Liquid Level

0.2

0 40

50

60

70

80

90

100

Temperature in °C

320

325

EP

315

heat source, the green line marks the working fluid temperature and the yellow line shows the averaged temperature of the wall material. The corresponding heat flow from the heat source to the wall and from the wall to the working fluid is also displayed in the bottom of the figure. In comparison to the former model, the added inertial factors caused two major effects. First, the added implementation of the wall showed that the plate is further heated and stays hot during the discharge phase even though the working fluid temperature decreases following the falling process pressure (right section of Figure 6b). Second, at the end of the discharge phase, the whole system is at a high temperature. When the cycle is then switched to the next filling phase, the evaporator is feed from the bottom with heat source medium that was already cooled down in previous phases of the other heat exchangers. Figure 7 shows a snapshot of the temperature profiles over the plate shortly after the beginning of the new filling phase. During those first seconds of the next filling phase, the hot wall is covered with the cold working fluid from the condenser. It can be seen, that the stored heat is rapidly transferred from the wall to the fresh working fluid and, in some parts, even to the heat source medium. Furthermore, due to the finite flow velocity of the heat source, there remains an amount of hot fluid within the evaporator. Therefore, the heat source leaving the heat exchanger is still hot for the duration of the residence time. As the evaporator in the filling phase is the end of the cascaded setup, the remaining exergy in this high temperature fraction is lost. The negative heat flow from the wall to the heat source and the high temperature outlet can also be seen in the heat flow profiles (left part of Figure 6b). To partly compensate these dynamic effects of the inertial mass, the optimization ended up at a higher cycle duration of 83 s instead of 58 s. This led to a lower power output (37.2 kW) and second law efficiency (34.1 %) even though the two models generated about the same amount of net energy Enet per cycle pass (about 3 MJ, see Table 5).

AC C

310

TE

D

Figure 7: Snapshot of the temperature profiles over one plate shortly after the switchover from the discharge to the next filling phase: standard model (1) shows effects of residence time and thermal inertia

14

ACCEPTED MANUSCRIPT

1000

100

800

60

600

40

400

20

200

0 20

40

60

800

60

600

40

400

20

0 0

80

200

0

80

0

0

Time in s

1000

RI PT

80

1200

Heat Flow in kW

120

Temperature Heat Source in Temperature Heat Source out Average Heat Source Temperature Temperature Working Fluid Average Wall Temperature Heat Flow to Working Fluid Heat Flow to Wall

20

40

60

SC

Temperature in °C

100

1200

Temperature in °C

120

Heat Flow in kW

Temperature Heat Source in Temperature Heat Source out Average Heat Source Temperature Temperature Working Fluid Average Wall Temperature Heat Flow to Working Fluid Heat Flow to Wall

80

Time in s

(a) current model (2), unheated discharge phase

(b) current model (3), unheated discharge and filling phase

330

4.2.2. Unheated Discharge (2) To overcome the disadvantages of the standard cycle setup (1), an unheated discharge phase was implemented. During the last period, no hot water was fed to the heat exchanger and the prior heat source medium remained inside. Without the forced flow, the heat transfer switches to the moderate natural convection. Missing the permanent heating, the wall and the enclosed heat source are cooled down during the discharge phase (see Figure 8a). In addition, the hottest heat source could be used to reach even higher maximal boiling temperatures and pressures instead of destroying exergy during the discharge phase due to a large temperature difference. The cooling of the wall and the higher maximum pressure allow for a shorter cycle duration and increased the power output again. The optimized process parameters resulted in a net power of 46.8 kW and a second law efficiency of 42.9 %.

350

355

EP

345

4.2.3. Unheated Filling and Discharge (3) Figure 8a shows that even with the unheated discharge phase, the wall temperature does not exactly match the working fluid in the first seconds of the filling phase and there is still some heat stored in the wall material. In addition, the average temperature of the heat source during the discharge indicates that the heat source is not completely cooled down, too. In extension to the previous case, an unheated filling phase was therefore added. It can be seen in Figure 8b that the absence of heat input, and especially the lower heat transfer of the natural convection, soften the reaction during the first seconds of the filling phase. The heat flows are reduced to a moderate level, and the stored heat can be almost completely utilized. However, beyond that, the major part of the filling phase shows hardly any heat transfer at all. This period of the process does not contribute to the pressure build-up, and no improvement compared to the previous cases could be achieved. The slightly lower cycle time is outweighed by the reduced transferred energy. The resulting second law efficiency of 42.3 % (46.1 kW) was slightly lower than for case (2).

AC C

340

TE

D

335

M AN U

Figure 8: Temperature profiles over one cycle pass

4.2.4. Decreased Heat Capacity (4) The optimized parameters in Table 6 show that the optimal wave amplitude for cases (1) to (3) matches the lower boundary, whereas the plate thickness lies in-between the boundaries. This indicates that the energy, which is stored in the heat source within the wall gap, is the major influencing factor causing the thermal inertia. To further analyze this effect, the lower boundaries for both plate gap and wall thickness were lowered by the factor 10 in this case. 15

ACCEPTED MANUSCRIPT

Temperature Heat Source in Temperature Heat Source out Average Heat Source Temperature Temperature Working Fluid Average Wall Temperature Heat Flow to Working Fluid Heat Flow to Wall

1000

80

800

60

600

40

400

20

200

0

0 0

10

20

30

40

60

70

SC

Time in s

50

RI PT

Temperature in °C

100

1200

Heat Flow in kW

120

Figure 9: Temperature profiles over one cycle pass: current model (4), smaller geometry

5. Discussion

5.1. Comparison of the Different Cases It becomes obvious from case (1) that the thermal inertia initially prevents the cycle from reaching the system efficiency that was predicted by the simple NTU model. The optimized values in Table 6 show that the optimization algorithms minimized the thermal capacity by different means. First, it set the wave amplitude to the smallest possible value to reduce the mass and residence time of the heat source that is held in the heat exchanger. Second, it reduced the plate thickness for a smaller mass of the wall material. As these two measures did not suffice, the heat transfer ability was reduced. A smaller heat transfer area results in a smaller heat exchanger and a further reduced mass of heat source and wall material. However, due to the reduced heat transfer, the cycle needs more time to provide the heat for the evaporation, and the cycle duration increases. Thus, the optimal point of operation is a compromise of a small enough heat exchanger to prevent high thermal capacity but a big enough heat exchanger to provide the needed heat transfer. The unheated discharge in case (2) compensates these limitations. The reached efficiencies were almost the same as the ones predicted by the NTU model. A closer look at the optimized process parameters shows that the heat transfer ability was about 60 % higher than for case (1). The thermally relevant mass is already cooled down during the discharge phase and therefore the restrictions on the size of the heat exchanger could be eased. A temperature cross, as seen in Figure 7 can also be avoided. The enhanced heat transfer allows for a shorter cycle duration and an increased efficiency. In addition, the maximal process temperature increased from 82.5 ◦C to 92.5 ◦C. The corresponding boost of the maximal boiling pressure from 27.8 bar to 34.2 bar also contributed to the efficiency improvement. The changes in case (3) provided no obvious improvements. The cycle duration, the reached power output and the efficiencies were almost the same as for case (2). Only a closer look at the first seconds of the filling phase reveals that, without the unheated filling, the temperature gradients are very high. At the end of the previous discharge phase, the wall still has an average temperature of 50 ◦C. The average temperature of the enclosed heat source is even higher at about 70 ◦C. When the cold working fluid coming from the condenser at 30 ◦C is filled into the evaporator, the wall material is exposed to high thermal stress.

380

385

390

EP

375

AC C

370

TE

D

365

The overall cycle sequence of this setup (Figure 9) looks quite similar to the unheated discharge case (2). However, the lower wave amplitude results in a smaller wall gap in the heat exchanger. Due to the reduced mass of the enclosed heat source, the heat source outlet temperature always matches the working fluid temperature in this case. So, the reduced heat capacity allows for a shorter cycle duration of only 65 s, and the same amount of transferred energy (3.4 MJ) caused a power output of 52.6 kW (efficiency 48.3 %).

M AN U

360

16

ACCEPTED MANUSCRIPT

415

420

RI PT

SC

M AN U

410

5.2. Assessment and Future Considerations The four test cases showed that the limitations from the thermal inertia of the heat exchangers can be handled by the right process operation. The first essential part is the unheated discharge phase, which smooths the rapid transition from hot to cold operation points. The second influencing factor is the ratio of working fluid mass to the thermal capacity of the evaporator. It can be seen from the optimized process parameters that the algorithm always set the mass of the working fluid to the highest possible value and concurrently reduced the thermal capacity of the whole heat exchanger as far as possible. This stresses the importance of well-chosen boundaries for the optimization. In the purely thermodynamic context, only some of the free variables had an optimum within their range, while others always stuck to their upper or lower limit. With the widened boundaries in Case (4), the plate gap was set to the lower limit. At the same time, the heat exchanger area was sized to the upper limit as long as the thermal capacity was small enough. Extreme values result in high thermodynamic system efficiencies, but they are not always feasible once practical or economic aspects are considered. Therefore, the boundaries for the optimization algorithm were chosen in accordance with common dimensions of heat exchangers and reasonable process parameters in this model. Once the objective function of the optimization is set to economic aspects, extreme values for the free variables will lead to poor economic numbers. If the model is properly defined, the costs for a large amount of working fluid and the enormous pressure losses in a very small plate gap will restrict extreme configurations themselves. It is expected that the optimization will then find an optimum for every variable within a reasonable range.

D

405

TE

400

6. Conclusion

430

435

An extended dynamic model of a plate heat exchanger including the thermal capacity of the wall material and the residence time of the heat source medium was presented. Several test cases with regard to the transient Misselhorn Cycle were compared. The results of the standard case showed that the heat capacity and thermal inertia have a negative effect on the net power output of the cycle. However, adjustments to the operation mode, like an unheated discharge phase, can compensate this influence. The optimized models reached efficiencies in the same range as the previous, simple model predicted. The simulations prove the general feasibility of a practical implementation of the Misselhorn Cycle. With those initial doubts put aside, the advantage in system efficiency of over 40 % compared to the common ORC becomes more realistic. Further research regarding the Misselhorn Cycle should now focus on economic aspects. While there are still some open questions with regard to the more complicated setup of the Misselhorn Cycle, the high efficiency and the simplicity of the single components should compensate for these additional costs.

AC C

425

EP

395

Although the unheated filling phase cannot lower this wall temperature, the cooling of the enclosed heat source down to 50 ◦C diminishes the wall stress. While this effect is not relevant for purely thermodynamic analysis, it should be considered during the design and rating of the evaporators in more detailed studies (like economic feasibility). In case (4), the optimization algorithm fully used the lowered boundaries for the heat exchanger geometry. The set plate thickness was still at 0.9 mm, but with a = 0.4 mm, the wave amplitude was set close to the new lower limit of 0.15 mm. While this small distance between the plates is neither a common, commercial value nor is it expected to be applicable in practice, it provides a good idea of the detailed influence of thermal inertia. These new dimensions were small enough to reduce the thermal capacity to an appropriate amount. Therefore, there was no more limit for the size of the heat exchanger, and the heat transfer ability was set to its upper boundary of 40 kW/(m2 · K). The improved heat transfer allows for even shorter cycle durations, and the power output is further increased to 52.6 kW.

17

ACCEPTED MANUSCRIPT

7. Acknowledgment

RI PT

This work was performed and funded in the context of a research project with Maschinenwerk Misselhorn MWM GmbH. The MWM project was sponsored by the Bavarian Ministry of Economic Affairs and Media, Energy and Technology in the context of the BayINVENT program (funding program for innovative energy technologies and energy efficiency). Special thanks go to Manfred Moullion, MWM, for the close cooperation and the detailed information about the test cycle. 8. Declaration of Interest

445

The authors declare no conflict of interest. The founding sponsors had no role in the design of the study; in the collection, analyzes or interpretation of data; in the writing of the manuscript; nor in the decision to publish the results. 9. Abbreviations

m2 kJ/(kg · K) m m2 /s J kJ/kg kW kJ/kg kW W/(m2 · K) m kg kg/s bar kW kW/m2 kW kJ/(kg · K) s K m/s kJ m3 /kg m

EP

TE

D

area specific (isobaric) heat capacity plate thickness diffusion coefficient Energy specific exergy exergy flow rate specific enthalpy enthalpy flow rate spatial index global heat transfer coefficient height of plate mass time index quantity mass flow rate pressure power heat flux heat flow rate specific entropy time temperature velocity of wave propagation internal energy specific volume cell size

AC C

450

A c, (cp ) d D E ex ˙ Ex h H˙ j k L m n N m ˙ p P q˙ Q˙ s t T u U v ∆z

M AN U

The following abbreviations are used in this manuscript:

Greek Symbols α local heat transfer coefficient η efficiency λ heat conductivity ρ density τ residence time

SC

440

W/(m2 · K) % W/(m · K) kg/m3 s 18

ACCEPTED MANUSCRIPT

References

SC

455

RI PT

Subscripts 0 reference conditions av available avg average ext external HS heat source HX heat exchanger in input int internal out output WF working fluid

480

485

490

495

500

M AN U

D

475

TE

470

EP

465

AC C

460

[1] Buttler, A., Hentschel, J., Kahlert, S., Angerer, M.. Statusbericht Flexibilit¨ atsbedarf im Stromsektor: Eine Analyse der aktuellen marktwirtschaftlichen und technischen Herausforderungen an Speicher und Kraftwerke im Zuge der Energiewende. 2015. URL: https://www.es.mw.tum.de/fileadmin/w00bhq/www/pdf/Statusbericht_ Flexibilitaetsbedarf_2014_Final.pdf. [2] Nitsch, J., Pregger, T., Naegler, T., Heide, D., Luca de Tena, D., Trieb, F., et al. Langfristszenarien und Strategien f¨ ur den Ausbau der erneuerbaren Energien in Deutschland bei Ber¨ ucksichtigung der Entwicklung in Europa und global: BMU - FKZ 03MAP146. 2012. URL: http://elib.dlr.de/76043/. [3] Quoilin, S., Aumann, R., Grill, A., Schuster, A., Lemort, V., Spliethoff, H.. Dynamic modeling and optimal control strategy of waste heat recovery Organic Rankine Cycles. Applied Energy 2011;88(6):2183–2190. doi:10.1016/j.apenergy. 2011.01.015. [4] Orlandini, V., Pierobon, L., Schløer, S., de Pascale, A., Haglind, F.. Dynamic performance of a novel offshore power system integrated with a wind farm. Energy 2016;109:236–247. doi:10.1016/j.energy.2016.04.073. [5] Steffen, M., L¨ offler, M., Schaber, K.. Efficiency of a new Triangle Cycle with flash evaporation in a piston engine. Energy 2013;57:295–307. doi:10.1016/j.energy.2012.11.054. [6] Gleinser, M., Wieland, C.. The Misselhorn Cycle: Batch-Evaporation Process for Efficient Low-Temperature Waste Heat Recovery. Energies 2016;9(5):337. doi:10.3390/en9050337. [7] Bell, K.J.. Delaware method for shell side design. In: Kaka¸c, S., Bergles, A.E., Mayinger, F., editors. Heat exchangers. Washington and Auckland: Hemisphere Pub. Corp. and Distribution outside the U.S., McGraw-Hill International Book Co. ISBN 0-07-033284-3; 1981, p. 581–618. [8] Kern, D.Q.. Process heat transfer. McGraw-Hill classic textbook reissue series; New York: McGraw-Hill; 1950. ISBN 978-0070341906. [9] Desideri, A., Hernandez, A., Gusev, S., van den Broek, M., Lemort, V., Quoilin, S.. Steady-state and dynamic validation of a small-scale waste heat recovery system using the ThermoCycle Modelica library. Energy 2016;115:684–696. doi:10.1016/j.energy.2016.09.004. [10] Mazzi, N., Rech, S., Lazzaretto, A.. Off-design dynamic model of a real Organic Rankine Cycle system fuelled by exhaust gases from industrial processes. Energy 2015;90:537–551. doi:10.1016/j.energy.2015.07.083. [11] Yousefzadeh, M., Uzgoren, E.. Mass-conserving dynamic organic Rankine cycle model to investigate the link between mass distribution and system state. Energy 2015;93:1128–1139. doi:10.1016/j.energy.2015.09.102. [12] Bell, I.H., Quoilin, S., Georges, E., Braun, J.E., Groll, E.A., Horton, W.T., et al. A generalized moving-boundary algorithm to predict the heat transfer rate of counterflow heat exchangers for any phase configuration. Applied Thermal Engineering 2015;79:192–201. doi:10.1016/j.applthermaleng.2014.12.028. [13] Li, D., Luo, K., Dang, J.. A moving boundary model for two-phase flow heat exchanger incorporated with relative velocities between boundaries and fluid. International Journal of Heat and Mass Transfer 2016;95:35–44. doi:10.1016/j. ijheatmasstransfer.2015.11.095. [14] Lo Brano, V., Ciulla, G., Piacentino, A., Cardona, F.. Finite difference thermal model of a latent heat storage system coupled with a photovoltaic device: Description and experimental validation. Renewable Energy 2014;68:181–193. doi:10.1016/j.renene.2014.01.043. [15] Mandrusiak, G.D., Carey, V.P.. A finite difference computational model of annular film-flow boiling and two-phase flow in vertical channels with offset strip fins. International Journal of Multiphase Flow 1990;16(6):1071–1096. doi:10.1016/ 0301-9322(90)90107-T. [16] Quoilin, S., Bell, I., Desideri, A., Dewallef, P., Lemort, V.. Methods to Increase the Robustness of Finite-Volume Flow Models in Thermodynamic Systems. Energies 2014;7(3):1621–1640. doi:10.3390/en7031621. [17] MacArthur, J., Grald, E.. Unsteady compressible two-phase flow model for predicting cyclic heat pump performance and a comparison with experimental data. International Journal of Refrigeration 1989;12(1):29–41. doi:10.1016/0140-7007(89) 90009-1. [18] Jensen, J.M.. Dynamic Modeling of Thermo-Fluid Systems: With Focus on Evaporators for Refrigeration. Dissertation; Technical University of Denmark; Lyngby; 2003.

19

ACCEPTED MANUSCRIPT

530

535

540

RI PT

SC

M AN U

525

D

520

TE

515

EP

510

AC C

505

[19] Bendapudi, S., Braun, J.E., Groll, E.A.. A comparison of moving-boundary and finite-volume formulations for transients in centrifugal chillers. International Journal of Refrigeration 2008;31(8):1437–1452. doi:10.1016/j.ijrefrig.2008.03.006. [20] Feru, E., de Jager, B., Willems, F., Steinbuch, M.. Two-phase plate-fin heat exchanger modeling for waste heat recovery systems in diesel engines. Applied Energy 2014;133:183–196. doi:10.1016/j.apenergy.2014.07.073. [21] Kanno, H., Shikazono, N.. Experimental and modeling study on adiabatic two-phase expansion in a cylinder. International Journal of Heat and Mass Transfer 2015;86:755–763. doi:10.1016/j.ijheatmasstransfer.2015.02.059. [22] Fischer, J.. Comparison of trilateral cycles and organic Rankine cycles. Energy 2011;36(10):6208–6219. doi:10.1016/j. energy.2011.07.041. [23] Yari, M., Mehr, A.S., Zare, V., Mahmoudi, S., Rosen, M.A.. Exergoeconomic comparison of TLC (trilateral Rankine cycle), ORC (organic Rankine cycle) and Kalina cycle using a low grade heat source. Energy 2015;83:712–722. doi:10.1016/j.energy.2015.02.080. [24] Lemmon, E.W., Huber, M.L., McLinden, M.O.. NIST Standard Reference Database 23: Reference Fluid Thermodynamic and Transport Properties-REFPROP: Standard Reference Data Program. 2013. [25] Miyatake, O., Murakami, K., Kawata, Y., Fujii, T.. FUNDAMENTAL EXPERIMENTS WITH FLASH EVAPORATION. Heat Transfer - Japanese Research 1973;2(4):89–100. [26] Miyatake, O., Fujii, T., Tanaka, T., Nakaoka, T.. FLASH EVAPORATION PHENOMENA OF POOL WATER. Heat Transfer - Japanese Research 1977;6(2):13–24. [27] Kim, J.I., Lior, N.. Some critical transitions in pool flash evaporation. International Journal of Heat and Mass Transfer 1997;40(10):2363–2372. doi:10.1016/S0017-9310(96)00296-7. [28] Saury, D., Harmand, S., Siroux, M.. Flash evaporation from a water pool: Influence of the liquid height and of the depressurization rate. International Journal of Thermal Sciences 2005;44(10):953–965. doi:10.1016/j.ijthermalsci. 2005.03.005. [29] Baehr, H.D., Kabelac, S.. Thermodynamik: Grundlagen und technische Anwendungen. Springer-Lehrbuch; 14. aufl. ed.; Berlin, Heidelberg: Springer-Verlag Berlin Heidelberg; 2009. ISBN 978-3-642-00555-8. [30] Press, W.H.. Numerical recipes: The art of scientific computing. Cambridge [Cambridgeshire] and New York: Cambridge University Press; 1986. ISBN 0521308119. [31] Martin, H.. A theoretical approach to predict the performance of chevron-type plate heat exchangers. Chemical Engineering and Processing: Process Intensification 1996;35(4):301–310. doi:10.1016/0255-2701(95)04129-X. [32] Rohsenow, W.M.. A method of correlating heat transfer data for surface boiling of liquids. Trans ASME 1952;74:969. [33] Rohsenow, W.M.. Boiling. In: Rohsenow, W.M., Hartnett, J.P., Gani´ c, E.N., editors. Handbook of heat transfer fundamentals. New York and London and Paris: McGraw-Hill. ISBN 0-07-053554-x; 1985, p. Chapter 12. [34] Smith, R.A.. Vaporisers: Selection, design & operation. Harlow and Essex and England, New York: Longman Scientific & Technical and Wiley; 1986. ISBN 0470207094. [35] Churchill, S.W., Chu, H.H.. Correlating equations for laminar and turbulent free convection from a vertical plate. International Journal of Heat and Mass Transfer 1975;18(11):1323–1329. doi:10.1016/0017-9310(75)90243-4. [36] Shah, R.K., Sekuli´ c, D.P.. Fundamentals of heat exchanger design. Hoboken, NJ: John Wiley & Sons; 2003. ISBN 0-471-32171-0. [37] Incropera, F.P., DeWitt, D.P., Bergman, T.L., Lavine, A.S.. Fundamentals of heat and mass transfer. 6 ed.; Hoboken, NJ: Wiley; 2007. ISBN 978-0-471-45728-2. ¨ [38] Marek, R., Nitsche, K.. Praxis der W¨ arme¨ ubertragung: Grundlagen - Anwendungen - Ubungsaufgaben. 3., aktualisierte auflage ed.; M¨ unchen: Fachbuchverl. Leipzig im Hanser Verl.; 2012. ISBN 978-3-446-43241-3.

20

ACCEPTED MANUSCRIPT

Highlights

EP

TE D

M AN U

SC

RI PT

Dynamic heat exchanger model is extended by thermal inertia Thermal inertia initially constricts the performance of the transient power cycle Adjusted process operations compensate for the inertial effects Batch process shows advantages over continuous cycles despite transient character

AC C

   