
{Page 122}
School of Chemical Engineering and Materials Science, The University of Oklahoma, Norman, Oklahoma
Presented are the logic and methodology used to develop a computer program for simulation of the geothermal binary cycle energy conversion process using both pure fluids and mixtures as working fluids. The simulation was developed sequentially to three levels of description: (a) thermodynamic cycle, (b) equipment sizing, and (c) capital cost specification. Pertinent background information, formulas, and computation methods are discussed.
Geothermal binary cycles are presently receiving serious consideration for geothermal energy conversion when the resource geothermal fluid occurs in the liquid phase, particularly when the resource temperature is below 500°F and/or the geothermal fluid contains high concentrations of noncondensible gases. By definition, in a geothermal binary cycle, shown schematically in Figure 1, the working fluid in the power production cycle (e.g., Rankinetype cycle) receives energy by heat transfer with the geothermal fluid. To satisfactorily evaluate candidate working fluids for use in such cycles and to consider the effects of varying operating conditions on resource utilization for alternative cycles, a computer program capable of working fluid properties prediction, equipment sizing, and capital cost estimation has been developed. A summary of the methodology used in the simulator is presented herein.
The working fluid thermodynamic cycle simulation program, referred to as GEO 1, was developed by writing an executive program which utilizes subroutines from a previously developed thermodynamic behavior prediction program (2), referred to herein as HSGC. Because the HSGC computer program listing and summary documentation are available in the literature (2), only the correlation basis and computation capabilities are presented herein. The correlation basis is a generalized modified BWR (BenedictWebbRubin) equation of state correlation. Properties calculated using HSGC include liquid and vapor density, enthalpy, entropy, heat capacity, and component fugacity. Unit computations which can be performed using HSGC include flash, dew point, bubble point, isentropic, and isenthalpic search calculations.
The executive program of GEO 1 performs a set of sequenced thermodynamic calculations, given certain specified input data, in order to simulate the sequence of states undergone by the working fluid in the thermodynamic cycle. The input data to GEO 1 are (a) inlet and exit cooling water temperatures, (b) inlet brine temperature, (c) exit brine temperature (estimated), (d) minimum approach temperatures (minimum temperature differences between fluids in heat exchange) at the exchanger entrances and exits, (e) turbine and pump efficiencies, (f) average heat capacities for cooling water and brine, (g) working fluid temperature and pressure changes in pipe flow between equipment units, (h) temperature interval for energy balances through the brine heat exchanger (also referred to as the evaporator), (i) turbine inlet working fluid pressure, (j) working fluid dew point pressure (estimated) at the condenser inlet, (k) minimum allowable vapor mole fraction at the turbine outlet, and (l) parameters for the thermodynamic properties calculations: convergence criteria, number of components in the working fluid, printout parameters, working fluid component physical and thermodynamic parameters, and working fluid compositions.
The sequence of calculations in GEO 1 is as follows: 1. The brine heat exchanger working fluid exit temperature is determined from the brine inlet temperature
{Page 123}
and the exchanger minimum approach temperature. Countercurrent heat exchange is considered. A state point calculation is performed for the working fluid at the exit temperature and pressure. The turbine inlet temperature and pressure are then calculated using the input temperature and pressure changes in pipe flow between the brine heat exchanger and turbine inlet. A state point calculation is then performed to determine the entropy and other working fluid thermodynamic properties at the turbine inlet. 2. The initial estimate of the working fluid pressure at the condenser inlet is chosen to be the dew point pressure at the minimum possible temperature consistent with the input condenser minimum approach temperature constraint. 3. The bubble point pressure is then calculated at the condenser outlet at the temperature consistent with the condenser outlet minimum approach temperature constraint. 4. The larger of the pressures calculated in steps 2 and 3 above is then selected as the condenser operating pressure with an adjustment for the input condenser working fluid pressure difference between the inlet and outlet. Further adjustment is made for the input pressure or temperature difference between the turbine exit and condenser inlet. The adjusted pressure then becomes the turbine exit pressure. 5. The turbine work and outlet temperature are calculated by first performing an isentropic expansion from the turbine inlet temperature and pressure to the turbine exit pressure, then calculating the work and actual enthalpy change as the input turbine efficiency times the isentropic enthalpy change, followed by an isenthalpic calculation at the turbine exit enthalpy and pressure to determine the turbine outlet temperature. 6. The working fluid temperature at the condenser inlet is then updated using the turbine exit temperature and input temperature difference between turbine exit and condenser inlet; if the minimum approach temperature criterion is violated the turbine exit pressure is increased by a specified (input) increment (e.g., 1 psia) and step 5 is repeated. 7. After the working fluid pressure at the condenser exit is calculated, the working fluid temperature at the condenser exit is calculated to be the bubble point temperature at this exit pressure, and if the minimum approach temperature criterion is violated, the turbine exit pressure is increased and step 5 is repeated. 8. When the turbine exit presure has been finalized, the condenser exit temperature is chosen to represent 0.2° F subcooling (to avoid cavitation in the cycle pump) and the condenser average log mean temperature difference is calculated. 9. The cycle pump inlet temperature and pressure are calculated using the input temperature and pressure differences between the condenser exit and pump inlet. The cycle pump outlet pressure is calculated by adjusting the brine heat exchanger exit pressure with the input temperature and pressure differences between the cycle pump exit and evaporator inlet and the input pressure difference between the brine heat exchanger inlet and exit. 10. The pump work and outlet temperature are calculated by first performing an isentropic compression from the pump inlet temperature and pressure to the pump outlet pressure, then calculating the work and actual enthalpy change as the isentropic enthalpy change divided by the pump efficiency, followed by an isenthalpic calculation at the pump exit enthalpy and pressure to determine the pump exit temperature. 11. Adjustments, if specified by the input data, are made to the cycle pump exit temperature and pressure. 12. An overall energy balance for the brine heat exchanger using the input estimate of the exit brine temperature yields an initial estimate of the ratio of the brine flow rate to the working fluid flow rate. The input value of the heat capacity for the brine is assumed to be constant throughout the brine temperature range in the brine heat exchanger. The initial brine outlet temperature is taken as the working fluid inlet temperature added to the minimum approach temperature difference. The brine heat exchanger is then divided into sections according to the input working fluid temperature interval. The working fluid thermodynamic properties are calculated at each temperature interval. 13. The brine temperature at each interval is then calculated from an energy balance for each brine heat exchanger section using the brine heat capacity, the assumed brine flow rate, and the working fluid properties at that section. 14. If the brine temperature at any section is less than the "pinch point" temperature for that section, then the brine flow rate is increased and step 13 is repeated. 15. If the minimum approach temperature constraint
{Page 124}
is not violated, the log mean temperature difference (LMTD) is calculated for each section. The average LMTD for the brine heat exchanger is determined by dividing the total exchanger working fluid enthalpy change by the sum of the ratios of the section enthalpy changes to the section LMTD's.
The output data from GEO 1 include detailed thermodynamic results in molar and mass units for each state point in the cycle, and if desired, the temperatureenthalpyentropyLMTD profile for the working fluid in the brine heat exchanger, the brine temperature profile in the evaporator, the heat added to the working fluid, the brine heat exchanger average LMTD, the heat rejected to the cooling water, the condenser average LMTD, the flow rate ratio of brine to working fluid, the flow rate ratio of cooling water to working fluid, the gross turbine work, the net thermodynamic cycle work, the Rankine cycle efficiency, and the net thermodynamic cycle work per pound of brine. The output from GEO 1 allows plotting, as in Figure 2 for the working fluid isobutane, the thermodynamic states for the working fluid on the temperatureentropy diagram. The temperatures (but not entropies) of brine and cooling water corresponding to the working fluid temperature also are shown in Figure 2. Table 1 shows sample output of working fluid thermodynamic properties at cycle state points shown in Figure 1.
The GEO 1 simulator was tested by comparisons with separate thermodynamic calculations and hand calculations at numerous stages in its development. It has been demonstrated that the simulator can accurately calculate the properties for a thermodynamic Rankine cycle using either pure components or mixtures as working fluids in the binary cycle.
The major purposes of equipment sizing calculations in the geothermal binary cycle simulator are (a) to provide an overall conceptual design for the power plant in detail sufficient for equipment parameter and operating conditions selection and (b) to provide sufficient information for cost estimation. The levels of design calculations in the simulator vary from unit to unit, depending on the levels of detail required for equipment parameter and operating conditions selection. Brief summaries of pertinent formulas used in sizing major equipment units are presented below to provide the minimum information necessary for cost estimation, with references to appropriate literature sources for design calculations performed in greater detail within the simulator. More detailed equipment sizing information appears in a report (3) available through NTIS (National Technical Information Service.)
{Page 125}
Heat Exchanger Size
Both the brine heat exchanger and the condenser are considered to be singlepass carbon steel shellandtube triangularpitch heat exchangers which approach countercurrent flow behavior with brine and cooling water on the tube side and working fluid on the shell side. The design methods presented by Kern (4) constitute the basis for the heat exchanger design calculations.
The minimum items of information needed to estimate for heat exchanger cost are the tube and shell side pressures (determined by GEO 1 ) and the heat transfer surface area, A,
Eq. 1 
The heat exchanger duty, Q, and the log mean temperature difference, LMTD, are determined by GEO 1, so that only the calculation of the overall heat transfer coefficient, U, is required for determination of the heat transfer surface area, A, from Equation 1. The overall heat transfer coefficient is calculated using the wellknown equation (4)
Eq. 2 
In Equation 2, h, R, and d are, respectively, convective heat transfer coefficient, fouling resistance and diameter, the subscripts i and o refer to inside and outside of tubes, and k_{t} is the tube thermal conductivity; all are simulator inputs except h_{i} and h_{o}. Convective heat transfer coefficients for the singlephase region are calculated using the SiederTate correlation (5),
Eq. 3 
In Equation 3, d is the equivalent diameter for flow, i.e., square root of (four times the area normal to flow divided by p), k, µ, and Cr are, respectively, the average free stream values of the thermal conductivity, viscosity, and heat capacity of the fluid (working fluid, brine or cooling water), µ_{w} is the fluid viscosity at the tube wall, and G is the fluid mass velocity. For boiling heat transfer, Chen's modification (6) of Equation 3 is used for calculation of the heat transfer coefficient. Condensing coefficients were computed using a correlation (7) for film type condensation
Eq. 4 
In Equation 4, k, r, and µ are, respectively, the liquid phase thermal conductivity, density, and viscosity, H_{fg} is the enthalpy of condensation, L is equal to onehalf of the tube outside perimeter, T_{w} is the vaporliquid interface temperature, T_{w} is the tube wall temperature, and g_{c} is the gravitational constant.
Because local heat transfer coefficients vary significantly through the heat exchangers, the heat exchangers are treated by sections and the heat transfer surface area is calculated section by section and then summed to determine the total heat transfer surface area.
Parasitic power requirements were obtained from the pressure drops computed for the tubeside and the shellside flows in the heat exchangers.
{Page 126}
Turbine Size
The simulator utilizes axial flow turbine sizing formulas (8) based on specific diameter, D_{sp}, and specific speed, N_{sp},
Eq. 5 
Eq. 6 
In Equations 5 and 6, V_{t} is the turbine wheel tip speed, C_{o} is the spouting velocity, C_{o} = 223Ö¯¯¯Dh_{i}, D h_{i} is the isentropic specific enthalpy drop, q is the volumetric flow rate, and N_{t} is the turbine angular velocity (i.e., RPM). The high efficiency region for axial turbines occurs in the range of specific speeds from 80 to 120, where the following relation can be used,
Eq. 7 
Thus, a value of N_{sp} is selected, Equation 7 is used to calculate V_{t} / C_{o}, and then Equation 5 is used to calculate D_{sp}. The turbine diameter, D_{t}, is then calculated using the relation,
Eq. 8 
Aside from the sizing of heat exchangers, turbine, and major piping, detailed calculations of this type were not performed. Pump sizes were not determined in the simulator because centrifugal pumps are offtheshelf items with cost specified by power rating. A similar situation exists for generators. Auxiliary equipment sizing was not performed because this cost can be estimated from the major equipment cost. Although the well system was not simulated, well costs were estimated in order to arrive at total capital investment.
The level of simulation of the geothermal binary cycle obtained by incorporating the equipment sizing calculations into the working fluid thermodynamic cycle simulator (GEO 1) will be referred to herein as GEO 2. The input to GEO 2 includes specification of equipment parameters such as heat exchanger tube size, pitch and thermal conductivity and fouling factors, as well as brine and cooling water velocities, in addition to the parameters input to GEO 1. Tables 2, 3, and 4 show simulator sample output summarizing equipment sizing cal
{Page 127}
calulations, parasitic power requirements, and net plant efficiencies for a 25 MWe isobutane geothermal binary cycle power plant. The thermodynamic states and properties in Figure 2 and Table 1 correspond to this cycle.
The cost estimation method used by the simulator is based primarily on the method used by Milora and Testor (9). This method requires the determination of the delivered equipment cost and geothermal well cost. The other items included in the total direct plant cost are then estimated as percentages of the delivered equipment and the well costs. The additional components of the capital investment are based on average percentages of the total direct plant cost, total direct and indirect plant costs, or total capital investment. This can be summarized in the following cost equation,
Eq. 9 
where C_{n} = net capital investment, S_{i}C_{Ei} = total major equipment, i.e., heat exchangers, condensers, turbine and condensate pump cost, S_{i}f_{i} = multiplying factors for piping, buildings and structures, instrumentation, equipment installation, etc., f_{I} = indirect cost factors, e.g., engineering fees, contingency and related costs, C_{w} = cost of geothermal production and/or reinjection well, n_{w} = number of geothermal wells required for a particular size power plant, f_{w} = fraction of C_{w} which accounts for the piping from the wellhead to the power plant, and f*_{i} = indirect cost factors, e.g., costs associated with discovery of geothermal field such as land acquisition, drilling exploratory holes, and contingencies. Table 5 outlines the numerical values of the direct multiplying factors (fractions of the total equipment and well costs) that were used. Table 6 presents the fractions of the total equipment investment indirectly associated with the capital investment.
In order to calculate the total cost of a completed geothermal power plant, the annual costs associated with the capital investment for plant construction (e.g., interest charged on borrowed money, taxes, depreciation, return on equity, and insurance, etc.) must be evaluated. Annual operating, maintenance and repair, and general expenses (administrative and research, etc.), and electrical transmission costs (from
{Page 128}
the power plant busbar to the regional distribution network), should be added to the total capital investment to compute the cost per kilowatthour.
Completed Geothermal Well Costs
The cost of the geothermal well is essentially dependent on the depth and diameter of the well and the type of rock. The depth of the well in turn depends on the thermal gradient (°F/1000 feet) of the geothermal site. A complete description of the well cost model has been provided by Tester and Milora (9). The following equation is a linear approximation to their correlation,
Eq. 10 
where W_{D} is well depth in feet.
Equation 10 is applicable only to softer rock (liquiddominated and/or hot dry rock systems). The performance and lifetime of a naturally or artificially stimulated reservoir is perhaps the most difficult aspect of a geothermal system to predict, since the detailed character of the aquifer is usually unknown. A well lifetime of 20 years or more has been assumed for the example calculation herein. For the liquiddominated aquifers the well flow rates range from 60 to 110 lb./sec. An equal number of production and reinjection wells has been assumed. The well cost equation and subsequent equations are based on 1976 dollars.
Cost of Heat Exchangers and Condensers
The brine heat exchanger and condenser costs can be estimated using the following equations,
(i) Tubeside pressure = 200300 psia
Eq. 11 
where P_{shell} = shellside pressure, psia
(ii) Tubeside pressure = 1000 psia
Eq. 12 
(iii) Tubeside Pressure = 2000 psia
Eq. 13 
Equations 11, 12, and 13 are valid only for carbon steel as the material of construction for both tubes and shells, and are based on 1976 dollars. The applicable heat exchanger sizes range from 20,000 ft^{2} to 35,000 ft^{2} surface area.
Cost of Turbine
The turbine cost estimate is based on a model developed by the BarberNichols Company (10),
{Page 129}
Eq. 14 
where N_{e} = number of exhaust ends, n_{s} = number of internal stages (Pr/stage 0.7), D_{T} = laststage pitch diameter (ft), h/D_{T} = laststage blade height to pitch diameter,f_{u} = cost multiplier for tip speed V_{T} (ft/sec), and f_{p} = cost multiplier for inlet pressure. The multiplier term (1.04 N_{e} 0.04 N_{e}^{2}) accounts for multiple exhaust ends on a common shaft and reflects the economy of scale associated with the elimination of some main shaftseals and bearings. It is considered to be valid for one to four exhaust ends. If multiple units with single exhaust ends are utilized this term is changed to N_{e}. The factor f_{p} corrects for pressure effects. The correlation for fp is shown in Equation 15. The (2485.8 n_{s}f_{u}D_{T}^{2.1}) term represents the cost of turbine stages. The f_{u} multiplier accounts for the increased costs for high tip speed turbines. The correlation for f_{u} has been presented in Equation 16. The (474.94 D_{T}^{2}) term covers the precision components which include main shaft seals, radial and thrust bearings, and turbine shaft. This equation is considered to be valid for h/D_{T} values up to 0.11. The cost multiplier f_{p} which corrects for the inlet pressure effects can be estimated from the following equation,
Eq. 15 
where P_{max} = maximum turbine pressure (psia). It should be noted that there are no temperature effects below 800°F. The factor which corrects for the blade tip speed, V_{T}, can be calculated from the following equation,
Eq. 16 
where V_{T} is expressed in ft/sec.
Generator Cost
The generator cost is estimated from the following equation,
Eq. 17 
where MW_{e} = net electric output of the unit in megawatts. Equation 17 is applicable for power levels of 1 MW_{e} to 100 MW_{e} .
Condensate Pump Cost
The condensate pump cost can be estimated from the following equation,
Eq. 18 
where MW_{e} = condensate pump power rating in megawatts. This equation is valid only for multistage centrifugal pumps and does not include the cost of the drive mechanism.
The simulator program including capital cost estimation, equipment sizing, and working fluid thermodynamic cycle simulation will be referred to as GEO 3. Table 7 shows GEO 3 sample output summarizing equipment cost, well cost, direct and indirect costs, and total capital investment for the 25 MW_{e} isobutane geothermal binary cycle power plant referred to previously in Figure 2 and Tables 1, 2, 3, and 4.
A description has been presented of the methodology used to develop a geothermal binary cycle simulator applicable to both pure fluid and mixture cycles. Other energy conversion simulators can be devel
{Page 130}
oped along similar lines, including simulation of geothermal, ocean thermal, bottoming and waste heat processes. The procedure is to develop the simulator at three successive levels: (a) thermodynamic cycle, (b) equipment sizing and (c) economics. The resultant simulator can be employed for many purposes, including use as a conceptual design tool, in selection of working fluid, for equipment parameter selection (e.g., heat exchanger tube diameter and pitch), in process operating conditions selection, in sensitivity analysis, and as a subprogram in process optimization.
This work was supported by the University of Oklahoma and the Energy Research and Development Administration, Contract No. E (401) 4944.
1. K. Z. IQBAL, L. W. FISH, and K. E. STARLING, Proc. Okla. Acad. Sci. 56: 1103 (1976).
2. K. E. STARLING, Fluid Thermodynamic Properties for Light Petroleum Systems, Gulf Publishing Co., Houston, Texas, 1973.
3. K. E. STARLING, L. W. FISH, K. Z. IQBAL, and D. YIEH, Report ORO49443. Resource Utilization Efficiency Improvement of Geothermal Binary Cycles  Phase I, prepared for the Energy Research and Development Administration, December 15, 1975 (available from the National Technical Information Service).
4. D. Q. KERN, Process Heat Transfer, McGrawHill Book Co., New York, N.Y., 1950.
5. J. H. PERRY and C. H. CHILTON, Chemical Engineer's Handbook, 5th ed., McGrawHill Book Co., New York, 1973.
6. L. S. TONG, Boiling Heat Transfer and TwoPhase Flow, John Wiley and Sons, New York, N.Y., 1965.
7. J. H. ANDERSON, Report No. NSF/RANN/ SE/GI34979/73/6, University of Massachusetts, Amherst, May, 1973.
8. J. H. ANDERSON, Report No. NSF/RANN/ SE/GI34979/TR/73/7, University of Massachusetts, Amherst, June, 1973.
9. S. L. MILORA and J. W. TESTER, Geothermal Energy as a Source of Electric Power, MIT Press, Cambridge, Massachusetts, 1976.