Nothing Special   »   [go: up one dir, main page]

Energy Conversion and Management: Javier Peña-Lamas, Juan Martinez-Gomez, Mariano Martín, José María Ponce-Ortega

Download as pdf or txt
Download as pdf or txt
You are on page 1of 11

Energy Conversion and Management 165 (2018) 172–182

Contents lists available at ScienceDirect

Energy Conversion and Management


journal homepage: www.elsevier.com/locate/enconman

Optimal production of power from mid-temperature geothermal sources: T


Scale and safety issues

Javier Peña-Lamasa, Juan Martinez-Gomezb,c, Mariano Martína, , José María Ponce-Ortegab
a
Departamento de Ingeniería Química, Universidad de Salamanca, Pza. Caídos 1-5, 37008 Salamanca, Spain
b
Chemical Engineering Department, Universidad Michoacana de San Nicolás de Hidalgo, Edificio V1, Ciudad Universitaria, Morelia, Michoacán 58060, Mexico
c
Departamento de Ingeniería y Ciencias Químicas, Universidad Iberoamericana, Prolongación Paseo de Reforma 880, Lomas de Santa Fe, C.P. 01219 Ciudad de México,
Mexico

A R T I C LE I N FO A B S T R A C T

Keywords: In this work we have optimized the operation of a geothermal facility. Benzene, toluene and cyclohexane are
Geothermal energy selected as candidate fluids based on previous studies. We formulated a superstructure optimization problem for
Safety the optimal use of geothermal energy using a binary cycle with either one or two turbine expansions. Surrogate
Rankine cycle correlations for the thermodynamic properties of the fluids have been developed (i.e. enthalpy, entropy) to
Environmental impact
model the operation of the turbine and the heat exchanger network. A single economic objective and a nor-
Mathematical optimization
malized multiobjective function accounting for economic, environmental and safety issues are used. The eco-
Scale
nomic optimization selects a two expansion cycle using toluene as organic fluid, producing 10.4 MW at 0.075 €/
kWh with an investment of 102 M€. Safety considerations slightly change the operating conditions, reducing the
pressures and temperatures, but not the selection of the working fluid. Sustainable and economic terms over-
come safety issues. The results are competitive with other renewable-based technologies for thermal power
production such as CSP or biomass.

1. Introduction up and evaporate a thermal fluid whose expansion in the turbine will
provide the power [2,3]. Therefore, dry-steam and flash systems are
Geothermal energy is a promising source from within the Earth’s widely used to produce electricity from high-temperature resources,
interior. In spite of its potential, it has not been extensively exploited and binary power plants (or Organic Rankine Cycle units; ORC) are the
yet. High-temperature reservoirs, with temperatures above 220 °C, are best energy conversion systems to exploit lower temperature reservoirs,
the ones most suitable for commercial production of electricity. both from a technical and environmental points of view. Thus, we have
However, the medium- and low-temperature water-dominated systems, focused on these ones.
with temperatures between 110 and 160 °C, are the most abundant [1]. The second issue related to the use of geothermal energy is to
So far most of the work focuses on evaluating the Organic Rankine evaluate the ORC and the proper fluid. There are a number of recent
Cycles (ORC), to determine their capacity of extracting energy from low studies in the literature [4–11]. Most of the studies are based on simple
temperature sources and the fluid that can be used. With regards to the models for the Rankine cycle, comparing a number of fluids [4,8].
structure of the cycle, there are a number of alternatives. We can find, Other works focused on the optimization of the cycle using different
either single- or double-flash geothermal power plants, dry-steam approaches. Franco and Villani [10] proposed a heuristic methodology
power plants and binary cycle power plants. Single-flash plants were to optimize geothermal plants comparing several working fluids, re-
the first type of geothermal facilities installed. The thermal fluid is frigerants and hydrocarbons. Roy and Misra [9] performed a sensitivity
evaporated one time only. Double-flash facilities include a second analysis for the optimization of the cycle for one single fluid, refrigerant
evaporation step to increase the power produced by 15–25%. Dry-steam R123. Another typical refrigerant, R134, was evaluated by Sun and Li
power plants consist of a simpler technology since no condensation is [6]. The authors developed a model based on first principles con-
expected in the expansion and typically only one expansion in the sidering thermal efficiency for the different units such as the expander,
turbine is allowed. Finally, binary cycle power plants are used when the the condenser, an air cooler, and the evaporator in order to optimize the
steam from the well has impurities that prevent from its direct use as power production. Simple thermodynamics is assumed as well as one
well as for low temperature sources. Steam or hot brine is used to heat single expansion. Wang et al. [7] compared a limited number of fluids,


Corresponding author.
E-mail address: mariano.m3@usal.es (M. Martín).

https://doi.org/10.1016/j.enconman.2018.03.048
Received 21 December 2017; Received in revised form 19 February 2018; Accepted 17 March 2018
Available online 24 March 2018
0196-8904/ © 2018 Elsevier Ltd. All rights reserved.
J. Peña-Lamas et al. Energy Conversion and Management 165 (2018) 172–182

Nomenclature Unit equipment from the flowsheet


