diff --git a/setup.py b/setup.py index e6989a316..b7a75997d 100755 --- a/setup.py +++ b/setup.py @@ -87,6 +87,9 @@ def read(*names, **kwargs): 'pylocker', 'nrel-pysam', 'tabulate', + 'tqdm', + #'sklearn', + 'scikit-learn', ], extras_require={ # eg: diff --git a/src/fgem/__init__.py b/src/fgem/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/src/fgem/data/powerplants/HighEnthalpyCLGWGPowerPlant.pkl b/src/fgem/data/powerplants/HighEnthalpyCLGWGPowerPlant.pkl new file mode 100755 index 000000000..7fea133f5 Binary files /dev/null and b/src/fgem/data/powerplants/HighEnthalpyCLGWGPowerPlant.pkl differ diff --git a/src/fgem/data/powerplants/NRELORCPowerPlant_30ckm.csv b/src/fgem/data/powerplants/NRELORCPowerPlant_30ckm.csv new file mode 100755 index 000000000..1c0e24fe1 --- /dev/null +++ b/src/fgem/data/powerplants/NRELORCPowerPlant_30ckm.csv @@ -0,0 +1,25 @@ +T_amb [deg C],T [deg C],m_turbine [kg/s],Net Output [MWh],T_inj [deg C],power_output [MWh/kg] +15,125,80,2.25,60,7.8125E-06 +10,125,80,2.34,60,0.000008125 +5,125,80,2.38,60,8.26389E-06 +0,125,80,2.4,60,8.33333E-06 +20,125,80,2.01,61.937,6.97917E-06 +25,125,80,1.76,64.648,6.11111E-06 +30,125,80,1.52,66.366,5.27778E-06 +15,130,80,2.47,60,8.57639E-06 +15,135,80,2.65,61.28,9.20139E-06 +15,140,80,2.81,62.288,9.75694E-06 +15,145,80,2.96,63.224,1.02778E-05 +15,120,80,2,59.64,6.94444E-06 +15,115,80,1.77,59.354,6.14583E-06 +15,110,80,1.54,59.166,5.34722E-06 +15,105,80,1.31,58.981,4.54861E-06 +15,100,80,1.06,56.704,3.68056E-06 +15,125,90,2.42,64.087,7.46914E-06 +15,125,100,2.56,67.754,7.11111E-06 +15,125,110,2.65,71.361,6.69192E-06 +15,125,120,2.74,74.997,6.34259E-06 +15,125,70,2.04,55.19,8.09524E-06 +15,125,60,1.82,49.334,8.42593E-06 +15,125,50,1.58,42.461,8.77778E-06 +15,125,40,1.26,35.186,0.00000875 diff --git a/src/fgem/data/powerplants/NRELORCPowerPlant_30ckm.pkl b/src/fgem/data/powerplants/NRELORCPowerPlant_30ckm.pkl new file mode 100755 index 000000000..814d76ffb Binary files /dev/null and b/src/fgem/data/powerplants/NRELORCPowerPlant_30ckm.pkl differ diff --git a/src/fgem/data/powerplants/NRELORCPowerPlant_60ckm.csv b/src/fgem/data/powerplants/NRELORCPowerPlant_60ckm.csv new file mode 100755 index 000000000..c07aa9b8b --- /dev/null +++ b/src/fgem/data/powerplants/NRELORCPowerPlant_60ckm.csv @@ -0,0 +1,25 @@ +T_amb [deg C],T [deg C],m_turbine [kg/s],Net Output [MWh],T_inj [deg C],power_output [MWh/kg] +15,210,80,8.6,60,2.98611E-05 +10,210,80,8.8,59.954,3.05556E-05 +5,210,80,8.9,58.242,3.09028E-05 +20,210,80,8.1,64.186,0.000028125 +25,210,80,7.5,67.919,2.60417E-05 +30,210,80,7.1,70.003,2.46528E-05 +15,220,80,9.3,66.737,3.22917E-05 +15,230,80,10,70.524,3.47222E-05 +15,240,80,10.6,77.805,3.68056E-05 +15,250,80,11.1,85.385,3.85417E-05 +15,260,80,11.6,93.226,4.02778E-05 +15,270,80,12.1,101.31,4.20139E-05 +15,280,80,12.6,109.66,0.00004375 +15,200,80,7.9,56.87,2.74306E-05 +15,190,80,7.1,53.809,2.46528E-05 +15,180,80,6.2,51.276,2.15278E-05 +15,210,90,9.4,70.479,2.90123E-05 +15,210,100,10,81.915,2.77778E-05 +15,210,110,10.4,91.677,2.62626E-05 +15,210,120,10.8,100.07,0.000025 +15,210,70,7.6,51.267,3.01587E-05 +15,210,60,6.5,40.629,3.00926E-05 +15,210,50,5.2,35.239,2.88889E-05 +15,210,40,3.9,42.934,2.70833E-05 diff --git a/src/fgem/data/powerplants/NRELORCPowerPlant_60ckm.pkl b/src/fgem/data/powerplants/NRELORCPowerPlant_60ckm.pkl new file mode 100755 index 000000000..dea73d665 Binary files /dev/null and b/src/fgem/data/powerplants/NRELORCPowerPlant_60ckm.pkl differ diff --git a/src/fgem/data/powerplants/NRELORCPowerPlant_normalized.csv b/src/fgem/data/powerplants/NRELORCPowerPlant_normalized.csv new file mode 100755 index 000000000..e653992da --- /dev/null +++ b/src/fgem/data/powerplants/NRELORCPowerPlant_normalized.csv @@ -0,0 +1,49 @@ +,T_amb [deg C],T [deg C],m_turbine [kg/s],Net Output [MWh],T_inj [deg C],power_output [MWh/kg] +0,0,125,80,1.0666666666666667,60.0,3.7037022222222225e-06 +1,5,125,80,1.0577777777777777,60.0,3.67284e-06 +2,5,210,80,1.0348837209302326,58.242,3.5933488372093023e-06 +3,10,125,80,1.04,60.0,3.6111111111111106e-06 +4,10,210,80,1.0232558139534884,59.954,3.552976744186047e-06 +5,15,125,80,1.0,60.0,3.4722222222222224e-06 +6,15,220,80,1.0813953488372094,66.737,3.7548488372093027e-06 +7,15,230,80,1.1627906976744187,70.524,4.037465116279069e-06 +8,15,240,80,1.2325581395348837,77.805,4.279720930232558e-06 +9,15,250,80,1.2906976744186047,85.385,4.4815930232558135e-06 +10,15,260,80,1.3488372093023255,93.226,4.6834651162790705e-06 +11,15,270,80,1.4069767441860466,101.31,4.885337209302326e-06 +12,15,280,80,1.4651162790697674,109.66,5.087209302325582e-06 +13,15,190,80,0.8255813953488372,53.809,2.8666046511627907e-06 +14,15,180,80,0.7209302325581396,51.276,2.503232558139535e-06 +15,15,210,90,1.0930232558139537,70.479,3.3735232558139535e-06 +16,15,210,100,1.1627906976744187,81.915,3.229976744186047e-06 +17,15,210,110,1.2093023255813955,91.677,3.053790697674419e-06 +18,15,210,120,1.2558139534883723,100.07,2.9069767441860468e-06 +19,15,210,70,0.8837209302325582,51.267,3.506825581395349e-06 +20,15,210,60,0.7558139534883721,40.629,3.499139534883721e-06 +21,15,200,80,0.9186046511627908,56.87,3.189604651162791e-06 +22,15,210,80,1.0,60.0,3.4722209302325583e-06 +23,15,125,40,0.56,35.186,3.888888888888889e-06 +24,15,125,50,0.7022222222222223,42.461,3.901235555555556e-06 +25,15,130,80,1.097777777777778,60.0,3.811728888888889e-06 +26,15,135,80,1.1777777777777778,61.28,4.089506666666667e-06 +27,15,140,80,1.248888888888889,62.288,4.336417777777778e-06 +28,15,145,80,1.3155555555555556,63.224,4.567911111111111e-06 +29,15,120,80,0.8888888888888888,59.64,3.086417777777778e-06 +30,15,115,80,0.7866666666666666,59.354,2.73148e-06 +31,15,110,80,0.6844444444444444,59.166,2.376542222222222e-06 +32,15,210,50,0.6046511627906977,35.239,3.359174418604651e-06 +33,15,100,80,0.47111111111111115,56.704,1.6358044444444444e-06 +34,15,105,80,0.5822222222222222,58.981,2.0216044444444444e-06 +35,15,125,100,1.1377777777777778,67.754,3.1604933333333336e-06 +36,15,125,110,1.1777777777777778,71.361,2.9741866666666666e-06 +37,15,125,120,1.2177777777777778,74.997,2.818928888888889e-06 +38,15,125,70,0.9066666666666667,55.19,3.5978844444444443e-06 +39,15,125,60,0.8088888888888889,49.334,3.744857777777778e-06 +40,15,125,90,1.0755555555555556,64.087,3.319617777777778e-06 +41,15,210,40,0.4534883720930233,42.934,3.149220930232558e-06 +42,20,210,80,0.9418604651162791,64.186,3.2703488372093026e-06 +43,20,125,80,0.8933333333333332,61.937,3.1018533333333332e-06 +44,25,210,80,0.872093023255814,67.919,3.0281046511627907e-06 +45,25,125,80,0.7822222222222223,64.648,2.716048888888889e-06 +46,30,210,80,0.8255813953488372,70.003,2.8666046511627907e-06 +47,30,125,80,0.6755555555555556,66.366,2.34568e-06 diff --git a/src/fgem/data/powerplants/NRELORCPowerPlant_normalized.pkl b/src/fgem/data/powerplants/NRELORCPowerPlant_normalized.pkl new file mode 100755 index 000000000..030a20e83 Binary files /dev/null and b/src/fgem/data/powerplants/NRELORCPowerPlant_normalized.pkl differ diff --git a/src/fgem/data/powerplants/ORCPowerPlant.pkl b/src/fgem/data/powerplants/ORCPowerPlant.pkl new file mode 100755 index 000000000..b18826aaf Binary files /dev/null and b/src/fgem/data/powerplants/ORCPowerPlant.pkl differ diff --git a/src/fgem/powerplant.py b/src/fgem/powerplant.py new file mode 100644 index 000000000..166981e9c --- /dev/null +++ b/src/fgem/powerplant.py @@ -0,0 +1,657 @@ +import math +import os +import pickle +from datetime import timedelta +from pathlib import Path + +import numpy as np +from scipy.spatial import cKDTree + +parent_path = Path(__file__).parent + + +class _Utils: + @staticmethod + def heatcapacitywater(Twater): + """Computationally efficient correlation for water specific heat capacity based on the GEOPHIRES. + + Args: + Twater (Union[ndarray, float]): water temperature in deg C + + Returns: + Union[ndarray, float]: water specific heat capacity in J/kg-K + """ + + # Based on GEOPHIRES correlations (more stable XSTEAM) + Twater = (Twater + 273.15) / 1000 + A = -203.6060 + B = 1523.290 + C = -3196.413 + D = 2474.455 + E = 3.855326 + cpwater = (A + B * Twater + C * Twater**2 + D * Twater**3 + E / (Twater**2)) / 18.02 * 1000 + return cpwater + + @staticmethod + def nonzero(x, thresh=1e-6): + """Ensure input is nonzero + + Args: + x (ndarray): input numbers + thresh (float, optional): threshold below which numbers are set to zero. Defaults to 1E-6. + + Returns: + ndarray: nonzero version of input array + """ + return np.maximum(thresh, x) + + +class CustomUnpickler(pickle.Unpickler): + + def find_class(self, module, name): + try: + return super().find_class(__name__, name) + except AttributeError: + return super().find_class(module, name) + + +class BasePowerPlant: + """Base class for defining a power plant.""" + + def __init__(self, ppc, Tres, cf, timestep=timedelta(hours=1)): + """Define attributes for the BasePowerPlant class. + + Args: + ppc (float): power plant nameplate capacity in MW. + Tres (float): reservoir temperature in deg C. + cf (float): capacity factor. + timestep (datetime, optional): simulation timestep size. Defaults to timedelta(hours=1). + """ + self.powerplant_capacity = ppc + self.cf = cf + + # default starting point + self.power_output_MWh_kg = 1e-5 + self.power_output_MWe = ppc + self.power_generation_MWh = ppc * 1 + self.T_mix = Tres + self.T_inj = Tres / 2 + self.timestep = timestep + + def compute_geofluid_consumption(self, T, T_amb, m_turbine): + """Compute geofluid consumption. + + Args: + T (float): power plant inlet temperature in deg C. + T_amb (float): ambient temperautre in deg C. + + Raises: + NotImplementedError: you must implemented this. + """ + raise NotImplementedError + + def compute_injection_temp(self, T, T_amb, m_turbine): + """Compute injection temperature. + + Args: + T (float): power plant inlet temperature in deg C. + T_amb (float): ambient temperautre in deg C. + + Raises: + NotImplementedError: you must implement this. + """ + raise NotImplementedError + + def compute_cplant(self, MaxProducedTemperature, min_cost=0): + """Compute power plant capex. + + Args: + MaxProducedTemperature (float): maximum produced geofluid temperature in deg C. + + Raises: + NotImplementedError: you must implement this. + """ + raise NotImplementedError + + def compute_thermalexergy(self, T, T_amb): + """Compute power plant thermal exergy for an inflowing geofluid at a given ambient temperature. + + Args: + T (float): power plant inlet temperature in deg C. + T_amb (float): ambient temperature in deg C. + + Returns: + float: thermal exergy in MJ/kg. + """ + + A = 4.041650 + B = -1.204e-2 + C = 1.60500e-5 + T_amb_k = T_amb + 273.15 # deg K + T_k = T + 273.15 # deg K + thermal_exergy = ( + ( + (A - B * T_amb_k) * (T_k - T_amb_k) + + (B - C * T_amb_k) / 2.0 * (T_k**2 - T_amb_k**2) + + C / 3.0 * (T_k**3 - T_amb_k**3) + - A * T_amb_k * np.log(T_k / T_amb_k) + ) + * 2.2046 + / 947.83 + ) # MJ/kg + + return thermal_exergy + + def compute_power_output(self, m_turbine): + return self.cf * min(self.powerplant_capacity, self.power_output_MWh_kg * (m_turbine * 3600)) + + def step( + self, + m_turbine, + T_prd_wh, + T_amb, + m_tes_out=0, + T_tes_out=100, + ): + """Step and update power plant. + + Args: + m_turbine (float): total mass flowing into the turbine in kg/s (including flow from both producers and thermal energy storage tank). + T_prd_wh (float): producer wellheat temperature in deg C. + T_amb (float): ambient temperature in deg C. + m_tes_out (float, optional): mass flow rate from the thermal energy storage tank in kg/s. Defaults to 0. + T_tes_out (float, optional): temperature of water of the thermal energy storage tank. Defaults to 100. + + """ + + # If we are using all geofluid to charge the tank, then there is no output + if m_turbine == 0: + self.power_output_MWe = 0 + self.power_generation_MWh = 0 + self.T_inj = self.compute_injection_temp(self.T_mix, T_amb, m_turbine) + return 0.0, 0.0, 0.0, 0.0, self.T_inj + else: + if m_tes_out > 0: + m_wh_to_turbine = m_turbine - m_tes_out + self.T_mix = ( + m_wh_to_turbine * _Utils.heatcapacitywater(T_prd_wh) * T_prd_wh + + m_tes_out * _Utils.heatcapacitywater(T_tes_out) * T_tes_out + ) / ( + m_wh_to_turbine * _Utils.heatcapacitywater(T_prd_wh) + + m_tes_out * _Utils.heatcapacitywater(T_tes_out) + + 1e-3 + ) + else: + self.T_mix = T_prd_wh + + self.power_output_MWh_kg = self.compute_geofluid_consumption(self.T_mix, T_amb, m_turbine) + self.power_output_MWe = self.compute_power_output(m_turbine) + self.power_generation_MWh = self.power_output_MWe * (self.timestep.total_seconds() / 3600) + self.T_inj = self.compute_injection_temp(self.T_mix, T_amb, m_turbine) + + +class GEOPHIRESORCPowerPlant(BasePowerPlant): + """ORC Binary power plant based on GEOPHIRES.""" + + def __init__(self, ppc, Tres, cf=1.0): + """Define attributes for the BasePowerPlant class. + + Args: + ppc (float): power plant nameplate capacity in MW. + Tres (float): reservoir temperature in deg C. + cf (float): capacity factor. Defualts to 1. + """ + super(GEOPHIRESORCPowerPlant, self).__init__(ppc, Tres, cf) + + def compute_cplant(self, MaxProducedTemperature, min_cost=0): + """Compute power plant capex. + + Args: + MaxProducedTemperature (float): maximum produced geofluid temperature in deg C. + + Returns: + float: cost of power plant in USD/kW + """ + if MaxProducedTemperature < 150.0: + C3 = -1.458333e-3 + C2 = 7.6875e-1 + C1 = -1.347917e2 + C0 = 1.0075e4 + CCAPP1 = C3 * MaxProducedTemperature**3 + C2 * MaxProducedTemperature**2 + C1 * MaxProducedTemperature + C0 + else: + CCAPP1 = 2231 - 2 * (MaxProducedTemperature - 150.0) + + Cplantcorrelation = ( + max(min_cost, CCAPP1 * math.pow(self.powerplant_capacity / 15.0, -0.06)) + * 1e-6 + * self.powerplant_capacity + * 1e3 + ) # $MM + return Cplantcorrelation + + def compute_geofluid_consumption(self, T, T_amb, m_turbine): + """Compute geofluid consumption. + + Args: + T (float): power plant inlet temperature in deg C. + T_amb (float): ambient temperautre in deg C. + + Returns: + float: power plant output in MWh/kg. + """ + + self.thermal_exergy = self.compute_thermalexergy(T, T_amb) + if T_amb < 15.0: + C1 = 2.746e-3 + C0 = -8.3806e-2 + D1 = 2.713e-3 + D0 = -9.1841e-2 + Tfraction = (T_amb - 5.0) / 10.0 + else: + C1 = 2.713e-3 + C0 = -9.1841e-2 + D1 = 2.676e-3 + D0 = -1.012e-1 + Tfraction = (T_amb - 15.0) / 10.0 + etaull = C1 * T + C0 + etauul = D1 * T + D0 + self.etau = (1 - Tfraction) * etaull + Tfraction * etauul + + power_output_MWh_kg = max(self.thermal_exergy * self.etau / 3600, 0.0) # MWh/kg + + return power_output_MWh_kg # MWh/kg + + def compute_injection_temp(self, T, T_amb, m_turbine): + """Compute injection temperature. + + Args: + T (float): power plant inlet temperature in deg C. + T_amb (float): ambient temperautre in deg C. + + Returns: + float: power plant condenser outlet temperature in deg C. + """ + + """Compute the exiting (reinjection) water temperature of an ORC binary power plant.""" + + if T_amb < 15.0: + C1 = 0.0894 + C0 = 55.6 + D1 = 0.0894 + D0 = 62.6 + Tfraction = (T_amb - 5.0) / 10.0 + else: + C1 = 0.0894 + C0 = 62.6 + D1 = 0.0894 + D0 = 69.6 + Tfraction = (T_amb - 15.0) / 10.0 + reinjtll = C1 * T + C0 + reinjtul = D1 * T + D0 + Tinj = max((1.0 - Tfraction) * reinjtll + Tfraction * reinjtul, 0.0) # deg C + + return Tinj # deg C + + +class GEOPHIRESFlashPowerPlant(BasePowerPlant): + """Single Flash Binary power plant.""" + + def __init__(self, ppc, Tres, cf=1.0): + """Define attributes for the BasePowerPlant class. + + Args: + ppc (float): power plant nameplate capacity in MW. + Tres (float): reservoir temperature in deg C. + cf (float): capacity factor. Defualts to 1. + """ + + super(GEOPHIRESFlashPowerPlant, self).__init__(ppc, Tres, cf) + + def compute_cplant(self, MaxProducedTemperature, min_cost=0): + """Compute power plant capex. + + Args: + MaxProducedTemperature (float): maximum produced geofluid temperature in deg C. + + Returns: + float: cost of power plant in USD/kW + """ + + if self.powerplant_capacity < 10: + C2 = 4.8472e-2 + C1 = -35.2186 + C0 = 8.4474e3 + D2 = 4.0604e-2 + D1 = -29.3817 + D0 = 6.9911e3 + PLL = 5.0 + PRL = 10.0 + elif self.powerplant_capacity < 25: + C2 = 4.0604e-2 + C1 = -29.3817 + C0 = 6.9911e3 + D2 = 3.2773e-2 + D1 = -23.5519 + D0 = 5.5263e3 + PLL = 10.0 + PRL = 25.0 + elif self.powerplant_capacity < 50: + C2 = 3.2773e-2 + C1 = -23.5519 + C0 = 5.5263e3 + D2 = 3.4716e-2 + D1 = -23.8139 + D0 = 5.1787e3 + PLL = 25.0 + PRL = 50.0 + elif self.powerplant_capacity < 75: + C2 = 3.4716e-2 + C1 = -23.8139 + C0 = 5.1787e3 + D2 = 3.5271e-2 + D1 = -24.3962 + D0 = 5.1972e3 + PLL = 50.0 + PRL = 75.0 + else: + C2 = 3.5271e-2 + C1 = -24.3962 + C0 = 5.1972e3 + D2 = 3.3908e-2 + D1 = -23.4890 + D0 = 5.0238e3 + PLL = 75.0 + PRL = 100.0 + CCAPPLL = C2 * MaxProducedTemperature**2 + C1 * MaxProducedTemperature + C0 + CCAPPRL = D2 * MaxProducedTemperature**2 + D1 * MaxProducedTemperature + D0 + b = math.log(CCAPPRL / CCAPPLL) / math.log(PRL / PLL) + a = CCAPPRL / PRL**b + Cplantcorrelation = ( + max(min_cost, 0.8 * a * math.pow(self.powerplant_capacity, b)) * self.powerplant_capacity * 1000.0 / 1e6 + ) # factor 0.75 to make double flash 25% more expansive than single flash + + return Cplantcorrelation + + def compute_geofluid_consumption(self, T, T_amb, m_turbine): + """Compute geofluid consumption. + + Args: + T (float): power plant inlet temperature in deg C. + T_amb (float): ambient temperautre in deg C. + + Returns: + float: power plant output in MWh/kg. + """ + + self.thermal_exergy = self.compute_thermalexergy(T, T_amb) + if T_amb < 15.0: + C2 = -4.27318e-7 + C1 = 8.65629e-4 + C0 = 1.78931e-1 + D2 = -5.85412e-7 + D1 = 9.68352e-4 + D0 = 1.58056e-1 + Tfraction = (T_amb - 5.0) / 10.0 + else: + C2 = -5.85412e-7 + C1 = 9.68352e-4 + C0 = 1.58056e-1 + D2 = -7.78996e-7 + D1 = 1.09230e-3 + D0 = 1.33708e-1 + Tfraction = (T_amb - 15.0) / 10.0 + etaull = C2 * T**2 + C1 * T + C0 + etauul = D2 * T**2 + D1 * T + D0 + self.etau = (1.0 - Tfraction) * etaull + Tfraction * etauul + + power_output_MWh_kg = max(self.thermal_exergy * self.etau / 3600, 0.0) # MWh/kg + + return power_output_MWh_kg # MWh/kg + + def compute_injection_temp(self, T, T_amb, m_turbine): + """Compute injection temperature. + + Args: + T (float): power plant inlet temperature in deg C. + T_amb (float): ambient temperautre in deg C. + + Returns: + float: power plant condenser outlet temperature in deg C. + """ + + """Compute the exiting (reinjection) water temperature of an ORC binary power plant.""" + + if T_amb < 15.0: + C2 = -1.11519e-3 + C1 = 7.79126e-1 + C0 = -10.2242 + D2 = -1.10232e-3 + D1 = 7.83893e-1 + D0 = -5.17039 + Tfraction = (T_amb - 5.0) / 10.0 + else: + C2 = -1.10232e-3 + C1 = 7.83893e-1 + C0 = -5.17039 + D2 = -1.08914e-3 + D1 = 7.88562e-1 + D0 = -1.89707e-1 + Tfraction = (T_amb - 15.0) / 10.0 + reinjtll = C2 * T**2 + C1 * T + C0 + reinjtul = D2 * T**2 + D1 * T + D0 + Tinj = max((1.0 - Tfraction) * reinjtll + Tfraction * reinjtul, 0.0) + + return Tinj # deg C + + +class HighEnthalpyCLGWGPowerPlant(BasePowerPlant): + """On/Off-Design Air-cooled ORC Binary power plant developed in Python by Aljubran et al. (2024).""" + + def __init__(self, Tres, Tamb, m_prd, num_prd, cf=1.0, high_resolution=False, k=2): + """Define attributes for the BasePowerPlant class. + + Args: + ppc (float): power plant nameplate capacity in MW. + Tres (float): reservoir temperature in deg C. + m_prd (float): mass flow rate of a single producer in the project in kg/s. + cf (float): capacity factor. Defualts to 1. + high_resolution (bool): if True, a higher resolution power plant model is loaded, but it is slower + """ + self.Tres_design = Tres + self.Tamb_design = Tamb + self.m_prd_design = m_prd + self.num_prd = num_prd + self.k = k + self.model = CustomUnpickler( + open(os.path.join(parent_path, 'data/powerplants', 'HighEnthalpyCLGWGPowerPlant.pkl'), 'rb') + ).load() + + model_input = np.array([[self.Tamb_design, self.Tres_design, self.m_prd_design]]) + model_output = self.model.predict(model_input)[0] + + ppc = model_output[0] * num_prd # MWe power plant capacity for all wells + super(HighEnthalpyCLGWGPowerPlant, self).__init__(ppc, Tres, cf) + + self.power_output_MWh_kg = model_output[0] / self.m_prd_design / 3600 # MWh/kg + self.T_inj = model_output[1] # deg C + + def compute_cplant(self, MaxProducedTemperature, min_cost=0): + """Compute power plant capex. + + Args: + MaxProducedTemperature (float): maximum produced geofluid temperature in deg C. + + Returns: + float: cost of power plant in USD/kW + """ + if MaxProducedTemperature < 150.0: + C3 = -1.458333e-3 + C2 = 7.6875e-1 + C1 = -1.347917e2 + C0 = 1.0075e4 + CCAPP1 = C3 * MaxProducedTemperature**3 + C2 * MaxProducedTemperature**2 + C1 * MaxProducedTemperature + C0 + else: + CCAPP1 = 2231 - 2 * (MaxProducedTemperature - 150.0) + Cplantcorrelation = ( + max(min_cost, CCAPP1 * math.pow(self.powerplant_capacity / 15.0, -0.06)) + * 1e-6 + * self.powerplant_capacity + * 1e3 + ) # $MM + return Cplantcorrelation + + def compute_geofluid_consumption(self, T, T_amb, m_turbine): + """Compute geofluid consumption. + + Args: + T (float): power plant inlet temperature in deg C. + T_amb (float): ambient temperautre in deg C. + + Returns: + float: power plant output in MWh/kg. + """ + + m_prd = m_turbine / self.num_prd # single well flow rate + + model_input = np.array([[T_amb, T, m_prd]]) + model_output = self.model.predict(model_input)[0] + + self.power_output_MWh_kg = max(model_output[0] / m_prd / 3600, 0) # MWh/kg + if self.power_output_MWh_kg == 0: + self.T_inj = T + else: + self.T_inj = model_output[1] # deg C + + return self.power_output_MWh_kg # MWh/kg + + def compute_injection_temp(self, T, T_amb, m_turbine): + """Compute injection temperature. + + Args: + T (float): power plant inlet temperature in deg C. + T_amb (float): ambient temperautre in deg C. + + Returns: + float: power plant condenser outlet temperature in deg C. + """ + + """Compute the exiting (reinjection) water temperature of an ORC binary power plant.""" + + return self.T_inj # deg C + + +class ORCPowerPlant(BasePowerPlant): + """On/Off-Design Air-cooled ORC Binary power plant developed in Python by Aljubran et al. (2024).""" + + def __init__(self, Tres, Tamb, m_prd, num_prd=None, ppc=None, cf=1.0, k=2, timestep=timedelta(hours=1)): + """Define attributes for the BasePowerPlant class. + + Args: + Tres (float): reservoir temperature in deg C. + Tamb (float): ambient temperature in deg C. + m_prd (float): mass flow rate of a single producer in the project in kg/s. + num_prd (int): number of producers (Note: only used if power plant capacity 'ppc' is not provided) + ppc (float): power plant nameplate capacity in MW. + cf (float): capacity factor. Defualts to 1. + """ + + self.Tres_design = Tres + self.Tamb_design = Tamb + self.m_prd_design = m_prd + self.num_prd = num_prd + self.k = k + self.model = CustomUnpickler( + open(os.path.join(parent_path, 'data/powerplants', 'ORCPowerPlant.pkl'), 'rb') + ).load() + + self.keys = np.array(list(self.model.keys())) + tree = cKDTree(self.keys) + distances, indices = tree.query(np.array([[self.Tamb_design, self.Tres_design, self.m_prd_design]]), k=self.k) + self.distances, self.indices = distances[0], indices[0] + if isinstance(self.distances, float): + self.distances, self.indices = np.array([self.distances]), np.array([self.indices]) + self.inv_distances = ( + 1 / (self.distances + 1e-3) ** 2 + ) # important choice on how much to interpolate ... more impactful when simualted data is sparse ... + self.inv_distances = self.inv_distances / self.inv_distances.sum() + self.inv_distances = self.inv_distances[:, None] + self.model_list = [self.model[tuple(self.keys[self.indices][i])] for i in range(self.k)] + + model_input = np.array([[self.Tamb_design, self.Tres_design, self.m_prd_design]]) + preds = np.vstack([m.predict(model_input) for m in self.model_list]) + model_output = (preds * self.inv_distances).sum(axis=0) + + if ppc is None: + ppc = max(model_output[0] * num_prd, 1) # MWe power plant capacity for all wells + super(ORCPowerPlant, self).__init__(ppc, Tres, cf, timestep=timestep) + + self.power_output_MWh_kg = np.clip( + model_output[0] / self.m_prd_design / 3600, a_min=0.0, a_max=2.5e-4 + ) # MWh/kg + self.T_inj = np.clip(model_output[1], a_min=self.Tamb_design, a_max=self.Tres_design) # deg C + + if self.num_prd is None: + m_g = ppc / _Utils.nonzero(self.power_output_MWh_kg) / 3600 # kg/s + self.num_prd = np.ceil(m_g / m_prd).astype(int) + + def compute_cplant(self, MaxProducedTemperature, min_cost=0): + """Compute power plant capex. + + Args: + MaxProducedTemperature (float): maximum produced geofluid temperature in deg C. + + Returns: + float: cost of power plant in USD/kW + """ + + if MaxProducedTemperature < 150.0: + C3 = -1.458333e-3 + C2 = 7.6875e-1 + C1 = -1.347917e2 + C0 = 1.0075e4 + CCAPP1 = C3 * MaxProducedTemperature**3 + C2 * MaxProducedTemperature**2 + C1 * MaxProducedTemperature + C0 + else: + CCAPP1 = 2231 - 2 * (MaxProducedTemperature - 150.0) + Cplantcorrelation = ( + max(min_cost, CCAPP1 * math.pow(self.powerplant_capacity / 15.0, -0.06)) + * 1e-6 + * self.powerplant_capacity + * 1e3 + ) # $MM + return Cplantcorrelation + + def compute_geofluid_consumption(self, T, T_amb, m_turbine): + """Compute geofluid consumption. + + Args: + T (float): power plant inlet temperature in deg C. + T_amb (float): ambient temperautre in deg C. + + Returns: + float: power plant output in MWh/kg. + """ + + m_prd = m_turbine / self.num_prd # single well flow rate + model_input = np.array([[T_amb, T, m_prd]]) + + preds = np.vstack([m.predict(model_input) for m in self.model_list]) + model_output = (preds * self.inv_distances).sum(axis=0) + self.power_output_MWh_kg = np.clip(model_output[0] / m_prd / 3600, a_min=0.0, a_max=2.5e-4) # MWh/kg + self.T_inj = np.clip(model_output[1], a_min=T_amb, a_max=T) # deg C + + return self.power_output_MWh_kg # MWh/kg + + def compute_injection_temp(self, T, T_amb, m_turbine): + """Compute injection temperature. + + Args: + T (float): power plant inlet temperature in deg C. + T_amb (float): ambient temperautre in deg C. + + Returns: + float: power plant condenser outlet temperature in deg C. + """ + + """Compute the exiting (reinjection) water temperature of an ORC binary power plant.""" + + return self.T_inj # deg C diff --git a/src/fgem/utils/__init__.py b/src/fgem/utils/__init__.py new file mode 100755 index 000000000..2900eda72 --- /dev/null +++ b/src/fgem/utils/__init__.py @@ -0,0 +1,3 @@ +from .config import * +from .constants import * +from .utils import * diff --git a/src/fgem/utils/config.py b/src/fgem/utils/config.py new file mode 100755 index 000000000..69659f3c8 --- /dev/null +++ b/src/fgem/utils/config.py @@ -0,0 +1,17 @@ +import json + + +def get_config_from_json(json_file): + """Parse the configurations from the config json file provided. + + Args: + json_file (str): file path to configuration json file + + Returns: + dict: configuration in json format + """ + + with open(json_file) as config_file: + config_dict = json.load(config_file) + + return config_dict diff --git a/src/fgem/utils/constants.py b/src/fgem/utils/constants.py new file mode 100755 index 000000000..c95c3c4c8 --- /dev/null +++ b/src/fgem/utils/constants.py @@ -0,0 +1,198 @@ +us_state_to_abbrev = { + 'Alabama': 'AL', + 'Alaska': 'AK', + 'Arizona': 'AZ', + 'Arkansas': 'AR', + 'California': 'CA', + 'Colorado': 'CO', + 'Connecticut': 'CT', + 'Delaware': 'DE', + 'Florida': 'FL', + 'Georgia': 'GA', + 'Hawaii': 'HI', + 'Idaho': 'ID', + 'Illinois': 'IL', + 'Indiana': 'IN', + 'Iowa': 'IA', + 'Kansas': 'KS', + 'Kentucky': 'KY', + 'Louisiana': 'LA', + 'Maine': 'ME', + 'Maryland': 'MD', + 'Massachusetts': 'MA', + 'Michigan': 'MI', + 'Minnesota': 'MN', + 'Mississippi': 'MS', + 'Missouri': 'MO', + 'Montana': 'MT', + 'Nebraska': 'NE', + 'Nevada': 'NV', + 'New Hampshire': 'NH', + 'New Jersey': 'NJ', + 'New Mexico': 'NM', + 'New York': 'NY', + 'North Carolina': 'NC', + 'North Dakota': 'ND', + 'Ohio': 'OH', + 'Oklahoma': 'OK', + 'Oregon': 'OR', + 'Pennsylvania': 'PA', + 'Rhode Island': 'RI', + 'South Carolina': 'SC', + 'South Dakota': 'SD', + 'Tennessee': 'TN', + 'Texas': 'TX', + 'Utah': 'UT', + 'Vermont': 'VT', + 'Virginia': 'VA', + 'Washington': 'WA', + 'West Virginia': 'WV', + 'Wisconsin': 'WI', + 'Wyoming': 'WY', + 'District of Columbia': 'DC', + 'American Samoa': 'AS', + 'Guam': 'GU', + 'Northern Mariana Islands': 'MP', + 'Puerto Rico': 'PR', + 'United States Minor Outlying Islands': 'UM', + 'U.S. Virgin Islands': 'VI', +} + +WESTERN_STATE_NAMES = [ + 'Arizona', + 'California', + 'Colorado', + 'Idaho', + 'Montana', + 'Nevada', + 'New Mexico', + 'Oregon', + 'Utah', + 'Washington', + 'Wyoming', +] + +regions_to_states = { + 'Southeast': [ + 'West Virginia', + 'District of Columbia', + 'Maryland', + 'Virginia', + 'Kentucky', + 'Tennessee', + 'North Carolina', + 'Mississippi', + 'Arkansas', + 'Louisiana', + 'Alabama', + 'Georgia', + 'South Carolina', + 'Florida', + 'Delaware', + ], + 'Southwest': ['Arizona', 'New Mexico', 'Oklahoma', 'Texas'], + 'West': [ + 'Washington', + 'Oregon', + 'California', + 'Nevada', + 'Idaho', + 'Montana', + 'Wyoming', + 'Utah', + 'Colorado', + 'Alaska', + 'Hawaii', + ], + 'Midwest': [ + 'North Dakota', + 'South Dakota', + 'Nebraska', + 'Kansas', + 'Minnesota', + 'Iowa', + 'Missouri', + 'Wisconsin', + 'Illinois', + 'Michigan', + 'Indiana', + 'Ohio', + ], + 'Northeast': [ + 'Maine', + 'Vermont', + 'New York', + 'New Hampshire', + 'Massachusetts', + 'Rhode Island', + 'Connecticut', + 'New Jersey', + 'Pennsylvania', + ], +} + +states_to_regions = { + 'Washington': 'West', + 'Oregon': 'West', + 'California': 'West', + 'Nevada': 'West', + 'Idaho': 'West', + 'Montana': 'West', + 'Wyoming': 'West', + 'Utah': 'West', + 'Colorado': 'West', + 'Alaska': 'West', + 'Hawaii': 'West', + 'Maine': 'Northeast', + 'Vermont': 'Northeast', + 'New York': 'Northeast', + 'New Hampshire': 'Northeast', + 'Massachusetts': 'Northeast', + 'Rhode Island': 'Northeast', + 'Connecticut': 'Northeast', + 'New Jersey': 'Northeast', + 'Pennsylvania': 'Northeast', + 'North Dakota': 'Midwest', + 'South Dakota': 'Midwest', + 'Nebraska': 'Midwest', + 'Kansas': 'Midwest', + 'Minnesota': 'Midwest', + 'Iowa': 'Midwest', + 'Missouri': 'Midwest', + 'Wisconsin': 'Midwest', + 'Illinois': 'Midwest', + 'Michigan': 'Midwest', + 'Indiana': 'Midwest', + 'Ohio': 'Midwest', + 'West Virginia': 'Southeast', + 'District of Columbia': 'Southeast', + 'Maryland': 'Southeast', + 'Virginia': 'Southeast', + 'Kentucky': 'Southeast', + 'Tennessee': 'Southeast', + 'North Carolina': 'Southeast', + 'Mississippi': 'Southeast', + 'Arkansas': 'Southeast', + 'Louisiana': 'Southeast', + 'Alabama': 'Southeast', + 'Georgia': 'Southeast', + 'South Carolina': 'Southeast', + 'Florida': 'Southeast', + 'Delaware': 'Southeast', + 'Arizona': 'Southwest', + 'New Mexico': 'Southwest', + 'Oklahoma': 'Southwest', + 'Texas': 'Southwest', +} + +UNIFIED_CRS = 'EPSG:3857' +LATLON_CRS = 'EPSG:4326' + +US_LATMIN, US_LATMAX = 24.521208, 49.382808 +US_LONMIN, US_LONMAX = -124.736342, -66.945392 + +SMALL_NUM = 1e-3 +SMALLER_NUM = 1e-6 + +# invert the dictionary +abbrev_to_us_state = dict(map(reversed, us_state_to_abbrev.items())) diff --git a/src/fgem/utils/geocluster_utils.py b/src/fgem/utils/geocluster_utils.py new file mode 100755 index 000000000..ea5a445fe --- /dev/null +++ b/src/fgem/utils/geocluster_utils.py @@ -0,0 +1,2451 @@ +import itertools as iter +import sys + +import h5py +import numpy as np +import scipy.io +from matplotlib import cm +from matplotlib import pyplot as plt +from scipy.interpolate import interpn + + +class data: + def __init__(self, fname, case, fluid): + + self.fluid = fluid + self.case = case + + with h5py.File(fname, 'r') as file: + fixed_loc = '/' + case + '/fixed_params/' + input_loc = '/' + case + '/' + fluid + '/input/' + output_loc = '/' + case + '/' + fluid + '/output/' + + # independent vars + self.mdot = file[input_loc + 'mdot'][:] # i0 + self.L2 = file[input_loc + 'L2'][:] # i1 + self.L1 = file[input_loc + 'L1'][:] # i2 + self.grad = file[input_loc + 'grad'][:] # i3 + self.D = file[input_loc + 'D'][:] # i4 + self.Tinj = file[input_loc + 'T_i'][:] # i5 + self.k = file[input_loc + 'k_rock'][:] # i6 + self.time = file[input_loc + 'time'][:] # i7 + self.ivars = (self.mdot, self.L2, self.L1, self.grad, self.D, self.Tinj, self.k, self.time) + + # fixed vars + self.Pinj = file[fixed_loc + 'Pinj'][()] + self.Tamb = file[fixed_loc + 'Tamb'][()] + + # dim = Mdot x L2 x L1 x grad x D x Tinj x k + self.Wt = file[output_loc + 'Wt'][:] # int mdot * dh dt + self.We = file[output_loc + 'We'][:] # int mdot * (dh - Too * ds) dt + + self.GWhr = 1e6 * 3600000.0 + + self.kWe_avg = self.We * self.GWhr / (1000.0 * self.time[-1] * 86400.0 * 365.0) + self.kWt_avg = self.Wt * self.GWhr / (1000.0 * self.time[-1] * 86400.0 * 365.0) + + # dim = Mdot x L2 x L1 x grad x D x Tinj x k x time + self.shape = ( + len(self.mdot), + len(self.L2), + len(self.L1), + len(self.grad), + len(self.D), + len(self.Tinj), + len(self.k), + len(self.time), + ) + self.Tout = self.__uncompress(file, output_loc, 'Tout') + self.Pout = self.__uncompress(file, output_loc, 'Pout') + + self.CP_fluid = 'CO2' + if fluid == 'H2O': + self.CP_fluid = 'H2O' + + def __uncompress(self, file, output_loc, state): + U = file[output_loc + state + '/' + 'U'][:] + sigma = file[output_loc + state + '/' + 'sigma'][:] + Vt = file[output_loc + state + '/' + 'Vt'][:] + M_k = np.dot(U, np.dot(np.diag(sigma), Vt)) + + shape = self.shape + valid_runs = np.argwhere(np.isfinite(self.We.flatten()))[:, 0] + M_k_full = np.full((shape[-1], np.prod(shape[:-1])), np.nan) + M_k_full[:, valid_runs] = M_k + return np.reshape(M_k_full.T, shape) + + def interp_outlet_states(self, point): + + points = list( + iter.product( + (point[0],), (point[1],), (point[2],), (point[3],), (point[4],), (point[5],), (point[6],), self.time + ) + ) + Tout = interpn(self.ivars, self.Tout, points) + Pout = interpn(self.ivars, self.Pout, points) + + return Tout, Pout + + def interp_kWe_avg(self, point): + ivars = self.ivars[:-1] + return self.GWhr * interpn(ivars, self.We, point) / (1000.0 * self.time[-1] * 86400.0 * 365.0) + + def interp_kWt_avg(self, point): + ivars = self.ivars[:-1] + return self.GWhr * interpn(ivars, self.Wt, point) / (1000.0 * self.time[-1] * 86400.0 * 365.0) + + +class TEA: + def __init__( + self, + Fluid, + End_use, + Configuration, + Flow_user, + Hor_length_user, + Depth_user, + Gradient_user, + Diameter_user, + Tin_user, + krock_user, + Drilling_cost_per_m, + O_and_M_cost_plant, + Discount_rate, + Pump_efficiency, + Lifetime, + Direct_use_heat_cost_per_kWth, + Electricity_rate, + Power_plant_cost_per_kWe, + T0, + P0, + Turbine_isentropic_efficiency, + Generator_efficiency, + Compressor_isentropic_efficiency, + Pre_Cooling_Delta_T, + Turbine_outlet_pressure, + ): + self.Fluid = Fluid + self.End_use = End_use + self.Configuration = Configuration + self.Flow_user = Flow_user + self.Hor_length_user = Hor_length_user + self.Depth_user = Depth_user + self.Gradient_user = Gradient_user + self.Diameter_user = Diameter_user + self.Tin_user = Tin_user + self.krock_user = krock_user + self.Drilling_cost_per_m = Drilling_cost_per_m + self.O_and_M_cost_plant = O_and_M_cost_plant + self.Discount_rate = Discount_rate + self.Pump_efficiency = Pump_efficiency + self.Lifetime = Lifetime + self.Direct_use_heat_cost_per_kWth = Direct_use_heat_cost_per_kWth + self.Electricity_rate = Electricity_rate + self.Power_plant_cost_per_kWe = Power_plant_cost_per_kWe + self.T0 = T0 + self.P0 = P0 + self.Turbine_isentropic_efficiency = Turbine_isentropic_efficiency + self.Generator_efficiency = Generator_efficiency + self.Compressor_isentropic_efficiency = Compressor_isentropic_efficiency + self.Pre_Cooling_Delta_T = Pre_Cooling_Delta_T + self.Turbine_outlet_pressure = Turbine_outlet_pressure + + # self.filename = 'clgs_results.h5' #Filename of h5 database with simulation results [-] + self.filename = 'clgs_results_final.h5' # Filename of h5 database with simulation results [-] + self.Number_of_points_per_year = 4 # Number of time steps per year in database [-] (must be 4) + + self.point = (Flow_user, Hor_length_user, Depth_user, Gradient_user, Diameter_user, Tin_user, krock_user) + + self.P_in = 2e7 # Constant Injection pressure [Pa] + self.T_in = Tin_user - 273.15 # Injection temperature [deg.C] + + def verify(self): # Verify inputs are within allowable bounds + self.error = 0 + if self.Fluid != 1 and self.Fluid != 2: + print('Error: Fluid must be 1 (H2O) or 2 (CO2). Simulation terminated.') + self.error = 1 + if self.End_use != 1 and self.End_use != 2: + print('Error: End_use must be 1 (Direct-Use) or 2 (Electricity). Simulation terminated.') + self.error = 1 + if self.Flow_user < 5 or self.Flow_user > 100: + print('Error: Flow rate must be between 5 and 100 kg/s. Simulation terminated.') + self.error = 1 + if self.Hor_length_user < 1000 or self.Hor_length_user > 20000: + print('Error: Horizontal length must be between 1,000 and 20,000 m. Simulation terminated.') + self.error = 1 + if self.Depth_user < 1000 or self.Depth_user > 5000: + print('Error: Vertical depth must be between 1,000 and 5,000 m. Simulation terminated.') + self.error = 1 + if self.Gradient_user < 0.03 or self.Gradient_user > 0.07: + print('Error: Geothermal gradient must be between 0.03 and 0.07 degrees C per m. Simulation terminated.') + self.error = 1 + if self.Diameter_user < 0.2159 or self.Diameter_user > 0.4445: + print('Error: Wellbore diameter must be between 0.2159 and 0.4445 m. Simulation terminated.') + self.error = 1 + if self.Tin_user < 303.15 or self.Tin_user > 333.15: + print('Error: Injection temperature must be between 303.15 and 333.15 K. Simulation terminated.') + self.error = 1 + if self.krock_user < 1.5 or self.krock_user > 4.5: + print('Error: Rock thermal conductivity must be between 1.5 and 4.5 W/m/K. Simulation terminated.') + self.error = 1 + if self.Drilling_cost_per_m < 0 or self.Drilling_cost_per_m > 10000: + print( + 'Error: Drilling costs per m of measured depth must be between 0 and 10,000 $/m. Simulation terminated.' + ) + self.error = 1 + if self.O_and_M_cost_plant < 0 or self.O_and_M_cost_plant > 0.2: + print( + 'Error: Operation & maintance cost of surface plant (expressed as fraction of total surface plant capital cost) must be between 0 and 0.2. Simulation terminated.' + ) + self.error = 1 + if self.Discount_rate < 0 or self.Discount_rate > 0.2: + print('Error: Discount rate must be between 0 and 0.2. Simulation terminated.') + self.error = 1 + if self.Pump_efficiency < 0.5 or self.Pump_efficiency > 1: + print('Error: Pump efficiency must be between 0.5 and 1. Simulation terminated.') + self.error = 1 + if self.Lifetime < 5 or self.Lifetime > 40: + print('Error: System lifetime must be between 5 and 40 years. Simulation terminated.') + self.error = 1 + if isinstance(self.Lifetime, int) == False: + print('Error: System lifetime must be integer. Simulation terminated.') + self.error = 1 + if self.End_use == 1: + if self.Direct_use_heat_cost_per_kWth < 0 or self.Direct_use_heat_cost_per_kWth > 10000: + print( + 'Error: Capital cost for direct-use surface plant must be between 0 and 10,000 $/kWth. Simulation terminated.' + ) + self.error = 1 + if self.Electricity_rate < 0 or self.Electricity_rate > 0.5: + print( + 'Error: Electricity rate in direct-use for pumping power must be between 0 and 0.5 $/kWh. Simulation terminated.' + ) + self.error = 1 + if self.End_use == 2: + if self.Power_plant_cost_per_kWe < 0 or self.Power_plant_cost_per_kWe > 10000: + print('Error: Power plant capital cost must be between 0 and 10,000 $/kWe. Simulation terminated.') + self.error = 1 + if self.T0 < 278.15 or self.T0 > 303.15: + print('Error: Dead-state temperature must be between 278.15 and 303.15 K. Simulation terminated.') + self.error = 1 + if self.P0 < 0.8e5 or self.P0 > 1.1e5: + print('Error: Dead state pressure must be between 0.8e5 and 1.1e5 Pa. Simulation terminated.') + self.error = 1 + if self.Fluid == 2 and self.End_use == 2: + if self.Turbine_isentropic_efficiency < 0.8 or self.Turbine_isentropic_efficiency > 1: + print('Error: Turbine isentropic efficiency must be between 0.8 and 1. Simulation terminated.') + self.error = 1 + if self.Generator_efficiency < 0.8 or self.Generator_efficiency > 1: + print('Error: Generator efficiency must be between 0.8 and 1. Simulation terminated.') + self.error = 1 + if self.Compressor_isentropic_efficiency < 0.8 or self.Compressor_isentropic_efficiency > 1: + print('Error: Compressor isentropic efficiency must be between 0.8 and 1. Simulation terminated.') + self.error = 1 + if self.Pre_Cooling_Delta_T < 0 or self.Pre_Cooling_Delta_T > 15: + print( + 'Error: CO2 temperature decline after turbine and before compressor must be between 0 and 15 degrees C. Simulation terminated.' + ) + self.error = 1 + if self.Turbine_outlet_pressure < 75 or self.Turbine_outlet_pressure > 200: + print('Error: CO2 turbine outlet pressure must be between 75 and 200 bar. Simulation terminated.') + self.error = 1 + return self.error + + def initialize(self): + + if self.Fluid == 1: + if self.Configuration == 1: + self.u_H2O = data(self.filename, 'utube', 'H2O') + elif self.Configuration == 2: + self.u_H2O = data(self.filename, 'coaxial', 'H2O') + self.timearray = self.u_H2O.time + self.FlowRateVector = self.u_H2O.mdot # length of 26 + self.HorizontalLengthVector = self.u_H2O.L2 # length of 20 + self.DepthVector = self.u_H2O.L1 # length of 9 + self.GradientVector = self.u_H2O.grad # length of 5 + self.DiameterVector = self.u_H2O.D # length of 3 + self.TinVector = self.u_H2O.Tinj # length of 3 + self.KrockVector = self.u_H2O.k # length of 3 + self.Fluid_name = 'Water' + elif self.Fluid == 2: + if self.Configuration == 1: + self.u_sCO2 = data(self.filename, 'utube', 'sCO2') + elif self.Configuration == 2: + self.u_sCO2 = data(self.filename, 'coaxial', 'sCO2') + self.timearray = self.u_sCO2.time + self.FlowRateVector = self.u_sCO2.mdot # length of 26 + self.HorizontalLengthVector = self.u_sCO2.L2 # length of 20 + self.DepthVector = self.u_sCO2.L1 # length of 9 + self.GradientVector = self.u_sCO2.grad # length of 5 + self.DiameterVector = self.u_sCO2.D # length of 3 + self.TinVector = self.u_sCO2.Tinj # length of 3 + self.KrockVector = self.u_sCO2.k # length of 3 + self.Fluid_name = 'CarbonDioxide' + + self.numberofcases = ( + len(self.FlowRateVector) + * len(self.HorizontalLengthVector) + * len(self.DepthVector) + * len(self.GradientVector) + * len(self.DiameterVector) + * len(self.TinVector) + * len(self.KrockVector) + ) + + self.Time_array = np.linspace( + 0, self.Lifetime * 365 * 24 * 3600, 1 + self.Lifetime * self.Number_of_points_per_year + ) # [s] + self.Linear_time_distribution = self.Time_array / 365 / 24 / 3600 + self.TNOP = self.Lifetime * self.Number_of_points_per_year + 1 # Total number of points for selected lifetime + # Find closests lifetime + closestlifetime = self.timearray.flat[np.abs(self.timearray - self.Lifetime).argmin()] + self.indexclosestlifetime = np.where(self.timearray == closestlifetime)[0][0] + + # load property data + if self.Fluid == 1: + mat = scipy.io.loadmat('properties_H2O.mat') + else: + mat = scipy.io.loadmat('properties_CO2v2.mat') + additional_mat = scipy.io.loadmat('additional_properties_CO2v2.mat') + self.Pvector = mat['Pvector'][0] + self.Tvector = mat['Tvector'][0] + self.density = mat['density'] + self.enthalpy = mat['enthalpy'] + self.entropy = mat['entropy'] + if self.Fluid == 2: + self.Pvector_ap = additional_mat['Pvector_ap'][0] + self.hvector_ap = additional_mat['hvector_ap'][0] + self.svector_ap = additional_mat['svector_ap'][0] + self.TPh = additional_mat['TPh'] + self.hPs = additional_mat['hPs'] + + # Define ORC power plant conversion efficiencies + self.Utilization_efficiency_correlation_temperatures = np.array( + [100, 200, 385] + ) # Linear correlation assumed here based on GEOPHIRES ORC correlation between 100 and 200 deg C [deg.C] plus plateaued above 200 deg. C + self.Utilization_efficiency_correlation_conversion = np.array( + [0.2, 0.45, 0.45] + ) # Efficiency of ORC conversion from production exergy to electricity based on GEOPHIRES correlation [-] + self.Heat_to_power_efficiency_correlation_temperatures = np.array( + [100, 200, 385] + ) # Linear correlation based on Chad Augustine's thesis [deg.C] plus plateaued above 200 deg. C + self.Heat_to_power_efficiency_correlation_conversion = np.array( + [0.05, 0.14, 0.14] + ) # Conversion from enthalpy to electricity [-] + + # Calculate dead-state enthalpy and entropy in case of electricity production + if self.End_use == 2: + self.h_0 = interpn((self.Pvector, self.Tvector), self.enthalpy, np.array([self.P0, self.T0]))[ + 0 + ] # dead-state enthalpy [J/kg] + self.s_0 = interpn((self.Pvector, self.Tvector), self.entropy, np.array([self.P0, self.T0]))[ + 0 + ] # dead-state entropy [J/kg/K] + + # Pre-populate specific heat capacity of air in case of electricity production + if self.End_use == 2: + self.Tair_for_cp_array = np.linspace(0, 100, num=10) + # self.cp_air_array = CP.PropsSI('C','P',self.P0,'T',self.Tair_for_cp_array+273.15,'air') + self.cp_air_array = np.array( + [ + 1005.65818063, + 1005.87727966, + 1006.19281999, + 1006.60616167, + 1007.11890862, + 1007.73265999, + 1008.44882744, + 1009.26850304, + 1010.19236691, + 1011.2206266, + ] + ) + + # Initialize heat/electricity arrays + self.Instantaneous_production_enthalpy = np.zeros(len(self.Time_array)) + self.Instantaneous_temperature_after_isenthalpic_throttling = np.zeros(len(self.Time_array)) + self.Instantaneous_heat_production = np.zeros(len(self.Time_array)) + self.Annual_heat_production = np.zeros(self.Lifetime) + self.Annual_pumping_power = np.zeros(self.Lifetime) + self.Average_fluid_density = np.zeros(len(self.Time_array)) + if self.End_use == 2: # electricity generation + self.Instantaneous_exergy_production = np.zeros( + len(self.Time_array) + ) # Produced exergy only (independent from injection conditions) + self.Instantaneous_exergy_extraction = np.zeros( + len(self.Time_array) + ) # Difference between produced exergy and injected exergy + self.Instantaneous_electricity_production_method_1 = np.zeros( + len(self.Time_array) + ) # based on exergy produced (only for water) + self.Instantaneous_electricity_production_method_2 = np.zeros( + len(self.Time_array) + ) # based on exergy extracted + self.Instantaneous_electricity_production_method_3 = np.zeros( + len(self.Time_array) + ) # based on thermal efficiency + self.Instantaneous_electricity_production_method_4 = np.zeros( + len(self.Time_array) + ) # based on direct turbine expansion (for CO2) + self.Instantaneous_utilization_efficiency_method_1 = np.zeros( + len(self.Time_array) + ) # conversion from produced exergy to electricity + self.Instantaneous_utilization_efficiency_method_2 = np.zeros( + len(self.Time_array) + ) # conversion from extracted exergy to electricity + self.Instantaneous_themal_efficiency = np.zeros( + len(self.Time_array) + ) # conversion from enthalpy to electricity + self.Annual_electricity_production = np.zeros(self.Lifetime) + if self.Fluid == 2: + self.Instantaneous_turbine_power = np.zeros( + len(self.Time_array) + ) # Direct turbine expansion considered for systems using sCO2 + + # Initialize error code + self.error_codes = np.zeros(0) # if error occurs, code will be assigned to this tag + + def getTandP(self): + + if self.Fluid == 1: + self.Tout, self.Pout = self.u_H2O.interp_outlet_states(self.point) + elif self.Fluid == 2: + self.Tout, self.Pout = self.u_sCO2.interp_outlet_states(self.point) + + # Initial time correction (Correct production temperature and pressure at time 0 (the value at time 0 [=initial condition] is not a good representation for the first few months) + self.Tout[0] = self.Tout[1] + self.Pout[0] = self.Pout[1] + + # Extract Tout and Pout over lifetime + self.InterpolatedTemperatureArray = self.Tout[0 : self.indexclosestlifetime + 1] - 273.15 + self.InterpolatedPressureArray = self.Pout[0 : self.indexclosestlifetime + 1] + + def calculateLC(self): + self.Linear_production_temperature = self.InterpolatedTemperatureArray + self.Linear_production_pressure = self.InterpolatedPressureArray + self.AveProductionTemperature = np.average(self.Linear_production_temperature) + self.AveProductionPressure = np.average(self.Linear_production_pressure) / 1e5 # [bar] + self.Flow_rate = self.Flow_user # Total flow rate [kg/s] + self.calculatedrillinglength() + # print(self.Linear_production_temperature) + if min(self.Linear_production_temperature) > self.T_in: + self.calculateheatproduction() + if self.End_use == 2: + self.calculateelectricityproduction() + self.calculatecapex() + self.calculatopex() + + Discount_vector = 1.0 / np.power(1 + self.Discount_rate, np.linspace(0, self.Lifetime - 1, self.Lifetime)) + # if self.End_use == 1: #direct-use heating + self.LCOH = ( + (self.TotalCAPEX + np.sum(self.OPEX_Plant * Discount_vector)) + * 1e6 + / np.sum(self.Annual_heat_production / 1e3 * Discount_vector) + ) # $/MWh + self.LCOH100 = ( + (self.TotalCAPEX100 + np.sum(self.OPEX_Plant * Discount_vector)) + * 1e6 + / np.sum(self.Annual_heat_production / 1e3 * Discount_vector) + ) # $/MWh + self.LCOH500 = ( + (self.TotalCAPEX500 + np.sum(self.OPEX_Plant * Discount_vector)) + * 1e6 + / np.sum(self.Annual_heat_production / 1e3 * Discount_vector) + ) # $/MWh + self.LCOH1500 = ( + (self.TotalCAPEX1500 + np.sum(self.OPEX_Plant * Discount_vector)) + * 1e6 + / np.sum(self.Annual_heat_production / 1e3 * Discount_vector) + ) # $/MWh + if self.LCOH < 0: + self.LCOH = 9999 + self.error_codes = np.append(self.error_codes, 5000) + # elif self.End_use == 2: #electricity production + if self.Average_electricity_production == 0: + self.LCOE = 9999 + self.LCOE100 = float('nan') + self.LCOE500 = float('nan') + self.LCOE1500 = float('nan') + self.error_codes = np.append(self.error_codes, 6000) + else: + self.LCOE = ( + (self.TotalCAPEX + np.sum(self.OPEX_Plant * Discount_vector)) + * 1e6 + / np.sum((self.Annual_electricity_production - self.Annual_pumping_power) / 1e3 * Discount_vector) + ) # $/MWh + self.LCOE100 = ( + (self.TotalCAPEX100 + np.sum(self.OPEX_Plant * Discount_vector)) + * 1e6 + / np.sum((self.Annual_electricity_production - self.Annual_pumping_power) / 1e3 * Discount_vector) + ) # $/MWh + self.LCOE500 = ( + (self.TotalCAPEX500 + np.sum(self.OPEX_Plant * Discount_vector)) + * 1e6 + / np.sum((self.Annual_electricity_production - self.Annual_pumping_power) / 1e3 * Discount_vector) + ) # $/MWh + self.LCOE1500 = ( + (self.TotalCAPEX1500 + np.sum(self.OPEX_Plant * Discount_vector)) + * 1e6 + / np.sum((self.Annual_electricity_production - self.Annual_pumping_power) / 1e3 * Discount_vector) + ) # $/MWh + if self.LCOE < 0: + self.LCOE = 9999 + self.LCOE100 = float('nan') + self.LCOE500 = float('nan') + self.LCOE1500 = float('nan') + self.error_codes = np.append(self.error_codes, 7000) + + else: # Production temperature went below injection temperature + self.LCOE = 9999 + self.LCOH = 9999 + self.error_codes = np.append(self.error_codes, 1000) + + def calculatedrillinglength(self): + if self.Configuration == 1: + self.Drilling_length = ( + self.Hor_length_user + 2 * self.Depth_user + ) # Total drilling depth of both wells and lateral in U-loop [m] + elif self.Configuration == 2: + self.Drilling_length = ( + self.Hor_length_user + self.Depth_user + ) # Total drilling depth of well and lateral in co-axial case [m] + + def calculateheatproduction(self): + # Calculate instantaneous heat production + self.Average_fluid_density = interpn( + (self.Pvector, self.Tvector), + self.density, + np.dstack( + ( + 0.5 * self.P_in + 0.5 * self.Linear_production_pressure, + 0.5 * self.T_in + 0.5 * self.Linear_production_temperature + 273.15, + ) + )[0], + ) + self.hprod = interpn( + (self.Pvector, self.Tvector), + self.enthalpy, + np.dstack((self.Linear_production_pressure, self.Linear_production_temperature + 273.15))[0], + ) + self.hinj = interpn((self.Pvector, self.Tvector), self.enthalpy, np.array([self.P_in, self.T_in + 273.15])) + self.Instantaneous_heat_production = ( + self.Flow_rate * (self.hprod - self.hinj) / 1000 + ) # Heat production based on produced minus injected enthalpy [kW] + + # Calculate annual heat production (kWh) + self.Annual_heat_production = ( + 8760 + / 5 + * ( + self.Instantaneous_heat_production[0::4][0:-1] + + self.Instantaneous_heat_production[1::4] + + self.Instantaneous_heat_production[2::4] + + self.Instantaneous_heat_production[3::4] + + self.Instantaneous_heat_production[4::4] + ) + ) + + # Calculate average heat production + self.AveAnnualHeatProduction = np.average(self.Annual_heat_production) # kWh + self.AveInstHeatProduction = np.average(self.Instantaneous_heat_production) # kWth + + # Calculate average heat production and first year heat production + self.Average_heat_production = np.average(self.Instantaneous_heat_production) # [kW] + # Average_production_temperature = np.average(Linear_production_temperature) #[deg.C] + self.FirstYearHeatProduction = self.Annual_heat_production[0] # kWh + + self.calculatepumpingpower() + + def calculateelectricityproduction(self): + + # Calculate instantaneous exergy production, exergy extraction, and electricity generation (MW) and annual electricity generation [kWh] + self.h_prod = self.hprod # produced enthalpy [J/kg] + self.h_inj = self.hinj # injected enthalpy [J/kg] + self.s_prod = interpn( + (self.Pvector, self.Tvector), + self.entropy, + np.dstack((self.Linear_production_pressure, self.Linear_production_temperature + 273.15))[0], + ) # produced entropy [J/kg/K] + self.s_inj = interpn( + (self.Pvector, self.Tvector), self.entropy, np.array([self.P_in, self.T_in + 273.15]) + ) # injected entropy [J/kg/K] + + self.Instantaneous_exergy_production = ( + self.Flow_rate * (self.h_prod - self.h_0 - self.T0 * (self.s_prod - self.s_0)) + ) / 1000 # [kW] + self.Instantaneous_exergy_extraction = ( + self.Flow_rate * (self.h_prod - self.h_inj - self.T0 * (self.s_prod - self.s_inj)) + ) / 1000 # [kW] + + self.AverageInstNetExergyProduction = np.average(self.Instantaneous_exergy_production) # [kW] + self.AverageInstNetExergyExtraction = np.average(self.Instantaneous_exergy_extraction) # [kW] + + if self.Fluid == 1: + + if ( + self.T_in >= 50 + and min(self.Linear_production_temperature) >= 100 + and max(self.Linear_production_temperature) <= 385 + ): + self.Instantaneous_utilization_efficiency_method_1 = np.interp( + self.Linear_production_temperature, + self.Utilization_efficiency_correlation_temperatures, + self.Utilization_efficiency_correlation_conversion, + left=0, + ) # Utilization efficiency based on conversion of produced exergy to electricity + self.Instantaneous_electricity_production_method_1 = ( + self.Instantaneous_exergy_production * self.Instantaneous_utilization_efficiency_method_1 + ) # [kW] + self.Instantaneous_themal_efficiency = np.interp( + self.Linear_production_temperature, + self.Heat_to_power_efficiency_correlation_temperatures, + self.Heat_to_power_efficiency_correlation_conversion, + left=0, + ) # Utilization efficiency based on conversion of produced exergy to electricity + self.Instantaneous_electricity_production_method_3 = ( + self.Instantaneous_heat_production * self.Instantaneous_themal_efficiency + ) # [kW] + + else: # Water injection temperature and/or production tempeature fall outside the range used in the correlations + self.error_codes = np.append(self.error_codes, 2000) + self.Instantaneous_utilization_efficiency_method_1 = np.zeros(len(self.Time_array)) + self.Instantaneous_electricity_production_method_1 = np.zeros(len(self.Time_array)) + self.Instantaneous_themal_efficiency = np.zeros(len(self.Time_array)) + self.Instantaneous_electricity_production_method_3 = np.zeros(len(self.Time_array)) + + # based on method 1 for now (could be 50-50) + self.Annual_electricity_production = ( + 8760 + / 5 + * ( + self.Instantaneous_electricity_production_method_1[0::4][0:-1] + + self.Instantaneous_electricity_production_method_1[1::4] + + self.Instantaneous_electricity_production_method_1[2::4] + + self.Instantaneous_electricity_production_method_1[3::4] + + self.Instantaneous_electricity_production_method_1[4::4] + ) + ) + self.Inst_electricity_production = self.Instantaneous_electricity_production_method_1 # [kW] + self.AveInstElectricityProduction = np.average(self.Instantaneous_electricity_production_method_1) # [kW] + + elif self.Fluid == 2: + T_prod = self.Linear_production_temperature # Production temperature [deg.C] + P_prod = self.Linear_production_pressure # Production pressure [Pa] + + h_turbine_out_ideal = interpn( + (self.Pvector_ap, self.svector_ap), + self.hPs, + np.dstack((np.ones(self.TNOP) * self.Turbine_outlet_pressure * 1e5, self.s_prod))[0], + ) + self.Instantaneous_turbine_power = ( + self.Flow_rate * (self.h_prod - h_turbine_out_ideal) * self.Turbine_isentropic_efficiency / 1000 + ) # Turbine output [kW] + h_turbine_out_actual = ( + self.h_prod - self.Instantaneous_turbine_power / self.Flow_rate * 1000 + ) # Actual fluid enthalpy at turbine outlet [J/kg] + self.T_turbine_out_actual = ( + interpn( + (self.Pvector_ap, self.hvector_ap), + self.TPh, + np.dstack((np.ones(self.TNOP) * self.Turbine_outlet_pressure * 1e5, h_turbine_out_actual))[0], + ) + - 273.15 + ) + + if min(self.T_turbine_out_actual) > 37 and self.T_in > 32: + + self.Pre_cooling_temperature = min(self.T_turbine_out_actual) - self.Pre_Cooling_Delta_T + Post_cooling = 2000 + valuefound = 0 + lastrun = 0 + # print('y') + while valuefound == 0: + # print('b') + Pre_compressor_h = interpn( + (self.Pvector, self.Tvector), + self.enthalpy, + np.array([self.Turbine_outlet_pressure * 1e5, self.Pre_cooling_temperature + 273.15]), + ) + + Pre_cooling = ( + self.Flow_rate * (h_turbine_out_actual - Pre_compressor_h) / 1e3 + ) # Pre-compressor cooling [kWth] + Pre_compressor_s = interpn( + (self.Pvector, self.Tvector), + self.entropy, + np.array([self.Turbine_outlet_pressure * 1e5, self.Pre_cooling_temperature + 273.15]), + ) + + Post_compressor_h_ideal = interpn( + (self.Pvector_ap, self.svector_ap), self.hPs, np.array([self.P_in, Pre_compressor_s[0]]) + ) + Post_compressor_h_actual = ( + Pre_compressor_h + + (Post_compressor_h_ideal - Pre_compressor_h) / self.Compressor_isentropic_efficiency + ) # Actual fluid enthalpy at compressor outlet [J/kg] + self.Post_compressor_T_actual = ( + interpn( + (self.Pvector_ap, self.hvector_ap), + self.TPh, + np.array([self.P_in, Post_compressor_h_actual[0]]), + ) + - 273.15 + ) + Compressor_Work = self.Flow_rate * (Post_compressor_h_actual - Pre_compressor_h) / 1e3 # [kWe] + Post_cooling = ( + self.Flow_rate * (Post_compressor_h_actual - self.h_inj) / 1e3 + ) # Fluid cooling after compression [kWth] + # print(str(Post_cooling)) + # print(str(self.Pre_cooling_temperature)) + + if lastrun == 0: + if self.Pre_cooling_temperature < 32: + lastrun = 1 + self.Pre_cooling_temperature = self.Pre_cooling_temperature + 0.5 + elif Post_cooling < 0: + self.Pre_cooling_temperature = self.Pre_cooling_temperature + 0.5 + lastrun = 1 + elif Post_cooling > 0: + self.Pre_cooling_temperature = self.Pre_cooling_temperature - 0.5 + + elif lastrun == 1: + valuefound = 1 + print(self.Pre_cooling_temperature - min(self.T_turbine_out_actual)) + + # print(min(self.T_turbine_out_actual)) + # print(self.Pre_cooling_temperature) + + if Post_cooling < 0: + ResistiveHeating = -Post_cooling + Post_cooling = 0 + else: + ResistiveHeating = 0 + + Total_cooling = Pre_cooling + Post_cooling # Total CO2 cooling requirements [kWth] + + T_air_in_pre_cooler = self.T0 - 273.15 + T_air_out_pre_cooler = ( + self.T_turbine_out_actual + self.Pre_cooling_temperature + ) / 2 # Air outlet temperature in pre-cooler [deg.C] + cp_air = np.interp( + 0.5 * T_air_in_pre_cooler + 0.5 * T_air_out_pre_cooler, self.Tair_for_cp_array, self.cp_air_array + ) # Air specific heat capacity in pre-cooler [J/kg/K] + m_air_pre_cooler = ( + Pre_cooling * 1000 / (cp_air * (T_air_out_pre_cooler - T_air_in_pre_cooler)) + ) # Air flow rate in pre-cooler [kg/s] + + T_air_in_post_cooler = self.T0 - 273.15 + T_air_out_post_cooler = ( + self.Post_compressor_T_actual + self.T_in + ) / 2 # Air outlet temperature in post-cooler [deg.C] + cp_air = np.interp( + 0.5 * T_air_in_post_cooler + 0.5 * T_air_out_post_cooler, self.Tair_for_cp_array, self.cp_air_array + ) # Air specific heat capacity in post-cooler [J/kg/K] + + m_air_post_cooler = ( + Post_cooling * 1000 / (cp_air * (T_air_out_post_cooler - T_air_in_post_cooler)) + ) # Air flow rate in post-cooler [kg/s] + + Air_cooling_power = ( + m_air_pre_cooler + m_air_post_cooler + ) * 0.25 # Electricity for air-cooling, assuming 0.25 kWe per kg/s [kWe] + # print(self.Instantaneous_turbine_power) + # print(Compressor_Work) + self.Instantaneous_electricity_production_method_4 = ( + self.Instantaneous_turbine_power * self.Generator_efficiency + - Compressor_Work + - Air_cooling_power + - ResistiveHeating + ) # Net electricity using CO2 direct turbine expansion cycle [kWe] + self.Inst_electricity_production = self.Instantaneous_electricity_production_method_4 # [kW] + self.Annual_electricity_production = ( + 8760 + / 5 + * ( + self.Instantaneous_electricity_production_method_4[0::4][0:-1] + + self.Instantaneous_electricity_production_method_4[1::4] + + self.Instantaneous_electricity_production_method_4[2::4] + + self.Instantaneous_electricity_production_method_4[3::4] + + self.Instantaneous_electricity_production_method_4[4::4] + ) + ) + self.AveInstElectricityProduction = np.average( + self.Instantaneous_electricity_production_method_4 + ) # [kW] + # check if negative + if min(self.Instantaneous_electricity_production_method_4) < 0: + self.error_codes = np.append( + self.error_codes, 5500 + ) # Calculated electricity generation is negative + + self.Annual_electricity_production = np.zeros(self.Lifetime) + self.Inst_electricity_production = np.zeros(self.TNOP) + self.AveInstElectricityProduction = 0 + else: # turbine outlet or reinjection temperature too low + if self.T_in <= 32: + self.error_codes = np.append(self.error_codes, 3000) + + if min(self.T_turbine_out_actual) <= 37: + self.error_codes = np.append(self.error_codes, 4000) + + self.Annual_electricity_production = np.zeros(self.Lifetime) + self.Inst_electricity_production = np.zeros(self.TNOP) + self.AveInstElectricityProduction = 0 + + self.calculatepumpingpower() + + self.Average_electricity_production = np.average(self.Annual_electricity_production) / 8760 # [kW] + self.AveAnnualElectricityProduction = np.average(self.Annual_electricity_production) # [kWh] + self.AveInstNetElectricityProduction = self.AveInstElectricityProduction - np.average(self.PumpingPower) # [kW] + self.AveAnnualNetElectricityProduction = self.AveAnnualElectricityProduction - np.average( + self.Annual_pumping_power + ) # kWh + self.FirstYearElectricityProduction = self.Annual_electricity_production[0] # kWh + self.Inst_Net_Electricity_production = self.Inst_electricity_production - self.PumpingPower # [kW] + + if self.AveInstNetElectricityProduction < 0: + self.AveInstNetElectricityProduction = 0 + + def calculatepumpingpower(self): + # Calculate pumping power + self.PumpingPower = ( + (self.P_in - self.Linear_production_pressure) + * self.Flow_rate + / self.Average_fluid_density + / self.Pump_efficiency + / 1e3 + ) # Pumping power [kW] + self.PumpingPower[self.PumpingPower < 0] = ( + 0 # Set negative values to zero (if the production pressure is above the injection pressure, we throttle the fluid) + ) + self.Annual_pumping_power = ( + 8760 + / 5 + * ( + self.PumpingPower[0::4][0:-1] + + self.PumpingPower[1::4] + + self.PumpingPower[2::4] + + self.PumpingPower[3::4] + + self.PumpingPower[4::4] + ) + ) # kWh + + def calculatecapex(self): + self.CAPEX_Drilling = self.Drilling_length * self.Drilling_cost_per_m / 1e6 # Drilling capital cost [M$] + if self.End_use == 1: # direct-use heating + self.CAPEX_Surface_Plant = ( + np.max(self.Instantaneous_heat_production) * self.Direct_use_heat_cost_per_kWth / 1e6 + ) # [M$] + elif self.End_use == 2: # electricity production + if self.Fluid == 1: + self.CAPEX_Surface_Plant = ( + np.max(self.Instantaneous_electricity_production_method_1) * self.Power_plant_cost_per_kWe / 1e6 + ) # [M$] + elif self.Fluid == 2: + self.CAPEX_Surface_Plant = ( + np.max(self.Instantaneous_electricity_production_method_4) * self.Power_plant_cost_per_kWe / 1e6 + ) # [M$] + + self.TotalCAPEX = ( + self.CAPEX_Drilling + self.CAPEX_Surface_Plant + ) # Total system capital cost (only includes drilling and surface plant cost) [M$] + self.TotalCAPEX100 = ( + self.Drilling_length * 100 / 1e6 + self.CAPEX_Surface_Plant + ) # Total system capital cost (only includes drilling and surface plant cost) [M$] + self.TotalCAPEX500 = ( + self.Drilling_length * 500 / 1e6 + self.CAPEX_Surface_Plant + ) # Total system capital cost (only includes drilling and surface plant cost) [M$] + self.TotalCAPEX1500 = ( + self.Drilling_length * 1500 / 1e6 + self.CAPEX_Surface_Plant + ) # Total system capital cost (only includes drilling and surface plant cost) [M$] + + def calculatopex(self): + # Calculate OPEX + if self.End_use == 1: # direct-use heating + self.OPEX_Plant = ( + self.O_and_M_cost_plant * self.CAPEX_Surface_Plant + + self.Annual_pumping_power * self.Electricity_rate / 1e6 + ) # Annual plant O&M cost [M$/year] + elif self.End_use == 2: # electricity production + self.OPEX_Plant = self.O_and_M_cost_plant * self.CAPEX_Surface_Plant # Annual plant O&M cost [M$/year] + self.AverageOPEX_Plant = np.average(self.OPEX_Plant) + + def printresults(self): + ### Print results to screen + print('##################') + print('Simulation Results') + print('##################') + print('Number of cases in database = ' + str(self.numberofcases)) + print(' ') + print('### Configuration ###') + if self.End_use == 1: + print('End-Use = Direct-Use') + elif self.End_use == 2: + print('End-Use = Electricity') + if self.Fluid == 1: + print('Fluid = water') + elif self.Fluid == 2: + print('Fluid = sCO2') + if self.Configuration == 1: + print('Design = U-loop') + elif self.Configuration == 2: + print('Design = Co-axial') + + # Print conditions + print('Flow rate = ' + f'{self.Flow_user:.1f}' + ' kg/s') + print('Lateral Length = ' + str(round(self.Hor_length_user)) + ' m') + print('Vertical Depth = ' + str(round(self.Depth_user)) + ' m') + print('Geothermal Gradient = ' + f'{self.Gradient_user*1000:.1f}' + ' deg.C/km') + print('Wellbore Diameter = ' + f'{self.Diameter_user:.4f}' + ' m') + print('Injection Temperature = ' + f'{self.Tin_user-273.15:.1f}' + ' deg.C') + print('Thermal Conductivity = ' + f'{self.krock_user:.2f}' + ' W/m/K') + + # check if error occured + if len(self.error_codes) > 0: + print(' ') + if np.in1d(1000, self.error_codes): # plot the temperature and pressure for these + print('Error: production temperature drops below injection temperature. Simulation terminated.\n') + + if np.in1d(2000, self.error_codes): # plot the temperature and pressure for these + print( + 'Error: Water injection temperature and/or production tempeature fall outside the range of the ORC correlations. These correlations require injection temperature larger than 50 deg.C and production temperature in the range of 100 to 385 deg.C. Electricity production set to 0. \n' + ) + + if np.in1d( + 3000, self.error_codes + ): # CO2 injection temperature cannot be below 32 degrees C (CO2 must be supercritical) + print('Error: Too low CO2 reinjection temperature. CO2 must remain supercritical.\n') + + if np.in1d(4000, self.error_codes): # Turbine outlet CO2 temperature dropped below 37 degrees C + print( + 'Error: Too low CO2 turbine outlet temperature. Turbine outlet temperature must be above 37 degrees C.\n' + ) + + if np.in1d(5000, self.error_codes): # Calculated LCOH was negative, set to $9999/MWh + print('Error: Calculated LCOH was negative and has been reset to $9,999/MWh.\n') + + if np.in1d(5500, self.error_codes): # Negative electricity calculated with CO2 cycle + print('Error: Calculated net electricity generation is negative and reset to 0.\n') + + if np.in1d(6000, self.error_codes): # Zero electricity production. LCOE set to $9999/MWh + print('Error: Calculated net electricity production was 0. LCOE reset to $9,999/MWh.\n') + + if np.in1d(7000, self.error_codes): # Calculated LCOE was negative, set to $9999/MWh + print('Error: Calculated LCOE was negative and has been reset to $9,999/MWh.\n') + + # Print results for heating + if self.End_use == 1: + print(' ') + print('### Reservoir Simulation Results ###') + print('Average Production Temperature = ' + f'{self.AveProductionTemperature:.1f}' + ' deg.C') + print('Average Production Pressure = ' + f'{self.AveProductionPressure:.1f}' + ' bar') + if np.in1d(1000, self.error_codes) == False: + print('Average Heat Production = ' + f'{self.AveInstHeatProduction:.1f}' + ' kWth') + print('First Year Heat Production = ' + f'{self.FirstYearHeatProduction/1e3:.1f}' + ' MWh') + print(' ') + print('### Cost Results ###') + print('Total CAPEX = ' + f'{self.TotalCAPEX:.1f}' + ' M$') + print('Drilling Cost = ' + f'{self.CAPEX_Drilling:.1f}' + ' M$') + print('Surface Plant Cost = ' + f'{self.CAPEX_Surface_Plant:.1f}' + ' M$') + print('OPEX = ' + f'{self.AverageOPEX_Plant*1000:.1f}' + ' k$/year') + print('LCOH = ' + f'{self.LCOH:.1f}' + ' $/MWh') + + if self.End_use == 2: + print(' ') + print('### Reservoir Simulation Results ###') + print('Average Production Temperature = ' + f'{self.AveProductionTemperature:.1f}' + ' deg.C') + print('Average Production Pressure = ' + f'{self.AveProductionPressure:.1f}' + ' bar') + if np.in1d(1000, self.error_codes) == False: + print('Average Heat Production = ' + f'{self.AveInstHeatProduction:.1f}' + ' kWth') + print('Average Net Electricity Production = ' + f'{self.AveInstNetElectricityProduction:.1f}' + ' kWe') + print('First Year Heat Production = ' + f'{self.FirstYearHeatProduction/1e3:.1f}' + ' MWh') + print( + 'First Year Electricity Production = ' + f'{self.FirstYearElectricityProduction/1e3:.1f}' + ' MWh' + ) + print(' ') + print('### Cost Results ###') + print('Total CAPEX = ' + f'{self.TotalCAPEX:.1f}' + ' M$') + print('Drilling Cost = ' + f'{self.CAPEX_Drilling:.1f}' + ' M$') + print('Surface Plant Cost = ' + f'{self.CAPEX_Surface_Plant:.1f}' + ' M$') + print('OPEX = ' + f'{self.AverageOPEX_Plant*1000:.1f}' + ' k$/year') + print('LCOE = ' + f'{self.LCOE:.1f}' + ' $/MWh') + + def plotresults(self): + # Plot injection and production temperature + plt.figure() + plt.title('Injection and Production Temperature') + plt.xlabel('Time [year]') + plt.ylabel(r'Temperature [$^\circ$C]') + plt.plot([0, self.Lifetime], [self.T_in, self.T_in], 'b-') + plt.plot(self.Linear_time_distribution, self.Linear_production_temperature, 'r-') + plt.legend(['Injection Temperature', 'Production Temperature']) + plt.axis([0, self.Lifetime, self.T_in - 10, max(self.Linear_production_temperature) + 10]) + plt.xticks(np.arange(0, self.Lifetime + 1, 2), np.arange(0, self.Lifetime + 1, 2)) + plt.show() + + # Plot injection and production pressure + plt.figure() + plt.title('Injection and Production Pressure') + plt.xlabel('Time [year]') + plt.ylabel('Pressure [bar]') + plt.plot([0, self.Lifetime], [self.P_in / 1e5, self.P_in / 1e5], 'b-') + plt.plot(self.Linear_time_distribution, self.Linear_production_pressure / 1e5, 'r-') + plt.legend(['Injection Pressure', 'Production Pressure']) + maxpressure = np.max([self.P_in / 1e5, np.max(self.Linear_production_pressure) / 1e5]) + minpressure = np.min([self.P_in / 1e5, np.min(self.Linear_production_pressure) / 1e5]) + plt.axis([0, self.Lifetime, minpressure - 10, maxpressure + 10]) + plt.xticks(np.arange(0, self.Lifetime + 1, 2), np.arange(0, self.Lifetime + 1, 2)) + plt.show() + + # Plot heat production + plt.figure() + plt.title('Heat Production') + plt.xlabel('Time [year]') + params = {'mathtext.default': 'regular'} + plt.rcParams.update(params) + plt.ylabel('Heat [$MW_{t}$]') + plt.plot(self.Linear_time_distribution, self.Instantaneous_heat_production / 1e3, 'r-') + plt.axis([0, self.Lifetime, 0, max(self.Instantaneous_heat_production / 1e3) + 1]) + plt.xticks(np.arange(0, self.Lifetime + 1, 2), np.arange(0, self.Lifetime + 1, 2)) + plt.show() + + plt.figure() + plt.bar(np.arange(1, self.Lifetime + 1), self.Annual_heat_production / 1e6) + plt.xlabel('Time [year]') + plt.ylabel('Heat Production [GWh]') + plt.title('Annual Heat Production') + plt.xticks(np.arange(0, self.Lifetime + 1, 2), np.arange(0, self.Lifetime + 1, 2)) + plt.show() + + # plot electricity over time if end-use is 2 + if self.End_use == 2: + plt.figure() + plt.title('Electricity Production') + plt.xlabel('Time [year]') + params = {'mathtext.default': 'regular'} + plt.rcParams.update(params) + plt.ylabel('Electricity [$MW_{e}$]') + plt.plot(self.Linear_time_distribution, self.Inst_Net_Electricity_production / 1e3, 'r-') + plt.axis([0, self.Lifetime, 0, max(self.Inst_Net_Electricity_production / 1e3) + 1]) + plt.xticks(np.arange(0, self.Lifetime + 1, 2), np.arange(0, self.Lifetime + 1, 2)) + plt.show() + + plt.figure() + plt.bar(np.arange(1, self.Lifetime + 1), self.Annual_electricity_production / 1e6) + plt.xlabel('Time [year]') + plt.ylabel('Electricity Production [GWh]') + plt.title('Annual Electricity Production') + plt.xticks(np.arange(0, self.Lifetime + 1, 2), np.arange(0, self.Lifetime + 1, 2)) + plt.show() + + # create Ts-diagram in case of CO2 power production + if ( + self.Fluid == 2 + and np.in1d(3000, self.error_codes) == False + and np.in1d(4000, self.error_codes) == False + ): + # CO2 isobaric lines for Ts diagram (do not change data below) + pvectorforTS = np.array([40, 50, 60, 70, 80, 90, 100, 120, 150, 190, 240, 300]) # bar + svectorforTS = np.arange(1000, 2000, 20) # J/kg/K + tmatrix = np.array( + [ + [ + 273.58779987, + 275.73964358, + 277.82952322, + 278.44972407, + 278.44972407, + 278.44972407, + 278.44972407, + 278.44972407, + 278.44972407, + 278.44972407, + 278.44972407, + 278.44972407, + 278.44972407, + 278.44972407, + 278.44972407, + 278.44972407, + 278.44972407, + 278.44972407, + 278.44972407, + 278.44972407, + 278.44972407, + 278.44972407, + 278.44972407, + 278.44972407, + 278.44972407, + 278.44972407, + 278.44972407, + 278.44972407, + 278.44972407, + 278.44972407, + 278.44972407, + 278.44972407, + 278.44972407, + 278.44972407, + 278.44972407, + 278.44972407, + 278.44972407, + 278.44972407, + 278.44972407, + 278.44972407, + 278.44972407, + 279.16739097, + 281.98680218, + 285.11251959, + 288.53183612, + 292.23286854, + 296.20421415, + 300.43497368, + 304.91492276, + 309.63463723, + ], + [ + 274.42077418, + 276.63242586, + 278.79031756, + 280.88470759, + 282.90420093, + 284.8353605, + 286.66228321, + 287.43392381, + 287.43392381, + 287.43392381, + 287.43392381, + 287.43392381, + 287.43392381, + 287.43392381, + 287.43392381, + 287.43392381, + 287.43392381, + 287.43392381, + 287.43392381, + 287.43392381, + 287.43392381, + 287.43392381, + 287.43392381, + 287.43392381, + 287.43392381, + 287.43392381, + 287.43392381, + 287.43392381, + 287.43392381, + 287.43392381, + 287.43392381, + 287.43392381, + 287.43392381, + 287.43392381, + 287.43392381, + 287.43392381, + 287.43392381, + 287.43392381, + 287.96651004, + 290.06555059, + 292.48128297, + 295.20562474, + 298.23009373, + 301.54561859, + 305.14278777, + 309.01219991, + 313.14471796, + 317.53155667, + 322.16425682, + 327.03466505, + ], + [ + 275.23288781, + 277.50029097, + 279.72101728, + 281.88672793, + 283.98779516, + 286.01303489, + 287.94946975, + 289.78218424, + 291.4942541, + 293.06653014, + 294.47679734, + 295.12790099, + 295.12790099, + 295.12790099, + 295.12790099, + 295.12790099, + 295.12790099, + 295.12790099, + 295.12790099, + 295.12790099, + 295.12790099, + 295.12790099, + 295.12790099, + 295.12790099, + 295.12790099, + 295.12790099, + 295.12790099, + 295.12790099, + 295.12790099, + 295.12790099, + 295.12790099, + 295.12790099, + 295.12790099, + 295.12790099, + 295.12790099, + 295.9354915, + 297.35537749, + 299.07212335, + 301.08795261, + 303.40261519, + 306.01350043, + 308.91613741, + 312.10478208, + 315.57290651, + 319.31351482, + 323.31931035, + 327.58280441, + 332.09645076, + 336.85281943, + 341.84475813, + ], + [ + 276.02607997, + 278.34583338, + 280.62513674, + 282.85677401, + 285.03249573, + 287.14282112, + 289.17689654, + 291.12248905, + 292.96614008, + 294.69336153, + 296.28858454, + 297.73458483, + 299.011625, + 300.09758466, + 300.97073532, + 301.61423442, + 301.8325153, + 301.8325153, + 301.8325153, + 301.8325153, + 301.8325153, + 301.8325153, + 301.8325153, + 301.8325153, + 301.8325153, + 301.8325153, + 301.8325153, + 301.8325153, + 301.8325153, + 301.8325153, + 302.09542982, + 302.58879462, + 303.2793647, + 304.18893716, + 305.33653776, + 306.73773627, + 308.40431208, + 310.34433077, + 312.56253551, + 315.06089701, + 317.83917492, + 320.89539075, + 324.2261797, + 327.82705562, + 331.69265609, + 335.81701249, + 340.19383373, + 344.81675258, + 349.67949196, + 354.77594888, + ], + [ + 276.80189842, + 279.17107055, + 281.50533372, + 283.79838944, + 286.04309117, + 288.23128859, + 290.3537384, + 292.40016137, + 294.35947993, + 296.22015087, + 297.97036237, + 299.59784329, + 301.08934154, + 302.43049789, + 303.60731254, + 304.60958484, + 305.43472959, + 306.09016372, + 306.59420747, + 306.97336927, + 307.25532028, + 307.47069845, + 307.65403034, + 307.82936031, + 308.0076719, + 308.20853737, + 308.45497452, + 308.76827506, + 309.16916382, + 309.67981016, + 310.32386748, + 311.12522831, + 312.10669659, + 313.28900162, + 314.69018476, + 316.32529927, + 318.20635339, + 320.34242373, + 322.73986477, + 325.4025539, + 328.33214361, + 331.52832983, + 334.98916153, + 338.71139868, + 342.69088804, + 346.92290167, + 351.40239399, + 356.12416833, + 361.08297566, + 366.27357671, + ], + [ + 277.56162433, + 279.97764204, + 282.36372714, + 284.71434274, + 287.02324203, + 289.28333988, + 291.48665471, + 293.6243954, + 295.68722534, + 297.66562966, + 299.55018853, + 301.33153232, + 302.9999648, + 304.54520748, + 305.95710135, + 307.22771897, + 308.35404538, + 309.33970039, + 310.19516516, + 310.935894, + 311.57873857, + 312.143772, + 312.65796099, + 313.14933519, + 313.64039782, + 314.1467086, + 314.6901326, + 315.29582886, + 315.98700684, + 316.78469962, + 317.70901078, + 318.78001191, + 320.01762104, + 321.44086311, + 323.06709059, + 324.91143497, + 326.98650416, + 329.30226115, + 331.86603447, + 334.68264398, + 337.75463976, + 341.08263885, + 344.66571751, + 348.50179783, + 352.58797567, + 356.92076582, + 361.49627286, + 366.31031192, + 371.3585025, + 376.63634824, + ], + [ + 278.30634695, + 280.76692289, + 283.20206879, + 285.60688807, + 287.97587967, + 290.30282905, + 292.58077049, + 294.80208911, + 296.95878792, + 299.04284976, + 301.04652031, + 302.96231809, + 304.78273895, + 306.49997116, + 308.10624416, + 309.59524733, + 310.96417741, + 312.21517992, + 313.35532245, + 314.39511557, + 315.34669749, + 316.22485477, + 317.04972335, + 317.84500607, + 318.63352456, + 319.43547086, + 320.26866588, + 321.15491544, + 322.11728356, + 323.17735523, + 324.3552376, + 325.67016149, + 327.14083776, + 328.78527246, + 330.6202218, + 332.6605961, + 334.91901704, + 337.40559839, + 340.12794395, + 343.09132309, + 346.29896643, + 349.7524136, + 353.45185087, + 357.39639806, + 361.58433378, + 366.01326948, + 370.68028972, + 375.58207254, + 380.71499591, + 386.07523169, + ], + [ + 279.75444778, + 282.29817944, + 284.82433455, + 287.32902789, + 289.80791774, + 292.25612472, + 294.66822243, + 297.0383542, + 299.36048536, + 301.62872141, + 303.83755132, + 305.98187696, + 308.05680937, + 310.05743003, + 311.97889683, + 313.81720263, + 315.57044276, + 317.23986013, + 318.82984388, + 320.3468679, + 321.79853724, + 323.19410279, + 324.5458906, + 325.8693588, + 327.18153051, + 328.49985457, + 329.84235722, + 331.22807789, + 332.67693979, + 334.20933551, + 335.84483513, + 337.60187632, + 339.49779226, + 341.54879485, + 343.76981952, + 346.17429027, + 348.77391483, + 351.57858366, + 354.59638447, + 357.83370121, + 361.29535632, + 364.98476646, + 368.9040978, + 373.05441794, + 377.43584441, + 382.04768885, + 386.88859397, + 391.9566601, + 397.24955922, + 402.76463545, + ], + [ + 281.83443065, + 284.49054627, + 287.13869075, + 289.77604486, + 292.3994679, + 295.00544884, + 297.59012963, + 300.14943366, + 302.67928723, + 305.17586241, + 307.63573556, + 310.05588235, + 312.43352359, + 314.76595749, + 317.05058504, + 319.28527574, + 321.46901173, + 323.60249084, + 325.68827882, + 327.73038804, + 329.73373947, + 331.70423537, + 333.64955593, + 335.5798509, + 337.50755558, + 339.44657727, + 341.4116138, + 343.41789639, + 345.48108684, + 347.61704308, + 349.84146682, + 352.169554, + 354.61583061, + 357.19410582, + 359.91740391, + 362.79786554, + 365.84662864, + 369.07371357, + 372.48793503, + 376.09685442, + 379.90677609, + 383.92278397, + 388.14881086, + 392.58773063, + 397.24146404, + 402.11109043, + 407.19695909, + 412.49879649, + 418.01580673, + 423.74676402, + ], + [ + 284.46010614, + 287.24808249, + 290.03809343, + 292.82824725, + 295.61643858, + 298.40033711, + 301.17744368, + 303.94522119, + 306.70126758, + 309.44346417, + 312.1700312, + 314.87946451, + 317.57039544, + 320.24147407, + 322.89138567, + 325.51905575, + 328.12400395, + 330.70672155, + 333.2689251, + 335.81360214, + 338.34489257, + 340.86796442, + 343.38903245, + 345.91549948, + 348.45602389, + 351.02034174, + 353.61888872, + 356.26243929, + 358.96193526, + 361.72849507, + 364.5734711, + 367.50843626, + 370.54506747, + 373.69496089, + 376.9694292, + 380.37931828, + 383.93485464, + 387.64552617, + 391.51999386, + 395.56603159, + 399.79049255, + 404.19930077, + 408.79746559, + 413.58911618, + 418.57755178, + 423.76530373, + 429.15420499, + 434.74546357, + 440.53973699, + 446.53720545, + ], + [ + 287.54370173, + 290.47477034, + 293.41743992, + 296.37055692, + 299.33284669, + 302.30293816, + 305.27943869, + 308.26104366, + 311.24664098, + 314.23536353, + 317.22656351, + 320.21972203, + 323.21434555, + 326.20991425, + 329.20592964, + 332.20206785, + 335.19840422, + 338.19565227, + 341.19536186, + 344.20003798, + 347.21316276, + 350.23912622, + 353.28309478, + 356.35086517, + 359.44874828, + 362.58349711, + 365.76225958, + 368.99253217, + 372.28211442, + 375.63908334, + 379.07179464, + 382.58888913, + 386.19926861, + 389.9120183, + 393.7362794, + 397.68109447, + 401.75525243, + 405.96715298, + 410.32469917, + 414.8352191, + 419.50541343, + 424.34132465, + 429.34832416, + 434.5311145, + 439.89374434, + 445.43963403, + 451.17161001, + 457.09194579, + 463.2024078, + 469.50430417, + ], + [ + 291.00408259, + 294.08333456, + 297.18293362, + 300.3023086, + 303.44085231, + 306.59796523, + 309.77312115, + 312.96593133, + 316.17617866, + 319.40380478, + 322.64885629, + 325.91141812, + 329.19156995, + 332.48939224, + 335.80502725, + 339.13878159, + 342.4912472, + 345.86342063, + 349.25680717, + 352.67350208, + 356.11624126, + 359.58841294, + 363.09402366, + 366.63762115, + 370.22418891, + 373.85903781, + 377.54771937, + 381.29597333, + 385.10970762, + 388.99499966, + 392.95810832, + 397.00548825, + 401.14379812, + 405.37989264, + 409.72078858, + 414.17360117, + 418.74545536, + 423.4433829, + 428.2742184, + 433.24450557, + 438.36042029, + 443.62771334, + 449.05167197, + 454.63709839, + 460.38830223, + 466.30910437, + 472.40284992, + 478.67242829, + 485.12029884, + 491.74852062, + ], + ] + ) + + # dome data for Ts diagram + sdome = np.array( + [ + 1001.37898699, + 1046.06882698, + 1088.21416035, + 1128.86481841, + 1169.09596369, + 1210.24526707, + 1254.66210399, + 1309.32653659, + 1323.58218729, + 1340.74189756, + 1364.61655928, + 1400.86587985, + 1468.75018408, + 1516.27446567, + 1546.66849166, + 1567.63482054, + 1584.41961907, + 1643.56442582, + 1686.22125974, + 1722.06785103, + 1754.37017763, + 1784.81415615, + 1814.51645442, + 1844.39848637, + ] + ) + Tdome = np.array( + [ + 273.31081538, + 278.44972407, + 283.13044128, + 287.43392381, + 291.41872472, + 295.12790099, + 298.59249754, + 301.8325153, + 302.45422394, + 303.06687319, + 303.6699029, + 304.08533585, + 304.08533585, + 303.6699029, + 303.06687319, + 302.45422394, + 301.8325153, + 298.59249754, + 295.12790099, + 291.41872472, + 287.43392381, + 283.13044128, + 278.44972407, + 273.31081538, + ] + ) + + # Ts at turbine inlet + turbinletT = self.Linear_production_temperature[-1] + turbinletS = self.s_prod[-1] + + # Ts at turbine outlet + turbineoutletT = self.T_turbine_out_actual[-1] + turbineoutletS = interpn( + (self.Pvector, self.Tvector), + self.entropy, + np.array([self.Turbine_outlet_pressure * 1e5, turbineoutletT + 273.15]), + )[0] + + # Ts at compressor inlet + compressorinletT = self.Pre_cooling_temperature + compressorinletS = interpn( + (self.Pvector, self.Tvector), + self.entropy, + np.array([self.Turbine_outlet_pressure * 1e5, compressorinletT + 273.15]), + )[0] + + # Ts at compressor outlet + compressoroutletT = self.Post_compressor_T_actual[0] + compressoroutletS = interpn( + (self.Pvector, self.Tvector), self.entropy, np.array([self.P_in, compressoroutletT + 273.15]) + )[0] + + # Ts at re-injection + reinjectionT = self.T_in + reinjectionS = interpn( + (self.Pvector, self.Tvector), self.entropy, np.array([self.P_in, reinjectionT + 273.15]) + )[0] + + cycleTarray = np.array( + [turbinletT, turbineoutletT, compressorinletT, compressoroutletT, reinjectionT, turbinletT] + ) + cycleSarray = np.array( + [turbinletS, turbineoutletS, compressorinletS, compressoroutletS, reinjectionS, turbinletS] + ) + + # + plt.figure() + plt.title('CO2 TS-diagram') + plt.xlabel('Entropy [J/kg/K]') + plt.ylabel('Temperature [K]') + for i in range(len(pvectorforTS)): + plt.plot(svectorforTS, tmatrix[i, :], 'b-') + plt.plot(sdome, Tdome, 'r-') + plt.plot(cycleSarray, cycleTarray + 273.15, 'k-+') + # plt.axis([0, 200,0,50]) + plt.show() + + def lineplot(self, parameterrange, parameterindex): + print(' ') + print('### Line Plot Calculation ###') + + # create arrays to store result + lcoerange = np.zeros(len(parameterrange)) + avekwthrange = np.zeros(len(parameterrange)) + avekwerange = np.zeros(len(parameterrange)) + + # 0 = Flow rate; 1 = Horizontal length; 2 = Depth; 3 = Gradient; 4 = Diameter; 5 = Injection temperature; 6 = Rock thermal conductivity + # store originals + flow_original = self.Flow_user + hor_length_original = self.Hor_length_user + depth_original = self.Depth_user + gradient_original = self.Gradient_user + diameter_original = self.Diameter_user + Tin_original = self.Tin_user + krock_original = self.krock_user + + # verify input + if parameterindex == 0: + if min(parameterrange) < 5 or max(parameterrange) > 100: + print('Error: Flow rate outise allowable range [5 to 100 kg/s]') + sys.exit() + elif parameterindex == 1: + if min(parameterrange) < 1000 or max(parameterrange) > 20000: + print('Error: Horizontal length outise allowable range [1,000 to 20,000 m]') + sys.exit() + elif parameterindex == 2: + if min(parameterrange) < 1000 or max(parameterrange) > 5000: + print('Error: Depth outise allowable range [1,000 to 5,000 m]') + sys.exit() + elif parameterindex == 3: + if min(parameterrange) < 0.03 or max(parameterrange) > 0.07: + print('Error: Gradient outside allowable range [0.03 to 0.07 deg.C/m]') + sys.exit() + elif parameterindex == 4: + if min(parameterrange) < 0.2159 or max(parameterrange) > 0.4445: + print('Error: Diameter outise allowable range [0.2159 to 0.4445 m]') + sys.exit() + elif parameterindex == 5: + if min(parameterrange) < 303.15 or max(parameterrange) > 333.15: + print('Error: Injection temperature outside allowable range [303.15 to 333.15 K]') + sys.exit() + elif parameterindex == 6: + if min(parameterrange) < 1.5 or max(parameterrange) > 4.5: + print('Error: Rock thermal conductivity outise allowable range [1.5 to 4.5 W/m-K]') + sys.exit() + + # cycle through array to update object and recalculated TEA results, then store them in arrays + for i in range(len(parameterrange)): + print('Line plot calculation: ' + str(i + 1) + '/' + str(len(parameterrange))) + if parameterindex == 0: + self.Flow_user = parameterrange[i] + xlabel = 'Flow rate [kg/s]' + elif parameterindex == 1: + self.Hor_length_user = parameterrange[i] + xlabel = 'Horizontal length [m]' + elif parameterindex == 2: + self.Depth_user = parameterrange[i] + xlabel = 'Depth [m]' + elif parameterindex == 3: + self.Gradient_user = parameterrange[i] + xlabel = 'Geothermal gradient [deg.C/m]' + elif parameterindex == 4: + self.Diameter_user = parameterrange[i] + xlabel = 'Diameter [m]' + elif parameterindex == 5: + self.Tin_user = parameterrange[i] + self.T_in = self.Tin_user - 273.15 # Injection temperature [deg.C] + xlabel = 'Injection temperature [deg.C]' + elif parameterindex == 6: + self.krock_user = parameterrange[i] + xlabel = 'Rock thermal conductivity [W/m-K]' + + # update self.point + self.point = ( + self.Flow_user, + self.Hor_length_user, + self.Depth_user, + self.Gradient_user, + self.Diameter_user, + self.Tin_user, + self.krock_user, + ) + + self.getTandP() + self.calculateLC() + if np.in1d(1000, self.error_codes) == False: + avekwthrange[i] = self.AveInstHeatProduction + if self.End_use == 1: + lcoerange[i] = self.LCOH + else: + lcoerange[i] = self.LCOE + avekwerange[i] = self.AveInstNetElectricityProduction + else: + lcoerange[i] = 9999 + + # restore object to its original condition + self.Flow_user = flow_original + self.Hor_length_user = hor_length_original + self.Depth_user = depth_original + self.Gradient_user = gradient_original + self.Diameter_user = diameter_original + self.Tin_user = Tin_original + self.T_in = self.Tin_user - 273.15 + self.krock_user = krock_original + self.point = ( + self.Flow_user, + self.Hor_length_user, + self.Depth_user, + self.Gradient_user, + self.Diameter_user, + self.Tin_user, + self.krock_user, + ) + + # create 1D line plots + plt.figure() + plt.plot(parameterrange, lcoerange, color='blue', linestyle='-', linewidth='2') + plt.grid() + plt.xlim([parameterrange[0], parameterrange[-1]]) + plt.ylim([0, 1.05 * np.max(lcoerange)]) + plt.xlabel(xlabel) + if self.End_use == 1: + plt.ylabel('Levelized cost of heat [$/MWh]') + else: + plt.ylabel('Levelized cost of electricity [$/MWh]') + plt.show() + + plt.figure() + plt.plot(parameterrange, avekwthrange / 1000, color='blue', linestyle='-', linewidth='2') + plt.grid() + plt.xlim([parameterrange[0], parameterrange[-1]]) + plt.ylim([0, 1.05 * np.max(avekwthrange / 1000)]) + plt.xlabel(xlabel) + plt.ylabel('Average heat production [$MW_{t}$]') + plt.show() + + if self.End_use == 2: + plt.figure() + plt.plot(parameterrange, avekwerange, color='blue', linestyle='-', linewidth='2') + plt.grid() + plt.xlim([parameterrange[0], parameterrange[-1]]) + plt.ylim([0, 1.05 * np.max(avekwerange)]) + plt.xlabel(xlabel) + plt.ylabel('Average electricity production [$kW_{e}$]') + plt.show() + + def contourplot(self, parameter1range, parameter1index, parameter2range, parameter2index, newtitle): + print(' ') + print('### Contour Plot Calculation ###') + + # create arrays to store result + lcoerange = np.zeros((len(parameter1range), len(parameter2range))) + lcoerange100 = np.zeros((len(parameter1range), len(parameter2range))) + lcoerange500 = np.zeros((len(parameter1range), len(parameter2range))) + lcoerange1500 = np.zeros((len(parameter1range), len(parameter2range))) + avekwthrange = np.zeros((len(parameter1range), len(parameter2range))) + avekwerange = np.zeros((len(parameter1range), len(parameter2range))) + + # 0 = Flow rate; 1 = Horizontal length; 2 = Depth; 3 = Gradient; 4 = Diameter; 5 = Injection temperature; 6 = Rock thermal conductivity + # store originals + flow_original = self.Flow_user + hor_length_original = self.Hor_length_user + depth_original = self.Depth_user + gradient_original = self.Gradient_user + diameter_original = self.Diameter_user + Tin_original = self.Tin_user + krock_original = self.krock_user + + # verify input + if parameter1index == 0: + if min(parameter1range) < 5 or max(parameter1range) > 100: + print('Error: Flow rate outise allowable range [5 to 100 kg/s]') + sys.exit() + elif parameter1index == 1: + if min(parameter1range) < 1000 or max(parameter1range) > 20000: + print('Error: Horizontal length outise allowable range [1,000 to 20,000 m]') + sys.exit() + elif parameter1index == 2: + if min(parameter1range) < 1000 or max(parameter1range) > 5000: + print('Error: Depth outise allowable range [1,000 to 5,000 m]') + sys.exit() + elif parameter1index == 3: + if min(parameter1range) < 0.03 or max(parameter1range) > 0.07: + print('Error: Gradient outside allowable range [0.03 to 0.07 deg.C/m]') + sys.exit() + elif parameter1index == 4: + if min(parameter1range) < 0.2159 or max(parameter1range) > 0.4445: + print('Error: Diameter outise allowable range [0.2159 to 0.4445 m]') + sys.exit() + elif parameter1index == 5: + if min(parameter1range) < 303.15 or max(parameter1range) > 333.15: + print('Error: Injection temperature outside allowable range [303.15 to 333.15 K]') + sys.exit() + elif parameter1index == 6: + if min(parameter1range) < 1.5 or max(parameter1range) > 4.5: + print('Error: Rock thermal conductivity outise allowable range [1.5 to 4.5 W/m-K]') + sys.exit() + if parameter2index == 0: + if min(parameter2range) < 5 or max(parameter2range) > 100: + print('Error: Flow rate outise allowable range [5 to 100 kg/s]') + sys.exit() + elif parameter2index == 1: + if min(parameter2range) < 1000 or max(parameter2range) > 20000: + print('Error: Horizontal length outise allowable range [1,000 to 20,000 m]') + sys.exit() + elif parameter2index == 2: + if min(parameter2range) < 1000 or max(parameter2range) > 5000: + print('Error: Depth outise allowable range [1,000 to 5,000 m]') + sys.exit() + elif parameter2index == 3: + if min(parameter2range) < 0.03 or max(parameter2range) > 0.07: + print('Error: Gradient outside allowable range [0.03 to 0.07 deg.C/m]') + sys.exit() + elif parameter2index == 4: + if min(parameter2range) < 0.2159 or max(parameter2range) > 0.4445: + print('Error: Diameter outise allowable range [0.2159 to 0.4445 m]') + sys.exit() + elif parameter2index == 5: + if min(parameter2range) < 303.15 or max(parameter2range) > 333.15: + print('Error: Injection temperature outside allowable range [303.15 to 333.15 K]') + sys.exit() + elif parameter2index == 6: + if min(parameter2range) < 1.5 or max(parameter2range) > 4.5: + print('Error: Rock thermal conductivity outise allowable range [1.5 to 4.5 W/m-K]') + sys.exit() + + # cycle through double array to update object and recalculated TEA results, then store them in arrays + for i in range(len(parameter1range)): + if parameter1index == 0: + self.Flow_user = parameter1range[i] + xlabel = 'Flow rate [kg/s]' + elif parameter1index == 1: + self.Hor_length_user = parameter1range[i] + xlabel = 'Horizontal length [m]' + elif parameter1index == 2: + self.Depth_user = parameter1range[i] + xlabel = 'Depth [m]' + elif parameter1index == 3: + self.Gradient_user = parameter1range[i] + xlabel = 'Geothermal gradient [deg.C/m]' + elif parameter1index == 4: + self.Diameter_user = parameter1range[i] + xlabel = 'Diameter [m]' + elif parameter1index == 5: + self.Tin_user = parameter1range[i] + self.T_in = self.Tin_user - 273.15 # Injection temperature [deg.C] + xlabel = 'Injection temperature [deg.C]' + elif parameter1index == 6: + self.krock_user = parameter1range[i] + + for j in range(len(parameter2range)): + print( + 'Contour plot calculation: ' + + str(j + 1 + i * len(parameter2range)) + + '/' + + str(len(parameter1range) * len(parameter2range)) + ) + if parameter2index == 0: + self.Flow_user = parameter2range[j] + ylabel = 'Flow rate [kg/s]' + elif parameter2index == 1: + self.Hor_length_user = parameter2range[j] + ylabel = 'Horizontal length [m]' + elif parameter2index == 2: + self.Depth_user = parameter2range[j] + ylabel = 'Depth [m]' + elif parameter2index == 3: + self.Gradient_user = parameter2range[j] + ylabel = 'Geothermal gradient [deg.C/m]' + elif parameter2index == 4: + self.Diameter_user = parameter2range[j] + ylabel = 'Diameter [m]' + elif parameter2index == 5: + self.Tin_user = parameter2range[j] + self.T_in = self.Tin_user - 273.15 # Injection temperature [deg.C] + ylabel = 'Injection temperature [deg.C]' + elif parameter2index == 6: + self.krock_user = parameter2range[j] + + # update self.point + self.point = ( + self.Flow_user, + self.Hor_length_user, + self.Depth_user, + self.Gradient_user, + self.Diameter_user, + self.Tin_user, + self.krock_user, + ) + self.error_codes = np.zeros(0) + self.getTandP() + self.calculateLC() + if np.in1d(1000, self.error_codes) == False: + avekwthrange[(i, j)] = self.AveInstHeatProduction + + if self.End_use == 1: + lcoerange[(i, j)] = self.LCOH + lcoerange100[(i, j)] = self.LCOH100 + lcoerange500[(i, j)] = self.LCOH500 + lcoerange1500[(i, j)] = self.LCOH1500 + else: + lcoerange[(i, j)] = self.LCOE + lcoerange100[(i, j)] = self.LCOE100 + lcoerange500[(i, j)] = self.LCOE500 + lcoerange1500[(i, j)] = self.LCOE1500 + avekwerange[(i, j)] = self.AveInstNetElectricityProduction + else: + lcoerange[(i, j)] = float('nan') + lcoerange100[(i, j)] = float('nan') + lcoerange500[(i, j)] = float('nan') + lcoerange1500[(i, j)] = float('nan') + + lcoerange100[lcoerange100 > 1000] = float('nan') + lcoerange500[lcoerange500 > 2000] = float('nan') + lcoerange1500[lcoerange1500 > 5000] = float('nan') + + self.avekwthrange = avekwthrange + self.lcoerange = lcoerange + self.lcoerange100 = lcoerange100 + self.lcoerange500 = lcoerange500 + self.lcoerange1500 = lcoerange1500 + + # restore object to its original condition + self.Flow_user = flow_original + self.Hor_length_user = hor_length_original + self.Depth_user = depth_original + self.Gradient_user = gradient_original + self.Diameter_user = diameter_original + self.Tin_user = Tin_original + self.T_in = self.Tin_user - 273.15 + self.krock_user = krock_original + self.point = ( + self.Flow_user, + self.Hor_length_user, + self.Depth_user, + self.Gradient_user, + self.Diameter_user, + self.Tin_user, + self.krock_user, + ) + + # make plots + X, Y = np.meshgrid(parameter1range, parameter2range) + fig, ax = plt.subplots() + contour = plt.contourf(X, Y, np.transpose(avekwthrange), levels=32, cmap=cm.coolwarm) + ax.set(xlim=(X.min(), X.max()), ylim=(Y.min(), Y.max())) + ax.set_xlabel(xlabel) + ax.set_ylabel(ylabel) + plt.title('Average Heat Production') + cbar = plt.colorbar(contour) + cbar.set_label('Average Heat Production [kWth]') + plt.show() + + newtitle100 = newtitle + str(100) + '/m' + newtitle500 = newtitle + str(500) + '/m' + newtitle1500 = newtitle + str(1500) + '/m' + + # fig, ax = plt.subplots() + # contour = plt.contourf(X,Y, np.transpose(lcoerange), levels = 32, cmap=cm.coolwarm) + # ax.set(xlim=(X.min(), X.max()),ylim=(Y.min(), Y.max())) + # ax.set_xlabel(xlabel) + # ax.set_ylabel(ylabel) + # plt.title(newtitle) + # cbar = plt.colorbar(contour) + # if self.End_use == 1: + # cbar.set_label('LCOH [$/MWh]') + # elif self.End_use == 2: + # cbar.set_label('LCOE [$/MWh]') + # plt.show() + + fig, ax = plt.subplots() + contour = plt.contourf(X, Y, np.transpose(lcoerange100), levels=32, cmap=cm.coolwarm) + ax.set(xlim=(X.min(), X.max()), ylim=(Y.min(), Y.max())) + ax.set_xlabel(xlabel) + ax.set_ylabel(ylabel) + # if self.End_use == 1: + # newtile = 'U-loop | water | $'+str(self.Drilling_cost_per_m)+'/m' + # plt.title(newtile) + # elif self.End_use == 2: + plt.title(newtitle100) + cbar = plt.colorbar(contour) + if self.End_use == 1: + cbar.set_label('LCOH [$/MWh]') + elif self.End_use == 2: + cbar.set_label('LCOE [$/MWh]') + plt.show() + + fig, ax = plt.subplots() + contour = plt.contourf(X, Y, np.transpose(lcoerange500), levels=32, cmap=cm.coolwarm) + ax.set(xlim=(X.min(), X.max()), ylim=(Y.min(), Y.max())) + ax.set_xlabel(xlabel) + ax.set_ylabel(ylabel) + # if self.End_use == 1: + # newtile = 'U-loop | water | $'+str(self.Drilling_cost_per_m)+'/m' + # plt.title(newtile) + # elif self.End_use == 2: + plt.title(newtitle500) + cbar = plt.colorbar(contour) + if self.End_use == 1: + cbar.set_label('LCOH [$/MWh]') + elif self.End_use == 2: + cbar.set_label('LCOE [$/MWh]') + plt.show() + + fig, ax = plt.subplots() + contour = plt.contourf(X, Y, np.transpose(lcoerange1500), levels=32, cmap=cm.coolwarm) + ax.set(xlim=(X.min(), X.max()), ylim=(Y.min(), Y.max())) + ax.set_xlabel(xlabel) + ax.set_ylabel(ylabel) + # if self.End_use == 1: + # newtile = 'U-loop | water | $'+str(self.Drilling_cost_per_m)+'/m' + # plt.title(newtile) + # elif self.End_use == 2: + + plt.title(newtitle1500) + cbar = plt.colorbar(contour) + if self.End_use == 1: + cbar.set_label('LCOH [$/MWh]') + elif self.End_use == 2: + cbar.set_label('LCOE [$/MWh]') + plt.show() + + if self.End_use == 2: + fig, ax = plt.subplots() + contour = plt.contourf(X, Y, np.transpose(avekwerange), levels=32, cmap=cm.coolwarm) + ax.set(xlim=(X.min(), X.max()), ylim=(Y.min(), Y.max())) + ax.set_xlabel(xlabel) + ax.set_ylabel(ylabel) + plt.title('Average Electricity Production') + cbar = plt.colorbar(contour) + cbar.set_label('Average Electricity Production [kWe]') + plt.show() + + def optimalLCOH(self): + # parameters to vary: FlowRateVector,HorizontalLengthVector,DepthVector, , TinVector, KrockVector + + # store originals + flow_original = self.Flow_user + hor_length_original = self.Hor_length_user + depth_original = self.Depth_user + gradient_original = self.Gradient_user + diameter_original = self.Diameter_user + Tin_original = self.Tin_user + krock_original = self.krock_user + + # create optimum vectors + flowratevectorforoptimum = self.FlowRateVector + horizontallengthvectorforoptimum = self.HorizontalLengthVector + depthvectorforoptimum = self.DepthVector + tinvectorforoptimum = self.TinVector + + # create zero vectors to store results + avekwthrange = np.zeros( + [ + len(flowratevectorforoptimum), + len(horizontallengthvectorforoptimum), + len(depthvectorforoptimum), + len(tinvectorforoptimum), + ] + ) + lcoerange = np.zeros( + [ + len(flowratevectorforoptimum), + len(horizontallengthvectorforoptimum), + len(depthvectorforoptimum), + len(tinvectorforoptimum), + ] + ) + avekwerange = np.zeros( + [ + len(flowratevectorforoptimum), + len(horizontallengthvectorforoptimum), + len(depthvectorforoptimum), + len(tinvectorforoptimum), + ] + ) + + # cycle through array to update object and recalculated TEA results, then store them in arrays + for i in range(len(flowratevectorforoptimum)): + print('Optimum calculation: ' + str(i + 1) + '/' + str(len(flowratevectorforoptimum))) + self.Flow_user = flowratevectorforoptimum[i] + + for j in range(len(horizontallengthvectorforoptimum)): + self.Hor_length_user = horizontallengthvectorforoptimum[j] + + for k in range(len(depthvectorforoptimum)): + self.Depth_user = depthvectorforoptimum[k] + + for l in range(len(tinvectorforoptimum)): + self.Tin_user = tinvectorforoptimum[l] + self.T_in = self.Tin_user - 273.15 + + # elif parameterindex == 2: + # self.Depth_user = parameterrange[i] + # xlabel = 'Depth [m]' + # elif parameterindex == 3: + # self.Gradient_user = parameterrange[i] + # xlabel = 'Geothermal gradient [deg.C/m]' + # elif parameterindex == 4: + # self.Diameter_user = parameterrange[i] + # xlabel = 'Diameter [m]' + # elif parameterindex == 5: + # self.Tin_user = parameterrange[i] + # self.T_in = self.Tin_user-273.15 #Injection temperature [deg.C] + # xlabel = 'Injection temperature [deg.C]' + # elif parameterindex == 6: + # self.krock_user = parameterrange[i] + # xlabel = 'Rock thermal conductivity [W/m-K]' + + # update self.point + self.point = ( + self.Flow_user, + self.Hor_length_user, + self.Depth_user, + self.Gradient_user, + self.Diameter_user, + self.Tin_user, + self.krock_user, + ) + + self.getTandP() + self.calculateLC() + + if np.in1d(1000, self.error_codes) == False: + avekwthrange[i, j, k, l] = self.AveInstHeatProduction + if self.End_use == 1: + lcoerange[i, j, k, l] = self.LCOH + else: + lcoerange[i, j, k, l] = self.LCOE + avekwerange[i, j, k, l] = self.AveInstNetElectricityProduction + else: + lcoerange[i, j, k, l] = 9999 + + # store results + self.avekwthrange = avekwthrange + self.lcoerange = lcoerange + self.avekwerange = avekwerange + + # restore object to its original condition + self.Flow_user = flow_original + self.Hor_length_user = hor_length_original + self.Depth_user = depth_original + self.Gradient_user = gradient_original + self.Diameter_user = diameter_original + self.Tin_user = Tin_original + self.T_in = self.Tin_user - 273.15 + self.krock_user = krock_original + self.point = ( + self.Flow_user, + self.Hor_length_user, + self.Depth_user, + self.Gradient_user, + self.Diameter_user, + self.Tin_user, + self.krock_user, + ) + + def optimalLCOH2(self): + # parameters to vary: FlowRateVector,HorizontalLengthVector,DepthVector, , TinVector, KrockVector + + # store originals + flow_original = self.Flow_user + hor_length_original = self.Hor_length_user + depth_original = self.Depth_user + gradient_original = self.Gradient_user + diameter_original = self.Diameter_user + Tin_original = self.Tin_user + krock_original = self.krock_user + + # create optimum vectors + flowratevectorforoptimum = self.FlowRateVector + horizontallengthvectorforoptimum = self.HorizontalLengthVector + tinvectorforoptimum = self.TinVector + + # create zero vectors to store results + avekwthrange = np.zeros( + [len(flowratevectorforoptimum), len(horizontallengthvectorforoptimum), len(tinvectorforoptimum)] + ) + lcoerange = np.zeros( + [len(flowratevectorforoptimum), len(horizontallengthvectorforoptimum), len(tinvectorforoptimum)] + ) + avekwerange = np.zeros( + [len(flowratevectorforoptimum), len(horizontallengthvectorforoptimum), len(tinvectorforoptimum)] + ) + + # cycle through array to update object and recalculated TEA results, then store them in arrays + for i in range(len(flowratevectorforoptimum)): + print('Optimum calculation: ' + str(i + 1) + '/' + str(len(flowratevectorforoptimum))) + self.Flow_user = flowratevectorforoptimum[i] + + for j in range(len(horizontallengthvectorforoptimum)): + self.Hor_length_user = horizontallengthvectorforoptimum[j] + + for l in range(len(tinvectorforoptimum)): + self.Tin_user = tinvectorforoptimum[l] + self.T_in = self.Tin_user - 273.15 + + # elif parameterindex == 2: + # self.Depth_user = parameterrange[i] + # xlabel = 'Depth [m]' + # elif parameterindex == 3: + # self.Gradient_user = parameterrange[i] + # xlabel = 'Geothermal gradient [deg.C/m]' + # elif parameterindex == 4: + # self.Diameter_user = parameterrange[i] + # xlabel = 'Diameter [m]' + # elif parameterindex == 5: + # self.Tin_user = parameterrange[i] + # self.T_in = self.Tin_user-273.15 #Injection temperature [deg.C] + # xlabel = 'Injection temperature [deg.C]' + # elif parameterindex == 6: + # self.krock_user = parameterrange[i] + # xlabel = 'Rock thermal conductivity [W/m-K]' + + # update self.point + self.point = ( + self.Flow_user, + self.Hor_length_user, + self.Depth_user, + self.Gradient_user, + self.Diameter_user, + self.Tin_user, + self.krock_user, + ) + + self.getTandP() + self.calculateLC() + + if np.in1d(1000, self.error_codes) == False: + avekwthrange[i, j, l] = self.AveInstHeatProduction + if self.End_use == 1: + lcoerange[i, j, l] = self.LCOH + else: + lcoerange[i, j, l] = self.LCOE + avekwerange[i, j, l] = self.AveInstNetElectricityProduction + else: + lcoerange[i, j, l] = 9999 + + # store results + self.avekwthrange = avekwthrange + self.lcoerange = lcoerange + self.avekwerange = avekwerange + + # restore object to its original condition + self.Flow_user = flow_original + self.Hor_length_user = hor_length_original + self.Depth_user = depth_original + self.Gradient_user = gradient_original + self.Diameter_user = diameter_original + self.Tin_user = Tin_original + self.T_in = self.Tin_user - 273.15 + self.krock_user = krock_original + self.point = ( + self.Flow_user, + self.Hor_length_user, + self.Depth_user, + self.Gradient_user, + self.Diameter_user, + self.Tin_user, + self.krock_user, + ) diff --git a/src/fgem/utils/utils.py b/src/fgem/utils/utils.py new file mode 100755 index 000000000..6ed9ff189 --- /dev/null +++ b/src/fgem/utils/utils.py @@ -0,0 +1,920 @@ +import os +from datetime import datetime +from datetime import timedelta + +import matplotlib.pyplot as plt +import numpy as np +import pandas as pd +import seaborn as sns +from meteostat import Hourly +from meteostat import Stations +from pyXSteam.XSteam import XSteam +from scipy.optimize import curve_fit +from shapely.geometry import Polygon +from shapely.geometry.polygon import orient +from shapely.prepared import prep +from timezonefinder import TimezoneFinder + +steamTable = XSteam(XSteam.UNIT_SYSTEM_MKS) # m/kg/sec/°C/bar/W + +colors = 24 * sns.color_palette() + +linestyles = {0: 'solid', 1: 'dashed', 2: 'dashdot', 3: 'dotted'} + + +class FastXsteam: + """Faster corrleations for Steam Tables.""" + + def __init__(self, T_max=300, steps=30000): + """Initiate class for fast Xsteam tables. + + Args: + T_max (int, optional): maximum temperature. Defaults to 300. + steps (int, optional): number of steps. Defaults to 30000. + """ + + self.xdata = np.linspace(1, T_max, steps) + self.tc = 647.096 # K + + ydata = [steamTable.hL_t(x) for x in self.xdata] + self.popt_hl, pcov = curve_fit(self.func_hl, self.xdata, ydata) + ydata = [steamTable.hV_t(x) for x in self.xdata] + self.popt_hv, pcov = curve_fit(self.func_hv, self.xdata, ydata) + ydata = [steamTable.vL_t(x) for x in self.xdata] + self.popt_vl, pcov = curve_fit(self.func_vl, self.xdata, ydata) + ydata = [steamTable.vV_t(x) for x in self.xdata] + self.popt_vv, pcov = curve_fit(self.func_vv, self.xdata, ydata) + + def func_hl(self, t, a, b, c, d, e): + """Liquid enthalpy correlation. + + Args: + t (float): temperature in C + a (float): correlation coefficient. + b (float): correlation coefficient. + c (float): correlation coefficient. + d (float): correlation coefficient. + e (float): correlation coefficient. + + Returns: + float: liquid enthalpy in kJ/kg at the given temperature + """ + tr = (t + 273.15) / self.tc + return a + b * np.log(1 / tr) ** 0.35 + c / tr**2 + d / tr**3 + e / tr**4 + + def func_hv(self, t, a, b, c, d, e): + """Vapor enthalpy correlation. + + Args: + t (float): temperature in C + a (float): correlation coefficient. + b (float): correlation coefficient. + c (float): correlation coefficient. + d (float): correlation coefficient. + e (float): correlation coefficient. + + Returns: + float: vapor enthalpy in kJ/kg at the given temperature + """ + + tr = (t + 273.15) / self.tc + return (a + b * np.log(1 / tr) ** 0.35 + c / tr**2 + d / tr**3 + e / tr**4) ** (1 / 2) + + def func_vl(self, t, a, b, c, d, e): + """Liquid specific volume correlation. + + Args: + t (float): temperature in C + a (float): correlation coefficient. + b (float): correlation coefficient. + c (float): correlation coefficient. + d (float): correlation coefficient. + e (float): correlation coefficient. + + Returns: + float: liquid specific volume in m3/kg at the given temperature + """ + + tr = (t + 273.15) / self.tc + return a + b * np.log(1 / tr) ** 0.35 + c / tr**2 + d / tr**3 + e / tr**4 + + def func_vv(self, t, a, b, c, d, e): + """Vapor specific volume correlation. + + Args: + t (float): temperature in C + a (float): correlation coefficient. + b (float): correlation coefficient. + c (float): correlation coefficient. + d (float): correlation coefficient. + e (float): correlation coefficient. + + Returns: + float: vapor specific volume in m3/kg at the given temperature + """ + + tr = (t + 273.15) / self.tc + return np.exp(a + b * np.log(1 / tr) ** 0.35 + c / tr**2 + d / tr**3 + e / tr**4) + + def hL_t(self, t): + """Liquid enthalpy correlation. + + Args: + t (float): temperature in C + + Returns: + float: liquid enthalpy in kJ/kg at the given temperature + """ + t = t[0] if isinstance(t, np.ndarray) else t + return self.func_hl(t, *self.popt_hl) + + def hV_t(self, t): + """Vapor enthalpy correlation. + + Args: + t (float): temperature in C + + Returns: + float: vapor enthalpy in kJ/kg at the given temperature + """ + + t = t[0] if isinstance(t, np.ndarray) else t + return self.func_hv(t, *self.popt_hv) + + def vL_t(self, t): + """Liquid specific volume correlation. + + Args: + t (float): temperature in C + + Returns: + float: liquid specific volume in m3/kg at the given temperature + """ + + t = t[0] if isinstance(t, np.ndarray) else t + return self.func_vl(t, *self.popt_vl) + + def vV_t(self, t): + """Vapor specific volume correlation. + + Args: + t (float): temperature in C + + Returns: + float: vapor specific volume in m3/kg at the given temperature + """ + + t = t[0] if isinstance(t, np.ndarray) else t + return self.func_vv(t, *self.popt_vv) + + +def plot_ex(worlds, figsize=(10, 10), dpi=150, colors=colors, fontsize=10): + """Visualize expenditure of multiple projects. + + Args: + worlds (dict): projects or worlds to visualize + figsize (tuple, optional): figure size. Defaults to (10,10). + dpi (int, optional): figure dpi resolution. Defaults to 150. + colors (list, optional): list of colors to be used. Defaults to seaborn.color_palette(). + fontsize (int, optional): font size. Defaults to 10. + + Returns: + matplotlib.figure.Figure: figure object + """ + + fig, axes = plt.subplots(1, len(worlds), figsize=figsize, dpi=dpi) + axes = [axes] if len(worlds) == 1 else axes + + def func(pct, allvals): + absolute = int(np.round(pct / 100.0 * np.sum(allvals))) + # return "{:.0f}%\n({:d} $MM)".format(pct, absolute) + return f'{pct:.0f}%' + + for i, (title, present_per_unit) in enumerate(worlds.items()): + ex = [v for _, v in present_per_unit.items() if v > 0] + labels = [k for k, v in present_per_unit.items() if v > 0] + wedges, _, _ = axes[i].pie( + x=ex, + pctdistance=0.8, + labels=labels, + colors=colors, + autopct=lambda pct: func(pct, ex), + textprops=dict(color='w', weight='bold', fontsize=fontsize), + ) + + axes[i].legend(wedges, labels, title=title, loc='center left', bbox_to_anchor=(0.3, 0.6, 0.0, 1)) + + return fig + + +def plot_cols( + dfs, + span, + quantities, + figsize=(10, 10), + xlabel='World Time', + ylabels=None, + ylogscale=False, + use_title=True, + legend_loc='lower right', + manual_legends=None, + color_per_col=True, + use_linestyles=True, + blackout_first=False, + formattime=False, + dpi=100, +): + """Visualize columns of a dataframe. + + Args: + dfs (dict): dictionary of dataframes to plot with key as names. + span (Union[range, list]): range of indices to plot. + quantities (list): dataframe columns to plot. + figsize (tuple, optional): figure size. Defaults to (10,10). + xlabel (str, optional): label of the x axis. Defaults to "World Time". + ylabels (list, optional): labels of the y axes. Defaults to None. + ylogscale (bool, optional): whether or not to use log-scale on the y axis. Defaults to False. + use_title (bool, optional): whether or not to use a figure title. Defaults to True. + legend_loc (str, optional): legend location. Defaults to "lower right". + manual_legends (list, optional): list to override default legends. Defaults to None. + color_per_col (bool, optional): whether or not to assign one color to each column. Defaults to True. + use_linestyles (bool, optional): whether or not to use linestyles. If False, all quantities are plotted with solid linestyle. Defaults to True. + blackout_first (bool, optional): whether or not to plot the very first quantity in black color. Defaults to False. + formattime (bool, optional): whether or not to apply time formatting. Defaults to False. + dpi (int, optional): figure dpi resolution. Defaults to 100. + + Returns: + matplotlib.figure.Figure: figure object + + """ + + fig, axes = plt.subplots(len(quantities), 1, figsize=figsize, sharex=True, dpi=dpi) + + df_plots = [df.iloc[span].copy() for df in dfs.values()] + counter = 0 + + for k, df_plot in enumerate(df_plots): + for i, q in enumerate(quantities): + q = q if isinstance(q, list) else [q] + for col in q: + if all([qi in df_plot.columns for qi in q]): + # axes[i].plot(df_plot.index, df_plot[col].rolling(10).mean(), color=colors[counter], linestyle=linestyles[k if use_linestyles else 0]) + axes[i].plot( + df_plot.index, + df_plot[col], + color=colors[counter], + linestyle=linestyles[k if use_linestyles else 0], + ) + if color_per_col: + counter += 1 + if legend_loc: + axes[i].legend(q, loc=legend_loc) + if not color_per_col: + counter += 1 + else: + counter = 0 + + if use_title: + if len(dfs.keys()) > 1: + axes[0].set_title('\n'.join([f'{k}: {linestyles[i]}' for i, k in enumerate(dfs.keys())])) + else: + axes[0].set_title(list(dfs.keys())[0]) + + axes[i].set_xlabel(xlabel) + if formattime: + if 'time' in xlabel.lower() or 'date' in xlabel.lower(): + plt.gcf().autofmt_xdate() + + if ylabels: + for i, ylabel in enumerate(ylabels): + axes[i].set_ylabel(ylabel) + if ylogscale: + for i, log in enumerate(ylogscale): + if log: + axes[i].set_yscale('log') + if manual_legends: + axes[0].legend(list(dfs.keys()), loc='upper right', fontsize=10) + if blackout_first: + axes[0].plot(df_plot.index, df_plot[quantities[0]], color='black') + + plt.show() + + return fig + + +def viz_wholesale(filepath, span=range(2025, 2056, 2)): + """Visualize wholesale market data. + + Args: + filepath (str): file path to csv file. + span (range, optional): range or list of years to visualize. Defaults to range(2025,2056,2). + """ + + if isinstance(filepath, str): + df = pd.read_csv(filepath) + df.Date = pd.to_datetime(df.Date, format='%m/%d/%y %H:%M') + else: + df = filepath.copy() + + df['Year'], df['Hour'] = df.Date.dt.year, df.Date.dt.hour + df = pd.pivot_table(df, values='price', index=['Hour'], columns=['Year'], aggfunc=np.mean) + mappaple = df.groupby('Hour').mean().plot(y=span, cmap='jet') + plt.legend(ncol=4, loc='upper right') + plt.xlabel('Hour of Day') + plt.ylabel('Price [Nominal $/MWh]') + plt.show() + + +def preprocess_cambium_capacity( + scenario, + state, + infl, + base_year=2023, + span=range(2024, 2061), + plot=False, + dst_dir='.', + save=None, + adjustment=0, + states_with_no_capacity_market=['TX'], +): + """Preprocess and plot Cambium capacity market data for a chosen state and forecast scenario""" + + filepath = f'../../Cambium_2022/{scenario}/Cambium22_{scenario}_annual_state.csv' + df = pd.read_csv(filepath, header=5) + df = df.loc[df.state == state, ['t', 'capacity_cost_enduse']] + df = df.rename(columns={'t': 'Year', 'capacity_cost_enduse': 'capacity cost'}) + df['capacity cost'] *= 8760 / 1e3 # to convert $/MWh to $/kW-yr + + if state in states_with_no_capacity_market: + df['capacity cost'] = 0.0 + + df_full = pd.merge(pd.DataFrame({'Year': span}), df, how='outer', on='Year').interpolate(method='linear', order=1) + + if adjustment > 0: + df_full[ + 'capacity cost' + ] += adjustment # based on Fig 3-C of https://eta-publications.lbl.gov/sites/default/files/berkeley_lab_2021.11-_integrating_cambium_prices_into_electric-sector_decisions.pdf + + df_full['capacity cost'] = df_full['capacity cost'] * (1 + infl) ** (df_full.Year - base_year) + if plot: + # plt.scatter(df["Year"], df["capacity cost"], color="black") + plt.plot(df_full['Year'], df_full['capacity cost']) + plt.xlabel('Year') + plt.ylabel('Capacity Value [Nominal $/kW-yr]') + plt.show() + + if save: + df_full.to_csv(os.path.join(dst_dir, 'Capacity.csv'), index=False) + + return df_full + + +def preprocess_cambium_wholesale( + scenario, + state, + infl, + daily_fat_factor=1.0, + seasonal_fat_factor=1.0, + base_year=2022, + span=range(2024, 2061), + plot=False, + dst_dir='.', + save=False, +): + """Preprocess and plot Cambium wholesale market data for a chosen state and forecast scenario.""" + + base_dir = os.path.join('../../Cambium_2022/', scenario) + data_dir = os.path.join(base_dir, 'hourly_state') + filepaths = [os.path.join(data_dir, filename) for filename in os.listdir(data_dir) if f'_{state}_' in filename] + avail_years = np.sort([int(f.split('.csv')[0].split('_')[-1]) for f in filepaths]) + + df_recs = pd.read_csv(os.path.join(base_dir, 'RECs.csv')) + df_recs['Date'] = pd.to_datetime(df_recs['Date'], format='%m/%d/%y %H:%M') + df_recs = df_recs.loc[df_recs.Date.dt.year >= 2024] + + df = [] + for i, year in enumerate(avail_years): + filepath = os.path.join(data_dir, f'Cambium22_{scenario}_hourly_{state}_{year}.csv') + df_temp = pd.read_csv(filepath, header=5) + + if i == 0: + df_temp['Date'] = pd.to_datetime(df_temp['timestamp_local']) + df_temp['Month'] = df_temp.Date.dt.month + df_temp['Day'] = df_temp.Date.dt.day + df_temp['Hour'] = df_temp.Date.dt.hour + df_temp['DayofYear'] = df_temp.Date.dt.dayofyear + df_temp['HourofYear'] = np.arange(1, len(df_temp) + 1) + df_temp[year] = df_temp['energy_cost_enduse'].values + df = df_temp[['Month', 'Day', 'Hour', 'DayofYear', 'HourofYear', year]].copy() + else: + df[year] = df_temp['energy_cost_enduse'].values + + for i, year in enumerate(span): + + if year in avail_years: + continue + + if any(year < avail_years): + upper_years = avail_years[avail_years > year] + lower_years = avail_years[avail_years < year] + upper_year = upper_years[np.argmin(np.abs(year - upper_years))] + lower_year = lower_years[np.argmin(np.abs(lower_years - year))] + assert (year < upper_year) and (year > lower_year), 'Issue with lower/upper year allocation ...' + + dpdy = (df[upper_year] - df[lower_year]) / (upper_year - lower_year) + + df[year] = df[lower_year] + dpdy * (year - lower_year) + + for i, year in enumerate(span): + df[year] = df[year] * (1 + infl) ** (year - base_year) + + daily_price_means = np.zeros(len(df)) + for i in range(int(len(daily_price_means) / 24)): + daily_price_means[24 * i : 24 * i + 24] = df[year][24 * i : 24 * i + 24].mean() + + seasonal_price_means = np.zeros(len(df)) + for i in range(int(len(seasonal_price_means) / 8760)): + seasonal_price_means[8760 * i : 8760 * i + 8760] = df[year][8760 * i : 8760 * i + 8760].mean() + + df[year] += daily_fat_factor * (df[year] - daily_price_means) + seasonal_fat_factor * ( + df[year] - seasonal_price_means + ) + + df_melt = pd.melt(df, id_vars=['Hour', 'Day', 'Month'], value_vars=span, var_name='Year', value_name='price') + df_melt['Date'] = pd.to_datetime(df_melt[['Year', 'Month', 'Day', 'Hour']]) + df_final = pd.merge(df_melt, df_recs, how='inner', on='Date')[['Date', 'Hour', 'price', 'recs_price']] + + baseload_score = df_final.price.mean() # $/MWh + arbitrage_score = ( + df_final.groupby('Hour').mean(numeric_only=True)['price'].max() + - df_final.groupby('Hour').mean(numeric_only=True)['price'].min() + ) # $/MWh + scores = {'baseload': round(baseload_score, 2), 'arbitrage': round(arbitrage_score, 2)} + if plot: + viz_wholesale(df_final, span=span[::2]) + + df_final['Date'] = df_final.Date.dt.strftime('%m/%d/%y %H:%M') + if save: + df_final.to_csv(os.path.join(dst_dir, 'DA.csv'), index=False) + + return df_final, df, scores + + +def preprocess_cambium_elcc( + scenario, state, span=range(2024, 2061), plot=False, dst_dir='../../Data/dummy', save=False +): + """Preprocess and plot Cambium ELCC data for a chosen state and forecast scenario.""" + + filepath = f'../../Cambium_2022/{scenario}/Cambium22_{scenario}_annual_state.csv' + df = pd.read_csv(filepath, header=5) + df = df.loc[df.state == state, ['t', 'battery_MW']] + df = df.rename(columns={'t': 'Year', 'battery_MW': 'battery_GW'}) + df['battery_GW'] = df['battery_GW'] / 1e3 + df_full = pd.merge(pd.DataFrame({'Year': span}), df, how='outer', on='Year').interpolate(method='linear', order=1) + + # based on Figure 20 in https://www.ethree.com/wp-content/uploads/2019/06/E3_Long_Run_Resource_Adequacy_CA_Deep-Decarbonization_Final.pdf + d = {0: 1, 5: 0.80, 10: 0.60, 15: 0.40, 20: 0.35, 25: 0.30, 30: 0.20, 40: 0.11, 45: 0.10, 55: 0.06, 75: 0.05} + df_elcc = pd.DataFrame([d.keys(), d.values()]).T.rename(columns={0: 'battery_GW', 1: 'elcc'}) + df_elcc = pd.merge(pd.DataFrame({'battery_GW': range(101)}), df_elcc, on='battery_GW', how='outer').interpolate( + method='linear' + ) + + df_full['elcc'] = df_full['battery_GW'].apply( + lambda x: df_elcc.loc[df_elcc['battery_GW'] == int(x), 'elcc'].iloc[0] + ) + + bat_filepath = os.path.join(dst_dir, 'battery_costs.csv') + df_bat = pd.read_csv(bat_filepath) + df_bat.drop(columns='elcc', inplace=True) + df_bat = pd.merge(df_full, df_bat, on='Year', how='left') + df_bat.drop(columns=[i for i in df_bat.columns if 'battery_GW' in i], inplace=True) + + if save: + df_bat.to_csv(bat_filepath, index=False) + + if plot: + df_bat.plot(x='Year', y='elcc', legend=False) + plt.ylabel('ELCC [fraction]') + + return df_bat + + +def viz_trial_logs(trial_dir, quantities, span=2000, figsize=(10, 10)): + """Visualize logs produced from an Rllib training session.""" + + df = pd.read_json(os.path.join(trial_dir, 'result.json'), lines=True)[:span] + fig, ax = plt.subplots(len(quantities), 1, figsize=figsize, sharex=True) + + for i, (q, label) in enumerate(quantities.items()): + data = df['custom_metrics'].apply(lambda x: dict(x).get(q + '_mean')).values + ax[i].plot(range(len(df)), data, label=label, color=colors[i]) + ax[i].legend([q], loc='upper right') + if q in ['NPV', 'ROI']: + ax[i].set_ylim([0, None]) + plt.show() + + +def compute_f(Rewaterprod, well_diam): + """Compute Darcy friction coefficient. + + Args: + Rewaterprod (Union[ndarray, float]): Reynold's number + well_diam (Union[ndarray, float]): well diameter + + Returns: + Union[ndarray, float]: Darcy friction factor + """ + + relroughness = 1e-4 / well_diam + f = 1.0 / np.power(-2 * np.log10(relroughness / 3.7 + 5.74 / np.power(Rewaterprod, 0.9)), 2.0) + f = 1.0 / np.power((-2 * np.log10(relroughness / 3.7 + 2.51 / Rewaterprod / np.sqrt(f))), 2.0) + f = 1.0 / np.power((-2 * np.log10(relroughness / 3.7 + 2.51 / Rewaterprod / np.sqrt(f))), 2.0) + f = 1.0 / np.power((-2 * np.log10(relroughness / 3.7 + 2.51 / Rewaterprod / np.sqrt(f))), 2.0) + f = 1.0 / np.power((-2 * np.log10(relroughness / 3.7 + 2.51 / Rewaterprod / np.sqrt(f))), 2.0) + f = 1.0 / np.power( + (-2 * np.log10(relroughness / 3.7 + 2.51 / Rewaterprod / np.sqrt(f))), 2.0 + ) # 6 iterations to converge + + # if laminar + f = np.where(Rewaterprod < 2300.0, 64.0 / nonzero(Rewaterprod), f) + + return f + + +def densitywater(Twater): + """Computationally efficient correlation for water density based on GEOPHIRES. + + Args: + Twater (Union[ndarray, float]): water temperature in deg C + + Returns: + Union[ndarray, float]: water density in kg/m3 + """ + + T = Twater + 273.15 + rhowater = ( + 0.7983223 + (1.50896e-3 - 2.9104e-6 * T) * T + ) * 1e3 # water density correlation as used in Geophires v1.2 [kg/m3] + return rhowater + + +def viscositywater(Twater): + """Computationally efficient correlation for water viscosity based on the GEOPHIRES. + + Args: + Twater (Union[ndarray, float]): water temperature in deg C + + Returns: + Union[ndarray, float]: water viscosity in m2/s + """ + + muwater = 2.414e-5 * np.power( + 10, 247.8 / (Twater + 273.15 - 140) + ) # accurate to within 2.5% from 0 to 370 degrees C [Ns/m2] + return muwater + + +def heatcapacitywater(Twater): + """Computationally efficient correlation for water specific heat capacity based on the GEOPHIRES. + + Args: + Twater (Union[ndarray, float]): water temperature in deg C + + Returns: + Union[ndarray, float]: water specific heat capacity in J/kg-K + """ + + # Based on GEOPHIRES correlations (more stable XSTEAM) + Twater = (Twater + 273.15) / 1000 + A = -203.6060 + B = 1523.290 + C = -3196.413 + D = 2474.455 + E = 3.855326 + cpwater = (A + B * Twater + C * Twater**2 + D * Twater**3 + E / (Twater**2)) / 18.02 * 1000 + return cpwater + + +def vaporpressurewater(Twater): + """Computationally efficient correlation for vapor pressure based on the GEOPHIRES. + + Args: + Twater (Union[ndarray, float]): water temperature in deg C + + Returns: + Union[ndarray, float]: vapor pressure in kPa + """ + + return np.where( + Twater < 100, + 133.322 * (10 ** (8.07131 - 1730.63 / (233.426 + Twater))) / 1000, + 133.322 * (10 ** (8.14019 - 1810.94 / (244.485 + Twater))) / 1000, + ) + + +def compute_drilling_cost( + well_tvd, well_diam, lateral_length=0, numberoflaterals=1, total_drilling_length=None, drilling_cost=None +): + """Correlations for computing drilling costs.""" + well_md = well_tvd + lateral_length * numberoflaterals + + if total_drilling_length: + if not drilling_cost: + raise ValueError('You must specify drilling drilling_cost for the provided system design.') + else: + return total_drilling_length * drilling_cost / 1e6 + + if drilling_cost: + return well_md * drilling_cost / 1e6 + # Divide by 1.15 to remove the 15% contiengency already included by SNL (ref: GETEM page 85) + if well_diam > 0.3: # units in meters + if lateral_length > 0: + return (0.2553 * well_md**2 + 1716.7157 * well_md + 500867.0) * 1e-6 / 1.15 + else: + return (0.2818 * well_tvd**2 + 1275.5213 * well_tvd + 632315.0) * 1e-6 / 1.15 + else: + if lateral_length > 0: + return (0.2898 * well_md**2 + 822.1507 * well_md + 680563.0) * 1e-6 / 1.15 + else: + return (0.3021 * well_tvd**2 + 584.9112 * well_tvd + 751368.0) * 1e-6 / 1.15 + + +def compute_latlon_distance(lat1, lon1, lat2, lon2): + """Compute Euclidean distance in kilometers based on latlon cooredinates of two locations.""" + + return ( + np.arccos(np.sin(lat1) * np.sin(lat2) + np.cos(lat1) * np.cos(lat2) * np.cos(lon2 - lon1)) * 6371 + ) # kilometers + + +def egs_mw_from_temp(T): + """Find GETEM default for power plant size in MW given resource temperature.""" + + if T < 140: + return 10 + elif T < 175: + return 15 + elif T < 250: + return 25 + else: + return 30 + + +def gpd_geographic_area(geodf): + """Compute the geographic area of latlon polygons in a geopandas df.""" + + if not geodf.crs and geodf.crs.is_geographic: + raise TypeError('geodataframe should have geographic coordinate system') + + geod = geodf.crs.get_geod() + + def area_calc(geom): + if geom.geom_type not in ['MultiPolygon', 'Polygon']: + return np.nan + + # For MultiPolygon do each separately + if geom.geom_type == 'MultiPolygon': + return np.sum([area_calc(p) for p in geom.geoms]) + + # orient to ensure a counter-clockwise traversal. + # See https://pyproj4.github.io/pyproj/stable/api/geod.html + # geometry_area_perimeter returns (area, perimeter) + return geod.geometry_area_perimeter(orient(geom, 1))[0] # m2 + + return geodf.geometry.apply(area_calc) + + +def latlon_tres_to_depth( + df_maps, query_northing, query_easting, tres, MJth_per_km3=5.1e9, thickness=500, eta=0.1, L=25 +): + """Get depth required to reach a target temperature at a latlon location.""" + + all_depths = np.array(list(df_maps.keys())) + depths = all_depths + temps = [] + for depth in all_depths: + df_temp = df_maps[depth] + row = df_temp.iloc[ + [ + ( + np.sqrt((query_northing - df_temp['Northing']) ** 2 + (query_easting - df_temp['Easting']) ** 2) + ).argmin() + ] + ] + temp = row['T'].values[0] + temps.append(temp) + if temp >= tres: + depths = all_depths[: len(temps)] + break + + temps = np.array(sorted(temps)) + surface_temp = temps[0] + geothermal_gradient = np.mean(np.clip(np.diff(temps[[0, -1]]) / (np.diff(depths[[0, -1]]) / 1000), 0.0, np.inf)) + if geothermal_gradient <= 0.0: + print(temps, depths, tres, temp) + print('Non-Positive Geothermal Gradient ... !!!') + + if tres > max(temps): + z = np.polyfit(temps, depths, 1) + p = np.poly1d(z) + well_depth = p(tres) + else: + well_depth = np.interp(tres, temps, depths) + + A_r = gpd_geographic_area(row).values[0] / 1e6 # km2 + V_r = A_r * thickness / 1000 # km3 + electric_energy_MJe = MJth_per_km3 * V_r * eta + electric_power_MWe = electric_energy_MJe / (L * 365 * 24 * 3600) + + return well_depth, surface_temp, geothermal_gradient, electric_power_MWe, A_r, V_r + + +def latlon_depth_to_tres( + df_maps, query_northing, query_easting, well_depth, MJth_per_km3=5.1e9, thickness=1000, eta=0.1, L=25 +): + """Get depth required to reach a target temperature at a latlon location.""" + + all_depths = np.sort(list(df_maps.keys())) + ref_depth = all_depths[np.argmin(np.abs(np.array(list(df_maps.keys())) - well_depth))] + df_temp = df_maps[ref_depth] + ref_row = df_temp.iloc[ + [(np.sqrt((query_northing - df_temp['Northing']).abs() + (query_easting - df_temp['Easting']).abs())).argmin()] + ] + + depths = all_depths[all_depths <= well_depth] + temps = [] + for depth in depths: + df_temp = df_maps[depth] + row = df_temp.iloc[ + [ + ( + np.sqrt((query_northing - df_temp['Northing']).abs() + (query_easting - df_temp['Easting']).abs()) + ).argmin() + ] + ] + temp = row['T'].values[0] + temps.append(temp) + + temps = np.array(sorted(temps)) + surface_temp = temps[0] + geothermal_gradient = np.mean(np.clip(np.diff(temps[[0, -1]]) / (np.diff(depths[[0, -1]]) / 1000), 0.0, np.inf)) + if geothermal_gradient <= 0.0: + print(temps, depths, tres, temp) + print('Non-Positive Geothermal Gradient ... !!!') + + z = np.polyfit(depths[-2:], temps[-2:], 1) + p = np.poly1d(z) + tres = p(well_depth) + + A_r = gpd_geographic_area(ref_row).values[0] / 1e6 # km2 + V_r = A_r * thickness / 1000 # km3 + electric_energy_MJe = MJth_per_km3 * V_r * eta + electric_power_MWe = electric_energy_MJe / (L * 365 * 24 * 3600) + + return tres, surface_temp, geothermal_gradient, electric_power_MWe, A_r, V_r + + +def augustine_MWeperkm3(tres): + """Based on page22 of Augustine (2011): https://www.nrel.gov/docs/fy12osti/47459.pdf""" + if tres <= 200: + return 0.59 + elif tres <= 250: + return 0.76 + elif tres <= 300: + return 0.86 + elif tres <= 350: + return 0.97 + else: + return 1.19 + + +def augustine_eff(tres): + """Based on page22 of Augustine (2011): https://www.nrel.gov/docs/fy12osti/47459.pdf""" + if tres <= 200: + return 0.11 + elif tres <= 250: + return 0.14 + elif tres <= 300: + return 0.16 + elif tres <= 350: + return 0.18 + else: + return 0.22 + + +def nonzero(x, thresh=1e-6): + """Ensure input is nonzero + + Args: + x (ndarray): input numbers + thresh (float, optional): threshold below which numbers are set to zero. Defaults to 1E-6. + + Returns: + ndarray: nonzero version of input array + """ + return np.maximum(thresh, x) + + +def grid_bounds(geom, delta): + """Create a grid/mesh based on geometry.""" + minx, miny, maxx, maxy = geom.bounds + nx = int((maxx - minx) / delta) + ny = int((maxy - miny) / delta) + gx, gy = np.linspace(minx, maxx, nx), np.linspace(miny, maxy, ny) + grid = [] + for i in range(len(gx) - 1): + for j in range(len(gy) - 1): + poly_ij = Polygon([[gx[i], gy[j]], [gx[i], gy[j + 1]], [gx[i + 1], gy[j + 1]], [gx[i + 1], gy[j]]]) + grid.append(poly_ij) + return grid + + +def partition(geom, delta): + """Partition a geometry based on delta value.""" + prepared_geom = prep(geom) + grid = list(filter(prepared_geom.intersects, grid_bounds(geom, delta))) + return grid + + +def harrison(z): + """Harrison temperature correction. + + Args: + z (ndarray): true vertical depth in meters + + Returns: + ndarray: harrison correction in deg C + """ + return -2.3449e-6 * z**2 + 0.018268 * z - 16.512 + + +def forster(z): + """Forster temperature correction. + + Args: + z (ndarray): true vertical depth in meters + + Returns: + ndarray: forster correction in deg C + """ + return 0.017 * z - 6.58 + + +def retrieve_weather(lat, lon, year, station_idx_limit=10): + """Retrieve weather.""" + + # Set time period + findtz = TimezoneFinder() + findstations = Stations() + start = datetime(year, 1, 1) + end = datetime(year, 12, 31, 23, 59) + + data = [] + station_idx = 0 + stations = findstations.nearby(lat, lon).fetch(station_idx_limit + 1) + + while len(data) != 8760: + station = stations.iloc[[station_idx]] + data = Hourly( + station, start - timedelta(days=2), end + timedelta(days=2) + ) # get more days to account for timezone shifting + data = data.fetch() + data = data.bfill().ffill() + data.reset_index(inplace=True) + tz = findtz.timezone_at(lng=lon, lat=lat) + if len(data): + data['time'] = data.time.dt.tz_localize('utc').dt.tz_convert(tz).dt.tz_localize(None) + year_col = data.time.dt.year + data = data[year_col == year].reset_index(drop=True) + + if station_idx >= station_idx_limit: + print(state, f'##### {len(data)} #####') + break + station_idx += 1 + + return data + + +class FaissKNeighbors: + def __init__(self, k=5): + self.index = None + self.y = None + self.k = k + + def fit(self, X, y): + self.index = faiss.IndexFlatL2(X.shape[1]) + self.index.add(X.astype(np.float32)) + self.y = y + + def predict(self, X): + + distances, indices = self.index.search(X.astype(np.float32), k=self.k) + inv_distances = 1 / (distances + 1e-3) ** 4 + inv_distances = inv_distances / inv_distances.sum(axis=1, keepdims=True) + predictions = (inv_distances[:, :, None] * self.y[indices]).sum(axis=1) + + return predictions + + +if __name__ == '__main__': + pass diff --git a/src/geophires_x/Model.py b/src/geophires_x/Model.py index 211c26a78..72d492b95 100644 --- a/src/geophires_x/Model.py +++ b/src/geophires_x/Model.py @@ -8,6 +8,7 @@ from geophires_x.GeoPHIRESUtils import read_input_file from geophires_x.OutputsAddOns import OutputsAddOns from geophires_x.OutputsS_DAC_GT import OutputsS_DAC_GT +from geophires_x.SurfacePlantFgemORC import SurfacePlantFgemOrc from geophires_x.TDPReservoir import TDPReservoir from geophires_x.WellBores import WellBores from geophires_x.SurfacePlant import SurfacePlant @@ -176,6 +177,8 @@ def __init__(self, enable_geophires_logging_config=True, input_file=None): self.surfaceplant = SurfacePlantSUTRA(self) elif self.InputParameters['Power Plant Type'].sValue in ['9', 'Industrial']: self.surfaceplant = SurfacePlantIndustrialHeat(self) + elif self.InputParameters['Power Plant Type'].sValue in ['10', 'FGEM ORC']: + self.surfaceplant = SurfacePlantFgemOrc(self) # if we find out we have an add-ons, we need to instantiate it, then read for the parameters if 'AddOn Nickname 1' in self.InputParameters: @@ -237,6 +240,8 @@ def read_parameters(self, default_output_path: Path = None) -> None: self.surfaceplant = SurfacePlantSupercriticalOrc(self) elif self.surfaceplant.plant_type.value == PlantType.SINGLE_FLASH: self.surfaceplant = SurfacePlantSingleFlash(self) + elif self.surfaceplant.plant_type.value == PlantType.FGEM_ORC: + self.surfaceplant = SurfacePlantFgemOrc(self) else: # default is double flash self.surfaceplant = SurfacePlantDoubleFlash(self) diff --git a/src/geophires_x/OptionList.py b/src/geophires_x/OptionList.py index 6d0dfc6ce..e1e845f71 100644 --- a/src/geophires_x/OptionList.py +++ b/src/geophires_x/OptionList.py @@ -86,6 +86,7 @@ class PlantType(GeophiresInputEnum): DISTRICT_HEATING = 7, "District Heating" RTES = 8, "Reservoir Thermal Energy Storage" INDUSTRIAL = 9, "Industrial" + FGEM_ORC = 10, "FGEM ORC" @staticmethod def from_input_string(input_string: str): diff --git a/src/geophires_x/SurfacePlant.py b/src/geophires_x/SurfacePlant.py index ccf2ac9a6..add973a44 100644 --- a/src/geophires_x/SurfacePlant.py +++ b/src/geophires_x/SurfacePlant.py @@ -279,7 +279,7 @@ def __init__(self, model: Model): self.plant_type = self.ParameterDict[self.plant_type.Name] = intParameter( "Power Plant Type", DefaultValue=PlantType.SUB_CRITICAL_ORC.int_value, - AllowableRange=[1, 2, 3, 4, 5, 6, 7, 8, 9], + AllowableRange=[1, 2, 3, 4, 5, 6, 7, 8, 9, 10], ValuesEnum=PlantType, UnitType=Units.NONE, ErrMessage="assume default power plant type (1: subcritical ORC)", diff --git a/src/geophires_x/SurfacePlantFgemORC.py b/src/geophires_x/SurfacePlantFgemORC.py new file mode 100644 index 000000000..e40fd5964 --- /dev/null +++ b/src/geophires_x/SurfacePlantFgemORC.py @@ -0,0 +1,187 @@ +from datetime import timedelta +from pathlib import Path + +import numpy as np + +from fgem.powerplant import GEOPHIRESORCPowerPlant, ORCPowerPlant +from geophires_x.Parameter import HasQuantity +from geophires_x.SurfacePlant import SurfacePlant +import geophires_x.Model as Model + + +class SurfacePlantFgemOrc(SurfacePlant): + def __init__(self, model: Model): + """ + The __init__ function is called automatically when a class is instantiated. + It initializes the attributes of an object, and sets default values for certain arguments that can be overridden + by user input. + The __init__ function is used to set up all the parameters in the Surfaceplant. + :param model: The container class of the application, giving access to everything else, including the logger + :type model: :class:`~geophires_x.Model.Model` + :return: None + """ + + model.logger.info("Init " + self.__class__.__name__ + ": " + __name__) + super().__init__(model) # Initialize all the parameters in the superclass + + # Set up all the Parameters that will be predefined by this class using the different types of parameter classes. + # Setting up includes giving it a name, a default value, The Unit Type (length, volume, temperature, etc.) and + # Unit Name of that value, sets it as required (or not), sets allowable range, the error message if that range + # is exceeded, the ToolTip Text, and the name of teh class that created it. + # This includes setting up temporary variables that will be available to all the class but noy read in by user, + # or used for Output + # This also includes all Parameters that are calculated and then published using the Printouts function. + + # These dictionaries contain a list of all the parameters set in this object, stored as "Parameter" and + # "OutputParameter" Objects. This will allow us later to access them in a user interface and get that list, + # along with unit type, preferred units, etc. + + # local variable initialization + sclass = self.__class__.__name__ + self.MyClass = sclass + self.MyPath = Path(__file__).resolve() + + model.logger.info("Complete " + self.__class__.__name__ + ": " + __name__) + + def __str__(self): + return "SurfacePlantFgemORC" + + def read_parameters(self, model: Model) -> None: + """ + The read_parameters function reads in the parameters from a dictionary and stores them in the parameters. + It also handles special cases that need to be handled after a value has been read in and checked. + If you choose to subclass this master class, you can also choose to override this method (or not), and if you do + :param model: The container class of the application, giving access to everything else, including the logger + :return: None + """ + model.logger.info("Init " + self.__class__.__name__ + ": " + __name__) + super().read_parameters(model) # Initialize all the parameters in the superclass + model.logger.info("complete " + self.__class__.__name__ + ": " + __name__) + + def Calculate(self, model: Model) -> None: + """ + :param model: The container class of the application + :type model: :class:`~geophires_x.Model.Model` + """ + + model.logger.info("Init " + self.__class__.__name__ + ": " + __name__) + + # calculate power plant entering temperature + self.TenteringPP.value = SurfacePlant.power_plant_entering_temperature(self, self.enduse_option.value, + model.reserv.timevector.value, + self.T_chp_bottom.value, + model.wellbores.ProducedTemperature.value) + + + + # Availability water + self.Availability.value = SurfacePlant.availability_water(self, self.ambient_temperature.value, + self.TenteringPP.value, + self.ambient_temperature.value) + + # Subcritical ORC-specific values. + if self.ambient_temperature.value < 15.: + C21 = 0.0 + C11 = 2.746E-3 + C01 = -8.3806E-2 + D21 = 0.0 + D11 = 2.713E-3 + D01 = -9.1841E-2 + C22 = 0.0 + C12 = 0.0894 + C02 = 55.6 + D22 = 0.0 + D12 = 0.0894 + D02 = 62.6 + else: + C21 = 0.0 + C11 = 2.713E-3 + C01 = -9.1841E-2 + D21 = 0.0 + D11 = 2.676E-3 + D01 = -1.012E-1 + C22 = 0.0 + C12 = 0.0894 + C02 = 62.6 + D22 = 0.0 + D12 = 0.0894 + D02 = 69.6 + + model.wellbores.Tinj.value, ReinjTemp, etau = SurfacePlant.reinjection_temperature(self, model, + self.ambient_temperature.value, self.TenteringPP.value, model.wellbores.Tinj.value, + C01, C11, C21, D01, D11, D21, C02, C12, C22, D02, D12, D22) + + # calculate electricity & heat production + self.ElectricityProduced.value, self.HeatExtracted.value, self.HeatProduced.value, HeatExtractedTowardsElectricity = \ + SurfacePlant.electricity_heat_production(self, self.enduse_option.value, self.Availability.value, etau, + model.wellbores.nprod.value, + model.wellbores.prodwellflowrate.value, + model.reserv.cpwater.value, + model.wellbores.ProducedTemperature.value, + model.wellbores.Tinj.value, ReinjTemp, self.T_chp_bottom.value, + self.enduse_efficiency_factor.value, self.chp_fraction.value) + + # subtract pumping power for net electricity and calculate first law efficiency + self.NetElectricityProduced.value = self.ElectricityProduced.value - model.wellbores.PumpingPower.value + self.FirstLawEfficiency.value = self.NetElectricityProduced.value / HeatExtractedTowardsElectricity + + # Calculate annual electricity, pum;ping, and heat production + self.HeatkWhExtracted.value, self.PumpingkWh.value, self.TotalkWhProduced.value, self.NetkWhProduced.value, self.HeatkWhProduced.value = \ + SurfacePlant.annual_electricity_pumping_power(self, self.plant_lifetime.value, self.enduse_option.value, + self.HeatExtracted.value, + model.economics.timestepsperyear.value, + self.utilization_factor.value, + model.wellbores.PumpingPower.value, + self.ElectricityProduced.value, + self.NetElectricityProduced.value, self.HeatProduced.value) + + # calculate reservoir heat content + self.RemainingReservoirHeatContent.value = SurfacePlant.remaining_reservoir_heat_content( + self, model.reserv.InitialReservoirHeatContent.value, self.HeatkWhExtracted.value) + + self._calculate_fgem(model) + self._calculate_derived_outputs(model) + model.logger.info(f"complete {self.__class__.__name__}: {__name__}") + + def _calculate_fgem(self, model: Model): + + def degC(p: HasQuantity) -> float: + return p.quantity().to('degC').magnitude + + # fgem_TenteringPP = np.zeros(len(model.reserv.timevector.value)) + + # nameplate_capacity_MW = 10 # FIXME WIP + # self._fgem_geophires_orc_power_plant: GEOPHIRESORCPowerPlant = GEOPHIRESORCPowerPlant( + # nameplate_capacity_MW, + # model.reserv.Tresoutput.value + # ) # FIXME TEMP WIP + + fgem_timestep = timedelta(days=365 * (1 / model.economics.timestepsperyear.value)) + fgem_orc_power_plant = ORCPowerPlant( + # degC(model.reserv.Tresoutput)[0], + degC(model.wellbores.ProducedTemperature)[0], + degC(self.ambient_temperature), + model.wellbores.prodwellflowrate.value, + num_prd=model.wellbores.nprod.value, + # TODO capacity factor + cf=self.utilization_factor.value, + + timestep=fgem_timestep + ) + + fgem_Tinj = np.zeros(len(model.reserv.timevector.value)) + fgem_ElectricityProduced = np.zeros(len(model.reserv.timevector.value)) + for t in range(len(model.reserv.timevector.value)): + m_turbine = model.wellbores.prodwellflowrate.value * model.wellbores.nprod.value + fgem_orc_power_plant.step( + m_turbine, # total mass flowing into turbine + model.wellbores.ProducedTemperature.value[t], # T_prd_wh = producer wellhead temperature degC + self.ambient_temperature.value + ) + fgem_Tinj[t] = fgem_orc_power_plant.T_inj + fgem_ElectricityProduced[t] = fgem_orc_power_plant.power_output_MWe + + + fgem_calculated = True + + diff --git a/tests/fgem_tests/__init__.py b/tests/fgem_tests/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/tests/fgem_tests/test_fgem.py b/tests/fgem_tests/test_fgem.py new file mode 100644 index 000000000..527c9b35c --- /dev/null +++ b/tests/fgem_tests/test_fgem.py @@ -0,0 +1,35 @@ +from __future__ import annotations + + +from base_test_case import BaseTestCase +from geophires_x_client import GeophiresXResult, GeophiresXClient, GeophiresInputParameters + + +# ruff: noqa: I001 # Successful module initialization is dependent on this specific import order. + +# noinspection PyProtectedMember + +# noinspection PyProtectedMember + + +# noinspection SpellCheckingInspection +class FgemTestCase(BaseTestCase): + + def test_fgem(self): + r = self._get_result({}) + self.assertIsNotNone(r) + + def _egs_test_file_path(self) -> str: + # return self._get_test_file_path('generic-egs-case-2_sam-single-owner-ppa.txt') + return self._get_test_file_path('../examples/example4.txt') + + def _get_result(self, _params, file_path=None) -> GeophiresXResult: + if file_path is None: + file_path = self._egs_test_file_path() + + return GeophiresXClient().get_geophires_result( + GeophiresInputParameters( + from_file_path=file_path, + params={'Power Plant Type': 10, 'Time steps per year': 12, **_params}, + ) + )