diff --git a/src/geophires_x/Economics.py b/src/geophires_x/Economics.py index 562c7ac6..75b80963 100644 --- a/src/geophires_x/Economics.py +++ b/src/geophires_x/Economics.py @@ -1,6 +1,5 @@ import math import sys -import os import numpy as np import numpy_financial as npf import geophires_x.Model as Model @@ -8,49 +7,7 @@ from geophires_x.Parameter import intParameter, floatParameter, OutputParameter, ReadParameter, boolParameter, \ coerce_int_params_to_enum_values from geophires_x.Units import * - - -def calculate_total_drilling_lengths_m(Configuration, numnonverticalsections: int, nonvertical_length_km: float, - InputDepth_km: float, OutputDepth_km: float, nprod:int, ninj:int) -> tuple: - """ - returns the total length, vertical length, and non-vertical lengths, depending on the configuration - :param Configuration: Configuration of the well - :type Configuration: :class:`~geophires - :param numnonverticalsections: number of non-vertical sections - :type numnonverticalsections: int - :param nonvertical_length_km: length of non-vertical sections in km - :type nonvertical_length_km: float - :param InputDepth_km: depth of the well in km - :type InputDepth_km: float - :param OutputDepth_km: depth of the output end of the well in km, if U shaped, and not horizontal - :type OutputDepth_km: float - :param nprod: number of production wells - :type nprod: int - :param ninj: number of injection wells - :return: total length, vertical length, and horizontal lengths in meters - :rtype: tuple - """ - if Configuration == Configuration.ULOOP: - # Total drilling depth of both wells and laterals in U-loop [m] - vertical_pipe_length_m = (nprod * InputDepth_km * 1000.0) + (ninj * OutputDepth_km * 1000.0) - nonvertical_pipe_length_m = numnonverticalsections * nonvertical_length_km * 1000.0 - elif Configuration == Configuration.COAXIAL: - # Total drilling depth of well and lateral in co-axial case [m] - vertical_pipe_length_m = (nprod + ninj) * InputDepth_km * 1000.0 - nonvertical_pipe_length_m = numnonverticalsections * nonvertical_length_km * 1000.0 - elif Configuration == Configuration.VERTICAL: - # Total drilling depth of well in vertical case [m] - vertical_pipe_length_m = (nprod + ninj) * InputDepth_km * 1000.0 - nonvertical_pipe_length_m = 0.0 - elif Configuration == Configuration.L: - # Total drilling depth of well in L case [m] - vertical_pipe_length_m = (nprod + ninj) * InputDepth_km * 1000.0 - nonvertical_pipe_length_m = numnonverticalsections * nonvertical_length_km * 1000.0 - else: - raise ValueError(f'Invalid Configuration: {Configuration}') - - tot_pipe_length_m = vertical_pipe_length_m + nonvertical_pipe_length_m - return tot_pipe_length_m, vertical_pipe_length_m, nonvertical_pipe_length_m +from geophires_x.WellBores import calculate_total_drilling_lengths_m def calculate_cost_of_one_vertical_well(model: Model, depth_m: float, well_correlation: int, @@ -1809,8 +1766,14 @@ def __init__(self, model: Model): PreferredUnits=CurrencyUnit.MDOLLARS, CurrentUnits=CurrencyUnit.MDOLLARS ) - self.cost_nonvertical_section = self.OutputParameterDict[self.cost_nonvertical_section.Name] = OutputParameter( - Name="Cost of the non-vertical section of a well", + self.cost_lateral_section = self.OutputParameterDict[self.cost_lateral_section.Name] = OutputParameter( + Name="Cost of the entire (multi-) lateral section of a well", + UnitType=Units.CURRENCY, + PreferredUnits=CurrencyUnit.MDOLLARS, + CurrentUnits=CurrencyUnit.MDOLLARS + ) + self.cost_to_junction_section = self.OutputParameterDict[self.cost_to_junction_section.Name] = OutputParameter( + Name="Cost of the entire section of a well from bottom of vertical to junction with laterals", UnitType=Units.CURRENCY, PreferredUnits=CurrencyUnit.MDOLLARS, CurrentUnits=CurrencyUnit.MDOLLARS @@ -2220,7 +2183,7 @@ def Calculate(self, model: Model) -> None: (self.cost_one_injection_well.value * model.wellbores.ninj.value)) else: if hasattr(model.wellbores, 'numnonverticalsections') and model.wellbores.numnonverticalsections.Provided: - self.cost_nonvertical_section.value = 0.0 + self.cost_lateral_section.value = 0.0 if not model.wellbores.IsAGS.value: input_vert_depth_km = model.reserv.depth.quantity().to('km').magnitude output_vert_depth_km = 0.0 @@ -2229,13 +2192,13 @@ def Calculate(self, model: Model) -> None: output_vert_depth_km = model.reserv.OutputDepth.quantity().to('km').magnitude model.wellbores.injection_reservoir_depth.value = input_vert_depth_km - tot_m, tot_vert_m, tot_horiz_m = calculate_total_drilling_lengths_m(model.wellbores.Configuration.value, - model.wellbores.numnonverticalsections.value, - model.wellbores.Nonvertical_length.value / 1000.0, - input_vert_depth_km, - output_vert_depth_km, - model.wellbores.nprod.value, - model.wellbores.ninj.value) + tot_m, tot_vert_m, tot_horiz_m, _ = calculate_total_drilling_lengths_m(model.wellbores.Configuration.value, + model.wellbores.numnonverticalsections.value, + model.wellbores.Nonvertical_length.value / 1000.0, + input_vert_depth_km, + output_vert_depth_km, + model.wellbores.nprod.value, + model.wellbores.ninj.value) else: tot_m = tot_vert_m = model.reserv.depth.quantity().to('km').magnitude @@ -2261,7 +2224,7 @@ def Calculate(self, model: Model) -> None: self.injection_well_cost_adjustment_factor.value) if hasattr(model.wellbores, 'numnonverticalsections') and model.wellbores.numnonverticalsections.Provided: - self.cost_nonvertical_section.value = calculate_cost_of_non_vertical_section(model, tot_horiz_m, + self.cost_lateral_section.value = calculate_cost_of_non_vertical_section(model, tot_horiz_m, self.wellcorrelation.value, self.Nonvertical_drilling_cost_per_m.value, model.wellbores.numnonverticalsections.value, @@ -2269,12 +2232,12 @@ def Calculate(self, model: Model) -> None: model.wellbores.NonverticalsCased.value, self.production_well_cost_adjustment_factor.value) else: - self.cost_nonvertical_section.value = 0.0 + self.cost_lateral_section.value = 0.0 # cost of the well field # 1.05 for 5% indirect costs self.Cwell.value = 1.05 * ((self.cost_one_production_well.value * model.wellbores.nprod.value) + (self.cost_one_injection_well.value * model.wellbores.ninj.value) + - self.cost_nonvertical_section.value) + self.cost_lateral_section.value) # reservoir stimulation costs (M$/injection well). These are calculated whether totalcapcost.Valid = 1 if self.ccstimfixed.Valid: diff --git a/src/geophires_x/Model.py b/src/geophires_x/Model.py index 4852e1d0..96c463e3 100644 --- a/src/geophires_x/Model.py +++ b/src/geophires_x/Model.py @@ -11,6 +11,9 @@ from geophires_x.TDPReservoir import TDPReservoir from geophires_x.WellBores import WellBores from geophires_x.SurfacePlant import SurfacePlant +from geophires_x.SBTEconomics import SBTEconomics +from geophires_x.SBTWellbores import SBTWellbores +from geophires_x.SBTReservoir import SBTReservoir from geophires_x.SurfacePlantIndustrialHeat import SurfacePlantIndustrialHeat from geophires_x.SurfacePlantSubcriticalORC import SurfacePlantSubcriticalOrc from geophires_x.SurfacePlantSupercriticalORC import SurfacePlantSupercriticalOrc @@ -72,68 +75,107 @@ def __init__(self, enable_geophires_logging_config=True, input_file=None): if input_file is None and len(sys.argv) > 1: input_file = sys.argv[1] + # Key step - read the entire provided input file read_input_file(self.InputParameters, logger=self.logger, input_file_name=input_file) + # initiate the outputs object + output_file = 'HDR.out' + if len(sys.argv) > 2: + output_file = sys.argv[2] + self.outputs = Outputs(self, output_file=output_file) + + # Initiate the elements of the Model object + # this is where you can change what class get initiated - the superclass, or one of the subclasses + self.logger.info("Initiate the elements of the Model") + + # Assume that SDAC and add-ons are not used self.sdacgtoutputs = None self.sdacgteconomics = None self.addoutputs = None self.addeconomics = None - # Initiate the elements of the Model - # this is where you can change what class get initiated - the superclass, or one of the subclasses - self.logger.info("Initiate the elements of the Model") - # we need to decide which reservoir to instantiate based on the user input (InputParameters), - # which we just read above for the first time - # Default is Thermal drawdown percentage model (GETEM) - self.reserv: TDPReservoir = TDPReservoir(self) - if 'Reservoir Model' in self.InputParameters: - if self.InputParameters['Reservoir Model'].sValue == '0': - self.reserv: CylindricalReservoir = CylindricalReservoir(self) # Simple Cylindrical Reservoir - elif self.InputParameters['Reservoir Model'].sValue == '1': - self.reserv: MPFReservoir = MPFReservoir(self) # Multiple parallel fractures model (LANL) - elif self.InputParameters['Reservoir Model'].sValue == '2': - self.reserv: LHSReservoir = LHSReservoir(self) # Multiple parallel fractures model (LANL) - elif self.InputParameters['Reservoir Model'].sValue == '3': - self.reserv: SFReservoir = SFReservoir(self) # Drawdown parameter model (Tester) - elif self.InputParameters['Reservoir Model'].sValue == '5': - self.reserv: UPPReservoir = UPPReservoir(self) # Generic user-provided temperature profile - elif self.InputParameters['Reservoir Model'].sValue == '6': - self.reserv: TOUGH2Reservoir = TOUGH2Reservoir(self) # Tough2 is called - elif self.InputParameters['Reservoir Model'].sValue == '7': - self.reserv: SUTRAReservoir = SUTRAReservoir(self) # SUTRA output is created - # initialize the default objects + self.reserv: TDPReservoir = TDPReservoir(self) self.wellbores: WellBores = WellBores(self) - self.surfaceplant: SurfacePlant = SurfacePlant(self) + self.surfaceplant = SurfacePlantIndustrialHeat(self) # default is Industrial Heat self.economics: Economics = Economics(self) - output_file = 'HDR.out' - if len(sys.argv) > 2: - output_file = sys.argv[2] - - self.outputs = Outputs(self, output_file=output_file) + # Now we need to handle the creation of all the special case objects based on the user settings + # We have to access the user setting from the overall master table because the read_parameters functions + # have not been called on the specific objects yet, so the instantiated objects only contain default values + # not user set values. The user set value can only come from the master dictionary "InputParameters" + # based on the user input, which we just read above for the first time + # First, we need to decide which reservoir to instantiate + # For reservoirs, the default is Thermal drawdown percentage model (GETEM); see above where it is initialized. + # The user can change this in the input file, so check the values and initiate the appropriate reservoir object if 'Reservoir Model' in self.InputParameters: - if self.InputParameters['Reservoir Model'].sValue == '7': - # if we use SUTRA output for simulating reservoir thermal energy storage, we use a special wellbore object that can handle SUTRA data + if self.InputParameters['Reservoir Model'].sValue in ['0', 'Simple cylindrical']: + self.reserv: CylindricalReservoir = CylindricalReservoir(self) + elif self.InputParameters['Reservoir Model'].sValue in ['1', 'Multiple Parallel Fractures']: + self.reserv: MPFReservoir = MPFReservoir(self) + elif self.InputParameters['Reservoir Model'].sValue in ['2', '1-D Linear Heat Sweep']: + self.reserv: LHSReservoir = LHSReservoir(self) + elif self.InputParameters['Reservoir Model'].sValue in ['3', 'Single Fracture m/A Thermal Drawdown']: + self.reserv: SFReservoir = SFReservoir(self) + elif self.InputParameters['Reservoir Model'].sValue in ['5', 'User-Provided Temperature Profile']: + self.reserv: UPPReservoir = UPPReservoir(self) + elif self.InputParameters['Reservoir Model'].sValue in ['6', 'TOUGH2 Simulator']: + self.reserv: TOUGH2Reservoir = TOUGH2Reservoir(self) + elif self.InputParameters['Reservoir Model'].sValue in ['7', 'SUTRA']: + # if we use SUTRA output for simulating reservoir thermal energy storage, + # we use a special wellbore object that handles SUTRA data, and special Economics and Outputs objects + self.logger.info('Setup the SUTRA elements of the Model and instantiate new attributes as needed') + self.reserv: SUTRAReservoir = SUTRAReservoir(self) self.wellbores: WellBores = SUTRAWellBores(self) self.surfaceplant: SurfacePlantSUTRA = SurfacePlantSUTRA(self) self.economics: SUTRAEconomics = SUTRAEconomics(self) self.outputs: SUTRAOutputs = SUTRAOutputs(self, output_file=output_file) + elif self.InputParameters['Reservoir Model'].sValue in ['8', 'SBT']: + self.logger.info('Setup the SBT elements of the Model and instantiate new attributes as needed') + self.reserv: SBTReservoir = SBTReservoir(self) + self.wellbores: SBTWellBores = SBTWellbores(self) + self.economics: SBTEconomics = SBTEconomics(self) + # Now handle the special cases for all AGS cases (CLGS, SBT, or CLGS) if 'Is AGS' in self.InputParameters: if self.InputParameters['Is AGS'].sValue in ['True', 'true', 'TRUE', 'T', '1']: - self.logger.info("Setup the AGS elements of the Model and instantiate new attributes as needed") - # If we are doing AGS, we need to replace the various objects we with versions of the objects - # that have AGS functionality. - # that means importing them, initializing them, then reading their parameters - # use the simple cylindrical reservoir for all AGS systems. - self.reserv: CylindricalReservoir = CylindricalReservoir(self) - self.wellbores: WellBores = AGSWellBores(self) - self.surfaceplant: SurfacePlantAGS = SurfacePlantAGS(self) - self.economics: AGSEconomics = AGSEconomics(self) - self.outputs: AGSOutputs = AGSOutputs(self, output_file=output_file) + self.logger.info('Setup the AGS elements of the Model and instantiate new attributes as needed') self.wellbores.IsAGS.value = True + if not isinstance(self.reserv, SBTReservoir): + if self.InputParameters['Economic Model'].sValue not in ['4', 'Simple (CLGS)']: + # must be doing wangju approach, # so go back to using default objects + self.surfaceplant = SurfacePlant(self) + self.economics = Economics(self) + # Must be doing CLGS, so we need to instantiate the right objects + self.reserv: CylindricalReservoir = CylindricalReservoir(self) + self.wellbores: WellBores = AGSWellBores(self) + self.surfaceplant: SurfacePlantAGS = SurfacePlantAGS(self) + self.economics: AGSEconomics = AGSEconomics(self) + self.outputs: AGSOutputs = AGSOutputs(self, output_file=output_file) + + # initialize the right Power Plant Type + if 'Power Plant Type' in self.InputParameters: + # electricity + if self.InputParameters['Power Plant Type'].sValue in ['1', 'Subcritical ORC']: + self.surfaceplant = SurfacePlantSubcriticalOrc(self) + elif self.InputParameters['Power Plant Type'].sValue in ['2', 'Supercritical ORC']: + self.surfaceplant = SurfacePlantSupercriticalOrc(self) + elif self.InputParameters['Power Plant Type'].sValue in ['3', 'Single-Flash']: + self.surfaceplant = SurfacePlantSingleFlash(self) + elif self.InputParameters['Power Plant Type'].sValue in ['4', 'Double-Flash']: + self.surfaceplant = SurfacePlantDoubleFlash(self) + # Heat applications + elif self.InputParameters['Power Plant Type'].sValue in ['5', 'Absorption Chiller']: + self.surfaceplant = SurfacePlantAbsorptionChiller(self) + elif self.InputParameters['Power Plant Type'].sValue in ['6', 'Heat Pump']: + self.surfaceplant = SurfacePlantHeatPump(self) + elif self.InputParameters['Power Plant Type'].sValue in ['7', 'District Heating']: + self.surfaceplant = SurfacePlantDistrictHeating(self) + elif self.InputParameters['Power Plant Type'].sValue in ['8', 'Reservoir Thermal Energy Storage']: + self.surfaceplant = SurfacePlantSUTRA(self) + elif self.InputParameters['Power Plant Type'].sValue in ['9', 'Industrial']: + self.surfaceplant = SurfacePlantIndustrialHeat(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: @@ -150,6 +192,7 @@ def __init__(self, enable_geophires_logging_config=True, input_file=None): self.logger.info(f'Complete {__class__}: {__name__}') + def __str__(self): return "Model" diff --git a/src/geophires_x/OptionList.py b/src/geophires_x/OptionList.py index 8abf17a5..7b1e3cb9 100644 --- a/src/geophires_x/OptionList.py +++ b/src/geophires_x/OptionList.py @@ -114,6 +114,7 @@ class ReservoirModel(GeophiresInputEnum): USER_PROVIDED_PROFILE = 5, "User-Provided Temperature Profile" TOUGH2_SIMULATOR = 6, "TOUGH2 Simulator" SUTRA = 7, "SUTRA" + SBT = 8, "SBT" @staticmethod def get_reservoir_model_from_input_string(input_string:str): @@ -252,6 +253,7 @@ class Configuration(GeophiresInputEnum): COAXIAL = 2, "coaxial" VERTICAL = 3, "vertical" L = 4, "L" + EAVORLOOP = 5, "EavorLoop" @staticmethod def from_int(int_val): @@ -266,3 +268,40 @@ def from_input_string(input_string: str): return member raise ValueError(f'Unknown Configuration input value: {input_string}') + + +class FlowrateModel(GeophiresInputEnum): + USER_SUPPLIED = 1, "user supplied" + FILE_SUPPLIED = 2, "file supplied" + + @staticmethod + def from_int(int_val): + for member in __class__: + if member.int_value == int_val: + return member + + @staticmethod + def from_input_string(input_string: str): + for member in __class__: + if input_string == str(member.int_value): + return member + + raise ValueError(f'Unknown Flow Rate Model input value: {input_string}') + +class InjectionTemperatureModel(GeophiresInputEnum): + USER_SUPPLIED = 1, "user supplied" + FILE_SUPPLIED = 2, "file supplied" + + @staticmethod + def from_int(int_val): + for member in __class__: + if member.int_value == int_val: + return member + + @staticmethod + def from_input_string(input_string: str): + for member in __class__: + if input_string == str(member.int_value): + return member + + raise ValueError(f'Unknown Injection Temperature Model input value: {input_string}') diff --git a/src/geophires_x/Outputs.py b/src/geophires_x/Outputs.py index 80a1f547..006e0125 100644 --- a/src/geophires_x/Outputs.py +++ b/src/geophires_x/Outputs.py @@ -1788,10 +1788,10 @@ def PrintOutputs(self, model: Model): model.economics.cost_one_injection_well.value != -1: f.write(f' Drilling and completion costs per production well: {econ.cost_one_production_well.value:10.2f} ' + econ.cost_one_production_well.CurrentUnits.value + NL) f.write(f' Drilling and completion costs per injection well: {econ.cost_one_injection_well.value:10.2f} ' + econ.cost_one_injection_well.CurrentUnits.value + NL) - elif econ.cost_nonvertical_section.value > 0.0: + elif econ.cost_lateral_section.value > 0.0: f.write(f' Drilling and completion costs per vertical production well: {econ.cost_one_production_well.value:10.2f} ' + econ.cost_one_production_well.CurrentUnits.value + NL) f.write(f' Drilling and completion costs per vertical injection well: {econ.cost_one_injection_well.value:10.2f} ' + econ.cost_one_injection_well.CurrentUnits.value + NL) - f.write(f' Drilling and completion costs per non-vertical sections: {econ.cost_nonvertical_section.value:10.2f} ' + econ.cost_nonvertical_section.CurrentUnits.value + NL) + f.write(f' Drilling and completion costs per non-vertical sections: {econ.cost_lateral_section.value:10.2f} ' + econ.cost_lateral_section.CurrentUnits.value + NL) else: f.write(f' Drilling and completion costs per well: {model.economics.Cwell.value/(model.wellbores.nprod.value+model.wellbores.ninj.value):10.2f} ' + model.economics.Cwell.CurrentUnits.value + NL) f.write(f' Stimulation costs: {model.economics.Cstim.value:10.2f} ' + model.economics.Cstim.CurrentUnits.value + NL) diff --git a/src/geophires_x/Reservoir.py b/src/geophires_x/Reservoir.py index 30f709e2..bb1cca8b 100644 --- a/src/geophires_x/Reservoir.py +++ b/src/geophires_x/Reservoir.py @@ -50,7 +50,7 @@ def __init__(self, model: Model): self.resoption = self.ParameterDict[self.resoption.Name] = intParameter( "Reservoir Model", DefaultValue=ReservoirModel.ANNUAL_PERCENTAGE.int_value, - AllowableRange=[0, 1, 2, 3, 4, 5, 6, 7], + AllowableRange=[0, 1, 2, 3, 4, 5, 6, 7, 8], ValuesEnum=ReservoirModel, Required=True, ErrMessage="run default reservoir model (Thermal Drawdown Percentage Model)", diff --git a/src/geophires_x/SBTEconomics.py b/src/geophires_x/SBTEconomics.py new file mode 100644 index 00000000..60721b99 --- /dev/null +++ b/src/geophires_x/SBTEconomics.py @@ -0,0 +1,849 @@ +import sys, math +import numpy as np +import geophires_x.Model as Model +from .Economics import Economics, calculate_cost_of_one_vertical_well, BuildPTCModel, BuildPricingModel, \ + CalculateRevenue, CalculateFinancialPerformance, CalculateLCOELCOHLCOC +from .OptionList import Configuration, WellDrillingCostCorrelation, PlantType +from geophires_x.Parameter import floatParameter +from geophires_x.Units import * +from geophires_x.OptionList import WorkingFluid, EndUseOptions, EconomicModel + + +def calculate_cost_of_lateral_section(model: Model, length_m: float, well_correlation: int, + lateral_drilling_cost_per_m: float, + num_lateral_sections: int, + fixed_well_cost_name: str, NonverticalsCased: bool, + well_cost_adjustment_factor: float) -> float: + """ + calculate_cost_of_lateral_section calculates the cost of the lateral section of the well. + Assume that the cost per meter for drilling of the lateral section is the same as the vertical section, + except the casing cost is half, if it is uncased. + :param model: The model object + :type model: :class:`~geophires + :param length_m: The depth of the well in meters + :type length_m: float + :param well_correlation: The well cost correlation + :type well_correlation: int + :param lateral_drilling_cost_per_m: The lateral drilling cost per meter in $/m + :type lateral_drilling_cost_per_m: float + :param num_lateral_sections: The number of lateral sections + :type num_lateral_sections: int + :param fixed_well_cost_name: The fixed well cost name + :type fixed_well_cost_name: str + :param NonverticalsCased: Are the laterals cased? + :type NonverticalsCased: bool + :param well_cost_adjustment_factor: The well cost adjustment factor + :type well_cost_adjustment_factor: float + :return: cost_of_one_well: The cost of the lateral section in MUSD + :rtype: float + """ + + # if we are drilling a vertical well, the lateral cost is 0 + if model.wellbores.Configuration.value == Configuration.VERTICAL: + return 0.0 + + # Check if well length is out of standard bounds for cost correlation + length_per_section_m = length_m / num_lateral_sections + correlations_min_valid_length_m = 500. + correlations_max_valid_length_m = 8000. + cost_of_lateral_section = 0.0 + cost_per_section = 0.0 + if length_per_section_m < correlations_min_valid_length_m and not well_correlation is WellDrillingCostCorrelation.SIMPLE: + well_correlation = WellDrillingCostCorrelation.SIMPLE + model.logger.warning( + f'Invalid cost correlation specified ({well_correlation}) for per lateral section drilling length ' + f'<{correlations_min_valid_length_m}m ({length_per_section_m}m). ' + f'Falling back to simple user-specified cost ' + f'({lateral_drilling_cost_per_m} per meter)' + ) + + if length_per_section_m > correlations_max_valid_length_m and not well_correlation is WellDrillingCostCorrelation.SIMPLE: + model.logger.warning( + f'{well_correlation} may be invalid for per lateral section rilling length ' + f'>{correlations_max_valid_length_m}m ({length_per_section_m}m). ' + f'Consider using {WellDrillingCostCorrelation.SIMPLE} (per-meter cost) or ' + f'{fixed_well_cost_name} (fixed cost per well) instead.' + ) + + casing_factor = 1.0 + if not NonverticalsCased: + # assume that casing & cementing costs 50% of drilling costs + casing_factor = 0.5 + if model.economics.Nonvertical_drilling_cost_per_m.Provided or well_correlation is WellDrillingCostCorrelation.SIMPLE: + cost_of_lateral_section = casing_factor * ((num_lateral_sections * lateral_drilling_cost_per_m * length_per_section_m)) * 1E-6 + else: + cost_per_section = well_correlation.calculate_cost_MUSD(length_per_section_m) + cost_of_lateral_section = casing_factor * num_lateral_sections * cost_per_section + + # account for adjustment factor + cost_of_lateral_section = well_cost_adjustment_factor * cost_of_lateral_section + return cost_of_lateral_section + + +class SBTEconomics(Economics): + """ + SBTEconomics Child class of Economics; it is the same, but has advanced SBT closed-loop functionality + """ + + def __init__(self, model: Model): + """ + The __init__ function is the constructor for a class. It is called whenever an instance of the class is created. + The __init__ function can take arguments, but self is always the first one. Self refers to the instance of the + object that has already been created, and it's used to access variables that belong to that object + :param model: The container class of the application, giving access to everything else, including the logger + :type model: Model + :return: Nothing, and is used to initialize the class + """ + model.logger.info(f'Init {__class__!s}: {sys._getframe().f_code.co_name}') + + # Initialize the superclass first to gain access to those variables + super().__init__(model) + sclass = str(__class__).replace("", "") + self.MyPath = os.path.abspath(__file__) + #self.Electricity_rate = None + #self.Discount_rate = None + #self.error = 0 + #self.AverageOPEX_Plant = 0 + #self.OPEX_Plant = 0 + #self.TotalCAPEX = 0 + #self.CAPEX_Surface_Plant = 0 + #self.CAPEX_Drilling = 0 + + # 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 the 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. + # If you choose to subclass this master class, you can do so before or after you create your own parameters. + # If you do, you can also choose to call this method from you class, which will effectively add + # and set all these parameters to your class. + # NB: inputs we already have ("already have it") need to be set at ReadParameter time so values + # are set at the last possible time + + """ + self.O_and_M_cost_plant = self.ParameterDict[self.O_and_M_cost_plant.Name] = floatParameter( + "Operation & Maintenance Cost of Surface Plant", + DefaultValue=0.015, + Min=0.0, + Max=0.2, + UnitType=Units.PERCENT, + PreferredUnits=PercentUnit.TENTH, + CurrentUnits=PercentUnit.TENTH, + Required=True, + ErrMessage="assume default Operation & Maintenance cost of surface plant expressed as fraction of total surface plant capital cost (0.015)" + ) + self.Direct_use_heat_cost_per_kWth = self.ParameterDict[ + self.Direct_use_heat_cost_per_kWth.Name] = floatParameter( + "Capital Cost for Surface Plant for Direct-use System", + DefaultValue=100.0, + Min=0.0, + Max=10000.0, + UnitType=Units.ENERGYCOST, + PreferredUnits=EnergyCostUnit.DOLLARSPERKW, + CurrentUnits=EnergyCostUnit.DOLLARSPERKW, + Required=False, + ErrMessage="assume default Capital cost for surface plant for direct-use system (100 $/kWth)" + ) + self.Power_plant_cost_per_kWe = self.ParameterDict[self.Power_plant_cost_per_kWe.Name] = floatParameter( + "Capital Cost for Power Plant for Electricity Generation", + DefaultValue=3000.0, + Min=0.0, + Max=10000.0, + UnitType=Units.ENERGYCOST, + PreferredUnits=EnergyCostUnit.DOLLARSPERKW, + CurrentUnits=EnergyCostUnit.DOLLARSPERKW, + Required=True, + ErrMessage="assume default Power plant capital cost per kWe (3000 USD/kWe)" + ) + + # results are stored here and in the parent ProducedTemperature array + +""" + model.logger.info(f'complete {__class__!s}: {sys._getframe().f_code.co_name}') + + def __str__(self): + return 'SBTEconomics' + + 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 + :type model: :class:`~geophires_x.Model.Model` + :return: None + """ + model.logger.info(f'Init {__class__!s}: {sys._getframe().f_code.co_name}') + super().read_parameters(model) # read the default parameters + # if we call super, we don't need to deal with setting the parameters here, + # just deal with the special cases for the variables in this class + # because the call to the super.readparameters will set all the variables, + # including the ones that are specific to this class + + # inputs we already have - needs to be set at ReadParameter time so values set at the latest possible time + # self.Discount_rate = model.economics.discountrate.value # same units are GEOPHIRES + # self.Electricity_rate = model.surfaceplant.electricity_cost_to_buy.value # same units are GEOPHIRES + + model.logger.info(f'complete {__class__!s}: {sys._getframe().f_code.co_name}') + + def Calculate(self, model: Model) -> None: + """ + The Calculate function is where all the calculations are done. + This function can be called multiple times, and will only recalculate what has changed each time it is called. + This is where all the calculations are made using all the values that have been set. + If you subclass this class, you can choose to run these calculations before (or after) your calculations, + but that assumes you have set all the values that are required for these calculations + If you choose to subclass this master class, you can also choose to override this method (or not), + and if you do, do it before or after you call you own version of this method. If you do, + you can also choose to call this method from you class, which can effectively run the calculations + of the superclass, making all thr values available to your methods. but you had + better have set all the parameters! + :param model: The container class of the application, giving access to everything else, including the logger + :type model: :class:`~geophires_x.Model.Model` + :return: Nothing, but it does make calculations and set values in the model + """ + model.logger.info(f'Init {__class__!s}: {sys._getframe().f_code.co_name}') + + #if hasattr(model.wellbores, 'numnonverticalsections') and model.wellbores.numnonverticalsections.Provided: + #self.cost_lateral_section.value = 0.0 + + # capital costs + # well costs (using GeoVision drilling correlations). These are calculated whether totalcapcostvalid = 1 + # start with the cost of one well + # C1well is well drilling and completion cost in M$/well + if self.per_production_well_cost.Valid: + self.cost_one_production_well.value = self.per_production_well_cost.value + if not self.per_injection_well_cost.Provided: + self.cost_one_injection_well.value = self.per_production_well_cost.value + else: + self.cost_one_injection_well.value = self.per_injection_well_cost.value + self.Cwell.value = ((self.cost_one_production_well.value * model.wellbores.nprod.value) + + (self.cost_one_injection_well.value * model.wellbores.ninj.value)) + else: + # calculate the cost of one vertical production well + # 1.05 for 5% indirect costs + self.cost_one_production_well.value = 1.05 * calculate_cost_of_one_vertical_well(model, model.wellbores.vertical_section_length.value, + self.wellcorrelation.value, + self.Vertical_drilling_cost_per_m.value, + self.per_production_well_cost.Name, + self.production_well_cost_adjustment_factor.value) + + # If there is no injector well, then we assume we are doing a coaxial closed-loop. + if model.wellbores.ninj.value == 0: + self.cost_one_injection_well.value = 0.0 + else: + # Now calculate the cost of one vertical injection well + # assume the cost of the injector and producer is the same + self.cost_one_injection_well.value = self.cost_one_production_well.value + + if hasattr(model.wellbores, 'numnonverticalsections') and model.wellbores.numnonverticalsections.Provided: + # now calculate the costs if we have a lateral section + # 1.05 for 5% indirect costs + self.cost_lateral_section.value = 1.05 * calculate_cost_of_lateral_section(model, model.wellbores.tot_lateral_m.value, + self.wellcorrelation.value, + self.Nonvertical_drilling_cost_per_m.value, + model.wellbores.numnonverticalsections.value, + self.per_injection_well_cost.Name, + model.wellbores.NonverticalsCased.value, + self.production_well_cost_adjustment_factor.value) + + # If it is an EavorLoop, we need to calculate the cost of the section of the well from + # the bottom of the vertical to the junction with the laterals. + # This section is not vertical, but it is cased, so we will estimate the cost + # of this section as if it were a vertical section. + if model.wellbores.Configuration.value == Configuration.EAVORLOOP: + self.cost_to_junction_section.value = 1.05 * calculate_cost_of_one_vertical_well(model, + model.wellbores.tot_to_junction_m.value, + self.wellcorrelation.value, + self.Vertical_drilling_cost_per_m.value, + self.per_injection_well_cost.Name, + self.injection_well_cost_adjustment_factor.value) + else: + self.cost_lateral_section.value = 0.0 + self.cost_to_junction_section.value = 0.0 + + # cost of the well field + self.Cwell.value = ((self.cost_one_production_well.value * model.wellbores.nprod.value) + + (self.cost_one_injection_well.value * model.wellbores.ninj.value) + + self.cost_lateral_section.value + self.cost_to_junction_section.value) + + # reservoir stimulation costs (M$/injection well). These are calculated whether totalcapcost.Valid = 1 + if self.ccstimfixed.Valid: + self.Cstim.value = self.ccstimfixed.value + else: + self.Cstim.value = 1.05 * 1.15 * self.ccstimadjfactor.value * model.wellbores.ninj.value * 1.25 # 1.15 for 15% contingency and 1.05 for 5% indirect costs + + # field gathering system costs (M$) + if self.ccgathfixed.Valid: + self.Cgath.value = self.ccgathfixed.value + else: + self.Cgath.value = self.ccgathadjfactor.value * 50 - 6 * np.max( + model.surfaceplant.HeatExtracted.value) * 1000. # (GEOPHIRES v1 correlation) + if model.wellbores.impedancemodelused.value: + pumphp = np.max(model.wellbores.PumpingPower.value) * 1341 + numberofpumps = np.ceil(pumphp / 2000) # pump can be maximum 2,000 hp + if numberofpumps == 0: + self.Cpumps = 0.0 + else: + pumphpcorrected = pumphp / numberofpumps + self.Cpumps = numberofpumps * 1.5 * ( + (1750 * pumphpcorrected ** 0.7) * 3 * pumphpcorrected ** (-0.11)) + else: + if model.wellbores.productionwellpumping.value: + prodpumphp = np.max(model.wellbores.PumpingPowerProd.value) / model.wellbores.nprod.value * 1341 + Cpumpsprod = model.wellbores.nprod.value * 1.5 * (1750 * prodpumphp ** 0.7 + 5750 * + prodpumphp ** 0.2 + 10000 + np.max( + model.wellbores.pumpdepth.value) * 50 * 3.281) # see page 46 in user's manual assuming rental of rig for 1 day. + else: + Cpumpsprod = 0 + + injpumphp = np.max(model.wellbores.PumpingPowerInj.value) * 1341 + numberofinjpumps = np.ceil(injpumphp / 2000) # pump can be maximum 2,000 hp + if numberofinjpumps == 0: + Cpumpsinj = 0 + else: + injpumphpcorrected = injpumphp / numberofinjpumps + Cpumpsinj = numberofinjpumps * 1.5 * ( + 1750 * injpumphpcorrected ** 0.7) * 3 * injpumphpcorrected ** (-0.11) + self.Cpumps = Cpumpsinj + Cpumpsprod + + # Based on GETEM 2016 #1.15 for 15% contingency and 1.12 for 12% indirect costs + self.Cgath.value = 1.15 * self.ccgathadjfactor.value * 1.12 * ( + (model.wellbores.nprod.value + model.wellbores.ninj.value) * 750 * 500. + self.Cpumps) / 1E6 + + # plant costs + if (model.surfaceplant.enduse_option.value == EndUseOptions.HEAT + and model.surfaceplant.plant_type.value not in [PlantType.ABSORPTION_CHILLER, PlantType.HEAT_PUMP, PlantType.DISTRICT_HEATING]): # direct-use + if self.ccplantfixed.Valid: + self.Cplant.value = self.ccplantfixed.value + else: + self.Cplant.value = 1.12 * 1.15 * self.ccplantadjfactor.value * 250E-6 * np.max( + model.surfaceplant.HeatExtracted.value) * 1000. # 1.15 for 15% contingency and 1.12 for 12% indirect costs + + # absorption chiller + elif model.surfaceplant.enduse_option.value == EndUseOptions.HEAT and model.surfaceplant.plant_type.value == PlantType.ABSORPTION_CHILLER: # absorption chiller + if self.ccplantfixed.Valid: + self.Cplant.value = self.ccplantfixed.value + else: + # this is for the direct-use part all the way up to the absorption chiller + self.Cplant.value = 1.12 * 1.15 * self.ccplantadjfactor.value * 250E-6 * np.max( + model.surfaceplant.HeatExtracted.value) * 1000. # 1.15 for 15% contingency and 1.12 for 12% indirect costs + if self.chillercapex.value == -1: # no value provided by user, use built-in correlation ($2500/ton) + self.chillercapex.value = 1.12 * 1.15 * np.max( + model.surfaceplant.cooling_produced.value) * 1000 / 3.517 * 2500 / 1e6 # $2,500/ton of cooling. 1.15 for 15% contingency and 1.12 for 12% indirect costs + + # now add chiller cost to surface plant cost + self.Cplant.value += self.chillercapex.value + + # heat pump + elif model.surfaceplant.enduse_option.value == EndUseOptions.HEAT and model.surfaceplant.plant_type.value == PlantType.HEAT_PUMP: + if self.ccplantfixed.Valid: + self.Cplant.value = self.ccplantfixed.value + else: + # this is for the direct-use part all the way up to the heat pump + self.Cplant.value = 1.12 * 1.15 * self.ccplantadjfactor.value * 250E-6 * np.max( + model.surfaceplant.HeatExtracted.value) * 1000. # 1.15 for 15% contingency and 1.12 for 12% indirect costs + if self.heatpumpcapex.value == -1: # no value provided by user, use built-in correlation ($150/kWth) + self.heatpumpcapex.value = 1.12 * 1.15 * np.max( + model.surfaceplant.HeatProduced.value) * 1000 * 150 / 1e6 # $150/kW. 1.15 for 15% contingency and 1.12 for 12% indirect costs + + # now add heat pump cost to surface plant cost + self.Cplant.value += self.heatpumpcapex.value + + # district heating + elif model.surfaceplant.enduse_option.value == EndUseOptions.HEAT and model.surfaceplant.plant_type.value == PlantType.DISTRICT_HEATING: + if self.ccplantfixed.Valid: + self.Cplant.value = self.ccplantfixed.value + else: + self.Cplant.value = 1.12 * 1.15 * self.ccplantadjfactor.value * 250E-6 * np.max( + model.surfaceplant.HeatExtracted.value) * 1000. # 1.15 for 15% contingency and 1.12 for 12% indirect costs + self.peakingboilercost.value = 65 * model.surfaceplant.max_peaking_boiler_demand.value / 1000 # add 65$/KW for peaking boiler + self.Cplant.value += self.peakingboilercost.value # add peaking boiler cost to surface plant cost + + + else: # all other options have power plant + if model.surfaceplant.plant_type.value == PlantType.SUB_CRITICAL_ORC: + MaxProducedTemperature = np.max(model.surfaceplant.TenteringPP.value) + if MaxProducedTemperature < 150.: + 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.) + x = np.max(model.surfaceplant.ElectricityProduced.value) + y = np.max(model.surfaceplant.ElectricityProduced.value) + if y == 0.0: + y = 15.0 + z = math.pow(y / 15., -0.06) + self.Cplantcorrelation = CCAPP1 * z * x * 1000. / 1E6 + + elif model.surfaceplant.plant_type.value == PlantType.SUPER_CRITICAL_ORC: + MaxProducedTemperature = np.max(model.surfaceplant.TenteringPP.value) + if MaxProducedTemperature < 150.: + 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.) + # factor 1.1 to make supercritical 10% more expansive than subcritical + self.Cplantcorrelation = 1.1 * CCAPP1 * math.pow( + np.max(model.surfaceplant.ElectricityProduced.value) / 15., -0.06) * np.max( + model.surfaceplant.ElectricityProduced.value) * 1000. / 1E6 + + elif model.surfaceplant.plant_type.value == PlantType.SINGLE_FLASH: + if np.max(model.surfaceplant.ElectricityProduced.value) < 10.: + C2 = 4.8472E-2 + C1 = -35.2186 + C0 = 8.4474E3 + D2 = 4.0604E-2 + D1 = -29.3817 + D0 = 6.9911E3 + PLL = 5. + PRL = 10. + elif np.max(model.surfaceplant.ElectricityProduced.value) < 25.: + C2 = 4.0604E-2 + C1 = -29.3817 + C0 = 6.9911E3 + D2 = 3.2773E-2 + D1 = -23.5519 + D0 = 5.5263E3 + PLL = 10. + PRL = 25. + elif np.max(model.surfaceplant.ElectricityProduced.value) < 50.: + C2 = 3.2773E-2 + C1 = -23.5519 + C0 = 5.5263E3 + D2 = 3.4716E-2 + D1 = -23.8139 + D0 = 5.1787E3 + PLL = 25. + PRL = 50. + elif np.max(model.surfaceplant.ElectricityProduced.value) < 75.: + C2 = 3.4716E-2 + C1 = -23.8139 + C0 = 5.1787E3 + D2 = 3.5271E-2 + D1 = -24.3962 + D0 = 5.1972E3 + PLL = 50. + PRL = 75. + else: + C2 = 3.5271E-2 + C1 = -24.3962 + C0 = 5.1972E3 + D2 = 3.3908E-2 + D1 = -23.4890 + D0 = 5.0238E3 + PLL = 75. + PRL = 100. + maxProdTemp = np.max(model.surfaceplant.TenteringPP.value) + CCAPPLL = C2 * maxProdTemp ** 2 + C1 * maxProdTemp + C0 + CCAPPRL = D2 * maxProdTemp ** 2 + D1 * maxProdTemp + D0 + b = math.log(CCAPPRL / CCAPPLL) / math.log(PRL / PLL) + a = CCAPPRL / PRL ** b + # factor 0.75 to make double flash 25% more expansive than single flash + self.Cplantcorrelation = (0.8 * a * math.pow(np.max(model.surfaceplant.ElectricityProduced.value), b) * + np.max(model.surfaceplant.ElectricityProduced.value) * 1000. / 1E6) + + elif model.surfaceplant.plant_type.value == PlantType.DOUBLE_FLASH: + if np.max(model.surfaceplant.ElectricityProduced.value) < 10.: + C2 = 4.8472E-2 + C1 = -35.2186 + C0 = 8.4474E3 + D2 = 4.0604E-2 + D1 = -29.3817 + D0 = 6.9911E3 + PLL = 5. + PRL = 10. + elif np.max(model.surfaceplant.ElectricityProduced.value) < 25.: + C2 = 4.0604E-2 + C1 = -29.3817 + C0 = 6.9911E3 + D2 = 3.2773E-2 + D1 = -23.5519 + D0 = 5.5263E3 + PLL = 10. + PRL = 25. + elif np.max(model.surfaceplant.ElectricityProduced.value) < 50.: + C2 = 3.2773E-2 + C1 = -23.5519 + C0 = 5.5263E3 + D2 = 3.4716E-2 + D1 = -23.8139 + D0 = 5.1787E3 + PLL = 25. + PRL = 50. + elif np.max(model.surfaceplant.ElectricityProduced.value) < 75.: + C2 = 3.4716E-2 + C1 = -23.8139 + C0 = 5.1787E3 + D2 = 3.5271E-2 + D1 = -24.3962 + D0 = 5.1972E3 + PLL = 50. + PRL = 75. + else: + C2 = 3.5271E-2 + C1 = -24.3962 + C0 = 5.1972E3 + D2 = 3.3908E-2 + D1 = -23.4890 + D0 = 5.0238E3 + PLL = 75. + PRL = 100. + maxProdTemp = np.max(model.surfaceplant.TenteringPP.value) + CCAPPLL = C2 * maxProdTemp ** 2 + C1 * maxProdTemp + C0 + CCAPPRL = D2 * maxProdTemp ** 2 + D1 * maxProdTemp + D0 + b = math.log(CCAPPRL / CCAPPLL) / math.log(PRL / PLL) + a = CCAPPRL / PRL ** b + self.Cplantcorrelation = (a * math.pow(np.max(model.surfaceplant.ElectricityProduced.value), b) * + np.max(model.surfaceplant.ElectricityProduced.value) * 1000. / 1E6) + + if self.ccplantfixed.Valid: + self.Cplant.value = self.ccplantfixed.value + self.CAPEX_cost_electricity_plant = self.Cplant.value * self.CAPEX_heat_electricity_plant_ratio.value + self.CAPEX_cost_heat_plant = self.Cplant.value * (1.0 - self.CAPEX_heat_electricity_plant_ratio.value) + else: + # 1.02 to convert cost from 2012 to 2016 #factor 1.15 for 15% contingency and 1.12 for 12% indirect costs. factor 1.10 to convert from 2016 to 2022 + self.Cplant.value = 1.12 * 1.15 * self.ccplantadjfactor.value * self.Cplantcorrelation * 1.02 * 1.10 + self.CAPEX_cost_electricity_plant = self.Cplant.value + + # add direct-use plant cost of co-gen system to Cplant (only of no total Cplant was provided) + if not self.ccplantfixed.Valid: # 1.15 below for contingency and 1.12 for indirect costs + if model.surfaceplant.enduse_option.value in [EndUseOptions.COGENERATION_TOPPING_EXTRA_ELECTRICITY, + EndUseOptions.COGENERATION_TOPPING_EXTRA_HEAT]: # enduse_option = 3: cogen topping cycle + self.CAPEX_cost_heat_plant = 1.12 * 1.15 * self.ccplantadjfactor.value * 250E-6 * np.max( + model.surfaceplant.HeatProduced.value / model.surfaceplant.enduse_efficiency_factor.value) * 1000. + elif model.surfaceplant.enduse_option.value in [EndUseOptions.COGENERATION_BOTTOMING_EXTRA_HEAT, + EndUseOptions.COGENERATION_BOTTOMING_EXTRA_ELECTRICITY]: # enduse_option = 4: cogen bottoming cycle + self.CAPEX_cost_heat_plant = 1.12 * 1.15 * self.ccplantadjfactor.value * 250E-6 * np.max( + model.surfaceplant.HeatProduced.value / model.surfaceplant.enduse_efficiency_factor.value) * 1000. + elif model.surfaceplant.enduse_option.value in [EndUseOptions.COGENERATION_PARALLEL_EXTRA_ELECTRICITY, + EndUseOptions.COGENERATION_PARALLEL_EXTRA_HEAT]: # cogen parallel cycle + self.CAPEX_cost_heat_plant = 1.12 * 1.15 * self.ccplantadjfactor.value * 250E-6 * np.max( + model.surfaceplant.HeatProduced.value / model.surfaceplant.enduse_efficiency_factor.value) * 1000. + + self.Cplant.value = self.Cplant.value + self.CAPEX_cost_heat_plant + if not self.CAPEX_heat_electricity_plant_ratio.Provided: + self.CAPEX_heat_electricity_plant_ratio.value = self.CAPEX_cost_electricity_plant/self.Cplant.value + + if not self.totalcapcost.Valid: + # exploration costs (same as in Geophires v1.2) (M$) + if self.ccexplfixed.Valid: + self.Cexpl.value = self.ccexplfixed.value + else: + self.Cexpl.value = 1.15 * self.ccexpladjfactor.value * 1.12 * ( + 1. + self.cost_one_production_well.value * 0.6) # 1.15 for 15% contingency and 1.12 for 12% indirect costs + + # Surface Piping Length Costs (M$) #assumed $750k/km + self.Cpiping.value = 750 / 1000 * model.surfaceplant.piping_length.value + + # district heating network costs + if model.surfaceplant.plant_type.value == PlantType.DISTRICT_HEATING: # district heat + if self.dhtotaldistrictnetworkcost.Provided: + self.dhdistrictcost.value = self.dhtotaldistrictnetworkcost.value + elif self.dhpipinglength.Provided: + self.dhdistrictcost.value = self.dhpipinglength.value * self.dhpipingcostrate.value / 1000 # M$ + elif self.dhroadlength.Provided: # check if road length is provided to calculate cost + self.dhdistrictcost.value = self.dhroadlength.value * 0.75 * self.dhpipingcostrate.value / 1000 # M$ (assuming 75% of road length is used for district network piping) + else: # calculate district network cost based on population density + if self.dhlandarea.Provided == False: + model.logger.warning("District heating network cost calculated based on default district area") + if self.dhpopulation.Provided: + self.populationdensity.value = self.dhpopulation.value / self.dhlandarea.value + elif model.surfaceplant.dh_number_of_housing_units.Provided: + self.populationdensity.value = model.surfaceplant.dh_number_of_housing_units.value * 2.6 / self.dhlandarea.value # estimate population based on 2.6 number of people per household + else: + model.logger.warning( + "District heating network cost calculated based on default number of people in district") + self.populationdensity.value = self.dhpopulation.value / self.dhlandarea.value + + if self.populationdensity.value > 1000: + self.dhpipinglength.value = 7.5 * self.dhlandarea.value # using constant 7.5km of pipe per km^2 when population density is >1500 + else: + self.dhpipinglength.value = max( + self.populationdensity.value / 1000 * 7.5 * self.dhlandarea.value, + self.dhlandarea.value) # scale the piping length based on population density, but with a minimum of 1 km of piping per km^2 of area + self.dhdistrictcost.value = self.dhpipingcostrate.value * self.dhpipinglength.value / 1000 + + else: + self.dhdistrictcost.value = 0 + + self.CCap.value = self.Cexpl.value + self.Cwell.value + self.Cstim.value + self.Cgath.value + self.Cplant.value + self.Cpiping.value + self.dhdistrictcost.value + else: + self.CCap.value = self.totalcapcost.value + + # update the capitol costs, assuming the entire ITC is used to reduce the capitol costs + if self.RITC.Provided: + self.RITCValue.value = self.RITC.value * self.CCap.value + self.CCap.value = self.CCap.value - self.RITCValue.value + + # Add in the FlatLicenseEtc, OtherIncentives, & TotalGrant + self.CCap.value = self.CCap.value + self.FlatLicenseEtc.value - self.OtherIncentives.value - self.TotalGrant.value + + # O&M costs + # calculate first O&M costs independent of whether oamtotalfixed is provided or not + # additional electricity cost for heat pump as end-use + if model.surfaceplant.plant_type.value == PlantType.HEAT_PUMP: # heat pump: + self.averageannualheatpumpelectricitycost.value = np.average( + model.surfaceplant.heat_pump_electricity_kwh_used.value) * model.surfaceplant.electricity_cost_to_buy.value / 1E6 # M$/year + + # district heating peaking fuel annual cost + if model.surfaceplant.plant_type.value == PlantType.DISTRICT_HEATING: # district heating + self.annualngcost.value = model.surfaceplant.annual_ng_demand.value * self.ngprice.value / 1000 / self.peakingboilerefficiency.value # array with annual O&M cost for peaking fuel + self.averageannualngcost.value = np.average(self.annualngcost.value) + + # calculate average annual pumping costs in case no electricity is provided + if model.surfaceplant.plant_type.value in [PlantType.INDUSTRIAL, PlantType.ABSORPTION_CHILLER, PlantType.HEAT_PUMP, PlantType.DISTRICT_HEATING]: + self.averageannualpumpingcosts.value = np.average(model.surfaceplant.PumpingkWh.value) * model.surfaceplant.electricity_cost_to_buy.value / 1E6 # M$/year + + if not self.oamtotalfixed.Valid: + # labor cost + if model.surfaceplant.enduse_option.value == EndUseOptions.ELECTRICITY: # electricity + if np.max(model.surfaceplant.ElectricityProduced.value) < 2.5: + self.Claborcorrelation = 236. / 1E3 # M$/year + else: + self.Claborcorrelation = (589. * math.log( + np.max(model.surfaceplant.ElectricityProduced.value)) - 304.) / 1E3 # M$/year + else: + if np.max(model.surfaceplant.HeatExtracted.value) < 2.5 * 5.: + self.Claborcorrelation = 236. / 1E3 # M$/year + else: + self.Claborcorrelation = (589. * math.log( + np.max(model.surfaceplant.HeatExtracted.value) / 5.) - 304.) / 1E3 # M$/year + # * 1.1 to convert from 2012 to 2016$ with BLS employment cost index (for utilities in March) + self.Claborcorrelation = self.Claborcorrelation * 1.1 + + # plant O&M cost + if self.oamplantfixed.Valid: + self.Coamplant.value = self.oamplantfixed.value + else: + self.Coamplant.value = self.oamplantadjfactor.value * ( + 1.5 / 100. * self.Cplant.value + 0.75 * self.Claborcorrelation) + + # wellfield O&M cost + if self.oamwellfixed.Valid: + self.Coamwell.value = self.oamwellfixed.value + else: + self.Coamwell.value = self.oamwelladjfactor.value * ( + 1. / 100. * (self.Cwell.value + self.Cgath.value) + 0.25 * self.Claborcorrelation) + + # water O&M cost + if self.oamwaterfixed.Valid: + self.Coamwater.value = self.oamwaterfixed.value + else: + # here is assumed 1 l per kg maybe correct with real temp. (M$/year) 925$/ML = 3.5$/1,000 gallon + self.Coamwater.value = self.oamwateradjfactor.value * (model.wellbores.nprod.value * + model.wellbores.prodwellflowrate.value * + model.reserv.waterloss.value * model.surfaceplant.utilization_factor.value * + 365. * 24. * 3600. / 1E6 * 925. / 1E6) + + # additional O&M cost for absorption chiller if used + if model.surfaceplant.plant_type.value == PlantType.ABSORPTION_CHILLER: # absorption chiller: + if self.chilleropex.value == -1: + self.chilleropex.value = self.chillercapex.value * 2 / 100 # assumed annual O&M for chiller is 2% of investment cost + + # correct plant O&M cost as otherwise chiller opex would be counted double (subtract chiller capex from plant cost when calculating Coandmplant) + if self.oamplantfixed.Valid == False: + self.Coamplant.value = self.oamplantadjfactor.value * ( + 1.5 / 100. * (self.Cplant.value - self.chillercapex.value) + 0.75 * self.Claborcorrelation) + + else: + self.chilleropex.value = 0 + + # district heating O&M cost + if model.surfaceplant.plant_type.value == PlantType.DISTRICT_HEATING: # district heating + self.annualngcost.value = model.surfaceplant.annual_ng_demand.value * self.ngprice.value / 1000 # array with annual O&M cost for peaking fuel + + if self.dhoandmcost.Provided: + self.dhdistrictoandmcost.value = self.dhoandmcost.value # M$/yr + else: + self.dhdistrictoandmcost.value = 0.01 * self.dhdistrictcost.value + 0.02 * sum( + model.surfaceplant.daily_heating_demand.value) * model.surfaceplant.electricity_cost_to_buy.value / 1000 # [M$/year] we assume annual district OPEX equals 1% of district CAPEX and 2% of total heat demand for pumping costs + + else: + self.dhdistrictoandmcost.value = 0 + + self.Coam.value = self.Coamwell.value + self.Coamplant.value + self.Coamwater.value + self.chilleropex.value + self.dhdistrictoandmcost.value # total O&M cost (M$/year) + + else: + self.Coam.value = self.oamtotalfixed.value # total O&M cost (M$/year) + + if model.wellbores.redrill.value > 0: + # account for well redrilling + self.Coam.value = self.Coam.value + \ + (self.Cwell.value + self.Cstim.value) * model.wellbores.redrill.value / model.surfaceplant.plant_lifetime.value + + # Add in the AnnualLicenseEtc and TaxRelief + self.Coam.value = self.Coam.value + self.AnnualLicenseEtc.value - self.TaxRelief.value + + # partition the OPEX for CHP plants based on the CAPEX ratio + self.OPEX_cost_electricity_plant = self.Coam.value * self.CAPEX_heat_electricity_plant_ratio.value + self.OPEX_cost_heat_plant = self.Coam.value * (1.0 - self.CAPEX_heat_electricity_plant_ratio.value) + + # The Reservoir depth measure was arbitrarily changed to meters despite being defined in the docs as kilometers. + # For display consistency sake, we need to convert it back + if model.reserv.depth.value > 500: + model.reserv.depth.value = model.reserv.depth.value / 1000.0 + model.reserv.depth.CurrentUnits = LengthUnit.KILOMETERS + + # build the PTC price models + self.PTCElecPrice = [0.0] * model.surfaceplant.plant_lifetime.value + self.PTCHeatPrice = [0.0] * model.surfaceplant.plant_lifetime.value + self.PTCCoolingPrice = [0.0] * model.surfaceplant.plant_lifetime.value + self.PTCCarbonPrice = [0.0] * model.surfaceplant.plant_lifetime.value + if self.PTCElec.Provided: + self.PTCElecPrice = BuildPTCModel(model.surfaceplant.plant_lifetime.value, + self.PTCDuration.value, self.PTCElec.value, self.PTCInflationAdjusted.value, + self.RINFL.value) + if self.PTCHeat.Provided: + self.PTCHeatPrice = BuildPTCModel(model.surfaceplant.plant_lifetime.value, + self.PTCDuration.value, self.PTCHeat.value, self.PTCInflationAdjusted.value, + self.RINFL.value) + if self.PTCCooling.Provided: + self.PTCCoolingPrice = BuildPTCModel(model.surfaceplant.plant_lifetime.value, + self.PTCDuration.value, self.PTCCooling.value, self.PTCInflationAdjusted.value, + self.RINFL.value) + + # build the price models + self.ElecPrice.value = BuildPricingModel(model.surfaceplant.plant_lifetime.value, + self.ElecStartPrice.value, self.ElecEndPrice.value, + self.ElecEscalationStart.value, self.ElecEscalationRate.value, + self.PTCElecPrice) + self.HeatPrice.value = BuildPricingModel(model.surfaceplant.plant_lifetime.value, + self.HeatStartPrice.value, self.HeatEndPrice.value, + self.HeatEscalationStart.value, self.HeatEscalationRate.value, + self.PTCHeatPrice) + self.CoolingPrice.value = BuildPricingModel(model.surfaceplant.plant_lifetime.value, + self.CoolingStartPrice.value, self.CoolingEndPrice.value, + self.CoolingEscalationStart.value, self.CoolingEscalationRate.value, + self.PTCCoolingPrice) + self.CarbonPrice.value = BuildPricingModel(model.surfaceplant.plant_lifetime.value, + self.CarbonStartPrice.value, self.CarbonEndPrice.value, + self.CarbonEscalationStart.value, self.CarbonEscalationRate.value, + self.PTCCarbonPrice) + + # do the additional economic calculations first, if needed, so the summaries below work. + if self.DoAddOnCalculations.value: + model.addeconomics.Calculate(model) + if self.DoSDACGTCalculations.value: + model.sdacgteconomics.Calculate(model) + + # Calculate cashflow and cumulative cash flow + total_duration = model.surfaceplant.plant_lifetime.value + model.surfaceplant.construction_years.value + self.ElecRevenue.value = [0.0] * total_duration + self.ElecCummRevenue.value = [0.0] * total_duration + self.HeatRevenue.value = [0.0] * total_duration + self.HeatCummRevenue.value = [0.0] * total_duration + self.CoolingRevenue.value = [0.0] * total_duration + self.CoolingCummRevenue.value = [0.0] * total_duration + self.CarbonRevenue.value = [0.0] * total_duration + self.CarbonCummCashFlow.value = [0.0] * total_duration + self.TotalRevenue.value = [0.0] * total_duration + self.TotalCummRevenue.value = [0.0] * total_duration + self.CarbonThatWouldHaveBeenProducedTotal.value = 0.0 + + # Based on the style of the project, calculate the revenue & cumulative revenue + if model.surfaceplant.enduse_option.value == EndUseOptions.ELECTRICITY: + self.ElecRevenue.value, self.ElecCummRevenue.value = CalculateRevenue( + model.surfaceplant.plant_lifetime.value, model.surfaceplant.construction_years.value, + model.surfaceplant.NetkWhProduced.value, self.ElecPrice.value) + self.TotalRevenue.value = self.ElecRevenue.value + #self.TotalCummRevenue.value = self.ElecCummRevenue.value + elif model.surfaceplant.enduse_option.value == EndUseOptions.HEAT and model.surfaceplant.plant_type.value not in [PlantType.ABSORPTION_CHILLER]: + self.HeatRevenue.value, self.HeatCummRevenue.value = CalculateRevenue( + model.surfaceplant.plant_lifetime.value, model.surfaceplant.construction_years.value, + model.surfaceplant.HeatkWhProduced.value, self.HeatPrice.value) + self.TotalRevenue.value = self.HeatRevenue.value + #self.TotalCummRevenue.value = self.HeatCummRevenue.value + elif model.surfaceplant.enduse_option.value == EndUseOptions.HEAT and model.surfaceplant.plant_type.value in [PlantType.ABSORPTION_CHILLER]: + self.CoolingRevenue.value, self.CoolingCummRevenue.value = CalculateRevenue( + model.surfaceplant.plant_lifetime.value, model.surfaceplant.construction_years.value, + model.surfaceplant.cooling_kWh_Produced.value, self.CoolingPrice.value) + self.TotalRevenue.value = self.CoolingRevenue.value + #self.TotalCummRevenue.value = self.CoolingCummRevenue.value + elif model.surfaceplant.enduse_option.value in [EndUseOptions.COGENERATION_TOPPING_EXTRA_HEAT, + EndUseOptions.COGENERATION_TOPPING_EXTRA_ELECTRICITY, + EndUseOptions.COGENERATION_BOTTOMING_EXTRA_ELECTRICITY, + EndUseOptions.COGENERATION_BOTTOMING_EXTRA_HEAT, + EndUseOptions.COGENERATION_PARALLEL_EXTRA_HEAT, + EndUseOptions.COGENERATION_PARALLEL_EXTRA_ELECTRICITY]: # co-gen + # else: + self.ElecRevenue.value, self.ElecCummRevenue.value = CalculateRevenue( + model.surfaceplant.plant_lifetime.value, model.surfaceplant.construction_years.value, + model.surfaceplant.NetkWhProduced.value, self.ElecPrice.value) + self.HeatRevenue.value, self.HeatCummRevenue.value = CalculateRevenue( + model.surfaceplant.plant_lifetime.value, model.surfaceplant.construction_years.value, + model.surfaceplant.HeatkWhProduced.value, self.HeatPrice.value) + + for i in range(0, model.surfaceplant.plant_lifetime.value + model.surfaceplant.construction_years.value, 1): + self.TotalRevenue.value[i] = self.ElecRevenue.value[i] + self.HeatRevenue.value[i] + #if i > 0: + # self.TotalCummRevenue.value[i] = self.TotalCummRevenue.value[i - 1] + self.TotalRevenue.value[i] + + if self.DoCarbonCalculations.value: + self.CarbonRevenue.value, self.CarbonCummCashFlow.value, self.CarbonThatWouldHaveBeenProducedAnnually.value, \ + self.CarbonThatWouldHaveBeenProducedTotal.value = CalculateCarbonRevenue(model, + model.surfaceplant.plant_lifetime.value, model.surfaceplant.construction_years.value, + self.CarbonPrice.value, self.GridCO2Intensity.value, self.NaturalGasCO2Intensity.value, + model.surfaceplant.NetkWhProduced.value, model.surfaceplant.HeatkWhProduced.value) + for i in range(model.surfaceplant.construction_years.value, model.surfaceplant.plant_lifetime.value + model.surfaceplant.construction_years.value, 1): + self.TotalRevenue.value[i] = self.TotalRevenue.value[i] + self.CarbonRevenue.value[i] + #self.TotalCummRevenue.value[i] = self.TotalCummRevenue.value[i] + self.CarbonCummCashFlow.value[i] + + # for the sake of display, insert zeros at the beginning of the pricing arrays + for i in range(0, model.surfaceplant.construction_years.value, 1): + self.ElecPrice.value.insert(0, 0.0) + self.HeatPrice.value.insert(0, 0.0) + self.CoolingPrice.value.insert(0, 0.0) + self.CarbonPrice.value.insert(0, 0.0) + + # Insert the cost of construction into the front of the array that will be used to calculate NPV + # the convention is that the upfront CAPEX is negative + # This is the same for all projects + ProjectCAPEXPerConstructionYear = self.CCap.value / model.surfaceplant.construction_years.value + for i in range(0, model.surfaceplant.construction_years.value, 1): + self.TotalRevenue.value[i] = -1.0 * ProjectCAPEXPerConstructionYear + self.TotalCummRevenue.value[i] = -1.0 * ProjectCAPEXPerConstructionYear +# self.TotalRevenue.value, self.TotalCummRevenue.value = CalculateTotalRevenue( +# model.surfaceplant.plant_lifetime.value, model.surfaceplant.construction_years.value, self.CCap.value, +# self.Coam.value, self.TotalRevenue.value, self.TotalCummRevenue.value) + + # Do a one-time calculation that accounts for OPEX - no OPEX in the first year. + for i in range(model.surfaceplant.construction_years.value, + model.surfaceplant.plant_lifetime.value + model.surfaceplant.construction_years.value, 1): + self.TotalRevenue.value[i] = self.TotalRevenue.value[i] - self.Coam.value + + # Now do a one-time calculation that calculates the cumulative cash flow after everything else has been accounted for + for i in range(1, model.surfaceplant.plant_lifetime.value + model.surfaceplant.construction_years.value, 1): + self.TotalCummRevenue.value[i] = self.TotalCummRevenue.value[i-1] + self.TotalRevenue.value[i] + + # Calculate more financial values using numpy financials + self.ProjectNPV.value, self.ProjectIRR.value, self.ProjectVIR.value, self.ProjectMOIC.value = \ + CalculateFinancialPerformance(model.surfaceplant.plant_lifetime.value, self.FixedInternalRate.value, + self.TotalRevenue.value, self.TotalCummRevenue.value, self.CCap.value, + self.Coam.value) + + # Calculate the project payback period + self.ProjectPaybackPeriod.value = 0.0 # start by assuming the project never pays back + for i in range(0, len(self.TotalCummRevenue.value), 1): + # find out when the cumm cashflow goes from negative to positive + if self.TotalCummRevenue.value[i] > 0 >= self.TotalCummRevenue.value[i - 1]: + # we just crossed the threshold into positive project cummcashflow, so we can calculate payback period + dFullDiff = self.TotalCummRevenue.value[i] + math.fabs(self.TotalCummRevenue.value[(i - 1)]) + dPerc = math.fabs(self.TotalCummRevenue.value[(i - 1)]) / dFullDiff + self.ProjectPaybackPeriod.value = i + dPerc + + + # Calculate LCOE/LCOH + self.LCOE.value, self.LCOH.value, self.LCOC.value = CalculateLCOELCOHLCOC(self, model) + + model.logger.info(f'complete {__class__!s}: {sys._getframe().f_code.co_name}') + diff --git a/src/geophires_x/SBTReservoir.py b/src/geophires_x/SBTReservoir.py new file mode 100644 index 00000000..41856951 --- /dev/null +++ b/src/geophires_x/SBTReservoir.py @@ -0,0 +1,1776 @@ +import sys +from functools import lru_cache + +import numpy as np +import pandas as pd +from scipy.special import erf, erfc, jv, yv, exp1 +from scipy.interpolate import interp1d +import matplotlib.pyplot as plt + +import geophires_x.Model as Model +from .CylindricalReservoir import CylindricalReservoir +from .OptionList import FlowrateModel, InjectionTemperatureModel, Configuration +from .Parameter import intParameter, floatParameter, OutputParameter, ReadParameter, strParameter, boolParameter +from .Reservoir import Reservoir +from .Units import * +import sys +from functools import lru_cache + + +def interpolator(time_value: float, times: np.ndarray, values: np.ndarray) -> float: + """ + Interpolates values based on time + :param time_value: Time value + :type time_value: float + :param values: Values + :type values: float + :return: Interpolated value + :rtype: float + """ + for i in range(1, len(times)): + if times[i] == time_value: + return values[i] + if time_value < times[i]: + time_diff = times[i] - times[i - 1] + value_diff = values[i] - values[i - 1] + ratio = (time_value - times[i - 1]) / time_diff + return values[i - 1] + ratio * value_diff + + +def generate_wireframe_model(lateral_endpoint_depth: float, number_of_laterals: int, lateral_spacing: float, element_length: float, + junction_depth:float, vertical_section_depth: float, angle: float, + vertical_well_spacing: float, generate_graphics: bool = False) -> tuple: + + # Generate inj well profile + zinj = np.linspace(0, -vertical_section_depth, round(vertical_section_depth / element_length)) + yinj = np.zeros(len(zinj)) + xinj = np.zeros(len(zinj)) + + inclined_length = abs(-junction_depth - zinj[-1]) / np.cos(angle) + number_of_elements_inclined_length = round(inclined_length / element_length) + + zinj_inclined_length = np.linspace(zinj[-1], -junction_depth, number_of_elements_inclined_length) + yinj_inclined_length = np.linspace(yinj[-1], yinj[-1], number_of_elements_inclined_length) + xinj_inclined_length = np.linspace(xinj[-1], inclined_length * np.sin(angle), number_of_elements_inclined_length) + + zinj = np.concatenate((zinj, zinj_inclined_length[1:])) + xinj = np.concatenate((xinj, xinj_inclined_length[1:])) + yinj = np.concatenate((yinj, yinj_inclined_length[1:])) + + # Generate prod well profile + zprod = np.flip(zinj) + xprod = np.flip(xinj) + yprod = np.flip(yinj + vertical_well_spacing) + + # Generate laterals + x_laterals_inj = np.zeros(number_of_laterals) + y_laterals_inj = np.zeros(number_of_laterals) + z_laterals_inj = np.zeros(number_of_laterals) + + for i in range(number_of_laterals): + y_laterals_inj[i] = yinj[-1] - (lateral_spacing * (number_of_laterals - 1)) / 2 + i * lateral_spacing - ( + yinj[-1] - yprod[0]) / 2 + x_laterals_inj[i] = xinj[-1] + element_length * 3 * np.sin(angle) + z_laterals_inj[i] = zinj[-1] - element_length * 3 * np.cos(angle) + + # Generate template + lateral_length = (lateral_endpoint_depth - abs(z_laterals_inj[0])) / np.cos(angle) + number_of_elements_lateral = round(lateral_length / element_length) + z_template_lateral = np.linspace(z_laterals_inj[-1], -lateral_endpoint_depth, number_of_elements_lateral) + x_template_lateral = np.linspace(x_laterals_inj[-1], x_laterals_inj[-1] + lateral_length * np.sin(angle), + number_of_elements_lateral) + y_template_lateral = np.full(number_of_elements_lateral, y_laterals_inj[-1]) + + # Add section upwards + z_template_lateral = np.concatenate((z_template_lateral, [z_template_lateral[-1] + element_length, z_template_lateral[-1] + element_length * 2])) + x_template_lateral = np.concatenate((x_template_lateral, [x_template_lateral[-1], x_template_lateral[-1]])) + y_template_lateral = np.concatenate((y_template_lateral, [y_template_lateral[-1], y_template_lateral[-1]])) + + # Add loop back + z_template_lateral = np.concatenate((z_template_lateral, np.flip(z_template_lateral[1:-3]) + element_length * 2)) + x_template_lateral = np.concatenate((x_template_lateral, np.flip(x_template_lateral[1:-3]))) + y_template_lateral = np.concatenate((y_template_lateral, np.flip(y_template_lateral[1:-3]))) + + # Generate feedways + xlat = np.zeros((3, number_of_laterals)) + ylat = np.zeros((3, number_of_laterals)) + zlat = np.zeros((3, number_of_laterals)) + + for i in range(number_of_laterals): + xlat[0:3, i] = np.linspace(xinj[-1], x_laterals_inj[i], 3) + ylat[0:3, i] = np.linspace(yinj[-1], y_laterals_inj[i], 3) + zlat[0:3, i] = np.linspace(zinj[-1], z_laterals_inj[i], 3) + + xlat = np.vstack((xlat, np.zeros((len(x_template_lateral) - 1, number_of_laterals)))) + ylat = np.vstack((ylat, np.zeros((len(x_template_lateral) - 1, number_of_laterals)))) + zlat = np.vstack((zlat, np.zeros((len(x_template_lateral) - 1, number_of_laterals)))) + + # Add template + for i in range(number_of_laterals): + adjusted_template = np.vstack((x_template_lateral, y_template_lateral, z_template_lateral)).T + np.array( + [x_laterals_inj[i], y_laterals_inj[i], z_laterals_inj[i]]) - np.array( + [x_template_lateral[0], y_template_lateral[0], z_template_lateral[0]]) + xlat[3:, i] = adjusted_template[1:, 0] + ylat[3:, i] = adjusted_template[1:, 1] + zlat[3:, i] = adjusted_template[1:, 2] + + # Add sections back to production point + xreturn = np.zeros((3, number_of_laterals)) + yreturn = np.zeros((3, number_of_laterals)) + zreturn = np.zeros((3, number_of_laterals)) + + for i in range(number_of_laterals): + xreturn[:, i] = np.linspace(xlat[-1, i], xprod[0], 3) + yreturn[:, i] = np.linspace(ylat[-1, i], yprod[0], 3) + zreturn[:, i] = np.linspace(zlat[-1, i], zprod[0], 3) + + xlat = np.vstack((xlat, xreturn[1:, :])) + ylat = np.vstack((ylat, yreturn[1:, :])) + zlat = np.vstack((zlat, zreturn[1:, :])) + + if generate_graphics: + # Plot profile + fig = plt.figure() + ax = fig.add_subplot(111, projection='3d') + ax.plot(xinj, yinj, zinj, 'b-o', linewidth=2) + ax.plot(xprod, yprod, zprod, 'r-o', linewidth=2) + + for i in range(number_of_laterals): + ax.plot(xlat[:, i], ylat[:, i], zlat[:, i], 'k-o', linewidth=2) + + ax.set_facecolor('white') + ax.tick_params(axis='both', which='major', labelsize=22) + ax.set_xlabel('x (m)', fontsize=22) + ax.set_ylabel('y (m)', fontsize=22) + ax.set_zlabel('Depth (m)', fontsize=22) + # ax.legend(['Injection Well', 'Production Well', 'Lateral(s)'], fontsize=22) # Uncomment to add legend + + az, el = 71.5676, 10.4739 + ax.view_init(az, el) + plt.show() + + return xinj, yinj, zinj, xprod, yprod, zprod, xlat, ylat, zlat + + +class SBTReservoir(CylindricalReservoir): + """ + The following code calculates the temperature profile in a U-shaped geothermal well using the SBT algorithm. + """ + + 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 Reservoir. + :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(f'Init {str(__class__)}: {sys._getframe().f_code.co_name}') + super().__init__(model) + + # 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 the 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. + # If you choose to subclass this master class, you can do so before or after you create your own parameters. + # If you do, you can also choose to call this method from you class, which will effectively add and set all + # these parameters to your class. + + # 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. + + self.flow_rate_model = self.ParameterDict[self.flow_rate_model.Name] = intParameter( + "Flowrate Model", + DefaultValue=FlowrateModel.USER_SUPPLIED, + AllowableRange=[1, 2], + UnitType=Units.NONE, + ErrMessage="assume constant user-provided flowrate (1)", + ToolTipText="Must be 1 or 2. '1' means the user provides a constant mass flow rate. \ + '1' means the user provides an excel file with a mass flow rate profile." + ) + self.flow_rate_file = self.ParameterDict[self.flow_rate_file.Name] = strParameter( + "Flowrate File", + DefaultValue="", + UnitType=Units.NONE, + ErrMessage="assume no flowrate file", + ToolTipText="Excel file with a mass flow rate profile" + ) + self.injection_temperature_model = self.ParameterDict[self.injection_temperature_model.Name] = intParameter( + "Injection Temperature Model", + DefaultValue=InjectionTemperatureModel.USER_SUPPLIED, + AllowableRange=[1, 2], + UnitType=Units.NONE, + ErrMessage="assume constant user-provided injection temperature (1)", + ToolTipText="Must be 1 or 2. '1' means the user provides a constant injection temperature. \ + '1' means the user provides an excel file with an injection temperature profile." + ) + self.injection_temperature_file = self.ParameterDict[self.injection_temperature_file.Name] = strParameter( + "Injection Temperature File", + DefaultValue="", + UnitType=Units.NONE, + ErrMessage="assume no injection temperature file", + ToolTipText="Excel file with an injection temperature profile" + ) + self.SBTAccuracyDesired = self.ParameterDict[self.SBTAccuracyDesired.Name] = intParameter( + "SBT Accuracy Desired", + DefaultValue=1, + AllowableRange=[1, 5], + UnitType=Units.NONE, + ErrMessage="assume default SBT accuracy desired (1)", + ToolTipText="Must be 1, 2, 3, 4 or 5 with 1 lowest accuracy and 5 highest accuracy. \ + Lowest accuracy runs fastest. Accuracy level impacts number of discretizations for \ + numerical integration and decision tree thresholds in SBT algorithm." + ) + self.percent_implicit = self.ParameterDict[self.percent_implicit.Name] = floatParameter( + "SBT Percent Implicit Euler Scheme", + DefaultValue=1.0, + Min=0.0, + Max=1.0, + UnitType=Units.PERCENT, + PreferredUnits=PercentUnit.TENTH, + CurrentUnits=PercentUnit.TENTH, + ErrMessage="assume default percent implicit (1.0)", + ToolTipText="Should be between 0 and 1. Most stable is setting it to 1 which results in \ + a fully implicit Euler scheme when calculating the fluid temperature at each time step. \ + With a value of 0, the convective term is modelled using explicit Euler. \ + A value of 0.5 would model the convective term 50% explicit and 50% implicit, \ + which may be slightly more accurate than fully implicit." + ) + self.initial_timestep_count = self.ParameterDict[self.initial_timestep_count.Name] = intParameter( + 'SBT Initial Timestep Count', + DefaultValue=5, + AllowableRange = [1,150], + UnitType=Units.NONE, + ErrMessage='assume default for Initial Timestep Count (5)', + ToolTipText='The number of timesteps in the first ~3 hours of model' + ) + self.final_timestep_count = self.ParameterDict[self.final_timestep_count.Name] = floatParameter( + 'SBT Final Timestep Count', + DefaultValue=70, + Min=5, + Max=1000, + UnitType=Units.NONE, + ErrMessage='assume default for Final Timestep Count 70)', + ToolTipText='The number of timesteps after the first ~3 hours of model' + ) + self.initial_final_timestep_transition = self.ParameterDict[self.initial_final_timestep_transition.Name] = floatParameter( + 'SBT Initial to Final Timestep Transition', + DefaultValue=9900, + Min=1, + Max=40_000_000, + UnitType=Units.TIME, + PreferredUnits=TimeUnit.SECOND, + CurrentUnits=TimeUnit.SECOND, + ErrMessage='assume default for Initial to Final Timestep Transition (9900 seconds)', + ToolTipText='The time in secs at which the time arrays switches from closely spaced linear to logarithmic' + ) + self.generate_wireframe_graphics = self.ParameterDict[self.generate_wireframe_graphics.Name] = boolParameter( + 'SBT Generate Wireframe Graphics', + DefaultValue=False, + UnitType=Units.NONE, + ErrMessage='assume default for SBT Generate Wireframe Graphics (False)', + ToolTipText='Switch to control the generation of a wireframe drawing of a SBT wells configuration' + ) + + sclass = str(__class__).replace("", "") + self.MyPath = os.path.abspath(__file__) + + # Saved Outputs + + self.NonLinearTime_temperature = self.OutputParameterDict[self.NonLinearTime_temperature.Name] = OutputParameter( + "NonLinear Time vs Temperature", + UnitType=Units.NONE + ) + + model.logger.info(f'Complete {str(__class__)}: {sys._getframe().f_code.co_name}') + + def __str__(self): + return "SBTReservoir" + + def read_parameters(self, model: Model) -> None: + """ + The read_parameters function reads in the parameters from a dictionary created by reading the user-provided file + and updates the parameter values for this object. + The function reads in all the parameters that relate to this object, including those that are inherited from + other objects. It then updates any of these parameter values that have been changed by the user. + It also handles any special cases. + :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(f'Init {str(__class__)}: {sys._getframe().f_code.co_name}') + super().read_parameters(model) + + # Deal with all the parameter values that the user has provided. They should really only provide values + # that they want to change from the default values, but they can provide a value that is already set + # because it is a default value set in __init__. It will ignore those. + # This also deals with all the special cases that need to be taken care of + # 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, do it before or after you call you own version of this method. If you do, you can + # also choose to call this method from you class, which can effectively modify all these + # superclass parameters in your class. + + if len(model.InputParameters) > 0: + # loop through all the parameters that the user wishes to set, looking for parameters that match this object + for item in self.ParameterDict.items(): + ParameterToModify = item[1] + key = ParameterToModify.Name.strip() + if key in model.InputParameters: + ParameterReadIn = model.InputParameters[key] + + # Before we change the parameter, let's assume that the unit preferences will match - + # if they don't, the later code will fix this. + # TODO: refactor GEOPHIRES such that parameters are read in immutably and only accessed with + # explicit units, with conversion only occurring in the getter as necessary + + ParameterToModify.CurrentUnits = ParameterToModify.PreferredUnits + ReadParameter(ParameterReadIn, ParameterToModify, model) # this handles all non-special cases + + # handle special cases + if ParameterToModify.Name == "Flowrate Model": + if ParameterReadIn.sValue == '1': + ParameterToModify.value = FlowrateModel.USER_SUPPLIED + elif ParameterReadIn.sValue == '2': + ParameterToModify.value = FlowrateModel.FILE_SUPPLIED + + elif ParameterToModify.Name == 'Injection Temperature Model': + if ParameterReadIn.sValue == '1': + ParameterToModify.value = InjectionTemperatureModel.USER_SUPPLIED + elif ParameterReadIn.sValue == '2': + ParameterToModify.value = InjectionTemperatureModel.FILE_SUPPLIED + else: + model.logger.info("No parameters read because no content provided") + + model.logger.info(f'complete {str(__class__)}: {sys._getframe().f_code.co_name}') + + def Calculate(self, model): + """ + The Calculate function is the main function that is called to run the calculations for this object. + In this case, it just calls the appropriate function based on the configuration of the wellbores. + """ + model.logger.info(f'Init {str(__class__)}: {sys._getframe().f_code.co_name}') + + # calculate the reservoir depth because there is no one single depth as there are for other reservoirs. + # Make the depth at which this pressure is calculates to be the average of the junction box depth and + # the lateral endpoint depths + self.depth.value = (model.wellbores.junction_depth.quantity().to('km').magnitude + + model.wellbores.lateral_endpoint_depth.quantity().to('km').magnitude) / 2.0 + + if model.wellbores.Configuration.value in [Configuration.ULOOP, Configuration.EAVORLOOP]: + self.Calculate_Uloop(model) + elif model.wellbores.Configuration.value == Configuration.COAXIAL: + self.Calculate_Coaxial(model) + else: + return + model.logger.info(f'complete {str(__class__)}: {sys._getframe().f_code.co_name}') + + + @lru_cache(maxsize=1024) + #@profile + def Calculate_Coaxial(self, model): + """ + Calculate the coaxial version of the SBT model + """ + model.logger.info(f'Init {str(__class__)}: {sys._getframe().f_code.co_name}') + + # Clear all equivalent: Initialize variables and import necessary libraries + + # SBT v2 for co-axial heat exchanger with high-temperature capability + # last update: January 14th, 2024 + + # 1. Input + # Generally, the user should only make changes to this section + + fluid = 1 # Heat transfer fluid selection: 1 = H2O; 2 = CO2 + m = 20 # Total fluid mass flow rate [kg/s] + Tin = 50 # Fluid input temperature [°C] + Pin = 100 # Fluid input pressure [bar] + Tsurf = 20 # Surface temperature [°C] + GeoGradient = 100 / 1000 # Geothermal gradient [°C/m] + self.krock.value = 2.83 # Rock thermal conductivity [W/m·K] + k_m_boiler = self.krock.value * 6 # Thermal conductivity in boiler [W/m·K] + self.cprock.value = 825 # Rock specific heat capacity [J/kg·K] + self.rhorock.value = 2875 # Rock density [kg/m³] + radius = 0.2286 / 2 # Wellbore radius [m] + radiuscenterpipe = 0.127 / 2 # Inner radius of inner pipe [m] + thicknesscenterpipe = 0.0127 # Thickness of inner pipe [m] + k_center_pipe = 0.01 # Thermal conductivity of insulation of center pipe wall [W/m·K] + perform_interpipe_heat_exchange_iteration = 0 # Perform additional iteration loop if problem with converging heat exchange + coaxialflowtype = 1 # 1 = CXA (fluid injection in annulus); 2 = CXC (fluid injection in center pipe) + piperoughness = 10 ** -6 # Pipe roughness to calculate friction pressure drop [m] + times = np.concatenate( + (np.arange(0, 1000, 100), np.logspace(np.log10(100 * 100), np.log10(20 * 365 * 24 * 3600), 75))) + reltolerance = 1e-5 # Target maximum acceptable relative tolerance each time step + maxnumberofiterations = 25 # Maximum number of iterations each time step + variablefluidproperties = 1 # 0 means constant fluid properties, 1 means properties vary with temperature and pressure + + # Constant fluid properties if variablefluidproperties is 0 + if variablefluidproperties == 0: + rho_f = 996 # Density of the fluid [kg/m³] + cp_f = 4177 # Specific heat capacity of the fluid [J/kg·K] + k_f = 0.615 # Thermal conductivity of the fluid [W/m·K] + mu_f = 7.97 * 10 ** -4 # Dynamic viscosity of the fluid [Pa·s] + + initialtemperatureprofile = 1 # 0 to calculate initial temperature using Tsurf and GeoGradient, 1 to provide array initialtemperaturedata + initialtemperaturedata = np.array([ + [0, 20], + [-1499.9, 600], + [-1500, 1000], + [-2000, 1000] + ]) + + # Coordinates of the centerline of the co-axial heat exchanger + # MIR z = np.arange(0, -2050, -50) + z = np.arange(0, -2050, -50) + x = np.zeros(len(z)) + y = np.zeros(len(z)) + + # Specify boiler elements + boilerelements = np.arange(31, 41) + + # Make 3D figure of borehole geometry to make sure it looks correct + # fig = plt.figure() + # ax = fig.add_subplot(111, projection='3d') + # ax.plot3D(x, y, z, 'k-o', linewidth=2) + # ax.set_xlabel('x (m)') + # ax.set_ylabel('y (m)') + # x.set_zlabel('Depth (m)') + # x.set_title('Co-axial heat exchanger geometry') + # plt.grid() + # plt.show() + + # 2. Pre-Processing + print('Start pre-processing ...') + g = 9.81 # Gravitational acceleration [m/s²] + gamma = 0.577215665 # Euler's constant + alpha_m = self.krock.value / (self.rhorock.value * self.cprock.value) # Rock thermal diffusivity [m²/s] + alpha_m_boiler = self.krock.value_boiler / (self.rhorock.value * self.cprock.value) # Boiler rock thermal diffusivity [m²/s] + + outerradiuscenterpipe = radiuscenterpipe + thicknesscenterpipe # Outer radius of inner pipe [m] + A_flow_annulus = np.pi * (radius ** 2 - outerradiuscenterpipe ** 2) # Flow area of annulus pipe [m²] + A_flow_centerpipe = np.pi * radiuscenterpipe ** 2 # Flow area of center pipe [m²] + Dh_annulus = 2 * (radius - outerradiuscenterpipe) # Hydraulic diameter of annulus [m] + eps_annulus = Dh_annulus * piperoughness # Relative roughness annulus [-] + eps_centerpipe = 2 * radiuscenterpipe * piperoughness # Relative roughness inner pipe [-] + + # Variable fluid properties logic + if variablefluidproperties == 0: + Pvector = [1, 1e9] + Tvector = [1, 1e4] + density = np.full((2, 2), rho_f) + heatcapacity = np.full((2, 2), cp_f) + thermalconductivity = np.full((2, 2), k_f) + viscosity = np.full((2, 2), mu_f) + thermalexpansion = np.zeros((2, 2)) + else: + print('Loading fluid properties ...') + if fluid == 1: + # Load properties for water from pre-generated CoolProp data + properties = sio.loadmat('properties_H2O_HT_v3.mat') + print('Fluid properties for water loaded successfully') + elif fluid == 2: + # Load properties for CO2 from pre-generated CoolProp data + properties = sio.loadmat('properties_CO2.mat') + print('Fluid properties for CO2 loaded successfully') + else: + print('No valid fluid selected') + exit() + + # Length of each segment [m] + Deltaz = np.sqrt(np.diff(x) ** 2 + np.diff(y) ** 2 + np.diff(z) ** 2) + TotalLength = np.sum(Deltaz) # Total length of co-axial heat exchanger [m] + + # Geometry Quality Control + LoverR = Deltaz / radius + smallestLoverR = np.min(LoverR) + if smallestLoverR < 10: + print( + 'Warning: smallest ratio of segment length over radius is less than 10. Good practice is to keep this ratio larger than 10.') + + RelativeLengthChanges = np.diff(Deltaz) / Deltaz[:-1] + if np.max(np.abs(RelativeLengthChanges)) > 0.5: + print( + 'Warning: abrupt change(s) in segment length detected, which may cause numerical instabilities. Good practice is to avoid abrupt length changes to obtain smooth results.') + + if self.SBTAccuracyDesired.value == 1: + NoArgumentsFinitePipeCorrection = 25 + NoDiscrFinitePipeCorrection = 200 + NoArgumentsInfCylIntegration = 25 + NoDiscrInfCylIntegration = 200 + LimitPointSourceModel = 1.5 + LimitCylinderModelRequired = 25 + LimitInfiniteModel = 0.05 + LimitNPSpacingTime = 0.1 + LimitSoverL = 1.5 + M = 3 + elif self.SBTAccuracyDesired.value == 2: + NoArgumentsFinitePipeCorrection = 50 + NoDiscrFinitePipeCorrection = 400 + NoArgumentsInfCylIntegration = 50 + NoDiscrInfCylIntegration = 400 + LimitPointSourceModel = 2.5 + LimitCylinderModelRequired = 50 + LimitInfiniteModel = 0.01 + LimitNPSpacingTime = 0.04 + LimitSoverL = 2 + M = 4 + elif self.SBTAccuracyDesired.value == 3: + NoArgumentsFinitePipeCorrection = 100 + NoDiscrFinitePipeCorrection = 500 + NoArgumentsInfCylIntegration = 100 + NoDiscrInfCylIntegration = 500 + LimitPointSourceModel = 5 + LimitCylinderModelRequired = 100 + LimitInfiniteModel = 0.004 + LimitNPSpacingTime = 0.02 + LimitSoverL = 3 + M = 5 + elif self.SBTAccuracyDesired.value == 4: + NoArgumentsFinitePipeCorrection = 200 + NoDiscrFinitePipeCorrection = 1000 + NoArgumentsInfCylIntegration = 200 + NoDiscrInfCylIntegration = 1000 + LimitPointSourceModel = 10 + LimitCylinderModelRequired = 200 + LimitInfiniteModel = 0.002 + LimitNPSpacingTime = 0.01 + LimitSoverL = 5 + M = 10 + elif self.SBTAccuracyDesired.value == 5: + NoArgumentsFinitePipeCorrection = 400 + NoDiscrFinitePipeCorrection = 2000 + NoArgumentsInfCylIntegration = 400 + NoDiscrInfCylIntegration = 2000 + LimitPointSourceModel = 20 + LimitCylinderModelRequired = 400 + LimitInfiniteModel = 0.001 + LimitNPSpacingTime = 0.005 + LimitSoverL = 9 + M = 20 + + # Use alpha_m instead of alpha_m_boiler for more conservative estimates + timeforpointssource = np.max( + Deltaz) ** 2 / alpha_m * LimitPointSourceModel # Minimum time step size when point source model becomes applicable [s] + timeforlinesource = radius ** 2 / alpha_m * LimitCylinderModelRequired # Minimum time step size when line source model becomes applicable [s] + timeforfinitelinesource = np.max( + Deltaz) ** 2 / alpha_m * LimitInfiniteModel # Minimum time step size when finite line source model should be considered [s] + + print('Precalculate SBT distributions ...') + # Precalculate the thermal response with a line and cylindrical heat source + # Precalculate finite pipe correction + fpcminarg = min(Deltaz) ** 2 / (4 * alpha_m * times[-1]) + fpcmaxarg = max(Deltaz) ** 2 / (4 * alpha_m * min(np.diff(times))) + Amin1vector = np.logspace(np.log10(fpcminarg) - 0.1, np.log10(fpcmaxarg) + 0.1, NoArgumentsFinitePipeCorrection) + finitecorrectiony = np.zeros(NoArgumentsFinitePipeCorrection) + + for i in range(len(Amin1vector)): + Amin1 = Amin1vector[i] + Amax1 = 16 ** 2 + if Amin1 > Amax1: + Amax1 = 10 * Amin1 + Adomain1 = np.logspace(np.log10(Amin1), np.log10(Amax1), NoDiscrFinitePipeCorrection) + finitecorrectiony[i] = trapezoid(-1. / (Adomain1 * 4 * np.pi * self.krock.value) * (erfc(0.5 * Adomain1 ** 0.5)), + Adomain1) + + # Precalculate Bessel integration for infinite cylinder + # MIR besselminarg = alpha_m * min(np.diff(times)) / max(radius)**2 + # MIR besselmaxarg = alpha_m * timeforlinesource / min(radius)**2 + besselminarg = alpha_m * min(np.diff(times)) / radius ** 2 + besselmaxarg = alpha_m * timeforlinesource / radius ** 2 + deltazbessel = np.logspace(-10, 8, NoDiscrInfCylIntegration) + argumentbesselvec = np.logspace(np.log10(besselminarg) - 0.5, np.log10(besselmaxarg) + 0.5, + NoArgumentsInfCylIntegration) + besselcylinderresult = np.zeros(NoArgumentsInfCylIntegration) + + for i in range(len(argumentbesselvec)): + argumentbessel = argumentbesselvec[i] + besselcylinderresult[i] = (2 / self.krock.value / np.pi ** 3 * + trapezoid((1 - np.exp(-deltazbessel ** 2 * argumentbessel)) / + (deltazbessel ** 3 * ( + jv(1, deltazbessel) ** 2 + yv(1, deltazbessel) ** 2)), + deltazbessel)) + + # Optional: Uncomment the lines below to plot the results + # import matplotlib.pyplot as plt + # plt.loglog(Amin1vector, finitecorrectiony, 'k-') + # plt.loglog(argumentbesselvec, besselcylinderresult, 'r-') + # plt.show() + + # Precalculate finite pipe correction for boiler + fpcminarg = min(Deltaz) ** 2 / (4 * alpha_m_boiler * times[-1]) + fpcmaxarg = max(Deltaz) ** 2 / (4 * alpha_m_boiler * min(np.diff(times))) + Amin1vector_boiler = np.logspace(np.log10(fpcminarg) - 0.1, np.log10(fpcmaxarg) + 0.1, + NoArgumentsFinitePipeCorrection) + finitecorrectiony_boiler = np.zeros(NoArgumentsFinitePipeCorrection) + + for i in range(len(Amin1vector_boiler)): + Amin1 = Amin1vector_boiler[i] + Amax1 = 16 ** 2 + if Amin1 > Amax1: + Amax1 = 10 * Amin1 + Adomain1 = np.logspace(np.log10(Amin1), np.log10(Amax1), NoDiscrFinitePipeCorrection) + finitecorrectiony_boiler[i] = trapezoid( + -1. / (Adomain1 * 4 * np.pi * k_m_boiler) * (erfc(0.5 * Adomain1 ** 0.5)), Adomain1) + + # Precalculate Bessel integration for infinite cylinder for boiler + # MIR besselminarg = alpha_m_boiler * min(np.diff(times)) / max(radius)**2 + # MIR besselmaxarg = alpha_m_boiler * timeforlinesource / min(radius)**2 + besselminarg = alpha_m_boiler * min(np.diff(times)) / radius ** 2 + besselmaxarg = alpha_m_boiler * timeforlinesource / radius ** 2 + deltazbessel = np.logspace(-10, 8, NoDiscrInfCylIntegration) + argumentbesselvec_boiler = np.logspace(np.log10(besselminarg) - 0.5, np.log10(besselmaxarg) + 0.5, + NoArgumentsInfCylIntegration) + besselcylinderresult_boiler = np.zeros(NoArgumentsInfCylIntegration) + + for i in range(len(argumentbesselvec_boiler)): + argumentbessel = argumentbesselvec_boiler[i] + besselcylinderresult_boiler[i] = (2 / k_m_boiler / np.pi ** 3 * + trapezoid((1 - np.exp(-deltazbessel ** 2 * argumentbessel)) / + (deltazbessel ** 3 * ( + jv(1, deltazbessel) ** 2 + yv(1, deltazbessel) ** 2)), + deltazbessel)) + + # Optional: Uncomment the lines below to plot the results + # import matplotlib.pyplot as plt + # plt.loglog(Amin1vector_boiler, finitecorrectiony_boiler, 'k-') + # plt.loglog(argumentbesselvec_boiler, besselcylinderresult_boiler, 'r-') + # plt.show() + print('SBT distributions calculated successfully') + + # MIR N = len(x) - 1 + N = len(x) - 1 + Nboiler = len(boilerelements) + # MIR Nreg = N - Nboiler + Nreg = 1 + N - Nboiler + self.krock.value_vector = np.full(N, self.krock.value) + k_m_vector[Nreg:] = k_m_boiler + alpha_m_vector = k_m_vector / self.rhorock.value / self.cprock.value + + SMatrix = np.zeros((N, N)) + for i in range(N): + SMatrix[i, :] = np.sqrt((0.5 * x[i] + 0.5 * x[i + 1] - 0.5 * x[:-1] - 0.5 * x[1:]) ** 2 + + (0.5 * y[i] + 0.5 * y[i + 1] - 0.5 * y[:-1] - 0.5 * y[1:]) ** 2 + + (0.5 * z[i] + 0.5 * z[i + 1] - 0.5 * z[:-1] - 0.5 * z[1:]) ** 2) + SoverL = SMatrix / (np.ones((N, 1)) * Deltaz) + SMatrixSorted = np.sort(SMatrix, axis=1) + SortedIndices = np.argsort(SMatrix, axis=1) + SoverLSorted = SMatrixSorted / (np.ones((N, 1)) * Deltaz) + mindexNPCP = np.argmax(np.min(SoverLSorted, axis=1) < LimitSoverL) + + # MIR midpointsx = 0.5 * x[1:] + 0.5 * x[:-1] + # MIR midpointsy = 0.5 * y[1:] + 0.5 * y[:-1] + # MIR midpointsz = 0.5 * z[1:] + 0.5 * z[:-1] + midpointsx = 0.5 * x + 0.5 * x + midpointsy = 0.5 * y + 0.5 * y + midpointsz = 0.5 * z + 0.5 * z + verticalchange = np.diff(z) + # MIR + verticalchange = np.append(verticalchange, verticalchange[-1]) + + if initialtemperatureprofile == 0: + BBinitial = Tsurf - GeoGradient * midpointsz + Tfluidupnodes = Tsurf - GeoGradient * z + Tfluiddownnodes = Tsurf - GeoGradient * z + elif initialtemperatureprofile == 1: + BBinitial = np.interp(midpointsz, initialtemperaturedata[:, 0], initialtemperaturedata[:, 1]) + Tfluidupnodes = np.interp(z, initialtemperaturedata[:, 0], initialtemperaturedata[:, 1]) + Tfluiddownnodes = np.interp(z, initialtemperaturedata[:, 0], initialtemperaturedata[:, 1]) + + # MIR Tfluiddownmidpoints = 0.5 * Tfluiddownnodes[1:] + 0.5 * Tfluiddownnodes[:-1] + # MIR Tfluidupmidpoints = 0.5 * Tfluidupnodes[1:] + 0.5 * Tfluidupnodes[:-1] + Tfluiddownmidpoints = 0.5 * Tfluiddownnodes + 0.5 * Tfluiddownnodes + Tfluidupmidpoints = 0.5 * Tfluidupnodes + 0.5 * Tfluidupnodes + + MaxSMatrixSorted = np.max(SMatrixSorted, axis=1) + indicesyoucanneglectupfront = alpha_m * np.outer(np.ones(N - 1), times) / ( + MaxSMatrixSorted[1:, np.newaxis] * np.ones((1, len(times)))) ** 2 / LimitNPSpacingTime + indicesyoucanneglectupfront[indicesyoucanneglectupfront > 1] = 1 + + lastneighbourtoconsider = np.zeros(len(times), dtype=int) + for i in range(len(times)): + lntc = np.where(indicesyoucanneglectupfront[:, i] == 1)[0] + if len(lntc) == 0: + lastneighbourtoconsider[i] = 1 + else: + lastneighbourtoconsider[i] = max(2, lntc[-1]) + + distributionx = np.array([np.linspace(x[i], x[i + 1], M + 1) for i in range(N)]) + distributiony = np.array([np.linspace(y[i], y[i + 1], M + 1) for i in range(N)]) + distributionz = np.array([np.linspace(z[i], z[i + 1], M + 1) for i in range(N)]) + + # Calculate initial pressure distribution + if fluid == 1: # H2O + Pfluidupnodes = Pin * 1e5 - 1000 * g * z + Pfluiddownnodes = Pfluidupnodes + Pfluidupmidpoints = Pin * 1e5 - 1000 * g * midpointsz + Pfluiddownmidpoints = Pfluidupmidpoints + elif fluid == 2: # CO2 + Pfluidupnodes = Pin * 1e5 - 500 * g * z + Pfluiddownnodes = Pfluidupnodes + Pfluidupmidpoints = Pin * 1e5 - 500 * g * midpointsz + Pfluiddownmidpoints = Pfluidupmidpoints + + kk = 1 + maxrelativechange = 1 + print(f'Calculating initial pressure field ... | Iteration = 1') + while kk < maxnumberofiterations and maxrelativechange > reltolerance: + Pfluidupmidpointsold = Pfluidupmidpoints.copy() + Pfluiddownmidpointsold = Pfluiddownmidpoints.copy() + densityfluidupmidpoints = CP.PropsSI('D', 'P', Pfluidupmidpoints, 'T', Tin + 273.15, 'Water') + densityfluiddownmidpoints = densityfluidupmidpoints + Pfluiddownnodes = Pin * 1e5 - np.cumsum([0] + g * verticalchange * densityfluiddownmidpoints) + Pfluidupnodes = Pfluiddownnodes + # MIR Pfluiddownmidpoints = 0.5 * Pfluiddownnodes[1:] + 0.5 * Pfluiddownnodes[:-1] + Pfluiddownmidpoints = 0.5 * Pfluiddownnodes + 0.5 * Pfluiddownnodes + Pfluidupmidpoints = Pfluiddownmidpoints + maxrelativechange = np.max(np.abs((Pfluiddownmidpointsold - Pfluiddownmidpoints) / Pfluiddownmidpointsold)) + kk += 1 + print( + f'Calculating initial pressure field ... | Iteration = {kk} | Max. Rel. change = {maxrelativechange:.5f}') + + densityfluiddownnodes = CP.PropsSI('D', 'P', Pfluiddownnodes, 'T', Tfluiddownnodes + 273.15, 'Water') + densityfluidupnodes = densityfluiddownnodes + + # HT Contribution + Phasefluiddownnodes = CP.PropsSI('Phase', 'P', Pfluiddownnodes, 'T', Tfluiddownnodes + 273.15, 'Water') + Phasefluidupnodes = CP.PropsSI('Phase', 'P', Pfluidupnodes, 'T', Tfluidupnodes + 273.15, 'Water') + Qfluiddownnodes = np.where(Phasefluiddownnodes == 0, 0, 1) + Qfluidupnodes = np.where(Phasefluidupnodes == 0, 0, 1) + + if maxrelativechange < reltolerance: + print('Initial pressure field calculated successfully') + else: + print('Initial pressure field calculated but maximum relative tolerance not met') + + # Calculate velocity field right at start-up using initial density distribution + if coaxialflowtype == 1: # CXA + velocityfluiddownmidpoints = m / A_flow_annulus / densityfluiddownmidpoints + velocityfluidupmidpoints = m / A_flow_centerpipe / densityfluidupmidpoints + velocityfluiddownnodes = m / A_flow_annulus / densityfluiddownnodes + velocityfluidupnodes = m / A_flow_centerpipe / densityfluidupnodes + elif coaxialflowtype == 2: # CXC + velocityfluiddownmidpoints = m / A_flow_centerpipe / densityfluiddownmidpoints + velocityfluidupmidpoints = m / A_flow_annulus / densityfluidupmidpoints + velocityfluiddownnodes = m / A_flow_centerpipe / densityfluiddownnodes + velocityfluidupnodes = m / A_flow_annulus / densityfluidupnodes + + # Obtain initial viscosity distribution of fluid [Pa·s] + viscosityfluiddownmidpoints = CP.PropsSI('V', 'P', Pfluiddownmidpoints, 'T', Tfluiddownmidpoints + 273.15, + 'Water') + viscosityfluidupmidpoints = CP.PropsSI('V', 'P', Pfluidupmidpoints, 'T', Tfluidupmidpoints + 273.15, 'Water') + + # Obtain initial specific heat capacity distribution of fluid [J/kg·K] + heatcapacityfluiddownmidpoints = CP.PropsSI('C', 'P', Pfluiddownmidpoints, 'T', Tfluiddownmidpoints + 273.15, + 'Water') + heatcapacityfluidupmidpoints = CP.PropsSI('C', 'P', Pfluidupmidpoints, 'T', Tfluidupmidpoints + 273.15, 'Water') + + # Obtain initial thermal conductivity distribution of fluid [W/m·K] + thermalconductivityfluiddownmidpoints = CP.PropsSI('L', 'P', Pfluiddownmidpoints, 'T', + Tfluiddownmidpoints + 273.15, 'Water') + thermalconductivityfluidupmidpoints = CP.PropsSI('L', 'P', Pfluidupmidpoints, 'T', Tfluidupmidpoints + 273.15, + 'Water') + + # Obtain initial thermal diffusivity distribution of fluid [m²/s] + alphafluiddownmidpoints = thermalconductivityfluiddownmidpoints / densityfluiddownmidpoints / heatcapacityfluiddownmidpoints + alphafluidupmidpoints = thermalconductivityfluidupmidpoints / densityfluidupmidpoints / heatcapacityfluidupmidpoints + + # Obtain initial thermal expansion coefficient distribution of fluid [1/K] + thermalexpansionfluiddownmidpoints = CP.PropsSI('ISOBARIC_EXPANSION_COEFFICIENT', 'P', Pfluiddownmidpoints, 'T', + Tfluiddownmidpoints + 273.15, 'Water') + thermalexpansionfluidupmidpoints = CP.PropsSI('ISOBARIC_EXPANSION_COEFFICIENT', 'P', Pfluidupmidpoints, 'T', + Tfluidupmidpoints + 273.15, 'Water') + + # Obtain initial Prandtl number distribution of fluid [-] + Prandtlfluiddownmidpoints = viscosityfluiddownmidpoints / densityfluiddownmidpoints / alphafluiddownmidpoints + Prandtlfluidupmidpoints = viscosityfluidupmidpoints / densityfluidupmidpoints / alphafluidupmidpoints + + # Obtain initial Reynolds number distribution of fluid [-] + if coaxialflowtype == 1: # CXA (injection in annulus) + Refluiddownmidpoints = densityfluiddownmidpoints * velocityfluiddownmidpoints * Dh_annulus / viscosityfluiddownmidpoints + Refluidupmidpoints = densityfluidupmidpoints * velocityfluidupmidpoints * 2 * radiuscenterpipe / viscosityfluidupmidpoints + elif coaxialflowtype == 2: # CXC (injection in center pipe) + Refluiddownmidpoints = densityfluiddownmidpoints * velocityfluiddownmidpoints * 2 * radiuscenterpipe / viscosityfluiddownmidpoints + Refluidupmidpoints = densityfluidupmidpoints * velocityfluidupmidpoints * Dh_annulus / viscosityfluidupmidpoints + + # Initialize SBT algorithm linear system of equation matrices + L = np.zeros((4 * N, 4 * N)) + R = np.zeros((4 * N, 1)) + Q = np.zeros((N, len(times))) + self.Tresoutput.value = np.zeros(len(times)) + Poutput = np.zeros(len(times)) + self.Tresoutput.value[0] = Tsurf + Poutput[0] = Pin * 1e5 + + Tfluidupnodesstore = np.zeros((N + 1, len(times))) + Tfluiddownnodesstore = np.zeros((N + 1, len(times))) + # MIR Tfluidupmidpointsstore = np.zeros((N, len(times))) + Tfluidupmidpointsstore = np.zeros((N + 1, len(times))) + # MIR Tfluiddownmidpointsstore = np.zeros((N, len(times))) + Tfluiddownmidpointsstore = np.zeros((N + 1, len(times))) + Pfluidupnodesstore = np.zeros((N + 1, len(times))) + Pfluiddownnodesstore = np.zeros((N + 1, len(times))) + # MIR Pfluidupmidpointsstore = np.zeros((N, len(times))) + # MIR Pfluiddownmidpointsstore = np.zeros((N, len(times))) + Pfluidupmidpointsstore = np.zeros((N + 1, len(times))) + Pfluiddownmidpointsstore = np.zeros((N + 1, len(times))) + Qfluidupnodesstore = np.zeros((N + 1, len(times))) + Qfluiddownnodesstore = np.zeros((N + 1, len(times))) + Phasefluidupnodesstore = np.zeros((N + 1, len(times))) + Phasefluiddownnodesstore = np.zeros((N + 1, len(times))) + Hfluidupnodesstore = np.zeros((N + 1, len(times))) + Hfluiddownnodesstore = np.zeros((N + 1, len(times))) + Qinterexchangestore = np.zeros((N, len(times))) + QinterexchangeUp = np.zeros(N) + + # Store initial values + Tfluidupnodesstore[:, 0] = Tfluidupnodes + Tfluiddownnodesstore[:, 0] = Tfluiddownnodes + Tfluidupmidpointsstore[:, 0] = BBinitial + Tfluiddownmidpointsstore[:, 0] = BBinitial + Pfluidupnodesstore[:, 0] = Pfluidupnodes + Pfluiddownnodesstore[:, 0] = Pfluiddownnodes + Pfluidupmidpointsstore[:, 0] = Pfluidupmidpoints + Pfluiddownmidpointsstore[:, 0] = Pfluiddownmidpoints + Qfluiddownnodesstore[:, 0] = Qfluiddownnodes + Qfluidupnodesstore[:, 0] = Qfluidupnodes + Phasefluidupnodesstore[:, 0] = Phasefluidupnodes + Phasefluiddownnodesstore[:, 0] = Phasefluiddownnodes + Qinterexchangestore[:, 0] = np.zeros(N) + + print('Pre-processing completed successfully. Starting simulation ...') + + # 3. Calculating + import time + start_time = time.time() + + for i in range(1, len(times)): + Deltat = times[i] - times[i - 1] + + if k_center_pipe > 1: + if times[i] < 1000: + k_center_pipe_corr = 1 + elif times[i] < 10000: + k_center_pipe_corr = max(1, 0.25 * k_center_pipe) + else: + k_center_pipe_corr = k_center_pipe + else: + k_center_pipe_corr = k_center_pipe + + if alpha_m * Deltat / radius ** 2 > LimitCylinderModelRequired: + CPCP = np.full(N, 1 / (4 * np.pi * self.krock.value) * expi(radius ** 2 / (4 * alpha_m * Deltat))) + else: + # MIR CPCP = np.interp(alpha_m * Deltat / radius ** 2, argumentbesselvec, besselcylinderresult) + # Calculate the interpolated values + interp_cylinder = interp1d(argumentbesselvec, besselcylinderresult, kind='linear', + fill_value="extrapolate") + interp_cylinder_boiler = interp1d(argumentbesselvec_boiler, besselcylinderresult_boiler, kind='linear', + fill_value="extrapolate") + + # Interpolation results + interp_value_cylinder = interp_cylinder(alpha_m * Deltat / radius ** 2) + interp_value_cylinder_boiler = interp_cylinder_boiler(alpha_m_boiler * Deltat / radius ** 2) + + # Create arrays filled with the interpolated values + ones_cylinder = np.ones(Nreg) * interp_value_cylinder + ones_cylinder_boiler = np.ones(Nboiler) * interp_value_cylinder_boiler + + # Combine the two arrays + CPCP = np.concatenate((ones_cylinder, ones_cylinder_boiler)) + + if Deltat > timeforfinitelinesource: + CPCP += np.interp(Deltaz ** 2 / (4 * alpha_m * Deltat), Amin1vector, finitecorrectiony) + + if i > 1: + CPOP = np.zeros((N, i - 1)) + for j in range(i - 1): + t1 = times[i] - times[j] + t2 = times[i] - times[j + 1] + CPOP[:, j] = Deltaz / (4 * np.pi * np.sqrt(alpha_m * np.pi) * self.krock.value) * ( + 1 / np.sqrt(t1) - 1 / np.sqrt(t2)) + + indexpsstart = 1 + indexpsend = np.where(timeforpointssource < (times[i] - times[:i]))[0] + if len(indexpsend) == 0: + indexpsend = indexpsstart - 1 + else: + indexpsend = indexpsend[-1] - 1 + + indexlsstart = indexpsend + 1 + indexlsend = np.where(timeforlinesource < (times[i] - times[:i]))[0] + if len(indexlsend) == 0: + indexlsend = indexlsstart - 1 + else: + indexlsend = indexlsend[-1] - 1 + + indexcsstart = max(indexpsend, indexlsend) + 1 + indexcsend = i - 2 + + if indexcsstart <= indexcsend: + CPOP[:, indexcsstart:indexcsend] = np.interp( + alpha_m * (times[i] - times[indexcsstart:indexcsend]) / radius ** 2, argumentbesselvec, + besselcylinderresult) + + indexflsstart = indexpsend + 1 + indexflsend = np.where(timeforfinitelinesource < (times[i] - times[:i]))[0] + if len(indexflsend) == 0: + indexflsend = indexflsstart - 1 + else: + indexflsend = indexflsend[-1] - 1 + + if indexflsend >= indexflsstart: + CPOP[:, indexflsstart:indexflsend] += np.interp( + Deltaz ** 2 / (4 * alpha_m * (times[i] - times[indexflsstart:indexflsend])), Amin1vector, + finitecorrectiony) + + NPCP = np.zeros((N, N)) + NPCP += np.diag(CPCP) + + spacingtest = alpha_m * Deltat / SMatrixSorted[:, 1:] ** 2 / LimitNPSpacingTime + maxspacingtest = np.max(spacingtest, axis=1) + if maxspacingtest[0] < 1: + maxindextoconsider = 0 + else: + maxindextoconsider = np.where(maxspacingtest > 1)[0][-1] + + if mindexNPCP < maxindextoconsider + 1: + indicestocalculate = SortedIndices[:, mindexNPCP + 1:maxindextoconsider + 1] + NPCP[range(N), indicestocalculate] = Deltaz[indicestocalculate] / ( + 4 * np.pi * k_m_vector[indicestocalculate] * SMatrix[range(N), indicestocalculate]) * erfc( + SMatrix[range(N), indicestocalculate] / np.sqrt(4 * alpha_m_vector[indicestocalculate] * Deltat)) + + BB = np.zeros(N) + if i > 1 and lastneighbourtoconsider[i] > 0: + SMatrixRelevant = SMatrixSorted[:, 1:lastneighbourtoconsider[i] + 1] + SoverLRelevant = SoverLSorted[:, 1:lastneighbourtoconsider[i] + 1] + SortedIndicesRelevant = SortedIndices[:, 1:lastneighbourtoconsider[i] + 1] + + maxtimeindexmatrix = alpha_m * (times[i] - times[1:i]) / (SMatrixRelevant.flatten()[:, None] ** 2) + allindices = np.arange(N * lastneighbourtoconsider[i] * (i - 1)) + pipeheatcomesfrom = SortedIndicesRelevant.flatten()[:, None] + pipeheatgoesto = np.tile(np.arange(N), (lastneighbourtoconsider[i], 1)).flatten()[:, None] + indicestoneglect = np.where(maxtimeindexmatrix.flatten() < LimitNPSpacingTime)[0] + maxtimeindexmatrix = np.delete(maxtimeindexmatrix, indicestoneglect, axis=0) + allindices = np.delete(allindices, indicestoneglect) + + indicesFoSlargerthan = np.where(maxtimeindexmatrix.flatten() > 10)[0] + indicestotakeforpsFoS = allindices[indicesFoSlargerthan] + + allindices2 = np.delete(allindices, indicesFoSlargerthan) + SoverLinearized = SoverLRelevant.flatten()[:, None] + indicestotakeforpsSoverL = np.where(SoverLinearized.flatten() > LimitSoverL)[0] + overallindicestotake = np.unique(np.concatenate((indicestotakeforpsSoverL, indicestotakeforpsFoS))) + + npipesheatsource = pipeheatcomesfrom[overallindicestotake] + npipesheatreceiving = pipeheatgoesto[overallindicestotake] + BB += np.bincount(npipesheatreceiving.flatten(), + weights=Q[pipeheatcomesfrom[overallindicestotake], :-1].flatten() * Deltaz[ + npipesheatsource.flatten()] / ( + 4 * np.pi * k_m_vector[npipesheatsource.flatten()] * SMatrix[ + pipeheatgoesto[overallindicestotake], pipeheatcomesfrom[ + overallindicestotake]].flatten()) * erfc(SMatrix[pipeheatgoesto[ + overallindicestotake], pipeheatcomesfrom[ + overallindicestotake]].flatten() / np.sqrt( + 4 * alpha_m_vector[npipesheatsource.flatten()] * (times[i] - times[:-1]))), + minlength=N) + + if i > 1: + maxindextocalculate = lastneighbourtoconsider[i] + maxindex = lastneighbourtoconsider[i] * (i - 1) + maxtimeindexmatrix = alpha_m * (times[i] - times[1:i]) / ( + SMatrixSorted[:, 1:maxindextocalculate + 1].flatten()[:, None] ** 2) + allindices = np.arange(N * lastneighbourtoconsider[i] * (i - 1)) + allindices = np.delete(allindices, np.where(maxtimeindexmatrix.flatten() < LimitNPSpacingTime)[0]) + pipeheatcomesfrom = SortedIndices[:, 1:maxindextocalculate + 1].flatten()[:, None] + pipeheatgoesto = np.tile(np.arange(N), (lastneighbourtoconsider[i], 1)).flatten()[:, None] + NPCP[pipeheatgoesto.flatten()[allindices], pipeheatcomesfrom.flatten()[allindices]] += Deltaz[ + pipeheatcomesfrom.flatten()[ + allindices]] / ( + 4 * np.pi * + k_m_vector[ + pipeheatcomesfrom.flatten()[ + allindices]] * + SMatrix[ + pipeheatgoesto.flatten()[ + allindices], + pipeheatcomesfrom.flatten()[ + allindices]]) * erfc( + SMatrix[pipeheatgoesto.flatten()[allindices], pipeheatcomesfrom.flatten()[allindices]] / np.sqrt( + 4 * alpha_m_vector[pipeheatcomesfrom.flatten()[allindices]] * Deltat)) + + Vvector = 2 * np.pi * radius ** 2 * self.rhorock.value * self.cprock.value + + L[0:N, 0:N] = np.diag(Vvector / Deltat + np.sum(NPCP, axis=1)) + L[0:N, 3 * N:4 * N] = -1 * np.eye(N) + + L[N:2 * N, N:2 * N] = np.diag(Vvector / Deltat + np.sum(NPCP, axis=1)) + L[N:2 * N, 3 * N:4 * N] = -1 * np.eye(N) + + L[2 * N:3 * N, 0:N] = -1 * np.eye(N) + L[2 * N:3 * N, 3 * N:4 * N] = np.diag(velocityfluidupmidpoints) + L[2 * N:3 * N, 3 * N:4 * N] -= np.diag(velocityfluiddownmidpoints) + L[2 * N:3 * N, 2 * N:3 * N] = np.diag(densityfluidupmidpoints) - np.diag(densityfluiddownmidpoints) + + L[3 * N:4 * N, 2 * N:3 * N] = np.diag(heatcapacityfluidupmidpoints) + L[3 * N:4 * N, 0:N] = np.diag(densityfluidupmidpoints * heatcapacityfluidupmidpoints) + L[3 * N:4 * N, 3 * N:4 * N] = np.diag(np.pi * ( + 2 * radiuscenterpipe) ** 2 * thermalconductivityfluidupmidpoints * velocityfluidupmidpoints * Deltat) + + Q[:, i] = BB + + R[0:N] = Vvector * Tfluiddownmidpointsstore[:, i - 1] / Deltat + R[0:N] += Q[:, i] + R[N:2 * N] = Vvector * Tfluidupmidpointsstore[:, i - 1] / Deltat + R[N:2 * N] += Q[:, i] + + R[2 * N:3 * N] = Pfluidupmidpointsstore[:, i - 1] - Pfluiddownmidpointsstore[:, + i - 1] + velocityfluiddownmidpointsstore[:, + i - 1] * Pfluiddownmidpointsstore[:, + i - 1] * Deltat + + R[3 * N:4 * N] = heatcapacityfluidupmidpointsstore[:, i - 1] * Tfluidupmidpointsstore[:, i - 1] + + try: + solutions = np.linalg.solve(L, R) + except np.linalg.LinAlgError: + print(f'Simulation terminated prematurely due to linear algebra error at time step {i}.') + break + + BB = solutions[0:N] + Tfluidupmidpointsstore[:, i] = BB + Tfluidupmidpoints = BB + Tfluidupmidpoints = Tfluidupmidpoints + + BB = solutions[N:2 * N] + Tfluiddownmidpointsstore[:, i] = BB + Tfluiddownmidpoints = BB + Tfluiddownmidpoints = Tfluiddownmidpoints + + BB = solutions[2 * N:3 * N] + Pfluidupmidpointsstore[:, i] = BB + Pfluidupmidpoints = BB + Pfluidupmidpoints = Pfluidupmidpoints + + BB = solutions[3 * N:4 * N] + Pfluiddownmidpointsstore[:, i] = BB + Pfluiddownmidpoints = BB + Pfluiddownmidpoints = Pfluiddownmidpoints + + Pfluidupmidpoints = Pfluidupmidpointsstore[:, i - 1] + velocityfluidupmidpoints * Deltat + Pfluiddownmidpoints = Pfluiddownmidpointsstore[:, i - 1] + velocityfluiddownmidpoints * Deltat + + Pfluidupnodes = np.append(Pfluidupnodesstore[0, i], np.diff(Pfluidupmidpoints)) + Pfluiddownnodes = np.append(Pfluiddownnodesstore[0, i], np.diff(Pfluiddownmidpoints)) + + Pfluidupnodesstore[:, i] = Pfluidupnodes + Pfluiddownnodesstore[:, i] = Pfluiddownnodes + + Qinterexchangestore[:, i] = QinterexchangeUp + + if Tfluidupmidpointsstore[-1, i] >= 100: + print(f'Warning: fluid temperature exceeds boiling point at time step {i}. Simulation may be invalid.') + + if Tfluiddownmidpointsstore[-1, i] >= 100: + print(f'Warning: fluid temperature exceeds boiling point at time step {i}. Simulation may be invalid.') + + end_time = time.time() + print(f'Simulation completed successfully in {end_time - start_time:.2f} seconds.') + + # plt.figure() + # plt.plot(times / (24 * 3600), Tfluiddownmidpointsstore[-1, :], label='Tfluiddownmidpoints') + # plt.plot(times / (24 * 3600), Tfluidupmidpointsstore[-1, :], label='Tfluidupmidpoints') + # plt.xlabel('Time (days)') + # plt.ylabel('Temperature (°C)') + # lt.legend() + # plt.title('Temperature vs Time') + # plt.grid() + # plt.show() + + model.logger.info(f'complete {str(__class__)}: {sys._getframe().f_code.co_name}') + + @lru_cache(maxsize=1024) + #@profile + def Calculate_Uloop(self, model): + """ + Calculate the U-loop version of the SBT model + """ + model.logger.info(f'Init {str(__class__)}: {sys._getframe().f_code.co_name}') + self.averagegradient.value = np.average(self.gradient.value[0:self.numseg.value]) + self.Trock.value = self.Tsurf.value + (self.averagegradient.value * model.wellbores.lateral_endpoint_depth.value) + + lateralflowallocation = [] + for i in range(model.wellbores.numnonverticalsections.value): + lateralflowallocation.append(1 / model.wellbores.numnonverticalsections.value) + + # interpolate time steps - but consider ignoring the first step. + # simulation times [s] (must start with 0; to obtain smooth results, + # abrupt changes in time step size should be avoided. logarithmic spacing is recommended) + initial_times = np.linspace(0, self.initial_final_timestep_transition.value, self.initial_timestep_count.value) + initial_time_interval = initial_times[1] - initial_times[0] + final_start = self.initial_final_timestep_transition.value + initial_time_interval + final_times = np.logspace(np.log10(final_start), np.log10(model.surfaceplant.plant_lifetime.value * 365 * 24 * 3600), int(self.final_timestep_count.value)) + times = np.concatenate([initial_times, final_times]) + # Note 1: When providing a variable injection temperature or flow rate, a finer time grid should be considered. + # Below is one with long term time steps of about 36 days. + # times = [0] + list(range(100, 10000, 100)) + list(np.logspace(np.log10(100*100), np.log10(0.1*365*24*3600), 40)) + list(np.arange(0.2*365*24*3600, 20*365*24*3600, 0.1*365*24*3600)) + # Note 2: To capture the start-up effects, several small time steps need to be taken during + # the first 99000 seconds in the time vector considered. + # To speed up the simulation, this can be avoided with limited impact on the long-term results. + # For example, an alternative time vector would be: + # times = [0] + list(range(100, 1000, 100)) + list(range(1000, 10000, 1000)) + list(np.logspace(np.log10(100*100), np.log10(20*365*24*3600), 75)) + + # (x,y,z)-coordinates of centerline of injection well, production well and laterals + # The vectors storing the x-, y- and z-coordinates should be column vectors + # To obtain smooth results, abrupt changes in segment lengths should be avoided. + # Coordinates of injection well (coordinates are provided from top to bottom in the direction of flow) + xinj, yinj, zinj, xprod, yprod, zprod, xlat, ylat, zlat = generate_wireframe_model(model.wellbores.lateral_endpoint_depth.value, + model.wellbores.numnonverticalsections.value, + model.wellbores.lateral_spacing.value, + model.wellbores.element_length.value, + model.wellbores.junction_depth.value, + model.wellbores.vertical_section_length.value, + model.wellbores.lateral_inclination_angle.value * np.pi / 180.0, + model.wellbores.vertical_wellbore_spacing.value, + self.generate_wireframe_graphics.value) + + # Merge x-, y-, and z-coordinates + x = np.concatenate((xinj, xprod)) + y = np.concatenate((yinj, yprod)) + z = np.concatenate((zinj, zprod)) + + for i in range(model.wellbores.numnonverticalsections.value): + x = np.concatenate((x, xlat[:, i].reshape(-1, 1).flatten())) + y = np.concatenate((y, ylat[:, i].reshape(-1, 1).flatten())) + z = np.concatenate((z, zlat[:, i].reshape(-1, 1).flatten())) + + gamma = 0.577215665 # Euler's constant + alpha_f = model.surfaceplant.k_fluid.value / model.surfaceplant.rho_fluid.value / model.surfaceplant.cp_fluid.value # Fluid thermal diffusivity [m2/s] + Pr_f = model.surfaceplant.mu_fluid.value / model.surfaceplant.rho_fluid.value / alpha_f # Fluid Prandtl number [-] + alpha_m = self.krock.value / self.rhorock.value / self.cprock.value # Thermal diffusivity medium [m2/s] + + # interconnections lists the indices of interconnections between inj, prod, and laterals + # (this will be used to take care of the duplicate coordinates of the start and end points of the laterals) + interconnections = np.concatenate((np.array([len(xinj)],dtype=int), np.array([len(xprod)],dtype=int), (np.ones(model.wellbores.numnonverticalsections.value - 1, dtype=int) * len(xlat)))) + interconnections = np.cumsum(interconnections) + + # radiusvector stores radius of each element in a vector [m] + radiusvector = np.concatenate([np.ones(len(xinj) + len(xprod) - 2) * (model.wellbores.prodwelldiam.quantity().to('m').magnitude / 2), np.ones(model.wellbores.numnonverticalsections.value * len(xlat) - model.wellbores.numnonverticalsections.value) * (model.wellbores.nonverticalwellborediameter.quantity().to('m').magnitude / 2.0)]) + Dvector = radiusvector * 2 # Diameter of each element [m] + lateralflowallocation = lateralflowallocation / np.sum(lateralflowallocation) # Ensure the sum equals 1 + + Deltaz = np.sqrt((x[1:] - x[:-1]) ** 2 + (y[1:] - y[:-1]) ** 2 + (z[1:] - z[:-1]) ** 2) # Length of each segment [m] + Deltaz = np.delete(Deltaz, interconnections - 1) # Removes the phantom elements due to duplicate coordinates + TotalLength = np.sum(Deltaz) # Total length of all elements (for informational purposes only) [m] + + # Quality Control + LoverR = Deltaz / radiusvector # Ratio of pipe segment length to radius along the wellbore [-] + smallestLoverR = np.min(LoverR) # Smallest ratio of pipe segment length to pipe radius. This ratio should be larger than 10. [-] + + if smallestLoverR < 10: + msg = 'Warning: smallest ratio of segment length over radius is less than 10. Good practice is to keep this ratio larger than 10.' + print(f'{msg}') + model.logger.warning(msg) + + if model.wellbores.numnonverticalsections.value > 1: + DeltazOrdered = np.concatenate((Deltaz[0:(interconnections[0]-1)], Deltaz[(interconnections[1]-2):(interconnections[2]-3)], Deltaz[(interconnections[0]-1):(interconnections[1]-2)])) + else: + DeltazOrdered = np.concatenate((Deltaz[0:interconnections[0] - 1], Deltaz[interconnections[1] - 1:-1], Deltaz[interconnections[0]:interconnections[1] - 2])) + + RelativeLengthChanges = (DeltazOrdered[1:] - DeltazOrdered[:-1]) / DeltazOrdered[:-1] + + if max(abs(RelativeLengthChanges)) > 0.6: + msg = 'Warning: abrupt change(s) in segment length detected, which may cause numerical instabilities. Good practice is to avoid abrupt length changes to obtain smooth results.' + print(f'{msg}') + model.logger.warning(msg) + + for dd in range(1, model.wellbores.numnonverticalsections.value + 1): + if abs(xinj[-1] - xlat[0][dd - 1]) > 1e-12 or abs(yinj[-1] - ylat[0][dd - 1]) > 1e-12 or abs(zinj[-1] - zlat[0][dd - 1]) > 1e-12: + msg = f'Error: Coordinate mismatch between bottom of injection well and start of lateral #{dd}' + print(f'{msg}') + model.logger.fatal(msg) + model.logger.info(f'Complete {str(__name__)}: {sys._getframe().f_code.co_name}') + raise ValueError(msg) + + if abs(xprod[0] - xlat[-1][dd - 1]) > 1e-12 or abs(yprod[0] - ylat[-1][dd - 1]) > 1e-12 or abs(zprod[0] - zlat[-1][dd - 1]) > 1e-12: + msg = f'Error: Coordinate mismatch between bottom of production well and end of lateral #{dd}' + print(f'{msg}') + model.logger.fatal(msg) + model.logger.info(f'Complete {str(__name__)}: {sys._getframe().f_code.co_name}') + raise ValueError(msg) + + if len(lateralflowallocation) != model.wellbores.numnonverticalsections.value: + msg = 'Error: Length of array "lateralflowallocation" does not match the number of laterals' + print(f'{msg}') + model.logger.fatal(msg) + model.logger.info(f'Complete {str(__name__)}: {sys._getframe().f_code.co_name}') + raise ValueError(msg) + + # Read injection temperature profile if provided + Tinstore = [0] * len(times) + if self.injection_temperature_model == 2: + # User has provided injection temperature in an Excel spreadsheet. + num = pd.read_excel(self.injection_temperature_file.value) + Tintimearray = np.array(num.iloc[:, 0]) + Tintemperaturearray = np.array(num.iloc[:, 1]) + # Quality control + if len(Tintimearray) < 2: + msg = 'Error: Provided injection temperature profile should have at least 2 values' + print(f'{msg}') + model.logger.fatal(msg) + model.logger.info(f'Complete {str(__name__)}: {sys._getframe().f_code.co_name}') + raise ValueError(msg) + + if Tintimearray[0] != 0: + msg = 'Error: First time value in the user-provided injection temperature profile does not equal 0 s' + print(f'{msg}') + model.logger.fatal(msg) + model.logger.info(f'Complete {str(__name__)}: {sys._getframe().f_code.co_name}') + raise ValueError(msg) + + if abs(Tintimearray[-1] - times[-1]) > 10e-6: + msg = 'Error: Last time value in the user-provided injection temperature profile does not equal the final value in the "times" array' + print(f'{msg}') + model.logger.fatal(msg) + model.logger.info(f'Complete {str(__name__)}: {sys._getframe().f_code.co_name}') + raise ValueError(msg) + + else: + # Ensure final time values "exactly" match to prevent interpolation issues at the final time step + Tintimearray[-1] = times[-1] + Tinstore[0] = Tintemperaturearray[0] + else: + Tinstore[0] = model.wellbores.Tinj.value + + # The value for m used at each time step is stored in this array (is either constant or interpolated from a user-provided mass flow rate profile) + mstore = np.zeros(len(times)) + + # User has provided mass flow rate in an Excel spreadsheet. + if self.flow_rate_model.value == 2: + data = pd.read_excel(self.flow_rate_file.value) + mtimearray = data.iloc[:, 0].values # This array has the times provided by the user + mflowratearray = data.iloc[:, 1].values # This array has the injection temperatures provided by the user + + # Quality control + if len(mtimearray) < 2: + msg = 'Error: Provided flow rate profile should have at least 2 values' + print(f'{msg}') + model.logger.fatal(msg) + model.logger.info(f'Complete {str(__name__)}: {sys._getframe().f_code.co_name}') + raise ValueError(msg) + + if mtimearray[0] != 0: + msg = 'Error: First time value in user-provided flow rate profile does not equal to 0 s' + print(f'{msg}') + model.logger.fatal(msg) + model.logger.info(f'Complete {str(__name__)}: {sys._getframe().f_code.co_name}') + raise ValueError(msg) + + if abs(mtimearray[-1] - times[-1]) > 10e-6: + msg = 'Error: Last time value in user-provided flow rate profile does not equal to final value in "times" array' + print(f'{msg}') + model.logger.fatal(msg) + model.logger.info(f'Complete {str(__name__)}: {sys._getframe().f_code.co_name}') + raise ValueError(msg) + + else: + # Ensure final time values "exactly" match to prevent interpolation issues at the final time step + mtimearray[-1] = times[-1] + + mstore[0] = mflowratearray[0] + else: + mstore[0] = model.wellbores.prodwellflowrate.value + + # assume default values + NoArgumentsFinitePipeCorrection = 25 + NoDiscrFinitePipeCorrection = 200 + NoArgumentsInfCylIntegration = 25 + NoDiscrInfCylIntegration = 200 + LimitPointSourceModel = 1.5 + LimitCylinderModelRequired = 25 + LimitInfiniteModel = 0.05 + LimitNPSpacingTime = 0.1 + LimitSoverL = 1.5 + M = 3 + if self.SBTAccuracyDesired.value == 2: + NoArgumentsFinitePipeCorrection = 50 + NoDiscrFinitePipeCorrection = 400 + NoArgumentsInfCylIntegration = 50 + NoDiscrInfCylIntegration = 400 + LimitPointSourceModel = 2.5 + LimitCylinderModelRequired = 50 + LimitInfiniteModel = 0.01 + LimitNPSpacingTime = 0.04 + LimitSoverL = 2 + M = 4 + elif self.SBTAccuracyDesired.value == 3: + NoArgumentsFinitePipeCorrection = 100 + NoDiscrFinitePipeCorrection = 500 + NoArgumentsInfCylIntegration = 100 + NoDiscrInfCylIntegration = 500 + LimitPointSourceModel = 5 + LimitCylinderModelRequired = 100 + LimitInfiniteModel = 0.004 + LimitNPSpacingTime = 0.02 + LimitSoverL = 3 + M = 5 + elif self.SBTAccuracyDesired.value == 4: + NoArgumentsFinitePipeCorrection = 200 + NoDiscrFinitePipeCorrection = 1000 + NoArgumentsInfCylIntegration = 200 + NoDiscrInfCylIntegration = 1000 + LimitPointSourceModel = 10 + LimitCylinderModelRequired = 200 + LimitInfiniteModel = 0.002 + LimitNPSpacingTime = 0.01 + LimitSoverL = 5 + M = 10 + elif self.SBTAccuracyDesired.value == 5: + NoArgumentsFinitePipeCorrection = 400 + NoDiscrFinitePipeCorrection = 2000 + NoArgumentsInfCylIntegration = 400 + NoDiscrInfCylIntegration = 2000 + LimitPointSourceModel = 20 + LimitCylinderModelRequired = 400 + LimitInfiniteModel = 0.001 + LimitNPSpacingTime = 0.005 + LimitSoverL = 9 + M = 20 + + timeforpointssource = max(Deltaz)**2 / alpha_m * LimitPointSourceModel # Calculates minimum time step size when point source model becomes applicable [s] + timeforlinesource = max(radiusvector)**2 / alpha_m * LimitCylinderModelRequired # Calculates minimum time step size when line source model becomes applicable [s] + timeforfinitelinesource = max(Deltaz)**2 / alpha_m * LimitInfiniteModel # Calculates minimum time step size when finite line source model should be considered [s] + + fpcminarg = min(Deltaz)**2 / (4 * alpha_m * times[-1]) + fpcmaxarg = max(Deltaz)**2 / (4 * alpha_m * (min(times[1:] - times[:-1]))) + Amin1vector = np.logspace(np.log10(fpcminarg) - 0.1, np.log10(fpcmaxarg) + 0.1, NoArgumentsFinitePipeCorrection) + finitecorrectiony = np.zeros(NoArgumentsFinitePipeCorrection) + + for i, Amin1 in enumerate(Amin1vector): + Amax1 = (16)**2 + if Amin1 > Amax1: + Amax1 = 10 * Amin1 + Adomain1 = np.logspace(np.log10(Amin1), np.log10(Amax1), NoDiscrFinitePipeCorrection) + finitecorrectiony[i] = np.trapz(-1 / (Adomain1 * 4 * np.pi * self.krock.value) * erfc(1/2 * np.power(Adomain1, 1/2)), Adomain1) + + besselminarg = alpha_m * (min(times[1:] - times[:-1])) / max(radiusvector)**2 + besselmaxarg = alpha_m * timeforlinesource / min(radiusvector)**2 + deltazbessel = np.logspace(-10, 8, NoDiscrInfCylIntegration) + argumentbesselvec = np.logspace(np.log10(besselminarg) - 0.5, np.log10(besselmaxarg) + 0.5, NoArgumentsInfCylIntegration) + besselcylinderresult = np.zeros(NoArgumentsInfCylIntegration) + + for i, argumentbessel in enumerate(argumentbesselvec): + besselcylinderresult[i] = 2 / (self.krock.value * np.pi**3) * np.trapz((1 - np.exp(-deltazbessel**2 * argumentbessel)) / (deltazbessel**3 * (jv(1, deltazbessel)**2 + yv(1, deltazbessel)**2)), deltazbessel) + + N = len(Deltaz) # Number of elements + elementcenters = 0.5 * np.column_stack((x[1:], y[1:], z[1:])) + 0.5 * np.column_stack((x[:-1], y[:-1], z[:-1])) # Matrix that stores the mid point coordinates of each element + interconnections = interconnections - 1 + elementcenters = np.delete(elementcenters, interconnections.reshape(-1,1), axis=0) # Remove duplicate coordinates + SMatrix = np.zeros((N, N)) # Initializes the spacing matrix, which holds the distance between center points of each element [m] + + for i in range(N): + SMatrix[i, :] = np.sqrt((elementcenters[i, 0] - elementcenters[:, 0])**2 + (elementcenters[i, 1] - elementcenters[:, 1])**2 + (elementcenters[i, 2] - elementcenters[:, 2])**2) + + SoverL = np.zeros((N, N)) # Initializes the ratio of spacing to element length matrix + + for i in range(N): + SMatrix[i, :] = np.sqrt((elementcenters[i, 0] - elementcenters[:, 0])**2 + (elementcenters[i, 1] - elementcenters[:, 1])**2 + (elementcenters[i, 2] - elementcenters[:, 2])**2) + SoverL[i, :] = SMatrix[i, :] / Deltaz[i] + + SortedIndices = np.argsort(SMatrix, axis=1, kind = 'stable') # Getting the indices of the sorted elements + SMatrixSorted = np.take_along_axis(SMatrix, SortedIndices, axis=1) # Sorting the spacing matrix + + SoverLSorted = SMatrixSorted / Deltaz + + mindexNPCP = np.where(np.min(SoverLSorted, axis=0) < LimitSoverL)[0][-1] # Finding the index where the ratio is less than the limit + + midpointsx = elementcenters[:, 0] # x-coordinate of center of each element [m] + midpointsy = elementcenters[:, 1] # y-coordinate of center of each element [m] + midpointsz = elementcenters[:, 2] # z-coordinate of center of each element [m] + BBinitial = self.Tsurf.value - (self.gradient1.value / 1000.0) * midpointsz # Initial temperature at center of each element [degC] + + previouswaterelements = np.zeros(N) + previouswaterelements[0:] = np.arange(-1,N-1) + + for i in range(model.wellbores.numnonverticalsections.value ): + previouswaterelements[interconnections[i + 1] - i-1] = len(xinj) - 2 + + previouswaterelements[len(xinj) - 1] = 0 + + lateralendpoints = [] + for i in range(1,model.wellbores.numnonverticalsections.value +1): + lateralendpoints.append(len(xinj) - 2 + len(xprod) - 1 + i * ((xlat[:, 0]).size- 1)) + lateralendpoints = np.array(lateralendpoints) + + MaxSMatrixSorted = np.max(SMatrixSorted, axis=0) + + indicesyoucanneglectupfront = alpha_m * (np.ones((N-1, 1)) * times) / (MaxSMatrixSorted[1:].reshape(-1, 1) * np.ones((1, len(times))))**2 / LimitNPSpacingTime + indicesyoucanneglectupfront[indicesyoucanneglectupfront > 1] = 1 + + lastneighbourtoconsider = np.zeros(len(times)) + for i in range(len(times)): + lntc = np.where(indicesyoucanneglectupfront[:, i] == 1)[0] + if len(lntc) == 0: + lastneighbourtoconsider[i] = 1 + else: + lastneighbourtoconsider[i] = max(2, lntc[-1] + 1) + + distributionx = np.zeros((len(x) - 1, M + 1)) + distributiony = np.zeros((len(x) - 1, M + 1)) + distributionz = np.zeros((len(x) - 1, M + 1)) + + for i in range(len(x) - 1): + distributionx[i, :] = np.linspace(x[i], x[i + 1], M + 1).reshape(-1) + distributiony[i, :] = np.linspace(y[i], y[i + 1], M + 1).reshape(-1) + distributionz[i, :] = np.linspace(z[i], z[i + 1], M + 1).reshape(-1) + + # Remove duplicates + distributionx = np.delete(distributionx, interconnections, axis=0) + distributiony = np.delete(distributiony, interconnections, axis=0) + distributionz = np.delete(distributionz, interconnections, axis=0) + + # Initialize SBT algorithm linear system of equation matrices + L = np.zeros((3 * N, 3 * N)) # Will store the "left-hand side" of the system of equations + R = np.zeros((3 * N, 1)) # Will store the "right-hand side" of the system of equations + Q = np.zeros((N, len(times))) # Initializes the heat pulse matrix, i.e., the heat pulse emitted by each element at each time step + self.Tresoutput.value = np.zeros(len(times)) # Initializes the production temperatures array + Twprevious = BBinitial # At time zero, the initial fluid temperature corresponds to the initial local rock temperature + self.Tresoutput.value[0] = self.Tsurf.value # At time zero, the outlet temperature is the initial local fluid temperature at the surface, which corresponds to the surface temperature + TwMatrix = np.zeros((len(times), N)) # Initializes the matrix that holds the fluid temperature over time + TwMatrix[0, :] = Twprevious + + count = 0 + for i in range(1, len(times)): + count = count + 1 + Deltat = times[i] - times[i - 1] # Current time step size [s] + + # If the user has provided an injection temperature profile, current value of model.wellbores.Tinj.value is calculated + if self.injection_temperature_model == 2: + model.wellbores.Tinj.value = np.interp(times[i], Tintimearray, Tintemperaturearray) + Tinstore[i] = model.wellbores.Tinj.value # Value that is used for model.wellbores.Tinj.value at each time step gets stored for postprocessing purposes + + # If the user has provided a flow rate profile, current value of model.wellbores.prodwellflowrate.value is calculated + if self.flow_rate_model.value == 2: + model.wellbores.prodwellflowrate.value = np.interp(times[i], mtimearray, mflowratearray) + mstore[i] = model.wellbores.prodwellflowrate.value # Value that is used for model.wellbores.prodwellflowrate.value at each time step gets stored for postprocessing purposes + + # Velocities and thermal resistances are calculated each time step as the flow rate is allowed to vary each time step + uvertical = model.wellbores.prodwellflowrate.value / model.surfaceplant.rho_fluid.value / (np.pi * (model.wellbores.prodwelldiam.quantity().to('m').magnitude / 2) ** 2) # Fluid velocity in vertical injector and producer [m/s] + ulateral = model.wellbores.prodwellflowrate.value / model.surfaceplant.rho_fluid.value / (np.pi * (model.wellbores.nonverticalwellborediameter.quantity().to('m').magnitude / 2.0) ** 2) * lateralflowallocation # Fluid velocity in each lateral [m/s] + uvector = np.hstack((uvertical * np.ones(len(xinj) + len(xprod) - 2))) + + for dd in range(model.wellbores.numnonverticalsections.value ): + uvector = np.hstack((uvector, ulateral[dd] * np.ones(len(xlat[:, 0]) - 1))) + + if model.wellbores.prodwellflowrate.value > 0.1: + Revertical = model.surfaceplant.rho_fluid.value * uvertical * (2 * (model.wellbores.prodwelldiam.quantity().to('m').magnitude / 2)) / model.surfaceplant.mu_fluid.value # Fluid Reynolds number in injector and producer [-] + Nuvertical = 0.023 * Revertical ** (4 / 5) * Pr_f ** 0.4 # Nusselt Number in injector and producer (we assume turbulent flow) [-] + else: + Nuvertical = 1 # At low flow rates, we assume we are simulating the condition of well shut-in and set the Nusselt number to 1 (i.e., conduction only) [-] + + hvertical = Nuvertical * model.surfaceplant.k_fluid.value / (2 * (model.wellbores.prodwelldiam.quantity().to('m').magnitude / 2)) # Heat transfer coefficient in injector and producer [W/m2/K] + Rtvertical = 1 / (np.pi * hvertical * 2 * (model.wellbores.prodwelldiam.quantity().to('m').magnitude / 2)) # Thermal resistance in injector and producer (open-hole assumed) + + if model.wellbores.prodwellflowrate.value > 0.1: + Relateral = model.surfaceplant.rho_fluid.value * ulateral * (2 * (model.wellbores.nonverticalwellborediameter.quantity().to('m').magnitude / 2.0)) / model.surfaceplant.mu_fluid.value # Fluid Reynolds number in lateral [-] + Nulateral = 0.023 * Relateral ** (4 / 5) * Pr_f ** 0.4 # Nusselt Number in lateral (we assume turbulent flow) [-] + else: + Nulateral = np.ones(model.wellbores.numnonverticalsections.value) # At low flow rates, we assume we are simulating the condition of well shut-in and set the Nusselt number to 1 (i.e., conduction only) [-] + + hlateral = Nulateral * model.surfaceplant.k_fluid.value / (2 * (model.wellbores.nonverticalwellborediameter.quantity().to('m').magnitude / 2.0)) # Heat transfer coefficient in lateral [W/m2/K] + Rtlateral = 1 / (np.pi * hlateral * 2 * (model.wellbores.nonverticalwellborediameter.quantity().to('m').magnitude / 2.0)) # Thermal resistance in lateral (open-hole assumed) + + Rtvector = Rtvertical * np.ones(len(radiusvector)) # Store thermal resistance of each element in a vector + + for dd in range(1, model.wellbores.numnonverticalsections.value + 1): + if dd < model.wellbores.numnonverticalsections.value: + Rtvector[interconnections[dd] - dd : interconnections[dd + 1] - dd] = Rtlateral[dd - 1] * np.ones(len(xlat[:, 0])) + else: + Rtvector[interconnections[dd] - model.wellbores.numnonverticalsections.value :] = Rtlateral[dd - 1] * np.ones(len(xlat[:, 0]) - 1) + + # CPCP (= Current pipe, current pulses) + if alpha_m * Deltat / max(radiusvector)**2 > LimitCylinderModelRequired: + CPCP = np.ones(N) * 1 / (4 * np.pi * self.krock.value) * exp1(radiusvector**2 / (4 * alpha_m * Deltat)) # Use line source model if possible + else: + CPCP = np.ones(N) * np.interp(alpha_m * Deltat / radiusvector**2, argumentbesselvec, besselcylinderresult) # Use cylindrical source model if required + + if Deltat > timeforfinitelinesource: # For long time steps, the finite length correction should be applied + CPCP = CPCP + np.interp(Deltaz**2 / (4 * alpha_m * Deltat), Amin1vector, finitecorrectiony) + + # CPOP (= Current pipe, old pulses) + if i > 1: # After the second time step, we need to keep track of previous heat pulses + CPOP = np.zeros((N, i-1)) + indexpsstart = 0 + indexpsend = np.where(timeforpointssource < (times[i] - times[1:i]))[-1] + if indexpsend.size > 0: + indexpsend = indexpsend[-1] + 1 + else: + indexpsend = indexpsstart - 1 + if indexpsend >= indexpsstart: # Use point source model if allowed + + CPOP[:, 0:indexpsend] = Deltaz * np.ones((N, indexpsend)) / (4 * np.pi * np.sqrt(alpha_m * np.pi) * self.krock.value) * ( + np.ones(N) * (1 / np.sqrt(times[i] - times[indexpsstart + 1:indexpsend + 2]) - + 1 / np.sqrt(times[i] - times[indexpsstart:indexpsend+1]))) + indexlsstart = indexpsend + 1 + indexlsend = np.where(timeforlinesource < (times[i] - times[1:i]))[0] + if indexlsend.size == 0: + indexlsend = indexlsstart - 1 + else: + indexlsend = indexlsend[-1] + + if indexlsend >= indexlsstart: # Use line source model for more recent heat pulse events + + CPOP[:, indexlsstart:indexlsend+1] = np.ones((N,1)) * 1 / (4*np.pi*self.krock.value) * (exp1((radiusvector**2).reshape(len(radiusvector ** 2),1) / (4*alpha_m*(times[i]-times[indexlsstart:indexlsend+1])).reshape(1,len(4 * alpha_m * (times[i] - times[indexlsstart:indexlsend+1]))))-\ + exp1((radiusvector**2).reshape(len(radiusvector ** 2),1) / (4 * alpha_m * (times[i]-times[indexlsstart+1:indexlsend+2])).reshape(1,len(4 * alpha_m * (times[i] - times[indexlsstart+1:indexlsend+2]))))) + #pdb.set_trace() + indexcsstart = max(indexpsend, indexlsend) + 1 + indexcsend = i - 2 + + if indexcsstart <= indexcsend: # Use cylindrical source model for the most recent heat pulses + + CPOPPH = np.zeros((CPOP[:, indexcsstart:indexcsend+1].shape)) + CPOPdim = CPOP[:, indexcsstart:indexcsend+1].shape + CPOPPH = CPOPPH.T.ravel() + CPOPPH = (np.ones(N) * ( + np.interp(alpha_m * (times[i] - times[indexcsstart:indexcsend+1]).reshape(len(times[i] - times[indexcsstart:indexcsend+1]),1) / (radiusvector ** 2).reshape(len(radiusvector ** 2),1).T, argumentbesselvec, besselcylinderresult) - \ + np.interp(alpha_m * (times[i] - times[indexcsstart+1: indexcsend+2]).reshape(len(times[i] - times[indexcsstart+1:indexcsend+2]),1) / (radiusvector ** 2).reshape(len(radiusvector ** 2),1).T, argumentbesselvec, besselcylinderresult))).reshape(-1,1) + CPOPPH=CPOPPH.reshape((CPOPdim),order='F') + CPOP[:, indexcsstart:indexcsend+1] = CPOPPH + indexflsstart = indexpsend + 1 + indexflsend = np.where(timeforfinitelinesource < (times[i] - times[1:i]))[-1] + if indexflsend.size == 0: + indexflsend = indexflsstart - 1 + else: + indexflsend = indexflsend[-1] - 1 + + if indexflsend >= indexflsstart: # Perform finite length correction if needed + CPOP[:, indexflsstart:indexflsend+2] = CPOP[:, indexflsstart:indexflsend+2] + (np.interp(np.matmul((Deltaz.reshape(len(Deltaz),1) ** 2),np.ones((1,indexflsend-indexflsstart+2))) / np.matmul(np.ones((N,1)),(4 * alpha_m * (times[i] - times[indexflsstart:indexflsend+2]).reshape(len(times[i] - times[indexflsstart:indexflsend+2]),1)).T), Amin1vector, finitecorrectiony) - \ + np.interp(np.matmul((Deltaz.reshape(len(Deltaz),1) ** 2),np.ones((1,indexflsend-indexflsstart+2))) / np.matmul(np.ones((N,1)),(4 * alpha_m * (times[i] - times[indexflsstart+1:indexflsend+3]).reshape(len(times[i] - times[indexflsstart:indexflsend+2]),1)).T), Amin1vector, finitecorrectiony)) + + NPCP = np.zeros((N, N)) + np.fill_diagonal(NPCP, CPCP) + + spacingtest = alpha_m * Deltat / SMatrixSorted[:, 1:]**2 / LimitNPSpacingTime + maxspacingtest = np.max(spacingtest,axis=0) + + if maxspacingtest[0] < 1: + maxindextoconsider = 0 + else: + maxindextoconsider = np.where(maxspacingtest > 1)[0][-1]+1 + + if mindexNPCP < maxindextoconsider + 1: + indicestocalculate = SortedIndices[:, mindexNPCP + 1:maxindextoconsider + 1] + indicestocalculatetranspose = indicestocalculate.T + indicestocalculatelinear = indicestocalculate.ravel() + indicestostorematrix = (indicestocalculate - 1) * N + np.arange(1, N) * np.ones((1, maxindextoconsider - mindexNPCP + 1)) + indicestostorematrixtranspose = indicestostorematrix.T + indicestostorelinear = indicestostorematrix.ravel() + NPCP[indicestostorelinear] = Deltaz[indicestocalculatelinear] / (4 * np.pi * self.krock.value * SMatrix[indicestostorelinear]) * erf(SMatrix[indicestostorelinear] / np.sqrt(4 * alpha_m * Deltat)) + + # Calculate and store neighbouring pipes for current pulse as set of line sources + if mindexNPCP > 1 and maxindextoconsider > 0: + lastindexfls = min(mindexNPCP, maxindextoconsider + 1) + indicestocalculate = SortedIndices[:, 1:lastindexfls] + indicestocalculatetranspose = indicestocalculate.T + indicestocalculatelinear = indicestocalculate.ravel() + indicestostorematrix = (indicestocalculate) * N + np.arange(N).reshape(-1,1) * np.ones((1, lastindexfls - 1), dtype=int) + indicestostorematrixtranspose = indicestostorematrix.T + indicestostorelinear = indicestostorematrix.ravel() + midpointindices = np.matmul(np.ones((lastindexfls - 1, 1)), np.arange( N ).reshape(1,N)).T + midpointsindices = midpointindices.ravel().astype(int) + rultimate = np.sqrt(np.square((midpointsx[midpointsindices].reshape(len(midpointsindices),1)*( np.ones((1, M + 1))) - distributionx[indicestocalculatelinear,:])) + + np.square((midpointsy[midpointsindices].reshape(len(midpointsindices),1)*( np.ones((1, M + 1))) - distributiony[indicestocalculatelinear,:])) + + np.square((midpointsz[midpointsindices].reshape(len(midpointsindices),1)*( np.ones((1, M + 1))) - distributionz[indicestocalculatelinear,:]))) + + NPCP[np.unravel_index(indicestostorelinear, NPCP.shape, 'F')] = Deltaz[indicestocalculatelinear] / M * np.sum((1 - erf(rultimate / np.sqrt(4 * alpha_m * Deltat))) / (4 * np.pi * self.krock.value * rultimate) * np.matmul(np.ones((N*(lastindexfls-1),1)),np.concatenate((np.array([1/2]), np.ones(M-1), np.array([1/2]))).reshape(-1,1).T), axis=1) + + # NPOP (= Neighbouring pipes, old pulses) + BB = np.zeros((N, 1)) + if i > 1 and lastneighbourtoconsider[i] > 0: + SMatrixRelevant = SMatrixSorted[:, 1 : int(lastneighbourtoconsider[i] + 1)] + SoverLRelevant = SoverLSorted[:, 1 : int(lastneighbourtoconsider[i]) + 1] + SortedIndicesRelevant = SortedIndices[:, 1 : int(lastneighbourtoconsider[i]) + 1] + maxtimeindexmatrix = alpha_m * np.ones((N * int(lastneighbourtoconsider[i]), 1)) * (times[i] - times[1:i]) / (SMatrixRelevant.ravel().reshape(-1,1) * np.ones((1,i-1)))**2 + + allindices = np.arange(N * int(lastneighbourtoconsider[i]) * (i - 1)) + pipeheatcomesfrom = np.matmul(SortedIndicesRelevant.T.ravel().reshape(len(SortedIndicesRelevant.ravel()),1), np.ones((1,i - 1))) + pipeheatgoesto = np.arange(N).reshape(N,1) * np.ones((1, int(lastneighbourtoconsider[i]))) + pipeheatgoesto = pipeheatgoesto.transpose().ravel().reshape(len(pipeheatgoesto.ravel()),1) * np.ones((1, i - 1)) + # Delete everything smaller than LimitNPSpacingTime + indicestoneglect = np.where((maxtimeindexmatrix.transpose()).ravel() < LimitNPSpacingTime)[0] + + maxtimeindexmatrix = np.delete(maxtimeindexmatrix, indicestoneglect) + allindices = np.delete(allindices, indicestoneglect) + indicesFoSlargerthan = np.where(maxtimeindexmatrix.ravel() > 10)[0] + indicestotakeforpsFoS = allindices[indicesFoSlargerthan] + + allindices2 = allindices.copy() + allindices2[indicesFoSlargerthan] = [] + SoverLinearized = SoverLRelevant.ravel().reshape(len(SoverLRelevant.ravel()),1) * np.ones((1, i - 1)) + indicestotakeforpsSoverL = np.where(SoverLinearized.transpose().ravel()[allindices2] > LimitSoverL)[0] + overallindicestotakeforpsSoverL = allindices2[indicestotakeforpsSoverL] + remainingindices = allindices2.copy() + + remainingindices=np.delete(remainingindices,indicestotakeforpsSoverL) + + NPOP = np.zeros((N * int(lastneighbourtoconsider[i]), i - 1)) + + # Use point source model when FoS is very large + if len(indicestotakeforpsFoS) > 0: + deltatlinear1 = np.ones(N * int(lastneighbourtoconsider[i]), 1) * (times[i] - times[1:i-1]) + deltatlinear1 = deltatlinear1.ravel()[indicestotakeforpsFoS] + deltatlinear2 = np.ones((N * int(lastneighbourtoconsider[i]), 1)) * (times[i] - times[0:i-2]) + deltatlinear2 = deltatlinear2[indicestotakeforpsFoS] + deltazlinear = pipeheatcomesfrom[indicestotakeforpsFoS] + SMatrixlinear = SMatrixRelevant.flatten(order='F') + NPOPFoS = Deltaz[deltazlinear] / (4 * np.pi * self.krock.value * SMatrixlinear[indicestotakeforpsFoS]) * (erfc(SMatrixlinear[indicestotakeforpsFoS] / np.sqrt(4 * alpha_m * deltatlinear2)) - + erfc(SMatrixlinear[indicestotakeforpsFoS] / np.sqrt(4 * alpha_m * deltatlinear1))) + + NPOP[indicestotakeforpsFoS] = NPOPFoS + + # Use point source model when SoverL is very large + if len(overallindicestotakeforpsSoverL) > 0: + deltatlinear1 = np.ones((N * int(lastneighbourtoconsider[i]), 1)) * (times[i] - times[1:i-2]).ravel() + deltatlinear1 = deltatlinear1[overallindicestotakeforpsSoverL] + deltatlinear2 = np.ones((N * int(lastneighbourtoconsider[i]), 1)) * (times[i] - times[0:i-2]).ravel() + deltatlinear2 = deltatlinear2[overallindicestotakeforpsSoverL] + deltazlinear = pipeheatcomesfrom[overallindicestotakeforpsSoverL] + SMatrixlinear = SMatrixRelevant.flatten(order='F') + NPOPSoverL = Deltaz[deltazlinear] / (4 * np.pi * self.krock.value * SMatrixlinear[overallindicestotakeforpsSoverL]) * (erfc(SMatrixlinear[overallindicestotakeforpsSoverL] / np.srt(4 * alpha_m * deltatlinear2)) - + erfc(SMatrixlinear[overallindicestotakeforpsSoverL] / np.sqrt(4 * alpha_m * deltatlinear1))) + + NPOP[overallindicestotakeforpsSoverL] = NPOPSoverL + + # Use finite line source model for remaining pipe segments + if len(remainingindices) > 0: + deltatlinear1 = np.ones((N * int(lastneighbourtoconsider[i]), 1)) * (times[i] - times[1:i]) + deltatlinear1 = (deltatlinear1.transpose()).ravel()[remainingindices] + deltatlinear2 = np.ones((N * int(lastneighbourtoconsider[i]), 1)) * (times[i] - times[0:i-1]) + deltatlinear2 = (deltatlinear2.transpose()).ravel()[remainingindices] + deltazlinear = (pipeheatcomesfrom.T).ravel()[remainingindices] + midpointstuff = (pipeheatgoesto.transpose()).ravel()[remainingindices] + rultimate = np.sqrt(np.square((midpointsx[midpointstuff.astype(int)].reshape(len(midpointsx[midpointstuff.astype(int)]),1)*( np.ones((1, M + 1))) - distributionx[deltazlinear.astype(int),:])) + + np.square((midpointsy[midpointstuff.astype(int)].reshape(len(midpointsy[midpointstuff.astype(int)]),1)*( np.ones((1, M + 1))) - distributiony[deltazlinear.astype(int),:])) + + np.square((midpointsz[midpointstuff.astype(int)].reshape(len(midpointsz[midpointstuff.astype(int)]),1)*( np.ones((1, M + 1))) - distributionz[deltazlinear.astype(int),:]))) + NPOPfls = Deltaz[deltazlinear.astype(int)].reshape(len(Deltaz[deltazlinear.astype(int)]),1).T / M * np.sum((-erf(rultimate / np.sqrt(4 * alpha_m * np.ravel(deltatlinear2).reshape(len(np.ravel(deltatlinear2)),1)*np.ones((1, M + 1)))) + erf(rultimate / np.sqrt(4 * alpha_m * np.ravel(deltatlinear1).reshape(len(np.ravel(deltatlinear1)),1)*np.ones((1, M + 1))))) / (4 * np.pi * self.krock.value * rultimate) * np.matmul((np.ones((len(remainingindices),1))),(np.concatenate((np.array([1/2]),np.ones(M - 1),np.array([1/2])))).reshape(-1,1).T), axis=1) + NPOPfls = NPOPfls.T + dimensions = NPOP.shape + NPOP=NPOP.T.ravel() + NPOP[remainingindices.reshape((len(remainingindices),1))] = NPOPfls + NPOP = NPOP.reshape((dimensions[1],dimensions[0])).T + + # Put everything together and calculate BB (= impact of all previous heat pulses from old neighbouring elements on current element at current time) + Qindicestotake = SortedIndicesRelevant.ravel().reshape((N * int(lastneighbourtoconsider[i]), 1))*np.ones((1,i-1)) + \ + np.ones((N * int(lastneighbourtoconsider[i]), 1)) * N * np.arange(i - 1) + Qindicestotake = Qindicestotake.astype(int) + Qlinear = Q.T.ravel()[Qindicestotake] + BBPS = NPOP * Qlinear + BBPS = np.sum(BBPS, axis=1) + BBPSindicestotake = np.arange(N).reshape((N, 1)) + N * np.arange(int(lastneighbourtoconsider[i])).reshape((1, int(lastneighbourtoconsider[i]))) + BBPSMatrix = BBPS[BBPSindicestotake] + BB = np.sum(BBPSMatrix, axis=1) + + if i > 1: + BBCPOP = np.sum(CPOP * Q[:, 1:i], axis=1) + else: + BBCPOP = np.zeros(N) + + # Populate L and R for fluid heat balance for first element (which has the injection temperature specified) + L[0, 0] = 1 / Deltat + uvector[0] / Deltaz[0] * (self.percent_implicit.value) * 2 + L[0, 2] = -4 / np.pi / Dvector[0]**2 / model.surfaceplant.rho_fluid.value / model.surfaceplant.cp_fluid.value + R[0, 0] = 1 / Deltat * Twprevious[0] + uvector[0] / Deltaz[0] * model.wellbores.Tinj.value * 2 - uvector[0] / Deltaz[0] * Twprevious[0] * (1 - self.percent_implicit.value) * 2 + + # Populate L and R for rock temperature equation for first element + L[1, 0] = 1 + L[1, 1] = -1 + L[1, 2] = Rtvector[0] + R[1, 0] = 0 + # Populate L and R for SBT algorithm for first element + L[2, np.arange(2,3*N,3)] = NPCP[0,0:N] + L[2,1] = 1 + R[2, 0] = -BBCPOP[0].item() - BB[0].item() + BBinitial[0].item() + + for iiii in range(2, N+1): + # Heat balance equation + L[0+(iiii - 1) * 3, (iiii - 1) * 3] = 1 / Deltat + uvector[iiii-1] / Deltaz[iiii-1] / 2 * (self.percent_implicit.value) * 2 + L[0+(iiii - 1) * 3, 2 + (iiii - 1) * 3] = -4 / np.pi / Dvector[iiii-1] ** 2 / model.surfaceplant.rho_fluid.value / model.surfaceplant.cp_fluid.value + + if iiii == len(xinj): # Upcoming pipe has first element temperature sum of all incoming water temperatures + for j in range(len(lateralendpoints)): + L[0+ (iiii - 1) * 3, 0 + (lateralendpoints[j]) * 3] = -ulateral[j] / Deltaz[iiii-1] / 2 / (self.percent_implicit.value) * 2 + R[0+(iiii - 1) * 3, 0] = 1 / Deltat * Twprevious[iiii-1] + uvector[iiii-1] / Deltaz[iiii-1] * ( + -Twprevious[iiii-1] + np.sum(lateralflowallocation[j] * Twprevious[lateralendpoints[j]])) / 2 * ( + 1 - self.percent_implicit.value) * 2 + else: + L[0+(iiii-1) * 3, 0 + (int(previouswaterelements[iiii-1])) * 3] = -uvector[iiii-1] / Deltaz[iiii-1] / 2 * ( + self.percent_implicit.value) * 2 + R[0+(iiii-1) * 3, 0] = 1 / Deltat * Twprevious[iiii-1] + uvector[iiii-1] / Deltaz[iiii-1] * ( + -Twprevious[iiii-1] + Twprevious[int(previouswaterelements[iiii-1])]) / 2 * (1 - self.percent_implicit.value) * 2 + + # Rock temperature equation + L[1 + (iiii - 1) * 3, (iiii - 1) * 3] = 1 + L[1 + (iiii - 1) * 3, 1 + (iiii - 1) * 3] = -1 + L[1 + (iiii - 1) * 3, 2 + (iiii - 1) * 3] = Rtvector[iiii-1] + R[1 + (iiii - 1) * 3, 0] = 0 + + # SBT equation + L[2 + (iiii - 1) * 3, np.arange(2,3*N,3)] = NPCP[iiii-1, :N] + L[2 + (iiii - 1) * 3, 1 + (iiii - 1) * 3] = 1 + R[2 + (iiii - 1) * 3, 0] = -BBCPOP[iiii - 1].item() - BB[iiii - 1].item() + BBinitial[iiii - 1].item() + + # Solving the linear system of equations + Sol = np.linalg.solve(L, R) + + # Extracting Q array for current heat pulses + Q[:, i] = Sol.ravel()[2::3] + + # Extracting fluid temperature + TwMatrix[i, :] = Sol.ravel()[np.arange(0,3*N,3)] + + # Storing fluid temperature for the next time step + Twprevious = Sol.ravel()[np.arange(0,3*N,3)] + + # Calculating the fluid outlet temperature at the top of the first element + top_element_index = len(xinj) + len(xprod) - 3 + self.Tresoutput.value[i] = Twprevious[top_element_index] + (Twprevious[top_element_index] - Twprevious[top_element_index - 1]) * 0.5 + if False: + print(times[i] / 3.154e+7, ',', self.Tresoutput.value[i]) + + # Save the nonlinear Time results as 2D array Output Parameter + self.NonLinearTime_temperature.value = np.column_stack((times, self.Tresoutput.value)) + + # define the linear time array, in years + self.timevector.value = np.linspace(0, model.surfaceplant.plant_lifetime.value, model.economics.timestepsperyear.value * model.surfaceplant.plant_lifetime.value) + + # Now interpolate that non-linear time array into a linear array in time with associated values. + # times locally are in seconds, so convert to years to match the linear time array, which is in years + times_years = times / 3.154e+7 + times_years[0] = 0.0 # Ensure the first time step is 0.0 + times_years[-1] = model.surfaceplant.plant_lifetime.value # Ensure the last time step is the plant lifetime + + # Calculate the maximum temperature for the SBT output. + max_temperature = np.max(self.Tresoutput.value) + + # Replace the first year of values in the SBT output array + # with the maximum temperature for the first year of the SBT output. + # This moderates the behavior for the first few months of the SBT output. + i = 0 + while times_years[i] < 1.0: + self.Tresoutput.value[i] = max_temperature + i = i + 1 + + linear_values = [] + # interpolate the values of self.Tresoutput.value to the linear time array + for t in self.timevector.value: + linear_values.append(interpolator(t, times_years, self.Tresoutput.value)) + self.Tresoutput.value = linear_values + + # Calculate the Initial Reservoir Heat Content + self.InitialReservoirHeatContent.value = mstore[0] * model.surfaceplant.cp_fluid.value * (self.Tresoutput.value[0] - Tinstore[0]) / 1e6 # Calculates the heat production [MW] + + model.logger.info(f'complete {str(__class__)}: {sys._getframe().f_code.co_name}') diff --git a/src/geophires_x/SBTWellbores.py b/src/geophires_x/SBTWellbores.py new file mode 100644 index 00000000..6873597b --- /dev/null +++ b/src/geophires_x/SBTWellbores.py @@ -0,0 +1,288 @@ +import math +import sys +import numpy as np +from pint.facets.plain import PlainQuantity + +from geophires_x.WellBores import (WellBores, RameyCalc, WellPressureDrop, \ + ProdPressureDropsAndPumpingPowerUsingImpedenceModel, + ProdPressureDropAndPumpingPowerUsingIndexes, calculate_total_drilling_lengths_m) +from .Parameter import floatParameter, intParameter, boolParameter, OutputParameter, ReadParameter +from geophires_x.GeoPHIRESUtils import vapor_pressure_water_kPa, quantity, static_pressure_MPa, \ + heat_capacity_water_J_per_kg_per_K +from geophires_x.GeoPHIRESUtils import density_water_kg_per_m3 +from geophires_x.GeoPHIRESUtils import viscosity_water_Pa_sec +from .Units import * +import geophires_x.Model as Model +from .OptionList import ReservoirModel, Configuration, WorkingFluid + + +class SBTWellbores(WellBores): + """ + SBTWellbores Child class of AGSWellBores; it is the same, but has advanced SBT closed-loop functionality + """ + + def __init__(self, model: Model): + """ + The __init__ function is the constructor for a class. It is called whenever an instance of the class is created. + The __init__ function can take arguments, but self is always the first one. Self refers to the instance of the + object that has already been created, and it's used to access variables that belong to that object. + :param model: The container class of the application, giving access to everything else, including the logger + :type model: :class:`~geophires_x.Model.Model` + :return: Nothing, and is used to initialize the class + """ + model.logger.info(f'Init {__class__!s}: {sys._getframe().f_code.co_name}') + + # Initialize the superclass first to gain access to those variables + super().__init__(model) + sclass = str(__class__).replace("", "") + self.MyPath = os.path.abspath(__file__) + self.Tini = 0.0 + + # Set up 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. + # If you choose to subclass this master class, you can do so before or after you create your own parameters. + # If you do, you can also choose to call this method from you class, which will effectively add and set all + # these parameters to your class. + # NB: inputs we already have ("already have it") need to be set at ReadParameter time so values are set at the + # last possible time + + self.vertical_section_length = self.ParameterDict[self.vertical_section_length.Name] = floatParameter( + 'Vertical Section Length', + DefaultValue=2000.0, + Min=0.01, + Max=10000.0, + UnitType=Units.LENGTH, + PreferredUnits=LengthUnit.METERS, + CurrentUnits=LengthUnit.METERS, + ErrMessage='assume default for Vertical Section Length (2000.0 meters)', + ToolTipText='length/depth to the bottom of the vertical wellbores' + ) + self.vertical_wellbore_spacing = self.ParameterDict[self.vertical_wellbore_spacing.Name] = floatParameter( + 'Vertical Wellbore Spacing', + DefaultValue=100.0, + Min=0.01, + Max=10000.0, + UnitType=Units.LENGTH, + PreferredUnits=LengthUnit.METERS, + CurrentUnits=LengthUnit.METERS, + ErrMessage='assume default for Vertical Wellbore Spacing (100.0 meters)', + ToolTipText='Horizontal distance between vertical wellbores' + ) + self.lateral_spacing = self.ParameterDict[self.lateral_spacing.Name] = floatParameter( + 'Lateral Spacing', + DefaultValue=100.0, + Min=0.01, + Max=10000.0, + UnitType=Units.LENGTH, + PreferredUnits=LengthUnit.METERS, + CurrentUnits=LengthUnit.METERS, + ErrMessage='assume default for Lateral Spacing (100.0 meters)', + ToolTipText='Horizontal distance between laterals' + ) + self.lateral_inclination_angle = self.ParameterDict[self.lateral_inclination_angle.Name] = floatParameter( + 'Lateral Inclination Angle', + DefaultValue=20.0, + Min=0.0, + Max=89.999999, + UnitType=Units.ANGLE, + PreferredUnits=AngleUnit.DEGREES, + CurrentUnits=AngleUnit.DEGREES, + ErrMessage='assume default for Lateral Inclination Angle (20.0 degrees)', + ToolTipText='Inclination of the lateral section, where 0 degrees would mean vertical while 90 degrees is pure horizontal' + ) + self.element_length = self.ParameterDict[self.element_length.Name] = floatParameter( + 'Discretization Length', + DefaultValue=250.0, + Min=0.01, + Max=10000.0, + UnitType=Units.LENGTH, + PreferredUnits=LengthUnit.METERS, + CurrentUnits=LengthUnit.METERS, + ErrMessage='assume default for Discretization Length (250.0 meters)', + ToolTipText='distance between sample point along length of model' + ) + self.junction_depth = self.ParameterDict[self.junction_depth.Name] = floatParameter( + 'Junction Depth', + DefaultValue=4000.0, + Min=1000, + Max=15000.0, + UnitType=Units.LENGTH, + PreferredUnits=LengthUnit.METERS, + CurrentUnits=LengthUnit.METERS, + ErrMessage='assume default for Junction Depth (4000.0 meters)', + ToolTipText='vertical depth where the different laterals branch out (where the multilateral section starts, second deepest depth of model)' + ) + self.lateral_endpoint_depth = self.ParameterDict[self.lateral_endpoint_depth.Name] = floatParameter( + 'Lateral Endpoint Depth', + DefaultValue=7000.0, + Min=1000, + Max=15000.0, + UnitType=Units.LENGTH, + PreferredUnits=LengthUnit.METERS, + CurrentUnits=LengthUnit.METERS, + ErrMessage='assume default for Lateral Endpoint Depth (7000.0 meters)', + ToolTipText='vertical depth where the lateral section ends (tip of the multilateral section, deepest depth of model)' + ) + + model.logger.info(f'complete {__class__!s}: {sys._getframe().f_code.co_name}') + + def __str__(self): + return 'SBTWellbores' + + 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 + :type model: :class:`~geophires_x.Model.Model` + :return: None + """ + model.logger.info(f'Init {str(__class__)}: {sys._getframe().f_code.co_name}') + super().read_parameters(model) # read the default parameters + # if we call super, we don't need to deal with setting the parameters here, just deal with the special cases + # for the variables in this class because the call to the super.readparameters will set all the variables, + # including the ones that are specific to this class + + if len(model.InputParameters) > 0: + # loop through all the parameters that the user wishes to set, looking for parameters that match this object + for item in self.ParameterDict.items(): + ParameterToModify = item[1] + key = ParameterToModify.Name.strip() + if key in model.InputParameters: + ParameterReadIn = model.InputParameters[key] + # just handle special cases for this class - the call to super set all the values, + # including the value unique to this class + else: + model.logger.info("No parameters read because no content provided") + + model.logger.info(f'complete {str(__class__)}: {sys._getframe().f_code.co_name}') + + def Calculate(self, model: Model) -> None: + """ + The calculate function verifies, initializes, and extracts the values from the AGS model + :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(f'Init {__class__!s}: {sys._getframe().f_code.co_name}') + self.Tini = np.max(model.reserv.Tresoutput.value) # initial temperature of the reservoir + self.ProducedTemperature.value = np.array(model.reserv.Tresoutput.value) + + # Calculate the drilled lengths + self.tot_vert_m.value = self.vertical_section_length.quantity().to('km').magnitude + input_vert_depth_km = self.vertical_section_length.quantity().to('km').magnitude + output_vert_depth_km = self.lateral_endpoint_depth.quantity().to('km').magnitude + junction_depth_km = self.junction_depth.quantity().to('km').magnitude + angle_rad = ((90.0 * np.pi) / 180) - self.lateral_inclination_angle.quantity().to('radians').magnitude + + # Use the parents class implementation of calculate_total_drilling_lengths_m + self.total_drilled_length.value, self.tot_vert_m.value, self.tot_lateral_m.value, self.tot_to_junction_m.value = \ + calculate_total_drilling_lengths_m(self.Configuration.value, + self.numnonverticalsections.value, + self.Nonvertical_length.value / 1000.0, + input_vert_depth_km, + output_vert_depth_km, + self.nprod.value, + self.ninj.value, + junction_depth_km, + angle_rad) + self.total_drilled_length.value = self.total_drilled_length.value / 1000.0 # convert to km + + # Now use the parent's calculation to calculate the upgoing and downgoing pressure drops and pumping power + self.PumpingPower.value = [0.0] * len(self.ProducedTemperature.value) # initialize the array + if self.productionwellpumping.value: + self.rhowaterinj = density_water_kg_per_m3( + model.reserv.Tsurf.value, + pressure=model.reserv.hydrostatic_pressure() + ) * np.linspace(1, 1, + len(self.ProducedTemperature.value)) + + self.rhowaterprod = density_water_kg_per_m3( + model.reserv.Trock.value, + pressure=model.reserv.hydrostatic_pressure() + ) * np.linspace(1, 1, len(self.ProducedTemperature.value)) + + self.DPProdWell.value, f3, vprod, self.rhowaterprod = WellPressureDrop(model, + self.ProducedTemperature.value, + self.prodwellflowrate.value, + self.prodwelldiam.value, + self.impedancemodelused.value, + self.vertical_section_length.value) + if self.impedancemodelused.value: # assumed everything stays liquid throughout + self.DPOverall.value, UpgoingPumpingPower, self.DPProdWell.value, self.DPReserv.value, self.DPBouyancy.value = \ + ProdPressureDropsAndPumpingPowerUsingImpedenceModel( + f3, vprod, + self.rhowaterprod, self.rhowaterinj, self.rhowaterprod, + self.vertical_section_length.value, self.prodwellflowrate.value, + self.prodwelldiam.value, self.impedance.value, + self.nprod.value, model.reserv.waterloss.value, model.surfaceplant.pump_efficiency.value) + self.DPOverall.value, DowngoingPumpingPower, self.DPProdWell.value, self.DPReserv.value, self.DPBouyancy.value = \ + ProdPressureDropsAndPumpingPowerUsingImpedenceModel( + f3, vprod, + self.rhowaterinj, self.rhowaterprod, self.rhowaterprod, + self.vertical_section_length.value, self.prodwellflowrate.value, + self.injwelldiam.value, self.impedance.value, + self.ninj.value, model.reserv.waterloss.value, model.surfaceplant.pump_efficiency.value, 90.0, False) + self.DPOverall.value, NonverticalPumpingPower, self.DPProdWell.value, self.DPReserv.value, self.DPBouyancy.value = \ + ProdPressureDropsAndPumpingPowerUsingImpedenceModel( + f3, vprod, + self.rhowaterinj, self.rhowaterprod, self.rhowaterprod, + self.vertical_section_length.value, self.prodwellflowrate.value, + self.injwelldiam.value, self.impedance.value, + self.ninj.value, model.reserv.waterloss.value, model.surfaceplant.pump_efficiency.value, + self.lateral_inclination_angle.value, False) + + else: # PI is used for both the verticals + UpgoingPumpingPower, self.PumpingPowerProd.value, self.DPProdWell.value, self.Pprodwellhead.value = \ + ProdPressureDropAndPumpingPowerUsingIndexes( + model, self.productionwellpumping.value, + self.usebuiltinppwellheadcorrelation, + model.reserv.Trock.value, self.vertical_section_length.value, + self.ppwellhead.value, self.PI.value, + self.prodwellflowrate.value, f3, vprod, + self.prodwelldiam.value, self.nprod.value, model.surfaceplant.pump_efficiency.value, + self.rhowaterprod) + + DowngoingPumpingPower, ppp2, dppw, ppwh = ProdPressureDropAndPumpingPowerUsingIndexes( + model, self.productionwellpumping.value, + self.usebuiltinppwellheadcorrelation, + model.reserv.Trock.value, self.vertical_section_length.value, + self.ppwellhead.value, self.PI.value, + self.prodwellflowrate.value, f3, vprod, + self.injwelldiam.value, self.nprod.value, model.surfaceplant.pump_efficiency.value, + self.rhowaterinj) + + # Calculate Nonvertical Pressure Drop + # TODO assume no pressure drop in the non-vertical section for now + #NonverticalPumpingPower = [0.0] * len(DowngoingPumpingPower) # initialize the array + #self.NonverticalPressureDrop.value = [0.0] * len(DowngoingPumpingPower) # initialize the array + #NonverticalPumpingPower = [0.0] * len(DowngoingPumpingPower) # initialize the array + + # recalculate the pumping power by looking at the difference between the upgoing and downgoing and the nonvertical + self.PumpingPower.value = np.array(DowngoingPumpingPower) + np.array(NonverticalPumpingPower) + np.array(UpgoingPumpingPower) + self.PumpingPower.value = [max(x, 0.) for x in self.PumpingPower.value] # cannot be negative, so set to 0 + + # calculate water values based on initial temperature, Need this for surface plant output calculation + rho_water = density_water_kg_per_m3(np.average(self.ProducedTemperature.value), model.reserv.lithostatic_pressure()) + model.reserv.cpwater.value = heat_capacity_water_J_per_kg_per_K(np.average(self.ProducedTemperature.value), model.reserv.hydrostatic_pressure(),) + + # set pumping power to zero for all times, assuming that the thermosphere wil always + # make pumping of working fluid unnecessary + self.PumpingPower.value = [0.0] * (len(self.DPOverall.value)) + self.PumpingPower.value = self.DPOverall.value * self.prodwellflowrate.value / rho_water / model.surfaceplant.pump_efficiency.value / 1E3 + # in GEOPHIRES v1.2, negative pumping power values become zero + # (b/c we are not generating electricity) = thermosiphon is happening! + self.PumpingPower.value = [max(x, 0.) for x in self.PumpingPower.value] + + model.logger.info(f'complete {str(__class__)}: {sys._getframe().f_code.co_name}') + + def CalculateNonverticalPressureDrop(self, model, value, time_max, al): + pass + diff --git a/src/geophires_x/SurfacePlant.py b/src/geophires_x/SurfacePlant.py index eb850255..90b97b19 100644 --- a/src/geophires_x/SurfacePlant.py +++ b/src/geophires_x/SurfacePlant.py @@ -384,6 +384,50 @@ def __init__(self, model: Model): ErrMessage="assume default number of years in construction (1)", ToolTipText="Number of years spent in construction (assumes whole years, no fractions)" ) + self.cp_fluid = self.ParameterDict[self.cp_fluid.Name] = floatParameter( + "Working Fluid Heat Capacity", + UnitType=Units.HEAT_CAPACITY, + PreferredUnits=HeatCapacityUnit.JPERKGPERK, + CurrentUnits=HeatCapacityUnit.JPERKGPERK, + DefaultValue=4200.0, + Min=0.0, + Max=10_000.0, + ErrMessage="assume default working fluid heat capacity (4200)", + ToolTipText="Heat capacity of the working fluid" + ) + self.rho_fluid = self.ParameterDict[self.rho_fluid.Name] = floatParameter( + "Working Fluid Density", + UnitType=Units.DENSITY, + PreferredUnits=DensityUnit.KGPERMETERS3, + CurrentUnits=DensityUnit.KGPERMETERS3, + DefaultValue=1000.0, + Min=0.0, + Max=10_000.0, + ErrMessage="assume default working fluid density (1000)", + ToolTipText="Density of the working fluid" + ) + self.k_fluid = self.ParameterDict[self.k_fluid.Name] = floatParameter( + "Working Fluid Thermal Conductivity", + UnitType=Units.THERMAL_CONDUCTIVITY, + PreferredUnits=ThermalConductivityUnit.WPERMPERK, + CurrentUnits=ThermalConductivityUnit.WPERMPERK, + DefaultValue=0.68, + Min=0.0, + Max=10.0, + ErrMessage="assume default working fluid thermal conductivity (0.68)", + ToolTipText="Thermal conductivity of the working fluid" + ) + self.mu_fluid = self.ParameterDict[self.mu_fluid.Name] = floatParameter( + "Working Fluid Dynamic Viscosity", + UnitType=Units.DYNAMIC_VISCOSITY, + PreferredUnits=Dynamic_ViscosityUnit.PASCALSEC, + CurrentUnits=Dynamic_ViscosityUnit.PASCALSEC, + DefaultValue=600*10**-6, + Min=0.0, + Max=1, + ErrMessage="assume default fluid dynamic density (600*10**-6)", + ToolTipText="Dynamic viscosity of the working fluid" + ) # local variable initialization self.setinjectionpressurefixed = False diff --git a/src/geophires_x/Units.py b/src/geophires_x/Units.py index 89d0b155..e1c4f07b 100644 --- a/src/geophires_x/Units.py +++ b/src/geophires_x/Units.py @@ -19,6 +19,7 @@ class Units(IntEnum): NONE = auto() CHOICE = auto() LENGTH = auto() + ANGLE = auto() AREA = auto() VOLUME = auto() MASS = auto() @@ -60,6 +61,13 @@ class Units(IntEnum): POWERPERUNITVOLUME = auto() DECAY_RATE=auto() INFLATION_RATE=auto() + DYNAMIC_VISCOSITY = auto() + + +class AngleUnit(str, Enum): + """Angle Units""" + DEGREES = "degrees" + RADIANS = "radians" class TemperatureUnit(str, Enum): @@ -351,3 +359,8 @@ class Decay_RateUnit(str,Enum): class Inflation_RateUnit(str,Enum): """Decay rate Units""" KPASCALPERYEAR = "kPa/yr" + + +class Dynamic_ViscosityUnit(str,Enum): + """Dynamic Viscosity Units""" + PASCALSEC = "PaSec" diff --git a/src/geophires_x/WellBores.py b/src/geophires_x/WellBores.py index 80f9b928..4866f9f0 100644 --- a/src/geophires_x/WellBores.py +++ b/src/geophires_x/WellBores.py @@ -15,6 +15,71 @@ # code from Koenraad +def calculate_total_drilling_lengths_m(Configuration, numnonverticalsections: int, nonvertical_length_km: float, + InputDepth_km: float, OutputDepth_km: float, nprod:int, ninj:int, + junction_depth_km: float = 0.0, angle_rad: float = 0.0) -> tuple: + """ + returns the total length, vertical length, and non-vertical lengths, depending on the configuration + :param Configuration: configuration of the well + :type Configuration: :class:`~geophires + :param numnonverticalsections: number of non-vertical sections + :type numnonverticalsections: int + :param nonvertical_length_km: length of non-vertical sections in km + :type nonvertical_length_km: float + :param InputDepth_km: depth of the well in km + :type InputDepth_km: float + :param OutputDepth_km: depth of the output end of the well in km, if U shaped, and not horizontal + :type OutputDepth_km: float + :param nprod: number of production wells + :type nprod: int + :param ninj: number of injection wells + :param junction_depth_km: depth of the junction in km + :type junction_depth_km: float + :param angle_rad: angle of the well in radians, from horizontal + :type angle_rad: float + :return: total length, vertical length, lateral, and junction lengths in meters + :rtype: tuple + """ + tot_pipe_length_m = vertical_pipe_length_m = lateral_pipe_length_m = tot_to_junction_m = 0.0 + if Configuration is Configuration.ULOOP: + # Total drilling depth of both wells and laterals in U-loop [m] + vertical_pipe_length_m = (nprod * InputDepth_km * 1000.0) + (ninj * OutputDepth_km * 1000.0) + lateral_pipe_length_m = numnonverticalsections * nonvertical_length_km * 1000.0 + + elif Configuration is Configuration.COAXIAL: + # Total drilling depth of well and lateral in co-axial case [m] - is not necessarily only vertical + vertical_pipe_length_m = (nprod + ninj) * InputDepth_km * 1000.0 + lateral_pipe_length_m = numnonverticalsections * nonvertical_length_km * 1000.0 + + elif Configuration is Configuration.VERTICAL: + # Total drilling depth of well in vertical case [m] + vertical_pipe_length_m = (nprod + ninj) * InputDepth_km * 1000.0 + lateral_pipe_length_m = 0.0 + + elif Configuration is Configuration.L: + # Total drilling depth of well in L case [m] + vertical_pipe_length_m = (nprod + ninj) * InputDepth_km * 1000.0 + lateral_pipe_length_m = numnonverticalsections * nonvertical_length_km * 1000.0 + + elif Configuration is Configuration.EAVORLOOP: + # Total drilling length of well in EavorLoop [m] + vertical_pipe_length_m = (nprod + ninj) * InputDepth_km * 1000.0 + + # now calculate the distance from the bottom of the vertical to the junction of the laterals [m] + O1 = (junction_depth_km - InputDepth_km) * 1000.0 # in meters + tot_to_junction_m = (O1 / math.sin(angle_rad)) * 2 # there are two of these of each EavorLoop + + # now calculate the distance from the junction of the laterals to the end of the laterals [m] + O2 = (OutputDepth_km - junction_depth_km) * 1000.0 # in meters + lateral_pipe_length_m = (O2 / math.sin(angle_rad)) * 2 # there are two of these of each lateral of an EavorLoop + lateral_pipe_length_m = lateral_pipe_length_m * numnonverticalsections # there are numnonverticalsections of these + else: + raise ValueError(f'Invalid Configuration: {Configuration}') + + tot_pipe_length_m = vertical_pipe_length_m + lateral_pipe_length_m + tot_to_junction_m + return tot_pipe_length_m, vertical_pipe_length_m, lateral_pipe_length_m, tot_to_junction_m + + def InjectionReservoirPressurePredictor(project_lifetime_yr: int, timesteps_per_year: int, initial_pressure_kPa: float, inflation_rate: float) -> list: """ @@ -261,12 +326,12 @@ def InjectionWellPressureDrop(model: Model, Taverage: float, wellflowrate: float return DPWell, f1, v, rhowater - def ProdPressureDropsAndPumpingPowerUsingImpedenceModel(f3: float, vprod: float, rhowaterinj: float, rhowaterprod: float, rhowaterreservoir: float, depth: float, wellflowrate: float, prodwelldiam: float, impedance: float, nprod: int, waterloss: float, - pumpeff: float) -> tuple: + pumpeff: float, tilt: float = 90.0, + trim_neg_to_zero: bool = True) -> tuple: """ Calculate Pressure Drops and Pumping Power needed for the production well using the Impedance Model :param f3: friction factor [-] @@ -293,6 +358,10 @@ def ProdPressureDropsAndPumpingPowerUsingImpedenceModel(f3: float, vprod: float, :type waterloss: float :param pumpeff: pump efficiency [-] :type pumpeff: float + :param tilt: tilt of the well from the junction box to the bottom of the laterals [degrees] + :type tilt: float + :param trim_neg_to_zero: whether to trim negative values of pumping power to zero + :type trim_neg_to_zero: bool :return: tuple of DPOverall, PumpingPower, DPProdWell, DPReserv, DPBouyancy :rtype: tuple """ @@ -303,7 +372,8 @@ def ProdPressureDropsAndPumpingPowerUsingImpedenceModel(f3: float, vprod: float, DPReserv = impedance * nprod * wellflowrate * 1000. / rhowaterreservoir # buoyancy pressure drop [kPa] - DPBouyancy = (rhowaterprod - rhowaterinj) * depth * 9.81 / 1E3 # /1E3 to convert from Pa to kPa + gravity_tilt = 9.81 * np.sin(np.radians(tilt)) + DPBouyancy = (rhowaterprod - rhowaterinj) * depth * gravity_tilt / 1E3 # /1E3 to convert from Pa to kPa # overall pressure drop DPOverall = DPReserv + DPProdWell + DPBouyancy @@ -313,7 +383,8 @@ def ProdPressureDropsAndPumpingPowerUsingImpedenceModel(f3: float, vprod: float, PumpingPower = DPOverall * nprod * wellflowrate * (1 + waterloss) / rhowaterinj / pumpeff / 1E3 # in GEOPHIRES v1.2, negative pumping power values become zero (b/c we are not generating electricity) - PumpingPower = [0. if x < 0. else x for x in PumpingPower] + if trim_neg_to_zero: + PumpingPower = [0. if x < 0. else x for x in PumpingPower] return DPOverall, PumpingPower, DPProdWell, DPReserv, DPBouyancy @@ -889,7 +960,7 @@ def __init__(self, model: Model): self.Configuration = self.ParameterDict[self.Configuration.Name] = intParameter( "Closed-loop Configuration", DefaultValue=Configuration.VERTICAL.int_value, - AllowableRange=[1, 2, 3, 4], + AllowableRange=[1, 2, 3, 4, 5], ValuesEnum=Configuration, UnitType=Units.NONE, Required=True, @@ -900,7 +971,7 @@ def __init__(self, model: Model): self.Configuration = self.ParameterDict[self.Configuration.Name] = intParameter( "Well Geometry Configuration", DefaultValue=Configuration.VERTICAL.int_value, - AllowableRange=[1, 2, 3, 4], + AllowableRange=[1, 2, 3, 4, 5], ValuesEnum=Configuration, UnitType=Units.NONE, Required=True, @@ -1105,6 +1176,34 @@ def __init__(self, model: Model): PreferredUnits=PressureUnit.KPASCAL, CurrentUnits=PressureUnit.KPASCAL ) + self.total_drilled_length = self.OutputParameterDict[self.total_drilled_length.Name] = OutputParameter( + Name="Total length of all drilling", + value=-1.0, + UnitType=Units.LENGTH, + PreferredUnits=LengthUnit.KILOMETERS, + CurrentUnits=LengthUnit.KILOMETERS + ) + self.tot_vert_m = self.OutputParameterDict[self.tot_vert_m.Name] = OutputParameter( + Name="Total length of vertical drilling", + value=-1.0, + UnitType=Units.LENGTH, + PreferredUnits=LengthUnit.METERS, + CurrentUnits=LengthUnit.METERS + ) + self.tot_to_junction_m = self.OutputParameterDict[self.tot_to_junction_m.Name] = OutputParameter( + Name="Total length from lateral junction to base of vertical drilling", + value=-1.0, + UnitType=Units.LENGTH, + PreferredUnits=LengthUnit.METERS, + CurrentUnits=LengthUnit.METERS + ) + self.tot_lateral_m = self.OutputParameterDict[self.tot_lateral_m.Name] = OutputParameter( + Name="Total length of lateral drilling", + value=-1.0, + UnitType=Units.LENGTH, + PreferredUnits=LengthUnit.METERS, + CurrentUnits=LengthUnit.METERS + ) def __str__(self): return "WellBores" diff --git a/tests/examples/example_SBT_Hi_T.out b/tests/examples/example_SBT_Hi_T.out new file mode 100644 index 00000000..30fc221c --- /dev/null +++ b/tests/examples/example_SBT_Hi_T.out @@ -0,0 +1,224 @@ + ***************** + ***CASE REPORT*** + ***************** + +Simulation Metadata +---------------------- + GEOPHIRES Version: 3.5.2 + Simulation Date: 2024-09-02 + Simulation Time: 13:59 + Calculation Time: 14.837 sec + + ***SUMMARY OF RESULTS*** + + End-Use Option: Electricity + Average Net Electricity Production: 9.62 MW + Electricity breakeven price: 18.23 cents/kWh + Number of production wells: 1 + Number of injection wells: 1 + Flowrate per production well: 80.0 kg/sec + Well depth (or total length, if not vertical): 5.8 kilometer + Geothermal gradient: 60 degC/km + + + ***ECONOMIC PARAMETERS*** + + Economic Model = BICYCLE + Accrued financing during construction: 0.00 + Project lifetime: 30 yr + Capacity factor: 90.0 % + Project NPV: 0.96 MUSD + Project IRR: 6.32 % + Project VIR=PI=PIR: 1.01 + Project MOIC: 0.65 + Project Payback Period: 13.59 yr + Estimated Jobs Created: 0 + + ***ENGINEERING PARAMETERS*** + + Number of Production Wells: 1 + Number of Injection Wells: 1 + Well depth (or total length, if not vertical): 5.8 kilometer + Water loss rate: 0.0 + Pump efficiency: 75.0 + Injection temperature: 56.9 degC + Production Wellbore heat transmission calculated with Ramey's model + Average production well temperature drop: 0.0 degC + Flowrate per production well: 80.0 kg/sec + Injection well casing ID: 8.500 in + Production well casing ID: 8.500 in + Number of times redrilling: 0 + Power plant type: Supercritical ORC + + + ***RESOURCE CHARACTERISTICS*** + + Maximum reservoir temperature: 400.0 degC + Number of segments: 1 + Geothermal gradient: 60 degC/km + + + ***RESERVOIR PARAMETERS*** + +The AGS models contain an intrinsic reservoir model that doesn't expose values that can be used in extensive reporting. + + + ***RESERVOIR SIMULATION RESULTS*** + + Maximum Production Temperature: 313.6 degC + Average Production Temperature: 218.7 degC + Minimum Production Temperature: 205.9 degC + Initial Production Temperature: 313.6 degC +The AGS models contain an intrinsic reservoir model that doesn't expose values that can be used in extensive reporting. + + + ***CAPITAL COSTS (M$)*** + + Drilling and completion costs: 104.87 MUSD + Drilling and completion costs per vertical production well: 6.19 MUSD + Drilling and completion costs per vertical injection well: 6.19 MUSD + Drilling and completion costs per non-vertical sections: 92.48 MUSD + Stimulation costs: 0.00 MUSD + Surface power plant costs: 52.76 MUSD + Field gathering system costs: 0.97 MUSD + Total surface equipment costs: 53.72 MUSD + Exploration costs: 0.00 MUSD + Total capital costs: 158.59 MUSD + + + ***OPERATING AND MAINTENANCE COSTS (M$/yr)*** + + Wellfield maintenance costs: 1.44 MUSD/yr + Power plant maintenance costs: 1.93 MUSD/yr + Water costs: 0.00 MUSD/yr + Total operating and maintenance costs: 3.37 MUSD/yr + + + ***SURFACE EQUIPMENT SIMULATION RESULTS*** + + Initial geofluid availability: 0.46 MW/(kg/s) + Maximum Total Electricity Generation: 17.60 MW + Average Total Electricity Generation: 9.62 MW + Minimum Total Electricity Generation: 8.40 MW + Initial Total Electricity Generation: 17.60 MW + Maximum Net Electricity Generation: 17.60 MW + Average Net Electricity Generation: 9.62 MW + Minimum Net Electricity Generation: 8.40 MW + Initial Net Electricity Generation: 17.60 MW + Average Annual Total Electricity Generation: 75.03 GWh + Average Annual Net Electricity Generation: 75.03 GWh + Average Pumping Power: 0.00 MW + + ************************************************************ + * HEATING, COOLING AND/OR ELECTRICITY PRODUCTION PROFILE * + ************************************************************ + YEAR THERMAL GEOFLUID PUMP NET FIRST LAW + DRAWDOWN TEMPERATURE POWER POWER EFFICIENCY + (degC) (MW) (MW) (%) + 1 1.0000 313.57 0.0000 17.6004 19.8758 + 2 0.7703 241.55 0.0000 11.9582 18.7707 + 3 0.7466 234.12 0.0000 11.2166 18.3442 + 4 0.7332 229.89 0.0000 10.7920 18.0811 + 5 0.7237 226.92 0.0000 10.4929 17.8876 + 6 0.7163 224.62 0.0000 10.2615 17.7333 + 7 0.7104 222.76 0.0000 10.0753 17.6061 + 8 0.7054 221.20 0.0000 9.9182 17.4968 + 9 0.7012 219.86 0.0000 9.7845 17.4022 + 10 0.6975 218.71 0.0000 9.6694 17.3197 + 11 0.6942 217.67 0.0000 9.5653 17.2443 + 12 0.6911 216.72 0.0000 9.4701 17.1745 + 13 0.6883 215.83 0.0000 9.3813 17.1088 + 14 0.6856 214.98 0.0000 9.2967 17.0456 + 15 0.6832 214.23 0.0000 9.2222 16.9895 + 16 0.6809 213.52 0.0000 9.1514 16.9358 + 17 0.6788 212.85 0.0000 9.0857 16.8855 + 18 0.6768 212.22 0.0000 9.0230 16.8373 + 19 0.6748 211.60 0.0000 8.9614 16.7895 + 20 0.6729 211.01 0.0000 8.9031 16.7442 + 21 0.6712 210.46 0.0000 8.8491 16.7018 + 22 0.6695 209.94 0.0000 8.7974 16.6611 + 23 0.6679 209.43 0.0000 8.7474 16.6214 + 24 0.6663 208.92 0.0000 8.6975 16.5817 + 25 0.6647 208.42 0.0000 8.6478 16.5420 + 26 0.6631 207.94 0.0000 8.6005 16.5039 + 27 0.6617 207.50 0.0000 8.5570 16.4688 + 28 0.6603 207.05 0.0000 8.5136 16.4335 + 29 0.6590 206.64 0.0000 8.4736 16.4009 + 30 0.6577 206.24 0.0000 8.4338 16.3684 + + + ******************************************************************* + * ANNUAL HEATING, COOLING AND/OR ELECTRICITY PRODUCTION PROFILE * + ******************************************************************* + YEAR ELECTRICITY HEAT RESERVOIR PERCENTAGE OF + PROVIDED EXTRACTED HEAT CONTENT TOTAL HEAT MINED + (GWh/year) (GWh/year) (10^15 J) (%) + 1 133.2 673.7 82.77 2.85 + 2 91.1 491.1 81.01 4.92 + 3 86.7 476.0 79.29 6.93 + 4 83.8 466.3 77.61 8.90 + 5 81.8 459.2 75.96 10.84 + 6 80.2 453.6 74.33 12.76 + 7 78.8 449.0 72.71 14.66 + 8 77.7 445.1 71.11 16.54 + 9 76.7 441.7 69.52 18.40 + 10 75.8 438.7 67.94 20.26 + 11 75.0 436.0 66.37 22.10 + 12 74.3 433.5 64.81 23.93 + 13 73.6 431.1 63.26 25.75 + 14 73.0 429.0 61.71 27.57 + 15 72.4 427.0 60.18 29.37 + 16 71.9 425.1 58.65 31.17 + 17 71.4 423.4 57.12 32.96 + 18 70.9 421.7 55.60 34.74 + 19 70.4 420.0 54.09 36.51 + 20 70.0 418.5 52.58 38.28 + 21 69.6 417.0 51.08 40.04 + 22 69.2 415.6 49.59 41.80 + 23 68.8 414.2 48.10 43.55 + 24 68.4 412.8 46.61 45.29 + 25 68.0 411.5 45.13 47.03 + 26 67.6 410.2 43.65 48.77 + 27 67.3 409.0 42.18 50.49 + 28 67.0 407.9 40.71 52.22 + 29 66.6 406.8 39.25 53.94 + 30 49.8 304.4 38.15 55.22 + + + ******************************** + * REVENUE & CASHFLOW PROFILE * + ******************************** +Year Electricity | Heat | Cooling | Carbon | Project +Since Price Ann. Rev. Cumm. Rev. | Price Ann. Rev. Cumm. Rev. | Price Ann. Rev. Cumm. Rev. | Price Ann. Rev. Cumm. Rev. | OPEX Net Rev. Net Cashflow +Start (cents/kWh)(MUSD/yr) (MUSD) |(cents/kWh) (MUSD/yr) (MUSD) |(cents/kWh) (MUSD/yr) (MUSD) |(USD/tonne) (MUSD/yr) (MUSD) |(MUSD/yr) (MUSD/yr) (MUSD) +________________________________________________________________________________________________________________________________________________________________________________________ + 1 0.00 -158.59 0.00 | 0.00 0.00 0.00 | 0.00 0.00 0.00 | 0.00 0.00 0.00 | 0.00 -158.59 -158.59 + 2 19.00 21.93 25.31 | 2.50 0.00 0.00 | 2.50 0.00 0.00 | 0.00 0.00 0.00 | 3.37 21.93 -136.66 + 3 19.00 13.93 42.61 | 2.50 0.00 0.00 | 2.50 0.00 0.00 | 0.00 0.00 0.00 | 3.37 13.93 -122.73 + 4 19.00 13.09 59.07 | 2.50 0.00 0.00 | 2.50 0.00 0.00 | 0.00 0.00 0.00 | 3.37 13.09 -109.64 + 5 19.00 12.56 75.00 | 2.50 0.00 0.00 | 2.50 0.00 0.00 | 0.00 0.00 0.00 | 3.37 12.56 -97.08 + 6 19.00 12.16 90.54 | 2.50 0.00 0.00 | 2.50 0.00 0.00 | 0.00 0.00 0.00 | 3.37 12.16 -84.92 + 7 19.00 11.86 105.77 | 2.50 0.00 0.00 | 2.50 0.00 0.00 | 0.00 0.00 0.00 | 3.37 11.86 -73.06 + 8 19.00 11.60 120.74 | 2.50 0.00 0.00 | 2.50 0.00 0.00 | 0.00 0.00 0.00 | 3.37 11.60 -61.46 + 9 19.00 11.38 135.50 | 2.50 0.00 0.00 | 2.50 0.00 0.00 | 0.00 0.00 0.00 | 3.37 11.38 -50.08 + 10 19.00 11.20 150.07 | 2.50 0.00 0.00 | 2.50 0.00 0.00 | 0.00 0.00 0.00 | 3.37 11.20 -38.89 + 11 19.00 11.03 164.47 | 2.50 0.00 0.00 | 2.50 0.00 0.00 | 0.00 0.00 0.00 | 3.37 11.03 -27.85 + 12 19.00 10.88 178.73 | 2.50 0.00 0.00 | 2.50 0.00 0.00 | 0.00 0.00 0.00 | 3.37 10.88 -16.97 + 13 19.00 10.75 192.85 | 2.50 0.00 0.00 | 2.50 0.00 0.00 | 0.00 0.00 0.00 | 3.37 10.75 -6.23 + 14 19.00 10.61 206.83 | 2.50 0.00 0.00 | 2.50 0.00 0.00 | 0.00 0.00 0.00 | 3.37 10.61 4.39 + 15 19.00 10.50 220.70 | 2.50 0.00 0.00 | 2.50 0.00 0.00 | 0.00 0.00 0.00 | 3.37 10.50 14.88 + 16 19.00 10.39 234.46 | 2.50 0.00 0.00 | 2.50 0.00 0.00 | 0.00 0.00 0.00 | 3.37 10.39 25.27 + 17 19.00 10.28 248.12 | 2.50 0.00 0.00 | 2.50 0.00 0.00 | 0.00 0.00 0.00 | 3.37 10.28 35.56 + 18 19.00 10.19 261.69 | 2.50 0.00 0.00 | 2.50 0.00 0.00 | 0.00 0.00 0.00 | 3.37 10.19 45.75 + 19 19.00 10.10 275.16 | 2.50 0.00 0.00 | 2.50 0.00 0.00 | 0.00 0.00 0.00 | 3.37 10.10 55.84 + 20 19.00 10.01 288.53 | 2.50 0.00 0.00 | 2.50 0.00 0.00 | 0.00 0.00 0.00 | 3.37 10.01 65.85 + 21 19.00 9.92 301.83 | 2.50 0.00 0.00 | 2.50 0.00 0.00 | 0.00 0.00 0.00 | 3.37 9.92 75.77 + 22 19.00 9.84 315.05 | 2.50 0.00 0.00 | 2.50 0.00 0.00 | 0.00 0.00 0.00 | 3.37 9.84 85.61 + 23 19.00 9.77 328.19 | 2.50 0.00 0.00 | 2.50 0.00 0.00 | 0.00 0.00 0.00 | 3.37 9.77 95.38 + 24 19.00 9.69 341.25 | 2.50 0.00 0.00 | 2.50 0.00 0.00 | 0.00 0.00 0.00 | 3.37 9.69 105.07 + 25 19.00 9.62 354.24 | 2.50 0.00 0.00 | 2.50 0.00 0.00 | 0.00 0.00 0.00 | 3.37 9.62 114.69 + 26 19.00 9.54 367.16 | 2.50 0.00 0.00 | 2.50 0.00 0.00 | 0.00 0.00 0.00 | 3.37 9.54 124.23 + 27 19.00 9.48 380.01 | 2.50 0.00 0.00 | 2.50 0.00 0.00 | 0.00 0.00 0.00 | 3.37 9.48 133.71 + 28 19.00 9.41 392.80 | 2.50 0.00 0.00 | 2.50 0.00 0.00 | 0.00 0.00 0.00 | 3.37 9.41 143.12 + 29 19.00 9.35 405.52 | 2.50 0.00 0.00 | 2.50 0.00 0.00 | 0.00 0.00 0.00 | 3.37 9.35 152.47 + 30 19.00 9.29 418.18 | 2.50 0.00 0.00 | 2.50 0.00 0.00 | 0.00 0.00 0.00 | 3.37 9.29 161.76 diff --git a/tests/examples/example_SBT_Hi_T.txt b/tests/examples/example_SBT_Hi_T.txt new file mode 100644 index 00000000..0aaf7e92 --- /dev/null +++ b/tests/examples/example_SBT_Hi_T.txt @@ -0,0 +1,48 @@ +************reservoir************* +Reservoir Model, 8 +Reservoir Depth, 7.5 kilometer +Gradient 1, 60 +Reservoir Volume Option, 4 +Reservoir Volume, 8136407202.64, pi * 1 mi * 1 mi * 1000 +Reservoir Heat Capacity, 1112 +Reservoir Density, 2663 +Reservoir Thermal Conductivity, 2.5 + +Lateral Endpoint Depth, 7.5 kilometer +Lateral Inclination Angle, 30.5 +Junction Depth, 4.2 kilometer +Vertical Section Length, 4.2 kilometer +Number of Multilateral Sections, 12 +SBT Accuracy Desired, 5 +Lateral Spacing, 75 +Discretization Length, 250 +SBT Initial Timestep Count, 5 +SBT Initial to Final Timestep Transition, 10000, seconds = one year +SBT Final Timestep Count, 120, ~4 per year + +****************Wellbore*********** +Is AGS, True +Well Geometry Configuration, 5 +Number of Production Wells, 1 +Number of Injection Wells, 1 +Production Well Diameter, 8.5 +Injection Well Diameter, 8.5 +Nonvertical Wellbore Diameter, 0.216 +Production Flow Rate per Well, 80 +Reservoir Impedance, 1E-4, there is very little impedance in an AGS system +Multilaterals Cased, False + +****************Surface plant************ +End-Use Option, 1 +Power Plant Type, 2 +Plant Lifetime, 30 +Ambient Temperature, 10 +Surface Temperature, 10 +Injection Temperature, 60 + +***************Economics**************** +Economic Model, 3 +Reservoir Stimulation Capital Cost, 0 +Exploration Capital Cost, 0 +Starting Electricity Sale Price, 0.19 +Ending Electricity Sale Price, 0.19 diff --git a/tests/examples/example_SBT_Lo_T.out b/tests/examples/example_SBT_Lo_T.out new file mode 100644 index 00000000..1390f1bd --- /dev/null +++ b/tests/examples/example_SBT_Lo_T.out @@ -0,0 +1,224 @@ + ***************** + ***CASE REPORT*** + ***************** + +Simulation Metadata +---------------------- + GEOPHIRES Version: 3.5.2 + Simulation Date: 2024-09-02 + Simulation Time: 14:35 + Calculation Time: 45.642 sec + + ***SUMMARY OF RESULTS*** + + End-Use Option: Electricity + Average Net Electricity Production: 0.00 MW + Electricity breakeven price: 18339.44 cents/kWh + Number of production wells: 1 + Number of injection wells: 1 + Flowrate per production well: 6.0 kg/sec + Well depth (or total length, if not vertical): 2.5 kilometer + Geothermal gradient: 31.25 degC/km + + + ***ECONOMIC PARAMETERS*** + + Economic Model = BICYCLE + Accrued financing during construction: 0.00 + Project lifetime: 30 yr + Capacity factor: 90.0 % + Project NPV: -47.18 MUSD + Project IRR: 0.00 % + Project VIR=PI=PIR: -0.22 + Project MOIC: -1.00 + Project Payback Period: N/A + Estimated Jobs Created: 0 + + ***ENGINEERING PARAMETERS*** + + Number of Production Wells: 1 + Number of Injection Wells: 1 + Well depth (or total length, if not vertical): 2.5 kilometer + Water loss rate: 0.0 + Pump efficiency: 75.0 + Injection temperature: 24.0 degC + Production Wellbore heat transmission calculated with Ramey's model + Average production well temperature drop: 0.0 degC + Flowrate per production well: 6.0 kg/sec + Injection well casing ID: 8.500 in + Production well casing ID: 8.500 in + Number of times redrilling: 0 + Power plant type: Supercritical ORC + + + ***RESOURCE CHARACTERISTICS*** + + Maximum reservoir temperature: 400.0 degC + Number of segments: 1 + Geothermal gradient: 31.25 degC/km + + + ***RESERVOIR PARAMETERS*** + +The AGS models contain an intrinsic reservoir model that doesn't expose values that can be used in extensive reporting. + + + ***RESERVOIR SIMULATION RESULTS*** + + Maximum Production Temperature: 59.0 degC + Average Production Temperature: 58.9 degC + Minimum Production Temperature: 58.4 degC + Initial Production Temperature: 59.0 degC +The AGS models contain an intrinsic reservoir model that doesn't expose values that can be used in extensive reporting. + + + ***CAPITAL COSTS (M$)*** + + Drilling and completion costs: 37.58 MUSD + Drilling and completion costs per vertical production well: 3.28 MUSD + Drilling and completion costs per vertical injection well: 3.28 MUSD + Drilling and completion costs per non-vertical sections: 31.02 MUSD + Stimulation costs: 0.00 MUSD + Surface power plant costs: 0.03 MUSD + Field gathering system costs: 0.97 MUSD + Total surface equipment costs: 0.99 MUSD + Exploration costs: 0.00 MUSD + Total capital costs: 38.57 MUSD + + + ***OPERATING AND MAINTENANCE COSTS (M$/yr)*** + + Wellfield maintenance costs: 0.45 MUSD/yr + Power plant maintenance costs: 0.20 MUSD/yr + Water costs: 0.00 MUSD/yr + Total operating and maintenance costs: 0.65 MUSD/yr + + + ***SURFACE EQUIPMENT SIMULATION RESULTS*** + + Initial geofluid availability: 0.02 MW/(kg/s) + Maximum Total Electricity Generation: 0.00 MW + Average Total Electricity Generation: 0.00 MW + Minimum Total Electricity Generation: 0.00 MW + Initial Total Electricity Generation: 0.00 MW + Maximum Net Electricity Generation: 0.00 MW + Average Net Electricity Generation: 0.00 MW + Minimum Net Electricity Generation: 0.00 MW + Initial Net Electricity Generation: 0.00 MW + Average Annual Total Electricity Generation: 0.02 GWh + Average Annual Net Electricity Generation: 0.02 GWh + Average Pumping Power: 0.00 MW + + ************************************************************ + * HEATING, COOLING AND/OR ELECTRICITY PRODUCTION PROFILE * + ************************************************************ + YEAR THERMAL GEOFLUID PUMP NET FIRST LAW + DRAWDOWN TEMPERATURE POWER POWER EFFICIENCY + (degC) (MW) (MW) (%) + 1 1.0000 58.97 0.0000 0.0024 0.2754 + 2 0.9908 58.44 0.0000 0.0019 0.2282 + 3 0.9960 58.74 0.0000 0.0022 0.2546 + 4 0.9979 58.85 0.0000 0.0023 0.2647 + 5 0.9989 58.91 0.0000 0.0023 0.2698 + 6 0.9994 58.94 0.0000 0.0024 0.2725 + 7 0.9998 58.96 0.0000 0.0024 0.2742 + 8 0.9999 58.97 0.0000 0.0024 0.2749 + 9 1.0000 58.97 0.0000 0.0024 0.2752 + 10 1.0000 58.97 0.0000 0.0024 0.2754 + 11 1.0000 58.97 0.0000 0.0024 0.2754 + 12 1.0000 58.97 0.0000 0.0024 0.2752 + 13 0.9999 58.97 0.0000 0.0024 0.2749 + 14 0.9998 58.96 0.0000 0.0024 0.2743 + 15 0.9996 58.95 0.0000 0.0024 0.2732 + 16 0.9995 58.94 0.0000 0.0024 0.2726 + 17 0.9993 58.94 0.0000 0.0024 0.2720 + 18 0.9992 58.93 0.0000 0.0024 0.2713 + 19 0.9991 58.92 0.0000 0.0023 0.2706 + 20 0.9989 58.91 0.0000 0.0023 0.2698 + 21 0.9988 58.90 0.0000 0.0023 0.2690 + 22 0.9986 58.89 0.0000 0.0023 0.2681 + 23 0.9984 58.88 0.0000 0.0023 0.2673 + 24 0.9982 58.87 0.0000 0.0023 0.2664 + 25 0.9981 58.86 0.0000 0.0023 0.2654 + 26 0.9979 58.85 0.0000 0.0023 0.2645 + 27 0.9977 58.84 0.0000 0.0023 0.2635 + 28 0.9975 58.83 0.0000 0.0023 0.2625 + 29 0.9973 58.82 0.0000 0.0023 0.2615 + 30 0.9971 58.81 0.0000 0.0022 0.2606 + + + ******************************************************************* + * ANNUAL HEATING, COOLING AND/OR ELECTRICITY PRODUCTION PROFILE * + ******************************************************************* + YEAR ELECTRICITY HEAT RESERVOIR PERCENTAGE OF + PROVIDED EXTRACTED HEAT CONTENT TOTAL HEAT MINED + (GWh/year) (GWh/year) (10^15 J) (%) + 1 0.0 6.8 0.86 2.79 + 2 0.0 6.8 0.83 5.55 + 3 0.0 6.8 0.81 8.33 + 4 0.0 6.8 0.78 11.12 + 5 0.0 6.8 0.76 13.91 + 6 0.0 6.8 0.73 16.70 + 7 0.0 6.8 0.71 19.50 + 8 0.0 6.8 0.68 22.29 + 9 0.0 6.8 0.66 25.08 + 10 0.0 6.8 0.64 27.88 + 11 0.0 6.8 0.61 30.67 + 12 0.0 6.8 0.59 33.47 + 13 0.0 6.8 0.56 36.26 + 14 0.0 6.8 0.54 39.05 + 15 0.0 6.8 0.51 41.84 + 16 0.0 6.8 0.49 44.64 + 17 0.0 6.8 0.46 47.43 + 18 0.0 6.8 0.44 50.22 + 19 0.0 6.8 0.41 53.01 + 20 0.0 6.8 0.39 55.79 + 21 0.0 6.8 0.37 58.58 + 22 0.0 6.8 0.34 61.37 + 23 0.0 6.8 0.32 64.16 + 24 0.0 6.8 0.29 66.94 + 25 0.0 6.8 0.27 69.73 + 26 0.0 6.8 0.24 72.51 + 27 0.0 6.8 0.22 75.29 + 28 0.0 6.8 0.19 78.07 + 29 0.0 6.8 0.17 80.85 + 30 0.0 5.1 0.15 82.94 + + + ******************************** + * REVENUE & CASHFLOW PROFILE * + ******************************** +Year Electricity | Heat | Cooling | Carbon | Project +Since Price Ann. Rev. Cumm. Rev. | Price Ann. Rev. Cumm. Rev. | Price Ann. Rev. Cumm. Rev. | Price Ann. Rev. Cumm. Rev. | OPEX Net Rev. Net Cashflow +Start (cents/kWh)(MUSD/yr) (MUSD) |(cents/kWh) (MUSD/yr) (MUSD) |(cents/kWh) (MUSD/yr) (MUSD) |(USD/tonne) (MUSD/yr) (MUSD) |(MUSD/yr) (MUSD/yr) (MUSD) +________________________________________________________________________________________________________________________________________________________________________________________ + 1 0.00 -38.57 0.00 | 0.00 0.00 0.00 | 0.00 0.00 0.00 | 0.00 0.00 0.00 | 0.00 -38.57 -38.57 + 2 19.00 -0.64 0.00 | 2.50 0.00 0.00 | 2.50 0.00 0.00 | 0.00 0.00 0.00 | 0.65 -0.64 -39.21 + 3 19.00 -0.64 0.01 | 2.50 0.00 0.00 | 2.50 0.00 0.00 | 0.00 0.00 0.00 | 0.65 -0.64 -39.86 + 4 19.00 -0.64 0.01 | 2.50 0.00 0.00 | 2.50 0.00 0.00 | 0.00 0.00 0.00 | 0.65 -0.64 -40.50 + 5 19.00 -0.64 0.01 | 2.50 0.00 0.00 | 2.50 0.00 0.00 | 0.00 0.00 0.00 | 0.65 -0.64 -41.14 + 6 19.00 -0.64 0.02 | 2.50 0.00 0.00 | 2.50 0.00 0.00 | 0.00 0.00 0.00 | 0.65 -0.64 -41.78 + 7 19.00 -0.64 0.02 | 2.50 0.00 0.00 | 2.50 0.00 0.00 | 0.00 0.00 0.00 | 0.65 -0.64 -42.42 + 8 19.00 -0.64 0.02 | 2.50 0.00 0.00 | 2.50 0.00 0.00 | 0.00 0.00 0.00 | 0.65 -0.64 -43.07 + 9 19.00 -0.64 0.03 | 2.50 0.00 0.00 | 2.50 0.00 0.00 | 0.00 0.00 0.00 | 0.65 -0.64 -43.71 + 10 19.00 -0.64 0.03 | 2.50 0.00 0.00 | 2.50 0.00 0.00 | 0.00 0.00 0.00 | 0.65 -0.64 -44.35 + 11 19.00 -0.64 0.03 | 2.50 0.00 0.00 | 2.50 0.00 0.00 | 0.00 0.00 0.00 | 0.65 -0.64 -44.99 + 12 19.00 -0.64 0.04 | 2.50 0.00 0.00 | 2.50 0.00 0.00 | 0.00 0.00 0.00 | 0.65 -0.64 -45.63 + 13 19.00 -0.64 0.04 | 2.50 0.00 0.00 | 2.50 0.00 0.00 | 0.00 0.00 0.00 | 0.65 -0.64 -46.28 + 14 19.00 -0.64 0.05 | 2.50 0.00 0.00 | 2.50 0.00 0.00 | 0.00 0.00 0.00 | 0.65 -0.64 -46.92 + 15 19.00 -0.64 0.05 | 2.50 0.00 0.00 | 2.50 0.00 0.00 | 0.00 0.00 0.00 | 0.65 -0.64 -47.56 + 16 19.00 -0.64 0.05 | 2.50 0.00 0.00 | 2.50 0.00 0.00 | 0.00 0.00 0.00 | 0.65 -0.64 -48.20 + 17 19.00 -0.64 0.06 | 2.50 0.00 0.00 | 2.50 0.00 0.00 | 0.00 0.00 0.00 | 0.65 -0.64 -48.84 + 18 19.00 -0.64 0.06 | 2.50 0.00 0.00 | 2.50 0.00 0.00 | 0.00 0.00 0.00 | 0.65 -0.64 -49.48 + 19 19.00 -0.64 0.06 | 2.50 0.00 0.00 | 2.50 0.00 0.00 | 0.00 0.00 0.00 | 0.65 -0.64 -50.13 + 20 19.00 -0.64 0.07 | 2.50 0.00 0.00 | 2.50 0.00 0.00 | 0.00 0.00 0.00 | 0.65 -0.64 -50.77 + 21 19.00 -0.64 0.07 | 2.50 0.00 0.00 | 2.50 0.00 0.00 | 0.00 0.00 0.00 | 0.65 -0.64 -51.41 + 22 19.00 -0.64 0.07 | 2.50 0.00 0.00 | 2.50 0.00 0.00 | 0.00 0.00 0.00 | 0.65 -0.64 -52.05 + 23 19.00 -0.64 0.08 | 2.50 0.00 0.00 | 2.50 0.00 0.00 | 0.00 0.00 0.00 | 0.65 -0.64 -52.69 + 24 19.00 -0.64 0.08 | 2.50 0.00 0.00 | 2.50 0.00 0.00 | 0.00 0.00 0.00 | 0.65 -0.64 -53.34 + 25 19.00 -0.64 0.08 | 2.50 0.00 0.00 | 2.50 0.00 0.00 | 0.00 0.00 0.00 | 0.65 -0.64 -53.98 + 26 19.00 -0.64 0.09 | 2.50 0.00 0.00 | 2.50 0.00 0.00 | 0.00 0.00 0.00 | 0.65 -0.64 -54.62 + 27 19.00 -0.64 0.09 | 2.50 0.00 0.00 | 2.50 0.00 0.00 | 0.00 0.00 0.00 | 0.65 -0.64 -55.26 + 28 19.00 -0.64 0.09 | 2.50 0.00 0.00 | 2.50 0.00 0.00 | 0.00 0.00 0.00 | 0.65 -0.64 -55.90 + 29 19.00 -0.64 0.10 | 2.50 0.00 0.00 | 2.50 0.00 0.00 | 0.00 0.00 0.00 | 0.65 -0.64 -56.55 + 30 19.00 -0.64 0.10 | 2.50 0.00 0.00 | 2.50 0.00 0.00 | 0.00 0.00 0.00 | 0.65 -0.64 -57.19 diff --git a/tests/examples/example_SBT_Lo_T.txt b/tests/examples/example_SBT_Lo_T.txt new file mode 100644 index 00000000..949ff1c5 --- /dev/null +++ b/tests/examples/example_SBT_Lo_T.txt @@ -0,0 +1,50 @@ +************reservoir************* +Reservoir Model, 8 +Reservoir Depth, 2.4 kilometer +Gradient 1, 31.25 +Reservoir Volume Option, 4 +Reservoir Volume, 8136407202.64, pi * 1 mi * 1 mi * 1000 +Reservoir Heat Capacity, 1112 +Reservoir Density, 2663 +Reservoir Thermal Conductivity, 2.25 + +Lateral Endpoint Depth, 2.5 kilometer +Lateral Inclination Angle, 89 +Junction Depth, 2.4 kilometer +Vertical Section Length, 2.4 kilometer +Number of Multilateral Sections, 2 +SBT Accuracy Desired, 5 +Lateral Spacing, 75 +Discretization Length, 250 +SBT Initial Timestep Count, 5 +SBT Initial to Final Timestep Transition, 10000, seconds = one year +SBT Final Timestep Count, 120, ~4 per year + +****************Wellbore*********** +Is AGS, True +Well Geometry Configuration, 5 +Number of Production Wells, 1 +Number of Injection Wells, 1 +Production Well Diameter, 8.5 +Injection Well Diameter, 8.5 +Nonvertical Wellbore Diameter, 0.216 +Production Flow Rate per Well, 6 +Reservoir Impedance, 1E-4, there is very little impedance in an AGS system +Multilaterals Cased, False + +****************Surface plant************ +End-Use Option, 1 +Power Plant Type, 2 +Plant Lifetime, 30 +Ambient Temperature, 3 +Surface Temperature, 3 +Injection Temperature, 24 + +***************Economics**************** +Economic Model, 3 +Reservoir Stimulation Capital Cost, 0 +Exploration Capital Cost, 0 +Starting Electricity Sale Price, 0.19 +Ending Electricity Sale Price, 0.19 + +SBT Generate Wireframe Graphics, False