HX heat exchanger
cp brine heat capacity (kJ/kg °C) IW injection well
f(t) dimensionless time function for heat transfer Mix mixer
f friction factor Pump pumps
fcbrine brine flow rate (kg/s) PW production well
fc(J,unit,unit1) mass flow of component J from unit to unit1 (kg/s) Turb turbine
Hi enthalpy of stage I in the diagram (kJ/kg) Splt splitter
Hb,(unit,unit1) enthalpy of the stream at the state b from the stream SR sands removal unit
from unit to unit1 (kJ/kg) Str storage tank
K thermal conductivity (W·(m°C)−1) Valv valve
k permeability (mDarcys)
Lyac, depth of the source (m) Subscripts (Fig. 2)
Pyac pressure down the well (MPa)
PIN pressure at the production well entrance (MPa) Stage 1 saturated liquid exiting the cooling superstructure, and
Phead pressure at the production well head (MPa) reaching the storage tank Str1
Pturb pressure at the expansion (MPa) Stage 2,s compressed liquid exiting the pump PMP1, supposing
r internal radius of the well (m) ideal compression (efficiency of the process equal to
ryac influence radius of the reservoir (m) 100%)
r well: external radius of the well (m) 0.120 Stage 2 compressed liquid exiting the pump PMP1
r int external radius of the well (m) 0.1133 Stage 3 saturated liquid exiting the heat exchanger HX1 and
sb(unit,unit1) entropy the stream at the state b for the stream from unit reaching HX2, as well as the pump PMP2
to uni1 kJ/kgK Stage 4,s compressed liquid exiting the pump PMP2, supposing
T(unit,unit1) temperature of the stream from unit to unit 1 (°C) ideal compression (efficiency of the process equal to
Tbh temperature down the well (°C) 100%)
To temperature at the bottom when no fluid movement (°C) Stage 4 compressed liquid exiting the pump PMP2
Tsat saturation temperature (°C) Stage 5 saturated liquid exiting HX3
Total Risk risk of the design (yr−1) Stage 6 saturated steam exiting HX4, which constitutes the feed of
U brine velocity (m/s) the first body of the turbine
y vertical distance from the bottom of the well to the surface Stage 7,s overheated steam produced at the expansion of the first
(m) 3600 body of the turbine, supposing ideal expansion (efficiency
of the process equal to 100%)
Symbols Stage 7 overheated steam produced at the expansion of the first
body of the turbine
α temperature gradient (43.49E−3 °C/m) Stage 8 overheated steam produced at the mix of stages 7 and 9, at
∝ thermal diffusivity of the formation (10−6 m2s−1) the inlet of the second body of the turbine
μ brine viscosity (poises) 0.0015854567 Stage 9 saturated steam exiting HX2, which constitutes the feed of
ρyac density of the brine (kg/m3) 922 the second body of the turbine
Stage 10,s overheated exhausted steam produced at the expansion
Subindexes of the second body of the turbine, supposing ideal ex-
pansion (efficiency of the process equal to 100%)
Brine hot brine from the well Stage 10 overheated exhausted steam produced at the expansion of
s isentropic the second body of the turbine
I point in Fig. 2 Stage 11 saturated steam exiting HX5
COOL cooling unit
Fluid thermal fluid

mostly refrigerants, on the performance of a simplified Rankine cycle. Cycles for the recovery of energy from medium and low geothermal
The cycle is optimized following a pinch analysis approach embedded sources. We focus on the three fluids that our previous study identified
in a genetic algorithm. Yu et al. [12] presented a pinch-based method as the most promising from environmental, economic and safety points
for the design of efficient ORC’s to recover waste heat from different of view [13]. In particular, we develop a superstructure for the heat
sources within a refinery. Most of the works considering detailed exchanger network to use the hot brine from a geothermal well al-
thermodynamics have used process simulators [11], otherwise simpli- lowing double and single extractions to optimally design the best cycle
fied flowsheets for the Rankine cycle have been considered. Only lately [14–16], where the thermodynamics of the fluids is included through
Yu et al. [12] presented the design of a heat exchanger network for the surrogate models for the enthalpies and entropies instead of the more
optimal heat recovery form waste heat using detailed thermodynamics theoretical approach [17,18]. Two solutions are presented, the first one
by implementing equations of state within the optimization for one optimizes the energy output and a second one where economic, en-
fluid, n-pentane. vironmental and safety considerations are included in a normalized
While in a previous work a fast screening procedure of fluids for the objective function to evaluate their effect on the selection of the cycle
operation of Organic Rankine Cycles (ORC’s) was carried out based on a and on the operating conditions. Finally, the effect of the scale on the
multiobjective optimization approach considering economic, environ- investment and production costs is evaluated. The paper is organized as
mental and safety issues, but the process was oversimplified. Therefore, follows. Section 2 describes the process flowsheet superstructure. Sec-
in this work we use a superstructure optimization approach for the tion 3 comments on the modelling aspects of all the units involved in
mathematical optimal design and operation of Organic Rankine binary the flowsheet. Section 4 presents the optimization procedure including

173
J. Peña-Lamas et al. Energy Conversion and Management 165 (2018) 172–182

both objective functions. Section 5 summarizes the operating conditions entire flow of thermal fluid is expanded again. The exhaust steam can
and the economic study of the three fluids. Finally, Section 6 draws be condensed using cooling water in a system of up to two heat ex-
some conclusions. changers in parallel to improve the heat transfer efficiency, HX5 and
HX6. Cooling water is split into two streams to improve the temperature
gradient and reduce the area and HX cost, since the operating tem-
2. Process description peratures at HX5 and HX6 are low. Finally, a cooling system is used to
cool down the water used for the condensation of the exhaust steam.
The process consists of three circuits, that of the brine, the thermal Fig. 1 shows the process flowsheet.
fluid and finally the one corresponding to the cooling system. The brine
is extracted at a flow of 850 GPM per well [2]. The pipe is 3600 m long,
and 0.2266 m of internal diameter [19]. These features determine the 3. Modelling
pressure at the well and at the head as well as their temperatures. We
assume a temperature of the well of 170 °C based on the region of study, 3.1. Well modelling
Zaragoza (Spain) [19].
The brine line starts at the production well or wells. The number of In this section we compute the pressure and the temperature of the
production wells determines the total flow of hot brine used to generate brine at the well’s head as a function of the reservoir properties.
power. That stream is the hot source to produce steam to be expanded Assuming monophasic fluid, the temperature of the fluid as it ascends
at the turbine. Once the energy is used, the brine is reinjected into the along the well can be computed as a function of the depth of the well,
ground. A superstructure of heat exchangers is developed to evaluate the temperature gradient inside the Earth and the heat transfer rate as
the optimal heat recovery from the brine. One or two expansions are given by Eq. (1) [2]:
allowed. Therefore, single or dual pressure cycles are considered in the T = (Tbh−a·y ) + a·A·(1−e−y / A) + (T0−Tbh)·e−y / A (1)
superstructure. Four heat exchangers are allocated to heat up and
evaporate the thermal fluid in order of descending temperature, HX4 to where Tbh is the downhole reservoir temperature and To is the
HX1. After the last of the heat exchangers, the brine is reinjected back inflowing fluid temperature (this is referred to the temperature that is
into the Earth. The flow of the thermal fluid is in counter-current with characteristic of the brine if there is no extraction of it from the geo-
respect to the brine. After condensation, it is pumped to HX1. Next, the thermal production well or, in other words, when all the brine is lying
flow of organic fluid is split into two streams. One of the streams is at the bottom of the reservoir without ascending by the well). Usually
evaporated in HX2 and fed to the medium pressure section of the tur- both of them have the same value. a is the temperature gradient
bine, whereas the rest is fed to HX3 after being pumped up. This stream (43.49·10−3 °C/m), y is the vertical distance from the bottom of the well
is heated up again and fed to HX4, where it is evaporated. Finally, it is to the surface (assuming that the feed zone is at the bottom) [2]. A
expanded at the high pressure section of the turbine. After the expan- corresponds to the diffusion depth parameter that is function of time as
sion, it is mixed with the second stream coming from HX2 and the follows:

