diff --git a/GHEtool/VariableClasses/BaseClass.py b/GHEtool/VariableClasses/BaseClass.py index 7de5b9e0..b27600a0 100644 --- a/GHEtool/VariableClasses/BaseClass.py +++ b/GHEtool/VariableClasses/BaseClass.py @@ -16,146 +16,11 @@ class BaseClass: """ This class contains basic functionality of different classes within GHEtool. - It contains the code to generate a dictionary from the class (in order to be able to export to JSON), - to load a class based on a dictionary and to check whether or not all attributes differ from None. This class should only be altered whenever a highly general method should be implemented. """ - def to_dict(self) -> dict: - """ - This function converts the class variables to a dictionary so it can be saved in a JSON format. - Currently, it can handle np.ndarray, list, set, str, int, float, tuple, - pygfunction.Borehole and classes within GHEtool. - - Returns - ------- - dict - Dictionary with all the attributes of the class - """ - - # get all variables in class - if hasattr(self, "__slots__"): - variables: List[str] = list(self.__slots__) - else: - variables: List[str] = [attr for attr in dir(self) if not callable(getattr(self, attr)) and not attr.startswith("__")] - - # initiate dictionary - dictionary: dict = {"__module__": f"{type(self).__module__}", "__name__": f"{type(self).__name__}"} - - # populate dictionary - for key in variables: - value = getattr(self, key) - - if value is None: - dictionary[key] = "None" - continue - - if isinstance(value, (int, bool, float, str)): - dictionary[key] = value - continue - - if isinstance(value, tuple): - dictionary[key] = {"value": list(value), "type": "tuple"} - continue - - if isinstance(value, np.ndarray): - dictionary[key] = {"value": value.tolist(), "type": "np.ndarray"} - continue - - if isinstance(value, set): - dictionary[key] = {"value": list(value), "type": "set"} - continue - - if isinstance(value, list): - if np.any(value) and isinstance(value[0], Borehole): - # pygfunction object - dictionary[key] = {"value": [{key: value for key, value in borehole.__dict__.items() - if key != "_is_tilted"} - for borehole in value], - "type": "pygfunction.Borehole"} - continue - dictionary[key] = value - - if isinstance(value, dict): - # note that this can cause problems whenever self.key contains values that or not int, bool, float or str - dictionary[key] = value - continue - - # for all self-defined classes - if hasattr(value, "to_dict"): - dictionary[key] = value.to_dict() - - return dictionary - - def from_dict(self, dictionary: dict) -> None: - """ - This function converts the dictionary values to the class attributes. - Currently, it can handle np.ndarray, list, set, str, int, float, tuple, pygfunction.Borehole - and classes within GHEtool. - - Parameters - ---------- - dictionary - Dictionary with all the attributes of the class - - Returns - ------- - None - """ - if hasattr(self, "__slots__"): - variables: List[str] = list(self.__slots__) - else: - variables: List[str] = [key for key, value in dictionary.items() if not key.startswith("__")] - - for key, value in dictionary.items(): - - if key not in variables: - continue - - if value == "None": - setattr(self, key, None) - continue - - if isinstance(value, (int, bool, float, str, list)): - setattr(self, key, value) - continue - - if isinstance(value, dict): - # note that this can mean that the value is a dictionary, or it is a np.ndarray, set or tuple - keys = value.keys() - # for all self-defined classes - if "__module__" in keys: - class_dict = getattr(import_module(value["__module__"]), value["__name__"]) - setattr(self, key, class_dict.__new__(class_dict)) - getattr(self, key).from_dict(value) - continue - - if len(keys) == 2 and "type" in keys and "value" in keys: - var_type = value["type"] - _value = value["value"] - if var_type == "set": - setattr(self, key, set(_value)) - continue - if var_type == "np.ndarray": - setattr(self, key, np.array(_value)) - continue - if var_type == "tuple": - setattr(self, key, tuple(_value)) - continue - if var_type == "pygfunction.Borehole": - borefield = [Borehole(H=borehole["H"], - D=borehole["D"], - r_b=borehole["r_b"], - x=borehole["x"], - y=borehole["y"], - tilt=borehole["tilt"], - orientation=borehole["orientation"]) - for borehole in _value] - setattr(self, key, borefield) - continue - # normal dictionary - setattr(self, key, value) + __allow_none__ = [] def check_values(self) -> bool: """ @@ -172,11 +37,10 @@ def check_values(self) -> bool: else: variables: List[str] = [attr for attr in dir(self) if not callable(getattr(self, attr)) and not attr.startswith("__")] - temp = [getattr(self, var) is not None for var in variables] + temp = [getattr(self, var) is not None for var in variables if var not in self.__class__.__allow_none__] if all(temp): return True else: - # print(f'There is a problem with the {[var for idx, var in enumerate(variables) if not temp[idx]]} variables.') return False diff --git a/GHEtool/VariableClasses/Dynamic_borhole_model.py b/GHEtool/VariableClasses/Dynamic_borhole_model.py new file mode 100644 index 00000000..10ce6ccc --- /dev/null +++ b/GHEtool/VariableClasses/Dynamic_borhole_model.py @@ -0,0 +1,598 @@ +import warnings +from enum import auto, IntEnum +from math import exp, log, pi, sqrt + +import numpy as np +import matplotlib.pyplot as plt +from scipy.interpolate import interp1d +from scipy.linalg.lapack import dgtsv + +import pygfunction as gt + +print('test import done') + +class CellProps(IntEnum): + R_IN = 0 + R_CENTER = auto() + R_OUT = auto() + K = auto() + RHO_CP = auto() + TEMP = auto() + VOL = auto() + + +class DynamicsBH(object): + """ + X. Xu and Jeffrey D. Spitler. 2006. 'Modeling of Vertical Ground Loop Heat + Exchangers with Variable Convective Resistance and Thermal Mass of the + Fluid.' in Proceedings of the 10th International Conference on Thermal + Energy Storage-EcoStock. Pomona, NJ, May 31-June 2. + """ + + """ + def foo(self): + from GHEtool import PipeData + bar() + """ + + def __init__(self, time, gFunc, boreholes, alpha, ground_data, fluid_data, pipe_data, borefield): + + print('initialize Numerical model') + + + data_ground = None + data_fluid = None + data_pipe = None + self.resist_bh_effective = None + self.Rb_cst = None + self.x = None + self.u_tube = None + self.fluid_factor = None + self.rho_cp_grout = None + self.rho_cp_pipe = None + far_field_radius = None + self.num_soil_cells = None + + + from GHEtool.VariableClasses import GroundFluxTemperature, FluidData, PipeData, GroundData + + print('gFunc in RM', gFunc, len(gFunc)) + + self.boreholes = boreholes + self.ground_ghe = ground_data + self.fluid_ghe = fluid_data + self.pipes_ghe = pipe_data + self.borefield = borefield + self.pipes_gt = gt.pipes + + # borefield not ok + + # make function based on number of boreholes (1 or more) and half of the distance between boreholes + distance_between_boreholes = 20 #aanpassen!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! + far_field_radius = distance_between_boreholes/2 + self.num_soil_cells = int(far_field_radius/10*500) + + self.init_temp = self.ground_ghe.Tg + + self.fluid_factor = 2 # 2 by default, make function to calculate this + self.x = 1 # parameter to modify final time + self.u_tube = 1 # 1 for single U tube, 2 for dubble U tube (not possible yet) + self.rho_cp_grout = 3800000.0 # Sandbox + self.rho_cp_pipe = 1800000.0 # not aan te passen voor sandbox + + self.Rb_cst = 1 # 0 for false, 1 for true (rewrite such that get_Rb takes right Rb) + self.resist_bh_effective = 0.165 # get this value with get_Rb + + print('data is loaded') + + # "The one dimensional model has a fluid core, an equivalent convective + # resistance layer, a tube layer, a grout layer and is surrounded by the + # ground." + + # cell numbers + self.num_fluid_cells = 3 + self.num_conv_cells = 1 + self.num_pipe_cells = 4 + self.num_grout_cells = 27 + + self.num_cells = self.num_fluid_cells + self.num_conv_cells + self.num_pipe_cells + self.num_cells += self.num_grout_cells + self.num_soil_cells + + self.bh_wall_idx = self.num_fluid_cells + self.num_conv_cells + self.num_pipe_cells + self.num_grout_cells + + # Geometry and grid procedure + + # far-field radius is set to half of the distance between two boreholes or 10m when a single borehole is used; + # the soil region is represented by (500/10 * far-field radius) cells + + self.r_far_field = far_field_radius - self.boreholes[0].r_b + + # borehole radius is set to the actual radius of the borehole + self.r_borehole = self.boreholes[0].r_b + + # outer tube radius is set to sqrt(2) * r_p_o, tube region has 4 cells + self.r_out_tube = sqrt(2 * self.u_tube) * self.pipes_ghe.r_out + + # inner tube radius is set to r_out_tube-t_p + self.thickness_pipe = self.pipes_ghe.r_out - self.pipes_ghe.r_in + self.r_in_tube = self.r_out_tube - self.thickness_pipe + + # r_in_convection is set to r_in_tube - 1/4 * t + self.r_in_convection = self.r_in_tube - self.thickness_pipe / 4.0 + + # r_fluid is set to r_in_convection - 3/4 * t + self.r_fluid = self.r_in_convection - (3.0 / 4.0 * self.thickness_pipe) + + # Thicknesses of the grid regions + self.thickness_soil = (self.r_far_field - self.r_borehole) / self.num_soil_cells + self.thickness_grout = (self.r_borehole - self.r_out_tube) / self.num_grout_cells + # pipe thickness is equivalent to original tube thickness + self.thickness_conv = (self.r_in_tube - self.r_in_convection) / self.num_conv_cells + self.thickness_fluid = (self.r_in_convection - self.r_fluid) / self.num_fluid_cells + + + print('all parameters are set') + + # other + self.g = np.array([], dtype=np.double) + self.Tf = np.array([], dtype=np.double) + self.Tb = np.array([], dtype=np.double) + self.diff = np.array([], dtype=np.double) + self.g_bhw = np.array([], dtype=np.double) + self.lntts = np.array([], dtype=np.double) + self.c_0 = 2.0 * pi * self.ground_ghe.k_s() + soil_diffusivity = self.ground_ghe.k_s() / (self.ground_ghe.volumetric_heat_capacity()) # = alpha + self.t_s = self.boreholes[0].H ** 2 / (9 * self.ground_ghe.alpha()) + # default is at least 49 hours, or up to -8.6 log time + self.calc_time_in_sec = max([self.t_s * exp(-8.6), 49.0 * 3600.0]) + + print('ts charact time', self.t_s) + print('ts * exp ', self.t_s * exp(-8.6)) + + self.t_b = 5 * (self.boreholes[0].r_b) ** 2 / self.ground_ghe.alpha() + print('t_b time vanaf wanneer stationaire benadering mag', self.t_b) + self.final_time = self.x * self.t_b + #self.final_time = 3601 + print('final time', self.final_time) + self.g_sts = None + + def partial_init(self): + # TODO: unravel how to eliminate this. + # - It was calling the full class ctor "self.__init__()" which is just plain wrong... + # - Now we're calling a stripped down version with only the most essential + # variables which are required. + # - This is here partially because equivalent boreholes are generated. + + soil_diffusivity = self.ground_ghe.k_s() / self.ground_ghe.volumetric_heat_capacity() # kon k_s niet terugvinden dus soil.k genomen + self.t_s = self.boreholes[0].H ** 2 / ( + 9 * self.ground_ghe.alpha()) # self.t_s = single_u_tube.b.H ** 2 / (9 * soil_diffusivity) + self.calc_time_in_sec = max([self.t_s * exp(-8.6), 49.0 * 3600.0]) + + def fill_radial_cell(self, radial_cell, resist_p_eq, resist_f_eq, resist_tg_eq): + + num_fluid_cells = self.num_fluid_cells + num_conv_cells = self.num_conv_cells + num_pipe_cells = self.num_pipe_cells + num_grout_cells = self.num_grout_cells + num_soil_cells = self.num_soil_cells + + cell_summation = 0 + + # load fluid cells + for idx in range(cell_summation, num_fluid_cells + cell_summation): + center_radius = self.r_fluid + idx * self.thickness_fluid + + if idx == 0: + inner_radius = center_radius + else: + inner_radius = center_radius - self.thickness_fluid / 2.0 + + outer_radius = center_radius + self.thickness_fluid / 2.0 + + # The equivalent thermal mass of the fluid can be calculated from + # equation (2) + # pi (r_in_conv ** 2 - r_f **2) C_eq_f = 2pi r_p_in**2 * C_f + rho_cp_eq = (self.fluid_factor * self.u_tube * 2.0 + * (self.pipes_ghe.r_in ** 2) + * self.fluid_ghe.Cp * self.fluid_ghe.rho + ) / ((self.r_in_convection ** 2) - (self.r_fluid ** 2)) + + k_eq = rho_cp_eq / self.fluid_ghe.Cp + + volume = pi * (outer_radius ** 2 - inner_radius ** 2) + radial_cell[:, idx] = np.array( + [ + inner_radius, + center_radius, + outer_radius, + k_eq, + rho_cp_eq, + self.init_temp, + volume, + ], + dtype=np.double, + ) + cell_summation += num_fluid_cells + + # TODO: verify whether errors are possible here and raise exception if needed + # assert cell_summation == num_fluid_cells + + # load convection cells + for idx in range(cell_summation, num_conv_cells + cell_summation): + j = idx - cell_summation + inner_radius = self.r_in_convection + j * self.thickness_conv + center_radius = inner_radius + self.thickness_conv / 2.0 + outer_radius = inner_radius + self.thickness_conv + k_eq = log(self.r_in_tube / self.r_in_convection) / (2.0 * pi * resist_f_eq) + rho_cp = 1.0 + volume = pi * (outer_radius ** 2 - inner_radius ** 2) + radial_cell[:, idx] = np.array( + [ + inner_radius, + center_radius, + outer_radius, + k_eq, + rho_cp, + self.init_temp, + volume, + ], + dtype=np.double, + ) + cell_summation += num_conv_cells + + # TODO: verify whether errors are possible here and raise exception if needed + # assert cell_summation == (num_fluid_cells + num_conv_cells) + + # load pipe cells + for idx in range(cell_summation, num_pipe_cells + cell_summation): + j = idx - cell_summation + inner_radius = self.r_in_tube + j * self.thickness_pipe + center_radius = inner_radius + self.thickness_pipe / 2.0 + outer_radius = inner_radius + self.thickness_pipe + conductivity = log(self.r_borehole / self.r_in_tube) / (2.0 * pi * resist_p_eq) + # rho_cp = self.single_u_tube.pipe.rhoCp + rho_cp = self.rho_cp_pipe + volume = pi * (outer_radius ** 2 - inner_radius ** 2) + radial_cell[:, idx] = np.array( + [ + inner_radius, + center_radius, + outer_radius, + conductivity, + rho_cp, + self.init_temp, + volume, + ], + dtype=np.double, + ) + cell_summation += num_pipe_cells + + # TODO: verify whether errors are possible here and raise exception if needed + # assert cell_summation == (num_fluid_cells + num_conv_cells + num_pipe_cells) + + # load grout cells + for idx in range(cell_summation, num_grout_cells + cell_summation): + j = idx - cell_summation + inner_radius = self.r_out_tube + j * self.thickness_grout + center_radius = inner_radius + self.thickness_grout / 2.0 + outer_radius = inner_radius + self.thickness_grout + conductivity = log(self.r_borehole / self.r_in_tube) / (2.0 * pi * resist_tg_eq) + # rho_cp = self.single_u_tube.grout.rhoCp + rho_cp = self.rho_cp_grout + volume = pi * (outer_radius ** 2 - inner_radius ** 2) + radial_cell[:, idx] = np.array( + [ + inner_radius, + center_radius, + outer_radius, + conductivity, + rho_cp, + self.init_temp, + volume, + ], + dtype=np.double, + ) + cell_summation += num_grout_cells + + # TODO: verify whether errors are possible here and raise exception if needed + # assert cell_summation == (num_fluid_cells + num_conv_cells + num_pipe_cells + num_grout_cells) + + # load soil cells + for idx in range(cell_summation, num_soil_cells + cell_summation): + j = idx - cell_summation + inner_radius = self.r_borehole + j * self.thickness_soil + center_radius = inner_radius + self.thickness_soil / 2.0 + outer_radius = inner_radius + self.thickness_soil + conductivity = self.ground_ghe.k_s() + rho_cp = self.ground_ghe.volumetric_heat_capacity() + volume = pi * (outer_radius ** 2 - inner_radius ** 2) + radial_cell[:, idx] = np.array( + [ + inner_radius, + center_radius, + outer_radius, + conductivity, + rho_cp, + self.init_temp, + volume, + ], + dtype=np.double, + ) + cell_summation += num_soil_cells + + def calc_sts_g_functions(self, final_time=None) -> tuple: + print('in calc sts g functions') + + self.partial_init() + self.m_flow_borehole = self.fluid_ghe.mfr + final_time = self.final_time + + self.pipe_roughness = 1e-06 + self.bh = gt.boreholes + self.H = self.boreholes[0].H + self.r_b = self.r_borehole + self.D = self.boreholes[0].D + borehole = self.bh.Borehole(self.H, self.D, self.r_b, 0, 0) + self.h_f = self.pipes_gt.convective_heat_transfer_coefficient_circular_pipe(self.m_flow_borehole, + self.pipes_ghe.r_in, + self.fluid_ghe.mu, + self.fluid_ghe.rho, + self.fluid_ghe.k_f, + self.fluid_ghe.Cp, + self.pipe_roughness) + R_f = 1 / (self.h_f * 2 * pi * self.pipes_ghe.r_in) + R_p = self.pipes_gt.conduction_thermal_resistance_circular_pipe(self.pipes_ghe.r_in, self.pipes_ghe.r_out, + self.pipes_ghe.k_p) + + # wanneer Rb niet gekend is, maar dan moet je dit model niet gebruiken eigenlijk + if self.Rb_cst == 0: + + if self.u_tube == 2: + pipe = self.tubes_ghe.MultipleUTube(self.pipes_ghe.pos, self.pipes_ghe.r_in, self.pipes_ghe.r_out, + borehole, self.ground_ghe.k_s(), self.pipes_ghe.k_g, + R_p + R_f, 2) + + self.resist_bh_effective = pipe.effective_borehole_thermal_resistance(self.m_flow_borehole, self.fluid_ghe.Cp) + + elif self.u_tube == 1: + pipe = self.tubes_ghe.SingleUTube(self.pipes_ghe.pos, self.pipes_ghe.r_in, self.pipes_ghe.r_out, + borehole, self.ground_ghe.k_s(), self.pipes_ghe.k_g, + R_p + R_f, 2) + + self.resist_bh_effective = pipe.effective_borehole_thermal_resistance(self.m_flow_borehole, self.fluid_ghe.Cp) + + else: + #self.resist_bh_effective = self.ground_ghe.Rb + lone = 22 + + print('Rb* in NM',self.resist_bh_effective) + + resist_f_eq = R_f / 2 + resist_p_eq = R_p / 2 + resist_tg_eq = self.resist_bh_effective - resist_f_eq + + print('resist_f_eq', resist_f_eq) + print('resist_p_eq', resist_p_eq) + print('resist_tg_eq', resist_tg_eq) + + # Pass radial cell by reference and fill here so that it can be + # destroyed when this method returns + radial_cell = np.zeros(shape=(len(CellProps), self.num_cells), dtype=np.double) + self.fill_radial_cell(radial_cell, resist_p_eq, resist_f_eq, resist_tg_eq) + + self.t_b = 5 * (self.boreholes[0].r_b) ** 2 / self.ground_ghe.alpha() + + final_time = self.final_time + print('final time t_b in s', final_time) + + self.t_s = self.boreholes[0].H ** 2 / (9 * self.ground_ghe.alpha()) + print('final time lntts', log(final_time / self.t_s)) + + g = [] + g_bhw = [] + g_comb = [] + g_plot = [] + lntts = [] + plottime = [2e-12] + Tf = [self.init_temp] + Tb = [self.init_temp] + diff = [0] + qb=[] + + _dl = np.zeros(self.num_cells - 1) + _d = np.zeros(self.num_cells) + _du = np.zeros(self.num_cells - 1) + _b = np.zeros(self.num_cells) + + heat_flux = 1.0 + init_temp = self.init_temp + + time = 1e-12 - 120 + time_step = 120 + + _fe_1 = np.zeros(shape=(self.num_cells - 2), dtype=np.double) + _fe_2 = np.zeros_like(_fe_1) + _ae = np.zeros_like(_fe_2) + _fw_1 = np.zeros_like(_ae) + _fw_2 = np.zeros_like(_fw_1) + _aw = np.zeros_like(_fw_2) + _ad = np.zeros_like(_aw) + + _west_cell = radial_cell[:, 0: self.num_cells - 2] + _center_cell = radial_cell[:, 1: self.num_cells - 1] + _east_cell = radial_cell[:, 2: self.num_cells - 0] + + fe_1 = log(radial_cell[CellProps.R_OUT, 0] / radial_cell[CellProps.R_CENTER, 0]) + fe_1 /= (2.0 * pi * radial_cell[CellProps.K, 0]) + + fe_2 = log(radial_cell[CellProps.R_CENTER, 1] / radial_cell[CellProps.R_IN, 1]) + fe_2 /= (2.0 * pi * radial_cell[CellProps.K, 1]) + + ae = 1 / (fe_1 + fe_2) + ad = radial_cell[CellProps.RHO_CP, 0] * radial_cell[CellProps.VOL, 0] / time_step + _d[0] = -ae / ad - 1 + _du[0] = ae / ad + + def fill_f1(fx_1, cell): + fx_1[:] = np.log(cell[CellProps.R_OUT, :] / cell[CellProps.R_CENTER, :]) / (2.0 * pi * cell[CellProps.K, :]) + + def fill_f2(fx_2, cell): + fx_2[:] = np.log(cell[CellProps.R_CENTER, :] / cell[CellProps.R_IN, :]) / (2.0 * pi * cell[CellProps.K, :]) + + fill_f1(_fe_1, _center_cell) + fill_f2(_fe_2, _east_cell) + _ae[:] = 1.0 / (_fe_1 + _fe_2) + + fill_f1(_fw_1, _west_cell) + fill_f2(_fw_2, _center_cell) + _aw[:] = -1.0 / (_fw_1 + _fw_2) + + _ad[:] = (_center_cell[CellProps.RHO_CP, :] * _center_cell[CellProps.VOL, :] / time_step) + _dl[0: self.num_cells - 2] = -_aw / _ad + _d[1: self.num_cells - 1] = _aw / _ad - _ae / _ad - 1.0 + _du[1: self.num_cells - 1] = _ae / _ad + + while True: + + time += time_step + + # For the idx == 0 case: + + _b[0] = -radial_cell[CellProps.TEMP, 0] - heat_flux / ad + + # For the idx == n-1 case + + _dl[self.num_cells - 2] = 0.0 + _d[self.num_cells - 1] = 1.0 + _b[self.num_cells - 1] = radial_cell[CellProps.TEMP, self.num_cells - 1] + + # Now handle the 1 to n-2 cases with numpy slicing and vectorization + _b[1: self.num_cells - 1] = -radial_cell[CellProps.TEMP, 1: self.num_cells - 1] + + # Tri-diagonal matrix solver + # High level interface to LAPACK routine + # https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.lapack.dgtsv.html#scipy.linalg.lapack.dgtsv + dgtsv(_dl, _d, _du, _b, overwrite_b=1) # TODO: Do we really need lapack just to do a TDMA solution? + + radial_cell[CellProps.TEMP, :] = _b + + # compute standard g-functions + g.append(self.c_0 * ((radial_cell[CellProps.TEMP, 0] - init_temp) / heat_flux - self.resist_bh_effective)) + + g_comb.append(self.c_0 * ((radial_cell[CellProps.TEMP, 0] - radial_cell[CellProps.TEMP, self.bh_wall_idx]) / heat_flux - self.resist_bh_effective)) + + g_plot.append(self.c_0 * ((radial_cell[CellProps.TEMP, 0] - init_temp) / heat_flux)) + + + qb.append(1- (self.resist_bh_effective - (radial_cell[CellProps.TEMP, 0]-radial_cell[CellProps.TEMP, self.bh_wall_idx]))) + + print(qb) + + print('f1', radial_cell[CellProps.TEMP, 0], 'f2', radial_cell[CellProps.TEMP, 1], 'f3', radial_cell[CellProps.TEMP, 2], 'c1', radial_cell[CellProps.TEMP, 3], 'time', time) + #print('idx -1 ', radial_cell[CellProps.TEMP, self.bh_wall_idx - 1], 'idx bhw', radial_cell[CellProps.TEMP,self.bh_wall_idx], 'idx bhw + 1', radial_cell[CellProps.TEMP, self.bh_wall_idx + 1], 'time', time) + + T0 = radial_cell[CellProps.TEMP, 0] + TBH = radial_cell[CellProps.TEMP, self.bh_wall_idx] + d= T0-TBH + + Tf.append(T0) + Tb.append(TBH) + diff.append(d) + + # compute g-functions at bh wall + bh_wall_temp = radial_cell[CellProps.TEMP, self.bh_wall_idx] + + #g_bhw.append(self.c_0 * ((bh_wall_temp - init_temp) / heat_flux)) + + lntts.append(time) + + plottime.append(time) + + """ + if d >= resist_bh_effective: + print('Tf - Tb heeft Rb* bereikt') + break + """ + if time >= final_time - time_step: + print('in time if') + print(final_time) + break + + + # plot de short term and long term g function + plt.rc('font', size=9) + plt.rc('xtick', labelsize=9) + plt.rc('ytick', labelsize=9) + plt.rc('lines', lw=1.5, markersize=5.0) + plt.rc('savefig', dpi=500) + fig = plt.figure() + + ax = fig.add_subplot(111) + ax.set_xlabel(r'$tijd$ [s] ', fontsize=18) + #ax.set_ylabel(r'$|T_f - T_b|$', fontsize=18) + + plt.tight_layout() + + dummy = [] + for i in range(0,len(plottime)): + dummy.append(0.09) + + #ax.plot(plottime, Tf) + #ax.plot(plottime, Tb) + ax.plot(plottime, diff) + ax.plot(plottime, dummy) + + """ + #ax.legend([f'$T_f$', f'$T_b$'], fontsize=10) + ax.legend([ f'$|T_f - T_b|$', f'$R_b^*$'], fontsize=20) + + # plot de short term and long term g function + plt.rc('font', size=9) + plt.rc('xtick', labelsize=9) + plt.rc('ytick', labelsize=9) + plt.rc('lines', lw=1.5, markersize=5.0) + plt.rc('savefig', dpi=500) + fig = plt.figure() + + ax = fig.add_subplot(111) + ax.set_xlabel(r'$tijd$ [s] ', fontsize=18) + ax.set_ylabel(r'g-functie', fontsize=18) + plt.tight_layout() + ax.plot(lntts, g_plot) + """ + + # quickly chop down the total values to a more manageable set + num_intervals = int(self.x * 30) + g_tmp = interp1d(lntts, g) + + + uniform_lntts_vals = np.linspace(lntts[0], lntts[-1], num_intervals) + uniform_g_vals = g_tmp(uniform_lntts_vals) + + #g_bhw_tmp = interp1d(lntts, g_bhw) + g_comb_tmp = interp1d(lntts,g_comb) + qb_tmp =interp1d(lntts,qb) + g_plot_tmp = interp1d(lntts, g_plot) + + #uniform_g_bhw_vals = g_bhw_tmp(uniform_lntts_vals) + uniform_g_comb_vals = g_comb_tmp(uniform_lntts_vals) + uniform_qb_vals = qb_tmp(uniform_lntts_vals) + uniform_g_plot_vals = g_plot_tmp(uniform_lntts_vals) + + # set the final arrays and interpolator objects + self.lntts = np.array(uniform_lntts_vals) + self.g = np.array(uniform_g_vals) + #self.g_bhw = np.array(uniform_g_bhw_vals) + self.g_comb = np.array(uniform_g_comb_vals) + self.qb = np.array(uniform_qb_vals) + self.g_sts = interp1d(self.lntts, self.g) + self.g_plot = np.array(uniform_g_plot_vals) + + print('lntts', self.lntts) + print('g', self.g) + print('g_bhw', self.g_bhw) + print('g_sts', self.g_sts) + print('g_comb',self.g_comb) + print('qb', self.qb) + + return self.lntts, self.g, self.g_sts, self.g_comb, self.g_plot, self.qb \ No newline at end of file diff --git a/GHEtool/VariableClasses/FluidData.py b/GHEtool/VariableClasses/FluidData.py index e7756816..80431ca6 100644 --- a/GHEtool/VariableClasses/FluidData.py +++ b/GHEtool/VariableClasses/FluidData.py @@ -13,13 +13,15 @@ class FluidData(BaseClass): Contains information regarding the fluid data of the borefield. """ - __slots__ = 'k_f', 'rho', 'Cp', 'mu', 'mfr' + __slots__ = 'k_f', 'rho', 'Cp', 'mu', '_mfr', '_vfr' + __allow_none__ = ['_vfr', '_mfr'] def __init__(self, mfr: float = None, k_f: float = None, rho: float = None, Cp: float = None, - mu: float = None): + mu: float = None, + vfr: float = None): """ Parameters @@ -34,12 +36,82 @@ def __init__(self, mfr: float = None, Thermal capacity of the fluid [J/kgK] mu : float Dynamic viscosity of the fluid [Pa/s] + vfr : float + Volume flow rate per borehole [l/s] """ self.k_f: float | None = k_f # Thermal conductivity W/mK - self.mfr: float | None = mfr # Mass flow rate per borehole kg/s + self._mfr: float | None = mfr # Mass flow rate per borehole kg/s self.rho: float | None = rho # Density kg/m3 self.Cp: float | None = Cp # Thermal capacity J/kgK self.mu: float | None = mu # Dynamic viscosity Pa/s + self._vfr: float | None = vfr # Volume flow rate l/s + + if self._mfr is not None and self._vfr is not None: + raise ValueError('You cannot set both the mass flow rate and volume flow rate') + + + @property + def vfr(self) -> float: + """ + This function returns the volume flow rate. + + Returns + ------- + float + volume flow rate [l/s] + """ + return self._vfr + + @vfr.setter + def vfr(self, vfr: float) -> None: + """ + This function sets the volume flow rate. + The mass flow rate will be set to 0 + + Parameters + ---------- + vfr : float + Volume flow rate [l/s] + + Returns + ------- + None + """ + self._vfr = vfr + self._mfr = None + + @property + def mfr(self) -> float: + """ + This function returns the mass flow rate. Either based on a given mass flow rate, + or calculated based on a volume flow rate. + + Returns + ------- + float + mass flow rate [kg/s] + """ + if self._mfr is not None: + return self._mfr + return self.vfr / 1000 * self.rho + + @mfr.setter + def mfr(self, mfr: float) -> None: + """ + This function sets the mass flow rate. + The vfr will hence be set to 0. + + Parameters + ---------- + mfr : float + Mass flow rate [kg/s] + + Returns + ------- + None + """ + self._mfr = mfr + self._vfr = None def set_mass_flow_rate(self, mfr: float) -> None: """ diff --git a/GHEtool/VariableClasses/GFunction.py b/GHEtool/VariableClasses/GFunction.py index 6d56020a..6a76fb53 100644 --- a/GHEtool/VariableClasses/GFunction.py +++ b/GHEtool/VariableClasses/GFunction.py @@ -9,11 +9,14 @@ from .CustomGFunction import _time_values +from GHEtool.VariableClasses.Short_term_effects import update_pygfunction_short_term_effects from GHEtool.VariableClasses.cylindrical_correction import update_pygfunction -# add cylindrical correction to pygfunction -update_pygfunction() - +# add short-term effects to pygfunction +print('update_pygfunction_short_term_effects STARTING') +update_pygfunction_short_term_effects() +print('update_pygfunction_short_term_effects FINISHED') +#update_pygfunction() class FIFO: """ @@ -87,7 +90,7 @@ class GFunction: DEFAULT_TIMESTEPS: np.ndarray = _time_values() DEFAULT_NUMBER_OF_TIMESTEPS: int = DEFAULT_TIMESTEPS.size - DEFAULT_STORE_PREVIOUS_VALUES: bool = True + DEFAULT_STORE_PREVIOUS_VALUES: bool = False def __init__(self): self._store_previous_values: bool = GFunction.DEFAULT_STORE_PREVIOUS_VALUES @@ -184,6 +187,7 @@ def gvalues(time_values: np.ndarray, borefield: List[gt.boreholes.Borehole], alp gvalues : np.ndarray 1D array with all the requested gvalues """ + # check if the value is in the fifo_list # if the value is in self.depth_array, there is no problem, since the interpolation will be exact anyway if self.fifo_list.in_fifo_list(depth) and depth not in self.depth_array: @@ -218,7 +222,8 @@ def gvalues(time_values: np.ndarray, borefield: List[gt.boreholes.Borehole], alp return gfunc_interpolated # calculate the g-values for uniform borehole wall temperature - gfunc_calculated = gt.gfunction.gFunction(borefield, alpha, time_values, options=self.options, method=self.options['method']).gFunc + gfunc_calculated = gt.gfunction.gFunction(borefield, alpha, time_values, options=self.options, method=self.options['method']).gFunc + gfunc_calculated = np.array(gfunc_calculated) if np.any(gfunc_calculated < 0): warnings.warn('There are negative g-values. This can be caused by a large borehole radius.') if self.use_cyl_correction_when_negative: diff --git a/GHEtool/VariableClasses/LoadData/GeothermalLoad/HourlyGeothermalLoad.py b/GHEtool/VariableClasses/LoadData/GeothermalLoad/HourlyGeothermalLoad.py index 6b7eba76..25a6ef8a 100644 --- a/GHEtool/VariableClasses/LoadData/GeothermalLoad/HourlyGeothermalLoad.py +++ b/GHEtool/VariableClasses/LoadData/GeothermalLoad/HourlyGeothermalLoad.py @@ -1,5 +1,6 @@ from __future__ import annotations +import matplotlib.pyplot as plt import numpy as np import pandas as pd @@ -414,3 +415,43 @@ def correct_for_start_month(self, array: np.ndarray) -> np.ndarray: if self.start_month == 1: return array return np.concatenate((array[self._start_hour :], array[: self._start_hour])) + + def plot_load_duration(self, legend: bool = False) -> Tuple[plt.Figure, plt.Axes]: + """ + This function makes a load-duration curve from the hourly values. + + Parameters + ---------- + legend : bool + True if the figure should have a legend + + Returns + ---------- + Tuple + plt.Figure, plt.Axes + """ + # sort heating and cooling load + heating = self.hourly_heating_load.copy() + heating[::-1].sort() + + cooling = self.hourly_cooling_load.copy() + cooling[::-1].sort() + cooling = cooling * (-1) + # create new figure and axes if it not already exits otherwise clear it. + fig = plt.figure() + ax = fig.add_subplot(111) + # add sorted loads to plot + ax.step(np.arange(0, 8760, 1), heating, "r-", label="Heating") + ax.step(np.arange(0, 8760, 1), cooling, "b-", label="Cooling") + # create 0 line + ax.hlines(0, 0, 8759, color="black") + # add labels + ax.set_xlabel("Time [hours]") + ax.set_ylabel("Power [kW]") + # set x limits to 8760 + ax.set_xlim(0, 8760) + # plot legend if wanted + if legend: + ax.legend() # + plt.show() + return fig, ax diff --git a/GHEtool/VariableClasses/LoadData/_LoadData.py b/GHEtool/VariableClasses/LoadData/_LoadData.py index 85c2054c..b878c337 100644 --- a/GHEtool/VariableClasses/LoadData/_LoadData.py +++ b/GHEtool/VariableClasses/LoadData/_LoadData.py @@ -36,6 +36,8 @@ def __init__(self, hourly_resolution: bool, simulation_period: int = DEFAULT_SIM self._all_months_equal: bool = True # true if it is assumed that all months are of the same length self._dhw_yearly: float = 0. self._start_month: float = 1 + self.exclude_DHW_from_peak: bool = False # by default, the DHW increase the peak load. Set to false, + # if you only want the heating load to determine the peak in extraction @property def start_month(self) -> int: @@ -225,7 +227,7 @@ def yearly_cooling_load(self) -> float: @property def baseload_heating_simulation_period(self) -> np.ndarray: """ - This function returns the baseload heating in kWh/month for a whole simulation period. + This function returns the baseload heating in kWh/month for the whole simulation period. Returns ------- @@ -249,7 +251,7 @@ def baseload_cooling_simulation_period(self) -> np.ndarray: @property def peak_heating_simulation_period(self) -> np.ndarray: """ - This function returns the peak heating in kW/month for a whole simulation period. + This function returns the peak heating in kW/month for the whole simulation period. Returns ------- @@ -261,7 +263,7 @@ def peak_heating_simulation_period(self) -> np.ndarray: @property def peak_cooling_simulation_period(self) -> np.ndarray: """ - This function returns the peak cooling in kW/month for a whole simulation period. + This function returns the peak cooling in kW/month for the whole simulation period. Returns ------- @@ -273,7 +275,7 @@ def peak_cooling_simulation_period(self) -> np.ndarray: @property def baseload_heating_power_simulation_period(self) -> np.ndarray: """ - This function returns the average heating power in kW avg/month for a whole simulation period. + This function returns the average heating power in kW avg/month for the whole simulation period. Returns ------- @@ -285,7 +287,7 @@ def baseload_heating_power_simulation_period(self) -> np.ndarray: @property def baseload_cooling_power_simulation_period(self) -> np.ndarray: """ - This function returns the average cooling power in kW avg/month for a whole simulation period. + This function returns the average cooling power in kW avg/month for the whole simulation period. Returns ------- @@ -294,6 +296,55 @@ def baseload_cooling_power_simulation_period(self) -> np.ndarray: """ return np.tile(self.baseload_cooling_power, self.simulation_period) + @property + def yearly_heating_load_simulation_period(self) -> np.array: + """ + This function returns the yearly heating demand in kWh/year for the whole simulation period. + + Returns + ------- + yearly heating : np.ndarray + yearly heating for the whole simulation period + """ + return np.sum(np.reshape(self.baseload_heating_simulation_period, (self.simulation_period, 12)), axis=1) + + @property + def yearly_cooling_load_simulation_period(self) -> np.array: + """ + This function returns the yearly cooling demand in kWh/year for the whole simulation period. + + Returns + ------- + yearly heating : np.ndarray + yearly cooling for the whole simulation period + """ + return np.sum(np.reshape(self.baseload_cooling_simulation_period, (self.simulation_period, 12)), axis=1) + + @property + def yearly_heating_peak_simulation_period(self) -> np.array: + """ + This function returns the yearly heating peak in kW/year for the whole simulation period. + + Returns + ------- + yearly heating : np.ndarray + yearly heating for the whole simulation period + """ + return np.max(np.reshape(self.peak_heating_simulation_period, (self.simulation_period, 12)), axis=1) + + @property + def yearly_cooling_peak_simulation_period(self) -> np.array: + """ + This function returns the yearly cooling peak in kW/year for the whole simulation period. + + Returns + ------- + yearly heating : np.ndarray + yearly cooling for the whole simulation period + """ + return np.max(np.reshape(self.peak_cooling_simulation_period, (self.simulation_period, 12)), axis=1) + + @property def imbalance(self) -> float: """ @@ -663,6 +714,8 @@ def dhw_power(self) -> float: ------- dhw power : float """ + if self.exclude_DHW_from_peak: + return 0 return self._dhw_yearly / 8760 @abc.abstractmethod diff --git a/GHEtool/VariableClasses/Short_term_effects.py b/GHEtool/VariableClasses/Short_term_effects.py index e0369f70..7f85f58f 100644 --- a/GHEtool/VariableClasses/Short_term_effects.py +++ b/GHEtool/VariableClasses/Short_term_effects.py @@ -1 +1,322 @@ -print('test') \ No newline at end of file +""" +This document contains the code to add the short-term dynamic effects to the pygfunction package, +as was descripted by Meertens (2024). +Note that this is a temporary solution, until the code will move to pygfunction. +""" + +import pygfunction as gt +import numpy as np + +from math import pi + +from pygfunction.boreholes import Borehole, _EquivalentBorehole, find_duplicates +from pygfunction.heat_transfer import finite_line_source, finite_line_source_vectorized, \ + finite_line_source_equivalent_boreholes_vectorized +from pygfunction.networks import network_thermal_resistance + +from scipy.integrate import quad_vec +from scipy.special import j0, j1, y0, y1 + +from scipy.interpolate import interp1d as interp1d +from time import perf_counter + +#update pygfunction +def evaluate_g_function(self, time): + """ + Evaluate the g-function. + + Parameters + ---------- + time : float or array + Values of time (in seconds) for which the g-function is evaluated. + + Returns + ------- + gFunction : float or array + Values of the g-function + + """ + time = np.maximum(np.atleast_1d(time), 0.) + assert len(time) == 1 or np.all(time[:-1] <= time[1:]), \ + "Time values must be provided in increasing order." + # Save time values + self.time = time + if self.solver.disp: + print(60*'-') + print(f"Calculating g-function for boundary condition : " + f"'{self.boundary_condition}'".center(60)) + print(60*'-') + # Initialize chrono + tic = perf_counter() + + # Evaluate g-function + self.gFunc = self.solver.solve(time, self.alpha) + + print('in new evaluat g function') + + print('self.short_term_effects', self.solver.short_term_effects) + + if self.solver.short_term_effects: + self.gFunc = _ShortTermEffects(self, self.time, self.gFunc, self.boreholes, self.alpha, self.solver.ground_data, + self.solver.fluid_data, self.solver.pipe_data, self.solver.borefield) + toc = perf_counter() + else: + toc = perf_counter() + + if self.solver.disp: + print(f'Total time for g-function evaluation: ' + f'{toc - tic:.3f} sec') + print(60*'-') + return self.gFunc + + + +def __init__(self, boreholes, network, time, boundary_condition, + nSegments=8, segment_ratios=gt.utilities.segment_ratios, + approximate_FLS=False, mQuad=11, nFLS=10, + linear_threshold=None, cylindrical_correction=False, short_term_effects=True, + ground_data=None, fluid_data=None, pipe_data=None, borefield=None, + disp=False, profiles=False, kind='linear', dtype=np.double, + **other_options): + self.boreholes = boreholes + self.network = network + # Convert time to a 1d array + + print('new init is used in pygfunction') + + self.time = np.atleast_1d(time).flatten() + self.linear_threshold = linear_threshold + self.cylindrical_correction = cylindrical_correction + self.short_term_effects = short_term_effects + print('short_term_effects in init', short_term_effects) + self.ground_data = ground_data + self.fluid_data = fluid_data + self.pipe_data = pipe_data + self.borefield = borefield + self.r_b_max = np.max([b.r_b for b in self.boreholes]) + self.boundary_condition = boundary_condition + nBoreholes = len(self.boreholes) + # Format number of segments and segment ratios + if type(nSegments) is int: + self.nBoreSegments = [nSegments] * nBoreholes + else: + self.nBoreSegments = nSegments + if isinstance(segment_ratios, np.ndarray): + segment_ratios = [segment_ratios] * nBoreholes + elif segment_ratios is None: + segment_ratios = [np.full(n, 1./n) for n in self.nBoreSegments] + elif callable(segment_ratios): + segment_ratios = [segment_ratios(n) for n in self.nBoreSegments] + self.segment_ratios = segment_ratios + # Shortcut for segment_ratios comparisons + self._equal_segment_ratios = \ + (np.all(np.array(self.nBoreSegments, dtype=np.uint) == self.nBoreSegments[0]) + and np.all([np.allclose(segment_ratios, self.segment_ratios[0]) for segment_ratios in self.segment_ratios])) + # Boreholes with a uniform discretization + self._uniform_segment_ratios = [ + np.allclose(segment_ratios, + segment_ratios[0:1], + rtol=1e-6) + for segment_ratios in self.segment_ratios] + # Find indices of first and last segments along boreholes + self._i0Segments = [sum(self.nBoreSegments[0:i]) + for i in range(nBoreholes)] + self._i1Segments = [sum(self.nBoreSegments[0:(i + 1)]) + for i in range(nBoreholes)] + self.approximate_FLS = approximate_FLS + self.mQuad = mQuad + self.nFLS = nFLS + self.disp = disp + self.profiles = profiles + self.kind = kind + self.dtype = dtype + # Check the validity of inputs + self._check_inputs() + # Initialize the solver with solver-specific options + self.nSources = self.initialize(**other_options) + return + +def _ShortTermEffects(self, time, gFunc, boreholes, alpha, ground_data, fluid_data, pipe_data, borefield): + + print('gFunc pygfunction', gFunc) + + print('in def _ShortTermEffects') + + from .Dynamic_borhole_model import DynamicsBH + + print('Dynamic_borehole_model imported') + + dynamic_numerical = DynamicsBH(time, gFunc, boreholes, alpha, ground_data, fluid_data, pipe_data, borefield) + + dynamic_numerical.calc_sts_g_functions() + + g = combine_sts_lts( + time, + gFunc, + dynamic_numerical.lntts, + dynamic_numerical.g, + dynamic_numerical.g_plot, + boreholes, alpha, + ) + + return g + +def combine_sts_lts(log_time_lts: list, g_lts: list, log_time_sts: list, g_sts: list, g_plot: list, boreholes, alpha) -> interp1d: + # make sure the short time step doesn't overlap with the long time step + H = boreholes[0].H + ts = H ** 2 / (9. * alpha) # Bore field characteristic time + print('in combine sts lts') + print('time sts', log_time_sts) + print(' size log time sts', len(log_time_sts)) + print('time lst', log_time_lts) + print('size log time lts', len(log_time_lts)) + print('g_lts', g_lts) + print('g_sts', g_sts) + print('stop') + + # log_time_sts: [-48.95425508 -47.6269381 -46.29962112 -44.972304 --> 30 FOUT!! STAAT NU OOK IN s + # log_time_lts: [3.600000e+03 7.200000e+03 1.080000e+04 1.440000e+04 --> 76 + + t_sts_in_s = [] + for k in range(0, len(log_time_sts)): + t_sts_in_s.append(log_time_sts[k]) + + print('time sts in s', t_sts_in_s) + + """ + dummy_list =[] + for b in range(0,len(log_time_lts)): + dummy_list.append(np.log(log_time_lts[b] / ts)) + log_time_lts = dummy_list + + dummy_list = [] + for b in range(0, len(g_lts)): + dummy_list.append(g_lts[b]) + g_lts = dummy_list + + print('log time sts', log_time_sts) + print(' size log time sts', len(log_time_sts)) + print('log time lst', log_time_lts) + print('size log time lts', len(log_time_lts)) + """ + print('make sure the short time step does not overlap with the long time step, MAKING ONE LIST NOW') + + # het noemt wel allemaal log, maar alles staat in s (nodig als input voor GHEtool) + log_time_sts = t_sts_in_s + max_log_time_sts = max(log_time_sts) + min_log_time_lts = min(log_time_lts) + + if max_log_time_sts < min_log_time_lts: + + log_time = log_time_sts + log_time_lts + g = g_sts + g_lts + else: + # find where to stop in sts + i = 0 + value = max_log_time_sts + while value >= log_time_lts[i]: + i += 1 + + g_tmp = interp1d(log_time_sts, g_sts) + g_plot_tmp = interp1d(log_time_sts, g_plot) + uniform_g_sts = g_tmp(log_time_lts[0:i]) + uniform_g_plot = g_plot_tmp(log_time_lts[0:i]) + + print('log time dat van sts', log_time_lts[0:i]) + + print('uniform g sts', uniform_g_sts) + log_time = log_time_lts + print('log time', log_time) + print(len(log_time)) + + dummy = [] + for b in range(0, len(uniform_g_sts)): + dummy.append(uniform_g_sts[b]) + + uniform_g_sts = dummy + print('uniform g sts na dummy list', uniform_g_sts) + + dummy = [] + for b in range(0, len(uniform_g_plot)): + dummy.append(uniform_g_plot[b]) + + uniform_g_plot = dummy + + dummy = [] + for b in range(0, len(g_lts)): + dummy.append(g_lts[b]) + + g_lts = dummy + print('g_lts na dummy list', g_lts) + + g = uniform_g_sts + g_lts[i:] + + g_lst_plot = [] + g_lst_plot = [x + 0.165*2*pi*2.88 for x in g_lts] + + print('g totaal klaar', g) + print(len(g)) + print('lekker bezig') + + print('log time ', log_time) + print('size log time', len(log_time)) + print('g voor interp1d', g) + print('size g voor interp1d', len(g)) + + """ + # plot de short term and long term g function + plt.rc('font', size=9) + plt.rc('xtick', labelsize=9) + plt.rc('ytick', labelsize=9) + plt.rc('lines', lw=1.5, markersize=5.0) + plt.rc('savefig', dpi=500) + fig = plt.figure() + + ax = fig.add_subplot(111) + ax.set_xlabel(r'$tijd$ [s] ', fontsize=18) + ax.set_ylabel(r'$g$-functie', fontsize=18) + + plt.tight_layout() + + ax.plot(log_time_lts[0:i], uniform_g_plot) + ax.plot(log_time_lts[0:25], g_lst_plot[0:25]) + #ax.plot(log_time_lts[0:11], g[0:11], '^') + + ax.legend([f'Numeriek model ', + f'ELB model ', + f'Gecombineerd model met 76 tijdstappen (voor interpolatie)'], fontsize=10) + # ax.set_title(f'Field of {nBoreholes} boreholes') + + """ + """ + x = log_time_lts[0:i] + y = uniform_g_sts + print(x,y) + plt.plot(x, y) + + x = log_time_lts + y = g_lts + print(x,y) + plt.plot(x, y) + plt.show() + """ + return g + +def update_pygfunction_short_term_effects() -> None: + """ + This function updates pygfunction by adding the cylindrical correction methods to it. + + Returns + ------- + None + """ + print('IN update_pygfunctions_short_term_effects') + + gt.gfunction._ShortTermEffects = _ShortTermEffects + print('short_term_effects succesfully added to pygfunction') + gt.gfunction.gFunction.evaluate_g_function = evaluate_g_function + print('evaluate gfunction succesfully replaced') + gt.gfunction._BaseSolver.__init__ = __init__ + print('__init__ succesfully replaced') + + diff --git a/GHEtool/VariableClasses/cylindrical_correction.py b/GHEtool/VariableClasses/cylindrical_correction.py index c0da749d..1a78fa6f 100644 --- a/GHEtool/VariableClasses/cylindrical_correction.py +++ b/GHEtool/VariableClasses/cylindrical_correction.py @@ -144,6 +144,7 @@ def thermal_response_factors(self, time, alpha, kind='linear'): if self.disp: print('Calculating segment to segment response factors ...', end='') + # Number of time values nt = len(np.atleast_1d(time)) # Initialize chrono @@ -460,6 +461,56 @@ def __init__(self, boreholes, network, time, boundary_condition, self.nSources = self.initialize(**other_options) return +def evaluate_g_function(self, time): + """ + Evaluate the g-function. + + Parameters + ---------- + time : float or array + Values of time (in seconds) for which the g-function is evaluated. + + Returns + ------- + gFunction : float or array + Values of the g-function + + """ + time = np.maximum(np.atleast_1d(time), 0.) + assert len(time) == 1 or np.all(time[:-1] <= time[1:]), \ + "Time values must be provided in increasing order." + # Save time values + self.time = time + if self.solver.disp: + print(60*'-') + print(f"Calculating g-function for boundary condition : " + f"'{self.boundary_condition}'".center(60)) + print(60*'-') + # Initialize chrono + tic = perf_counter() + + # Evaluate g-function + self.gFunc = self.solver.solve(time, self.alpha) + + print('in new evaluat g function') + + print('self.short_term_effects', self.alpha) + + if self.cylindrical_correction: + print(self.cylindrical_correction) + self.gFunc = _ShortTermEffects(self, self.time, self.gFunc, self.boreholes, self.alpha, self.ground_data, + self.fluid_data, self.pipe_data, self.borefield) + toc = perf_counter() + else: + toc = perf_counter() + + if self.solver.disp: + print(f'Total time for g-function evaluation: ' + f'{toc - tic:.3f} sec') + print(60*'-') + return self.gFunc + + def update_pygfunction() -> None: """ @@ -471,6 +522,7 @@ def update_pygfunction() -> None: """ gt.heat_transfer.cylindrical_heat_source = cylindrical_heat_source gt.heat_transfer.infinite_line_source = infinite_line_source + gt.gfunction.gFunction.evaluate_g_function = evaluate_g_function gt.gfunction._Equivalent.thermal_response_factors = thermal_response_factors gt.gfunction._BaseSolver.solve = solve gt.gfunction._BaseSolver.__init__ = __init__ \ No newline at end of file