Fig. 1. Scheme of the superstructure of the process.

174
J. Peña-Lamas et al. Energy Conversion and Management 165 (2018) 172–182

fcBrine ·cp·f (t ) 2fρyac u2y


A (t ) = PIn−PHead = + ρyac gy
2π·K (2) r (5)

where fcbrine is the brine flow rate, cp is the heat capacity of the brine, K
is the thermal conductivity and f(t) is a dimensionless time function that
represents the transient heat transfer to the formation. For stable sys- 3.2. Thermal cycle modelling
tems (flowing times greater than 30 days more or less) f(t) can be ap-
proximated as follows: The operation of the cycle follows the termodynamic cycle shown in
Fig. 2. The brine flows from HX4 to HX3, down to HX2 and finally to
r ⎞ HX1. However, the flow of the thermal fluid is more complex. The or-
f (t ) = −ln ⎛ −0.290
⎝ ∝·t ⎠
2 (3) ganic fluid is recycled from the cooling system, (1) in Fig. 2, pumped up
at PMP1 and fed to HX1. Thus, the blue line from stages (1) to (3) in
where r is the internal radius of the well, ∝ the thermal diffusivity of the
Fig. 2 takes place at HX1. The organic fluid is heated up and split in
formation and t is assumed to be 30 days. In case that the brine flow
Splt2. The split fraction at Splt2 is modelled using Eq. (6),
would not be monophasic, the temperature would be the one that
corresponds to the saturated state, for the precise pressure that the fc(Fluid,HX 1. Splt 2) = fc(Fluid,Splt 2,Turb1) + fc(Fluid,Splt 2,HX 2) (6)
brine has at the production wellhead.
To determine the type of flow, either mono or biphasic, the algo- where fc(Fluid,unit,unit1) is the mass flow of the fluid in the stream that
rithm proposed by DiPippo [20] is typically used. This method com- joins both units. One fraction of the stream is heated up to reach sa-
putes the depth from the heat at which the brine turns biphasic and turated vapor at HX2 and fed as such to the second stage of the turbine.
estimates the pressure needed there so that there is no such a transition The line from stages (3) to (9) in Fig. 2 shows the evaporation of the
when the flow ascends. Alternatively, using the pressure at the bottom thermal fluid in HX2. The other fraction of the stream is pumped again
of the well and Eq. (4), if the pressure computed at the head is above to a higher pressure and it is fed to HX3 to be heated up and to HX4 to
the saturation point of the brine, there will not be a biphasic flow, and generate a vapor stream. The yellow line in Fig. 2 shows the path from
therefore there is no need for an inner pump. Splt2 to PMP2, following to HX3 to heat up the fluid, and up to the
Thus, we have computed the pressure as follows. The pressure drop evaporation stage taking place in HX4 is represented by stages (5) and
between the bottom of the reservoir and the entrance of the production (6) in Fig. 2. This vapor stream is fed to the high pressure section of the
well is computed as given by Eq. (4) turbine.
Other main decision variables are the pressures at each of the tur-
μ·fcbrine,well ·ln (ryac ./ rwell )
Pyac−PIN = bine stages, as well as the feed temperature to the first body of the
2π·k·L yac .·ρyac . (4) turbine and the pressure of the exhaust steam computed by Eq. (7)

where Pyac is the pressure down the well, PIN is the pressure at the shaft’s fc(Fluid,turb1,HX 5) = fc(Fluid,HX 4,Turb1) + fc(Fluid,HX 2,Turb1) (7)
entrance, μ is brine viscosity, fcbrine is the flow rate of brine and ryac is
the influence radius of the reservoir, r well is the external radius of the Each of the heat exchangers presented in the description of the cycle
well, k the permeability, Lyac, is the depth of the source and ρyac is the is modelled using an energy balance for the brine, another oner for the
density of the brine. Then, the pressure drop along the shaft is com- thermal fluid and providing bounds for the inlet and outlet tempera-
puted as Eq. (5), where Phead is the pressure at the well’s head, u is the tures of both fluids to avoid temperature cross in the heat exchanger,
brine velocity, and y is the well depth, 3600 m [19]. f is the friction Eq. (8):
factor of the brine along the shaft.

Fig. 2. T-S Diagram for the process (Toluene case).

175
J. Peña-Lamas et al. Energy Conversion and Management 165 (2018) 172–182

Q(HXi) = −fc (Brine) ·cp·(T(Brine,Out )−T(Brine,in) ) cooling tower must be designed and integrated with the system [23,25].
Q(HXi) = fc(Fluid) ·(H(Fluid,Out )−H(Fluid,in) ) Thus, in our case we assume that the cooling water is cooled down back
T(Fluid,Out ) + ΔTmin ⩽ T(Brine,in) to ambient temperature in a pond. We estimate the cost of the pond as a
function of the cooling load as given in the Appendix A.
T(Fluid,in) + ΔTmin ⩽ T(Brine,Out ) (8)
where Q(HXi) is the energy transferred in heat exchanger i, cp is the li-
quid heat capacity of the brine and H is the thermal fluid entalphy. The 4. Optimization procedure
heat transferred from the brine is computed using heat capacities.
However, to capture the thermodynamics of the system, correlations for We use a mathematical optimization approach to determine the
the enthalphy and entropy of the organic fluids are developed as a optimal structure of the cycle, the operating conditions and the best
function of the pressure and temperature, Eq. (9). The correlations can fluid considering economic, environmental and safety issues. To eval-
be seen in the Appendix A for all the fluids considered. These correla- uate the effect of safety in process design we optimize the production of
tions are fitted from tabulated data for specific ranges of the variables, geothermal energy using the best fluids identified in a previous work
pressure and temperature. The fitting was good as it could be seen in [13], benzene, cyclohexane and toluene, where the process was sim-
Fig. 2 where the points are the solution of the problem and the two- plified. A multiobjective optimization approach was used to simulta-
phase bell is drawn from experimental data. Furthermore, we make sure neously consider all three aspects, namely economic, environmental
that the solution of the optimization lies within the range of the vari- and social (safety), for the optimal selection of working fluids in organic
ables. Note that saturation stages are characterized by only one vari- Rankine Cycles. In this work two solutions are evaluated.
able, pressure or temperature, and the correlations are simpler. A ΔTmin First, we optimize the operation of such a cycle aiming at maximum
of 5 °C is imposed. power production. We formulate the problem consisting of the equa-
tions that model all the units involved in the flowsheet presented in
H(Fluid,i) = f (TFluid,PFluid ) (9)
Fig. 1 and described briefly in Section 3. It is a nonlinear and nonconvex
sFluid(i) = g (TFluid,PFluid ) formulation whose size is, per fluid, around 204 equations and 262
variables. The model is written in GAMS and solved in an Intel i7,
The two expansions at the turbine, from stages (6) to (7) and from Windows 10 machine taking around 10–15 seconds. The main decision
(8) to (10) in Fig. 2, are assumed to have an isentropic efficiency of 0.85 variables are the split fraction of the organic fluid to HX2 and to PMP2,
based on rules of thumb [21]. Therefore, the stream exiting the first which define whether a one expansion or a two expansion system is
body can be calculated using Eqs. (10), (11) selected, the operating pressures and temperatures are each heat ex-
HFluid(6)−HFluid,(7) changer as well as the split fraction in the cooling cycle. Finally, the
η=
HFluid,(6)−HFluid,(7,s) (10) turbine operation pressures and temperatures and the steam exhaust
conditions.
The efficiency determines the operating temperature computed as Next, additional constraints to account for safety considerations are
Eq. (11): included in the formulation. We developed correlations for each com-
sFluid(6) = f (TFluid,PFluid ) = sFluid(7,s) (11) ponent to compute the effect of pressure and temperature on the risks
involved in any one of the units that constitute the process flowsheet,
Antoine correlation [22] is used as a bound for the points of satu- see Fig. 1, such as continuous and instantaneous release of mass due to
rated vapor, Eq. (12), so that we ensure that a particular stream is not in leaks. These releases can evolve into a set of accidents related to fire
a two-phase flow. Saturated temperatures are also computed using Eq. and explosion (jet fire, pool fire and VCE). These accidents have an
(12) effect or generate radiation, overpressure in addition to toxic con-
⎛A (OF ) − B (OF) ⎞ centrations, which cause harm to people if they are exposed. For the
Pturb1·760 = e⎝ (C (OF ) + TSat) ⎠ (12) sake of the length of the manuscript and because of the fact that
where P is the fluid pressure at saturation. The second expansion is computing the total risk involves a large number of intermediate cor-
modelled similarly. The energy obtained from the expansions is given relations, we refer the reader to the supplementary material and pre-
by Eq. (13): vious work by Martínez-Gomez et al. [13]. By implementing these
correlations, together with the model of the flowsheet presented above,
W(Turb1) = fc(Tol,HX 4,Turb1) ·(H7−H8) + [fc(Tol,HX 4,Turb1) + fc(Tol,HX 2,Turb1) ]·(H8 the problem size increased up to 329 equations and 406 variables, per
−H10) (13) organic fluid. The objective function is also based on the one used in
previous work where economic, environmental and safety objectives
The exhaust steam, stage (10) in Fig. 2, is condensed and com- are normalized [13]. For that, we optimize each objective in-
pressed at PMP1 before restarting the cycle again. We compute the dependently to compute a reference value and, finally, we use Eq. (15)
power consumed by each of the pumps as follows, Eq. (14): as the objective function that involves the calculation of the energy
fc(Fluid,unit ,Pump) produced, the emission mitigated by using geothermal and no fossil
W(Pump) = ΔP·Q = ΔP· resources and the total risk. The model was written in GAMS and solved
ρ (14)
for about 5 min in an Intel i7 Windows 10 machine.
The refrigeration of the exhaust steam is carried out using cooling
water in a system consisting of two heat exchangers so that the cooling GVP ⎤ ⎡ ENVIRONMENTAL ⎤
NF = ⎡ −
water is split into two, HX5 and HX6, to improve the temperature ⎢
⎣ GVPMAX ⎥⎦ ⎢⎣ ENVIRONMENTALMAX ⎥

gradient while the steam is fed in series, first to HX5 and from that to TotalRiskMAX −TotalRisk ⎤
+ ⎡
HX6, before being recycled. Energy balances based on heat capacities ⎢ TotalRiskMAX ⎥
⎣ ⎦ (15)
for the cooling water and enthalpies, as in Eq. (8), are performed.
Finally a scale-up study of the plant is presented to evaluate the
3.3. Cooling stage wells to be drilled as well as the investment and production costs as-
sociated with the production capacity following a methodology as
Power plants use different designs of cooling towers to condense the presented in previous papers [23].
exhaust steam from the turbine [23,24]. However, the small size of the
plant evaluated in this work suggests the use of lagoons, Otherwise, a

176
J. Peña-Lamas et al. Energy Conversion and Management 165 (2018) 172–182

5. Results Table 2
Multi-objective optimization.
5.1. Operation of the plant
Variable Benzene Toluene Cyclohexane

The brine down the well is at 170 °C and 38 MPa. It ascends along Power 9.134 10.458 10.280
the production well shaft and reaches the chemical/thermal power P_max (MPa) 0.310 0.098 0.266
P_mid (MPa) 0.080 0.021 0.067
plant at 164.34 °C and 3.912 MPa. This is the geothermal fluid used to
P(exhaust) (MPa) 0.016 0.003 0.012
evaporate the organic working fluid that generates power in the ex- Flow Fluid 1st Body Turb1 (kg/s) 66.547 83.012 82.729
pansion at the turbine. For the same flow of hot brine, we compare the Flow Fluid 2nd Body Turb1 (kg/s) 60.223 58.781 59.127
facilities operating with the three fluids. Table 1 shows the main op- T Feed (°C) 121.365 109.449 116.708
erating parameters across the flowsheet for the three fluids considered, Flow Brine (kg/s) 146.579 146.579 146.579
Cooling Water HX5 (kg/s) 6.705 136.481 30.240
benzene, toluene and cyclohexane maximizing the power produced
Cooling Water HX6 (kg/s) 1981.850 1815.069 1951.456
alone. The main operating variables presented include the power gen- Total Risk (yr−1) 9.0858E−4 9.0594E−4 9.0802E−4
erated by each fluid, the two operating pressure levels at the turbine,
the exhaust steam pressure, the feed temperature to the turbine, the
flows of the thermal fluid across the turbine, the brine flow and the Furthermore, the split fraction of the fluid towards the turbine is 60% to
cooling needs at HX5 and HX6. We see the differences in the operating the high-pressure section and the rest to the mid-pressure section for
pressures. The T-S diagram for Toluene is presented in Fig. 2 where the toluene and cyclohexane while benzene shows a split closer to 50%.
dots are computed using the surrogate models for the enthalpy and Among the three fluids, toluene shows the largest power production,
entropy. We see that the surrogate models are validated. The safer fluid 10.436 MW, followed by cyclohexane and benzene. It should be noted
identified in previous work [13], toluene, shows lower operating tem- that this is inversely related to the operating pressures. Furthermore,
peratures and pressures across the flowsheet. This is due to the features the pressure ratio at each expansion is around 4 and 7 in all cases. These
of the liquid that shows a higher volatility and therefore operates at results are not far from typical operating values presented in the lit-
lower pressures so as to better recover the energy within the brine. erature for the fluids under consideration [24].
Fig. 3 shows the T-Q diagram of the heat exchanger network used to
Table 1 recover the energy from the brine in the case of using toluene, the fluid
Main operating results for the three fluids. with better performance. We see that there are two pinch points, with a
ΔT of 5 °C, corresponding to the temperatures of saturated vapor fed to
Variable Benzene Toluene Cyclohexane
the turbine. Thus, the optimization finds the optimal heat exchange
Power (MW) 9.134 10.436 10.280 since both curves cannot came closer without temperature crossing.
P_max (MPa) 0.310 0.102 0.267 Similar results are obtained for the other two fluids.
P_mid (MPa) 0.080 0.025 0.067 By adding the safety issues into the model, Table 2 shows the op-
P(exhaust) (MPa) 0.016 0.003 0.012
Flow Fluid 1st body Turb1 (kg/s) 66.465 80.651 82.479
erating conditions. We see that the flows of organic fluid, as well as the
Flow Fluid 2nd body Turb1 (kg/s) 60.306 56.077 59.267 operating pressures change slightly. Since the power production does
T Feed (°C) 121.417 111.003 116.854 not only play a role in the economic terms of the optimization but also
Flow Brine (kg/s) 146.579 146.579 146.579 in avoiding CO2 emissions, larger CO2 emissions are avoided when
Cooling HX5 (kg/s) 8.701 19.130 30.233
larger power is produced from non fossil based resources, it results in
Cooling HX6 (kg/s) 2209.458 2165.822 1926.735
the fact that power production is similar to the previous case, see

Fig. 3. T-Q Curve for the heat integration leading to geothermal energy recovery.

177
J. Peña-Lamas et al. Energy Conversion and Management 165 (2018) 172–182

Table 3
Safety optimization. 90% of maximum power production.

Variable Benzene Toluene Cyclohexane

Power 8.221 9.392 9.252


P_max (MPa) 0.197 0.037 0.133
P_mid (MPa) 0.080 0.010 0.045
P(exhaust) (MPa) 0.020 0.003 0.012
Flow Fluid 1st Body Turb1 (kg/s) 93.823 127.072 125.444
Flow Fluid 2nd Body Turb1 (kg/s) 34.199 40.076 24.556
T Feed (°C) 103.417 78.385 90.049
Flow Brine (kg/s) 146.579 146.579 146.579
Cooling Water HX5 (kg/s) 3.447 2.508 28.240
Cooling Water HX6 (kg/s) 1963.632 2097.492 2143.452
Total risk (yr−1) 9.0749E−4 9.0525E−4 9.0652E−4

Fig. 5. Effect of total risk on the operating pressures of the turbine.


Table 1. Due to the additional weight associated to mitigating emis-
sions, the power production when using toluene is higher, but within production to 95% of the value, but later on this operating pressure
numerical error. Similarly, the operating conditions are either close to remains constant, see Fig. 5.
the ones reported when only power was maximized or just below, to
reduce the risk. On the other hand, the thermal fluid flows are similar
5.2. Economic evaluation
compared to the ones obtained when only power was maximized.
Furthermore, the safest design among the three fluids is the one that
We divide this section into investment and production cost esti-
uses toluene, the total risk values for the three options are also shown in
mations. Investment cost is computed using the factorial method [22]
Table 2.
that relies on the equipment cost. We consider the example of the
If only the total risk term of the multiobjective problem is opti-
maximum power production for toluene. First, the equipment cost is
mized, the system results in a large decrease of power production. Thus,
estimated by performing a short-cut sizing of the units involved such as
we evaluate, for a maximum acceptable loss in power production, the
heat exchangers, turbine, the cooling pond and wells [26–28]. For the
safest operating conditions. In Table 3, we have assumed 10% loss in
factorial method, we do not include the cost of the well, but we add it
power production to evaluate the safety issues by optimizing plant
afterwards. Equipment cost is estimated using the correlations devel-
safety. We see that, comparing to the results presented in Table 1, the
oped in Martín and Grossmann [27] and updated by Almena and Martín
operating pressures and temperatures decrease resulting in safer de-
[28]. The cooling pond cost estimation correlation can be found in the
signs, see risk values. However, the organic fluid flow to the high
Appendix A below. The installation of the units represents 1.5 times the
pressure turbine increases, while the flow to the second expansion of
equipment cost but for the well, whose cost is considered as installed as
the turbine decreases. Another alternative would be to select an ac-
given by Eq. (14).
ceptable level of risk and determine the power that can be obtained
The cost of the well is estimated by fitting the data from Mansure
from Organic Rankine Cycles.
and Blakenship [29] as a function of the well depth, Eq. (14). The cost is
Therefore, for the case of the toluene, Figs. 4 and 5 show the Pareto
later updated and the currency exchanged into Euros. Three production
curve between power production and total risk and the effect of redu-
wells are built, plus 2 reinjection wells as shown in Fig. 1.
cing the total risk on the operating pressures at the turbine respectively.
Since temperatures and pressures are related due to the vapor–liquid −04 (Depth(m))
Cost(M$) = 0.982355·e 4.128·10 (14)
equilibrium, we do not show the temperature profiles. We see that
power production remains above 95% in a large range of total risk Piping, isolation, instrumentation and utilities represent 20%, 15%,
values. However, total risk shows a steeper descent beyond 90% of the 20% and 10% of the units but the wells. Land and buildings cost is
maximum power production. This is directly related to the decrease in estimated to be 8 M€. We use the same distribution of costs as presented
the operating pressure at the high pressure section of the turbine, see in the work that evaluated a concentrated solar power plant [23,24] for
Fig. 4, where the highest pressure of the organic fluid is the one that the sake of a comparison. The fees represent 3% of the fix cost, other
suffers a larger reduction for a safer design. The exhaust pressure re- administrative expenses and overheads and the plant layout represent
mains constant in all cases, while the pressure at the intermediate feed 10% of the direct costs (fees plus fix capital) and 5% of the fix cost
to the turbine decreases when we go from 100% of maximum power respectively. The plant start-up cost represents 15% of the investment.
Table 3 shows the results for each fluid. Fig. 6 presents the breakdown
of the cost. We see that the wells represent around 45–50% of the
system followed by the heat exchanger network, up to 45% and the
turbine (12–14%). The distribution is similar for all the fluids with
small differences. The rest is the cooling pond.
Furthermore, we estimate the production cost of the electricity. For
the average annual cost, we consider the labour costs (0.5% of invest-
ment), equipment maintenance (2.5% of fix costs), amortization (linear
with time in 20 years), taxes (1% investment), overheads (1% invest-
ment) and administration (5% of labour, equipment maintenance,
amortization, taxes and overheads) [23]. Table 4 summarizes the pro-
duction costs of the three fluids.
Table 5 shows the summary of the economic evaluation when safety
considerations are included in the optimization. No major differences
can be seen.
To compare geothermal with CSP facilities [23,24], for a production
Fig. 4. Pareto curve between power production and total Risk.
plant of the same production capacity, around 10 MW, we scale down

178
J. Peña-Lamas et al. Energy Conversion and Management 165 (2018) 172–182

Fig. 6. Units cost share in geothermal plants. (a) Benzene; (b) Toluene; (c) Cyclohexane.

Table 4 developed unit cost correlations [27,28] so that, for instance, the cost of
Summary of economic parameters of the geothermal plants. a heat exchanger is based on its area that is related to the fluids oper-
ating, assuming average global heat transfer coefficients, and the ΔT
Benzene Toluene Cyclohexane
from the results of the flowsheet optimization. Consequently, the larger
Investment (M€) 103 102 112 the heat load, the larger the area and the cost is recomputed. Therefore,
Production cost (€/kWh) 0.083 0.072 0.081 scaling-up or down a facility only requires computing the flows of en-
€/kW 11,100 9800 10,900
ergy and mass involved proportionally to the base case and, with those,
the individual cost of the units. By increasing the production capacity,
the investment cost increases and it is possible to determine the re-
Table 5
Summary of economic parameters of the geothermal plants including safety considera- lationship between both. Thus, using the factorial method the cost of
tions. the entire facility can be recomputed. Note that, the unit cost is not
proportional to its size in most cases. Moreover, when the units reach a
Benzene Toluene Cyclohexane
certain size, depending on the standard, we cannot longer order a
Investment (M€) 103 102 112 bigger unit, but several units in parallel must be allocated. A linear
Production cost (€/kWh) 0.083 0.072 0.081 relationship is found based on the integer costs for any new well drilled
€/kW 11,100 9800 10,900 and the contribution of the heat exchanger network and their sizes. The
correlation for the investment cost is given by Eq. (15). Note that a
constant depth of 3600 m for all the wells is used for estimating their
the investment and production costs using the correlations developed in
cost.
previous work [23,24]. The results show that the investment cost for
geothermal energy is almost half (188 M€ for CSP), and the investment Investment(M€) = 5.0252·Power(MW) + 58.606 (15)
per kW is 50% lower showing the potential of using this energy versus
CSP. The main difference is basically the cost of the heliostats, which is Scaling-up production costs is somehow easier since most of the
far more expensive than well’s construction. Another advantage is the costs are proportional to the production capacity such as chemicals and
lack of daily variation in the production of power. However, the energy raw materials. Furthermore, the amortization is computed from the
available from the sun is quite large and therefore, an effort towards a scaled-up investment cost computed in the previous stage. In terms of
more economic production of power from solar plants is needed, fur- production cost of geothermal electricity, the correlation does not fit
thermore, the heat source is better since the temperatures of operation the data as good as for the investment cost. Fig. 7 shows the scaled-up
are higher. The production cost for CSP power is around 0.17 €/kWh by of the power production cost and Eq. (16) the correlation. The costs are
scaling down previous works [23,24], resulting in the fact that geo- promising for production capacities larger than 25 MW. When ponds
thermal power can be competitive with CSP. are substituted by cooling towers, at the cost of 0.06 €/kWh, a

5.3. Scale-up

Using the results of the best fluid, toluene, for the case of the op-
timization of the production of power, we evaluate the scale-up costs of
geothermal energy based on the fact that a new well is to be drilled per
48.86 kg/s of brine extracted. Note that over 25 MW power instead of a
pond, a wet cooling tower is used. Using the procedure described in
previous work [30], the effect of the production capacity on the in-
vestment and production costs is evaluated. In short, the investment of
the plant is estimated based on the factorial method that relies on the
units cost. The investment costs of the units are a function of the mass
or energy flow rate involved. While pressure and temperature are in-
tensive variables, constant with the production capacity, mass and flow
Fig. 7. Power production cost scale-up.
rates are proportional to the scale-up or down of the facility. We

179
J. Peña-Lamas et al. Energy Conversion and Management 165 (2018) 172–182

discontinuity can be seen in the curve. alone is optimized, it turns out that toluene is the most efficient fluid
from the economic point of view, capable of recovering energy at
Power production cost(€/kWh) = 0.1084(Power(MW))−0.163 (16)
competitive prices for electricity. The operating conditions of the three
The number of wells required increases linearly with the power fluids are similar, in terms of flows used whereas toluene shows the
production capacity, Eq. (17). Therefore, the availability of energy of a largest power production and lowest cost. In addition to this, it operates
particular reservoir will limit the number of wells drilled and, conse- at lower pressures and similar temperatures compared to the other
quently, the power that can be produced. Geological limitations of the fluids, which is also a safety advantage. If safety and environmental
location must be carefully considered to determine the capacity of a issues are included in the objective function, typically lower operating
particular reservoir. pressures and temperatures are chosen, but the emissions mitigated by
Wells = 0.2868·Power(MW) + 0.1964 (17) the production of power result in that two of the three terms of the
objective function are in favour of maximum power production. No
difference in the fluid selection is suggested.
6. Conclusions Scale-up costs depend on the availability and feasibility of ex-
tracting a larger flow of hot brine. However, if possible, the investment
The operation of a geothermal power plant based on Organic cost scale up linearly and promising production costs for power can be
Ranking Cycles has been optimized considering not only power pro- found over 25 MW.
duction but also environmental and safety objectives. We have con-
sidered the best three fluids determined by a previous study, benzene,
toluene and cyclohexane, and performed detailed process optimization Acknowledgement
modelling the thermodynamic operation of the turbine and the heat
exchanger network. Finally, an economic evaluation to compute in- The authors acknowledge the financial support from the Thematic
vestment and production costs has been also performed, including Network of Bioenergy (REMBIO) and Mexican Council for Science and
scale-up analysis of the facility. Technology (CONACyT) as well as Salamanca Research for Software
The optimization selects a double expansion scheme. When power licenses and MICINN project DPI2015-67341-C2-1-R.

Appendix A

Thermodynamic correlations
S (J/mol·K)

H (kJ/mol)

T (K)

P (bar)
Toluene

Saturated liquid
H = 2.413·10−13T 6−5.986·10−10T 5 + 6.123·10−7T 4−0.0003303T 3 + 0.09922T 2−15.6T + 971.6;

S = 4.583·10−9T 4−6.919·10−6T 3 + 0.003667T 2−0.3103T −129.2;

Saturated steam
H = −3.15523·10−7T 3 + 4.51006·10−4T 2−8.2643·10−2T −16.3670;

S = −8.01·10−11T 5 + 1.683·10−7T 4−0.0001422T 3 + 0.06072T 2−13.06T + 1208;

Overheated steam

Pressure 0.02–1.4 bar


H = 7.797 + 0.0007286T −1.743P + 0.000174T 2 + 0.003393PT −0.01305P 2;
Pressure 0.02 to 0.3 bar
S = −4.367 + 0.3199T −323.9P + 4.004·10−5T 2 + 0.2449PT + 1156P 2−0.0003046T 2P−0.03606TP 2−1774P 3;
Pressure 0.2 to 0.5 bar

S = 4.690·10−6P 3T 2−2.152·10−3P 3T −74.800P 3 + 3.165·10−5P 2T 2−2.482·10−2P 2T + 120.891P 2−1.004·10−4PT 2 + 0.0859PT −96.419P−5.595·10−7T 2


+ 3.5006·10−1T −21.865;
Pressure 0.3 to 1.4 bar

S = 80.06−0.7629T −53.46P + 0.004495T 2−0.2038TP + 119.7P 2−8.21·10−6T 3 + 0.001046T 2P−0.329TP 2−32.58P 3 + 5.706·10−9T 4−1.34·10−6T 3P
+ 0.0004433T 2P 2−0.01034TP 3 + 9.053P 4;

Liquid

180
J. Peña-Lamas et al. Energy Conversion and Management 165 (2018) 172–182

H = href + (140140−152.3·TAV + 0.695·(TAV )2)·10−6·(T −Tref );

With:
kJ T + Tref
href = −14.55425 ; TAV = ; Tref = 298.15K ;
mol 2
S = −228 + 0.1355P + 0.6739T −0.002778P 2−3.308·10−5PT −0.0002137T 2;
Antoine Eq.

3096.52 ⎞
ln (Psat ·750.064) = 16.0137−⎛ ⎜ ;⎟

⎝ Tsat −53.67 ⎠
Benzene

Saturated liquid
H = −5.3293023·10−14T 6 + 1.3528752·10−10T 5−1.4069533·10−7T 4 + 7.6948020·10−5T 3−2.3241987·10−2T 2 + 3.7956329T −2.8223102·102;

S = 2.9008962·10−11T 5−5.6525126·10−8T 4 + 4.4218808·10−5T 3−1.7507080·10−2T 2 + 3.9400717T −4.3576804·102;

Saturated steam
H = −5.9672469·10−12T 5 + 1.0752354·10−8T 4−8.0309120·10−6T 3 + 3.1637092·10−3T 2−5.6945064·10−1T + 56.601956;

S = 1.75589·10−15T 6−3.20075·10−11T 5 + 6.31984·10−8T 4−5.42456·10−5T 3 + 2.46331·10−2T 2−5.79659T + 640.791;

Overheated steam
H = 9.847 + 0.01923T −0.7586P + 0.0001165T 2 + 0.001124TP + 0.01733P 2;
Pressure 0–0.8 bar
S = 6.739 + 0.3097T −54.14P−3.761·10−5T 2 + 0.003373TP + 32.55P 2;
Pressure 0.8–2.0 bar
S = −10.06 + 0.3255T −14.02P−5.284·10−5T 2 + 0.001901TP + 2.361P 2;

Liquid
H = href + (162940−344.94·TAV + 0.85562·(TAV )2)·10−6·(T −Tref );

With:
kJ T + Tref
href = −7.6977 ; TAV = ; Tref = 298.15K ;
mol 2
S = −173.8 + 0.5649T −0.02493P−0.0002062T 2 + 8.083·10−6PT + 0.001357P 2;
Antoine Eq.

2788.51 ⎞
ln (Psat ·750.064) = 15.9008−⎛ ⎜ ; ⎟

⎝ sat −52.36 ⎠
T

Cyclohexane

Saturated liquid
H = 5.64877·10−11·T 5−1.10858·10−7·T 4 + 8.61667·10−5·T 3−0.0329151·T 2 + 6.32556·T −510.124;

S = 1.22384·10−8·T 4−1.92074·10−5·T 3 + 0.0110945·T 2−2.28934·T + 80.0381;

Saturated steam
H = −4.47359·10−13·T 6 + 1.06898·10−9·T 5−1.05657·10−6·T 4 + 0.000552199·T 3−0.160689·T 2 + 24.7614·T −1565.25;

S = −7.58895·10−13·T 6 + 1.80099·10−9·T 5−1.76164·10−6·T 4 + 0.000905964·T 3−0.256920·T 2 + 37.9053·T −2181.58;

Overheated steam
H = 5.066 + 0.007424T −0.8625P + 0.0001808T 2 + 0.001191PT + 0.02255P 2;

S = −19.97 + 0.3682T −94.6P + 7.489·10−6T 2 + 0.05242PT + 126.04P 2−4.627·10−5T 2P−0.02955P 2T −100.4P 3 + 2.569·10−5P 2T 2 + 0.008158TP 3
+ 40.39P 4−6.232·10−6T 2P 3−0.0003521TP 4−6.408P 5;

181
J. Peña-Lamas et al. Energy Conversion and Management 165 (2018) 172–182

Liquid
H = href + (−220600 + 3118.3·TAV −9.4216·(TAV )2 + 0.010687·(TAV )3)·10−6·(T −Tref );

With:
kJ T + Tref
href = −9.3247 ; TAV = ; Tref = 298.15K ;
mol 2
S = −184.9 + 0.5331T + 0.1245P−2.991·10−5T 2−0.0004279PT −0.0008751P 2;
Antoine Eq.

2766.63 ⎞
ln (Psat ·750.064) = 15.7527−⎛ ⎜ ;

⎝ Tsat −50.50 ⎠
Cost correlations

Cost (Cooling pond €) = 11.325 · (Cooling kW) + 6.406.

References Energy 2014;76:175–86.


[15] Manente G, Lazzaretoo A, Bonamico E. Design guidelines for the choice between
single and dual pressure layouts in organic Rankine cycle (ORC) systems. Energy
[1] Barbier E. Geothermal energy technology and current status: an overview. Renew 2017;123:413–31.
Sust Energy Rev 2002;6(1–2):3–65. [16] Fiaschi D, Lifshitz A, Manfrida G, Tempesti D. An innovative ORC power plant
[2] DiPippo R. Geothermal power plants. Principles, applications, case studies and layout for heat and power generation from medium- to low-temperature geothermal
environmental impact. Elsevier; 2008. ISBN: 978-00-8055-476-1. resources. Energy Convers Manage 2014;88:883–93.
[3] Moon H, Zarrouk SJ. Efficiency of geothermal power plants: a worldwide review. [17] Santos-Rodriguez MM, Flores-Tlacuahuac A, Zavala VM. A stochastic optimization
New Zealand geothermal workshop 2012 proceedings 19–21 november 2012 approach for the design of organic fluid mixtures for low-temperature heat re-
Auckland, New Zealand. covery. App Energy 2017;198:145–59.
[4] Borsukiewicz-Gozdur A, Nowak W. Comparative analysis of natural and synthetic [18] Yu H, Eason J, Biegler LT, Feng X. Process integration and superstructure optimi-
refrigerants in application to low temperature Clausius-Rankine cycle. Energy zation of OrganicRankine Cycles (ORCs) with heat exchanger network synthesis.
2007;32:344–52. Comput Chem Eng 2017. http://dx.doi.org/10.1016/j.compchemeng.2017.05.013.
[5] Chen H, Goswami DY, Stefanakos EK. A review of thermodynamic cycles and [19] Peña-Lamas J. Geothermal power plant (In Spanish) MSci Thesis University of
working fluids for the conversion of low grade heat. Renew Sustain Energy Rev Salamanca; 2016.
2010;14:2059–3067. [20] DiPippo R. Geothermal power generation. Developments and innovation.
[6] Sun J, Li W. Operation optimization of an organic Rankine cycle (ORC) heat re- Woodhead Publishing978-00-8100-337-4; 2016.
covery plant. Appl Therm Eng 2011;31:2032–41. [21] Wallas SM. Chemical process equipment – selection and design. Buttherworth
[7] Wang J, Yan Z, Wang M, Ma S, Dai Y. Thermodynamic analysis and optimization of Heinemman; 1990.
an organic Rankine Cycle (Orc) using low grade heat source. Energy [22] Sinnot R, Coulson K, Richardson. Chemical engineering. 3ª ed. Butterworth
2013;49:356–65. Heinemann, Singapore; 1999.
[8] Delgado Torres AM, García Rodríguez L. Analysis and optimization of the low [23] Martín L, Martín M. Optimal year-round operation of a concentrated solar Energy
temperature solar organic Rankin cycle (ORC). Energy Convers Manage plant in the south of Europe. Appl Thermal Eng 2013;59:627–33.
2010;5:2846–56. [24] Martin M. Optimal annual operation of the dry cooling system of a concentrated
[9] Roy JP, Misra A. Parametric optimization and performance analysis of a re- solar energy plant in the south of Spain. Energy 2015;84:774–82.
generative organic Rankine cycle using R 123 for waste heat recovery. Energy [25] Martin M, Martín M. Cooling limitations in power plants: optimal multiperiod de-
2012;39:227–35. sign of natural draft cooling towers. Energy 2017;135:625–36.
[10] Franco A, Villani M. Optimal design of binary cycle power plants for water-domi- [26] Girgin L, Ezgi C. Design and thermodynamic and thermoeconomic analysis of an
nated, medium-temperature geothermal fields. Geothermics 2009;38:379–91. organic Rankine cycle for naval surface ship applications. Energy Convers Manage
[11] Ghasemi H, Paci M, Tizzanini A, Mitsos A. Modeling and optimization of a binary 2017;148:623–34.
geothermal power plant. Energy 2013;50:412–28. [27] Martín M, Grossmann IE. Energy optimization of bioethanol production via gasifi-
[12] Yu H, Feng X, Wang Y, Biegler LT, Eason J. A systematic method to customize an cation of Switchgrass. AIChE J 2011;57(12):3408–28.
efficient organic Rankine cycle (ORC) to recover waste heat in refineries. Appl [28] Almena A, Martín M. Techno-economic analysis of the production of epiclorhidrin
Energy 2016;179:302–15. from glycerol. Ind Eng Chem Res 2016;55(12):3226–38.
[13] Martínez-Gómez J, Peña-Lamas J, Martín M, Ponce–Ortega JM, et al. A multi-ob- [29] Mansure AJ, Blakenship DA. Geothermal well cost analyses; 2008. SAND2008-
jective optimization approach for the selection of working fluids of geothermal 3807C.
facilities: economic, environmental and social aspects. J Environ Manage [30] Martín M. Artificial versus natural reuse of CO2 for DME production: are we any
2017;203(3). closer? Engineering 2017;3:166–70.
[14] Guzovic Z, Raskovic P, Blararic Z. The comparision of a basic and a dual-pressure
ORC (Organic Rankine Cycle): geothermal power plant velika ciglena case study.

182

You might also like