1096 lines
43 KiB
Python
1096 lines
43 KiB
Python
|
"""
|
|||
|
|
V5.0_beta: Komplette Neuprogrammierung der Potenzialberechnung
|
|||
|
|
Vorgegeben werden die Positionen des Erdsondenfeldes (beliebige Koordinaten), die Volllaststunden für Heizen und Kühlen, sowie das gewünschte Leistungsverhältnis (PHZ/PKL) auf der kalten Seite der Wärmepumpe
|
|||
|
|
Wenn keine Vorgabe der Volllaststunden vorhanden, werden die Norm_Betriebsstunden verwendet und als Leistungsverhältnis das deltaT von ungestörter Untergrundtemperatur zu den Heiz-Kühllimits
|
|||
|
|
Ergebnis bei Leistungsvorgabe: Sondenleistung pro Laufmeter Sonde für Heizen und Kühlen, unter Einhaltung des Leistungsverhältnisses und BEtriebsstunden sowie Deckungsgrad
|
|||
|
|
Ergebnis ohne Leistungsvorgabe: Sondenleistung pro Laufmeter Sonde für Heizen und Kühlen, unter Einhaltung Normbetriebsstunden
|
|||
|
|
zusätzlich wird das Leistungspotenzial bei ausgeglichener Jahresenergiebilanz berechnet
|
|||
|
|
zusätzlich werden die Ergebnisse grafisch ausgegeben!
|
|||
|
|
|
|||
|
|
Calculation of g-functions using uniform heat extraction rates.
|
|||
|
|
|
|||
|
|
Programm zur Erstdimensionierung von Erdwärmesonden für ein Grundstück, nach Vorgabe der Positionen des Erdsondenfeldes und Optional die Leistungen und Volllaststunden für Heizen und Kühlen
|
|||
|
|
|
|||
|
|
beta5: Grafikausgabe DEUTSCH
|
|||
|
|
COP und EER Vorgabe möglich - vorerst fix COP=5 und EER=0
|
|||
|
|
option für Berechnung der NORMBS
|
|||
|
|
option für Definition Sondenfeld als Rechteck oder mit Koordinaten
|
|||
|
|
Berechnung der durchschnittlichen Entfernung der nächsten Nachbarsonden, bei Vorgabe mit Koordinaten
|
|||
|
|
SA_FF ersetzt durch nBoreholes
|
|||
|
|
Berücksichtige COP, SCOP, EER und SEER bei der Leistungsvorgabe
|
|||
|
|
|
|||
|
|
beta6: Vorgabe der Vorlauftemperatur des Heizsystems zur Berechnung des COP (neuer parameter "T_radiator")
|
|||
|
|
Berücksichtige Steigerung der Unergrundtemperatur bei Sonden tiefer als 100 m: (neuer parameter "gradient", "GT_calc")
|
|||
|
|
deaktiviere 3D Bohrgrafik
|
|||
|
|
Berechne Steigerung der Leistung bei BALANCED mode (neuer parameter: "cover_rise")
|
|||
|
|
Enable/Disable DEBUG or RESULT print messages
|
|||
|
|
|
|||
|
|
"""
|
|||
|
|
|
|||
|
|
from __future__ import absolute_import, division, print_function
|
|||
|
|
from multiprocessing import freeze_support
|
|||
|
|
|
|||
|
|
import sys
|
|||
|
|
import os
|
|||
|
|
|
|||
|
|
import numpy as np
|
|||
|
|
|
|||
|
|
import pygfunction as gt
|
|||
|
|
from scipy.constants import pi
|
|||
|
|
import matplotlib.pyplot as plt
|
|||
|
|
|
|||
|
|
import base64
|
|||
|
|
import io
|
|||
|
|
|
|||
|
|
|
|||
|
|
def main():
|
|||
|
|
# -------------------------------------------------------------------------
|
|||
|
|
# Simulation parameters
|
|||
|
|
# -------------------------------------------------------------------------
|
|||
|
|
|
|||
|
|
# Find relative path of application for File I/O
|
|||
|
|
if getattr(sys, 'frozen', False):
|
|||
|
|
application_path = os.path.dirname(sys.argv[0])
|
|||
|
|
print("frozen")
|
|||
|
|
else:
|
|||
|
|
application_path = os.path.dirname(os.path.abspath(__file__))
|
|||
|
|
|
|||
|
|
PNGPath = os.path.join(application_path, 'Results_userdefined.png')
|
|||
|
|
PNGPath_bal = os.path.join(application_path, 'Results_autobalanced.png')
|
|||
|
|
PNGPath_borefield = os.path.join(application_path, 'borefield.png')
|
|||
|
|
|
|||
|
|
# Enable or Disable DEBUG print messages and RESULT print messages
|
|||
|
|
DEBUG = False
|
|||
|
|
RESULT = False
|
|||
|
|
|
|||
|
|
def log(s):
|
|||
|
|
if DEBUG:
|
|||
|
|
print(s)
|
|||
|
|
|
|||
|
|
# DEFINE FIXED PARAMETERS (FOR EXPERT USE)
|
|||
|
|
# Define coefficient of performance for heat pump, value "0" is direct heating or cooling
|
|||
|
|
# COP = float("5")
|
|||
|
|
# EER = float("0")
|
|||
|
|
|
|||
|
|
Tmin = float("-1.5") # minimum Fluid Temperature for Heating (°C)
|
|||
|
|
Tmax = float("28") # Maximum Fluid Temperature for Cooling (°C)
|
|||
|
|
years = int("20") # operating time
|
|||
|
|
D = float("1.") # Borehole buried depth (m)
|
|||
|
|
r_b = float("0.075") # Borehole radius (m)
|
|||
|
|
flowrate = float("0.4") # flow rate per borehole (kg/s)
|
|||
|
|
cv = float("2.2") # volumetric heat capacity of the earth (MJ/m³/K)
|
|||
|
|
# ground temperature gradient (K/m) for depth dependent ground temperature
|
|||
|
|
gradient = float("0.03")
|
|||
|
|
|
|||
|
|
# calculate scenario for balanced load if BALANCED > 0
|
|||
|
|
BALANCED = 1
|
|||
|
|
|
|||
|
|
# DEFINE SITE-DEPENDENT PARAMETERS
|
|||
|
|
# import mean surface temperature °C from map
|
|||
|
|
BT = float(sys.argv[1])
|
|||
|
|
GT = float(sys.argv[2]) # mean underground temperature °C
|
|||
|
|
# heat conductivity of the earth (W/m/K)
|
|||
|
|
lamda = float(sys.argv[3])
|
|||
|
|
# typical operational hours for heating (h/yr)
|
|||
|
|
BS_HZ_Norm = float(sys.argv[4])
|
|||
|
|
# typical operational hours for cooling (h/yr)
|
|||
|
|
BS_KL_Norm = float(sys.argv[5])
|
|||
|
|
# known operational hours for heating (h/yr)
|
|||
|
|
BS_HZ = float(sys.argv[6])
|
|||
|
|
BS_KL = float(sys.argv[7]) # known operational hours for cooling
|
|||
|
|
P_HZ = float(sys.argv[8]) # known heating power of heat pump
|
|||
|
|
# known cooling power (of heat pump), has to be negative
|
|||
|
|
P_KL = -float(sys.argv[9])
|
|||
|
|
H = float(sys.argv[10]) # Borehole length (m)
|
|||
|
|
T_radiator = float(sys.argv[11])
|
|||
|
|
|
|||
|
|
# GT = float("11.6") # mean underground temperature (°C)
|
|||
|
|
# lamda = float("2.0") # heat conductivity of the earth (W/m/K)
|
|||
|
|
# BS_HZ = float("0") # operational hours for heating (h/yr)
|
|||
|
|
# BS_KL = float("0") # operational hours for cooling (h/yr)
|
|||
|
|
# P_HZ = float("17") # heating power of building (kW)
|
|||
|
|
# P_KL = float("-11.5") # cooling power of building, use NEGATIVE VALUES! (kW)
|
|||
|
|
# H = float("250") # Borehole length (m)
|
|||
|
|
# T_radiator = float("35") # inlet temperature of heating system (30 - 60 °C) for COP calculation
|
|||
|
|
|
|||
|
|
# # option 1.1:
|
|||
|
|
# # calculate typical operation hours from BT:
|
|||
|
|
# BT = float("11.8") # mean surface temperature °C
|
|||
|
|
# BS_HZ_Norm = np.round(200. + 3380. * np.power(BT,-0.276),0)
|
|||
|
|
# BS_KL_Norm = np.round(max((BT-8.)*180.,0.),0)
|
|||
|
|
|
|||
|
|
# # option 1.2:
|
|||
|
|
# # import typical operation hours from maps
|
|||
|
|
# BS_HZ_Norm = float(sys.argv[5]) # typical operational hours for heating (h/yr)
|
|||
|
|
# BS_KL_Norm = float(sys.argv[6]) # typical operational hours for cooling (h/yr)
|
|||
|
|
# BS_HZ_Norm = float("1902") # typical operational hours for heating (h/yr)
|
|||
|
|
# BS_KL_Norm = float("1000") # typical operational hours for cooling (h/yr)
|
|||
|
|
|
|||
|
|
# DEFINE BHE Field
|
|||
|
|
# option 2.1: define rectangular field
|
|||
|
|
# if N_1 or N_2 zero, then borefiled is defined from coordinates of bore_position
|
|||
|
|
N_1 = int("0")
|
|||
|
|
# if N_1 and N_2 >0, then borefiled ignored and rectangular field is defined
|
|||
|
|
N_2 = int("0")
|
|||
|
|
# Borehole spacing for unbalanced load (m), if bore_position is ignored
|
|||
|
|
BB = float("5")
|
|||
|
|
|
|||
|
|
# option 2.2: define coordinates of bore_positions, can be irregular, if ractangular field is not defined
|
|||
|
|
# bore_position=np.array([[ 0., 0.], [ 0., 5.], [0., 10.], [ 5., 0.], [ 5., 5.]])
|
|||
|
|
l = list(map(lambda x: x.strip(" []"), sys.argv[12].split(",")))
|
|||
|
|
bore_position_temp = []
|
|||
|
|
for i in range(0, len(l) - 1, 2):
|
|||
|
|
bore_position_temp.append([float(l[i]), float(l[i + 1])])
|
|||
|
|
bore_position = np.array(bore_position_temp)
|
|||
|
|
|
|||
|
|
log('\n')
|
|||
|
|
log("=================================================================")
|
|||
|
|
|
|||
|
|
# calculate depth dependent underground temperature
|
|||
|
|
if (H > 100.):
|
|||
|
|
GTcalc = GT + (H-100)/2 * gradient
|
|||
|
|
else:
|
|||
|
|
GTcalc = GT
|
|||
|
|
# calculate maximum deltaT
|
|||
|
|
deltaT_HZ = GTcalc - Tmin
|
|||
|
|
deltaT_KL = GTcalc - Tmax # negative
|
|||
|
|
|
|||
|
|
# if the field is defined by retangular array N_1 x N_2
|
|||
|
|
log(f"calculate bore field")
|
|||
|
|
if (N_1 > 0 and N_2 > 0):
|
|||
|
|
if ((N_1*N_2) == 1):
|
|||
|
|
BB_L3 = (f"of infinite")
|
|||
|
|
else:
|
|||
|
|
BB_L3 = (f"of {BB} m")
|
|||
|
|
boreField = gt.boreholes.rectangle_field(N_1, N_2, BB, BB, H, D, r_b)
|
|||
|
|
bore_position = np.zeros([N_1*N_2, 2], dtype=float)
|
|||
|
|
nBoreholes = len(boreField)
|
|||
|
|
BB_mean = BB
|
|||
|
|
else:
|
|||
|
|
# for bore field with coordinates: create boreField with xy coordinates
|
|||
|
|
b = []
|
|||
|
|
for i in (range(len(bore_position))):
|
|||
|
|
b.append(gt.boreholes.Borehole(H=H, D=D, r_b=r_b,
|
|||
|
|
x=bore_position[i, 0], y=bore_position[i, 1]))
|
|||
|
|
# remove duplicates:
|
|||
|
|
b = gt.boreholes.remove_duplicates(b, disp=DEBUG)
|
|||
|
|
|
|||
|
|
nBoreholes = len(b)
|
|||
|
|
|
|||
|
|
if nBoreholes > 1:
|
|||
|
|
# calculate BHE distances and direction
|
|||
|
|
np.set_printoptions(threshold=np.inf, formatter={
|
|||
|
|
'float': '{: 0.3f}'.format})
|
|||
|
|
MATdistance = np.zeros((len(b), len(b)))
|
|||
|
|
MATquadrant = np.zeros((len(b), len(b)), dtype=int)
|
|||
|
|
NEAREST = np.tile(np.nan, (len(b), 4))
|
|||
|
|
for i in (range(len(b))):
|
|||
|
|
b1 = b[i]
|
|||
|
|
xy1 = b[i].position()
|
|||
|
|
for j in (range(len(b))):
|
|||
|
|
b2 = b[j]
|
|||
|
|
xy2 = b[j].position()
|
|||
|
|
MATdistance[i, j] = b1.distance(b2)
|
|||
|
|
XX, YY = np.subtract(xy2, xy1)
|
|||
|
|
# caculate Direction (Quadrant)
|
|||
|
|
if (XX > 0 and YY >= 0):
|
|||
|
|
MATquadrant[i, j] = 1 # in 1st quadrant or pos x-axis
|
|||
|
|
elif (XX <= 0 and YY > 0):
|
|||
|
|
MATquadrant[i, j] = 2 # in 2nd quadrant or pos y-axis
|
|||
|
|
elif (XX < 0 and YY <= 0):
|
|||
|
|
MATquadrant[i, j] = 3 # in 3rd quadrant or neg x-axis
|
|||
|
|
elif (XX >= 0 and YY < 0):
|
|||
|
|
MATquadrant[i, j] = 4 # in 4th quadrant or neg y-axis
|
|||
|
|
else:
|
|||
|
|
MATquadrant[i, j] = 5 # ident
|
|||
|
|
MATdistance[i, i] = 0. # distance b1 to themselve = zero
|
|||
|
|
|
|||
|
|
# calculate nearest distance in each quadrant for all BHEs
|
|||
|
|
for k in range(4):
|
|||
|
|
a = (MATdistance[i, :][MATquadrant[i, :] == k+1])
|
|||
|
|
if (np.size(a) > 0):
|
|||
|
|
NEAREST[i, k] = min(a)
|
|||
|
|
# log(MATdistance)
|
|||
|
|
# log(MATquadrant)
|
|||
|
|
# log(NEAREST)
|
|||
|
|
BB_mean = np.nanmean(NEAREST)
|
|||
|
|
BB_mean = np.round(np.nanmean(NEAREST), 1)
|
|||
|
|
log(f"mean distance to other BHEs: {BB_mean}")
|
|||
|
|
else:
|
|||
|
|
BB_mean = 0.
|
|||
|
|
|
|||
|
|
boreField = b
|
|||
|
|
BB_L3 = ("given by coordinates")
|
|||
|
|
|
|||
|
|
divider = nBoreholes * H
|
|||
|
|
|
|||
|
|
if (T_radiator > 29. and T_radiator < 90.):
|
|||
|
|
Qc, Qe, Pel = get_COP_potential(T_radiator, -3., 1.)
|
|||
|
|
COP = Qc/Pel
|
|||
|
|
Qc, Qe, Pel = get_COP_potential(30., 18., 1.)
|
|||
|
|
EER = Qe/Pel
|
|||
|
|
log(f"COP calculation with {T_radiator:.1f} °C and -3°C: {COP:.2f}")
|
|||
|
|
log(f"EER calculation with 30 °C and 18°C: {EER:.2f}")
|
|||
|
|
|
|||
|
|
# clarify heating/cooling cases
|
|||
|
|
if (BS_HZ == 0):
|
|||
|
|
P_HZ = 0.
|
|||
|
|
if (BS_KL == 0):
|
|||
|
|
P_KL = 0.
|
|||
|
|
if (P_HZ == 0):
|
|||
|
|
BS_HZ = 0.
|
|||
|
|
if (P_KL == 0):
|
|||
|
|
BS_KL = 0.
|
|||
|
|
|
|||
|
|
# if defined, use given building demand
|
|||
|
|
if (abs(P_HZ*BS_HZ) > 0. or abs(P_KL*BS_KL) > 0.):
|
|||
|
|
BS_HZ_L3 = BS_HZ
|
|||
|
|
BS_KL_L3 = BS_KL
|
|||
|
|
|
|||
|
|
if (COP > 0.):
|
|||
|
|
P_HZ = P_HZ * (1-1/COP)
|
|||
|
|
if (EER > 0.):
|
|||
|
|
P_KL = P_KL * (1+1/EER)
|
|||
|
|
|
|||
|
|
if (P_KL == 0.):
|
|||
|
|
Pfactor = 1000.
|
|||
|
|
elif (P_HZ == 0.):
|
|||
|
|
Pfactor = 1/1000.
|
|||
|
|
else:
|
|||
|
|
Pfactor = P_HZ / -P_KL
|
|||
|
|
|
|||
|
|
if (BS_HZ_L3 == 0):
|
|||
|
|
Efactor = 0.
|
|||
|
|
elif (BS_KL_L3 > 0 and Pfactor < 100.):
|
|||
|
|
Efactor = Pfactor * BS_HZ_L3 / BS_KL_L3
|
|||
|
|
else:
|
|||
|
|
Efactor, BS_KL_L3 = 1000., 0.
|
|||
|
|
case = "user"
|
|||
|
|
# if NOT defined, use typical operation hours
|
|||
|
|
else:
|
|||
|
|
BS_HZ_L3 = BS_HZ_Norm
|
|||
|
|
BS_KL_L3 = BS_KL_Norm
|
|||
|
|
Pfactor = deltaT_HZ/-deltaT_KL
|
|||
|
|
if (BS_KL_L3 > 0):
|
|||
|
|
Efactor = Pfactor * BS_HZ_L3 / BS_KL_L3
|
|||
|
|
else:
|
|||
|
|
Efactor = 1000.
|
|||
|
|
|
|||
|
|
Efactor = Pfactor*BS_HZ_Norm/BS_KL_Norm
|
|||
|
|
case = "norm"
|
|||
|
|
|
|||
|
|
if (Efactor > 0.9 and Efactor < 1.1):
|
|||
|
|
BALANCED = 0
|
|||
|
|
log("BALANCED MODE DEAKTIVATED")
|
|||
|
|
|
|||
|
|
# predefinitions
|
|||
|
|
Efactor_bal = 1.0
|
|||
|
|
Pfactor_bal = 1.0
|
|||
|
|
PHZ_L3_bal, PKL_L3_bal = 0., 0.
|
|||
|
|
PHZ_L3, PKL_L3 = 0., 0.
|
|||
|
|
coverHZ, coverKL, cover_user = np.nan, np.nan, np.nan
|
|||
|
|
coverHZ_bal, coverKL_bal, cover_bal, cover_rise = np.nan, np.nan, np.nan, np.nan
|
|||
|
|
Eel_heatpump_user, Eel_chiller_user = 0., 0.
|
|||
|
|
Eel_heatpump_bal, Eel_chiller_bal = 0., 0.
|
|||
|
|
Pel_heatpump_user, Pel_chiller_user = 0., 0.
|
|||
|
|
Pel_heatpump_bal, Pel_chiller_bal = 0., 0.
|
|||
|
|
SCOP, SEER, SCOP_bal, SEER_bal = 0., 0., 0., 0.
|
|||
|
|
|
|||
|
|
# if BHEs are defined, calculate POTENTIAL with user-defined or standard power ratio and operational hours
|
|||
|
|
image_hash = ""
|
|||
|
|
image_hash_bal = ""
|
|||
|
|
image_hash_borefield = ""
|
|||
|
|
if (nBoreholes > 0):
|
|||
|
|
PHZ_L3, PKL_L3, Efactor, SCOP, SEER, image_hash, image_hash_borefield = calculateL3(
|
|||
|
|
GTcalc, lamda, BS_HZ_L3, BS_KL_L3, Tmin, Tmax, years, D, H, r_b, cv, Pfactor, flowrate, boreField, PNGPath, PNGPath_borefield, T_radiator, DEBUG, RESULT)
|
|||
|
|
# if BHEs are defined, calculate POTENTIAL with automatic balanced power ratio and operational hours
|
|||
|
|
BS_HZ_bal = 0
|
|||
|
|
BS_KL_bal = 0
|
|||
|
|
if (nBoreholes > 0 and BALANCED > 0):
|
|||
|
|
Pfactor_bal = deltaT_HZ/-deltaT_KL
|
|||
|
|
if (BS_HZ_L3 >= BS_KL_L3):
|
|||
|
|
BS_HZ_bal = max(BS_HZ_L3, BS_KL_L3)
|
|||
|
|
BS_KL_bal = np.floor(min(4000., BS_HZ_bal*Pfactor_bal/Efactor_bal))
|
|||
|
|
else:
|
|||
|
|
BS_KL_bal = max(BS_HZ_L3, BS_KL_L3)
|
|||
|
|
BS_HZ_bal = np.floor(min(4000., BS_KL_bal/Pfactor_bal*Efactor_bal))
|
|||
|
|
PHZ_L3_bal, PKL_L3_bal, Efactor_bal, SCOP_bal, SEER_bal, image_hash_bal, _ = calculateL3(
|
|||
|
|
GTcalc, lamda, BS_HZ_bal, BS_KL_bal, Tmin, Tmax, years, D, H, r_b, cv, Pfactor_bal, flowrate, boreField, PNGPath_bal, "", T_radiator, DEBUG, RESULT)
|
|||
|
|
|
|||
|
|
if (RESULT):
|
|||
|
|
printdemand(case, P_HZ, P_KL, BS_HZ, BS_KL,
|
|||
|
|
COP, EER, T_radiator, GTcalc, lamda)
|
|||
|
|
if (RESULT):
|
|||
|
|
printratio(Efactor, Pfactor)
|
|||
|
|
|
|||
|
|
# return user defined results:
|
|||
|
|
P_HZ_user = PHZ_L3*divider/1000. # potential from BHEs in kW
|
|||
|
|
P_KL_user = PKL_L3*divider/1000. # potential from BHEs in kW
|
|||
|
|
E_HZ_user = P_HZ_user * BS_HZ_L3 # potential from BHEs in kWh
|
|||
|
|
E_KL_user = P_KL_user * BS_KL_L3 # potential from BHEs in kWh
|
|||
|
|
Efactor_user = Efactor
|
|||
|
|
if (COP > 0):
|
|||
|
|
# electric power for heating kW
|
|||
|
|
Pel_heatpump_user = P_HZ_user / (COP - 1)
|
|||
|
|
if (EER > 0):
|
|||
|
|
# electric power for cooling kW
|
|||
|
|
Pel_chiller_user = P_KL_user / (EER + 1)
|
|||
|
|
if (SCOP > 0):
|
|||
|
|
# electric power for heating kWh
|
|||
|
|
Eel_heatpump_user = E_HZ_user / (SCOP - 1)
|
|||
|
|
if (SEER > 0):
|
|||
|
|
# electric power for cooling kWh
|
|||
|
|
Eel_chiller_user = E_KL_user / (SEER + 1)
|
|||
|
|
if (P_HZ > 0):
|
|||
|
|
coverHZ = 100.0 / (P_HZ*BS_HZ) * E_HZ_user # %
|
|||
|
|
if (P_KL < 0):
|
|||
|
|
coverKL = 100.0 / (P_KL*BS_KL) * E_KL_user # %
|
|||
|
|
if (np.isfinite(coverHZ) or np.isfinite(coverKL)):
|
|||
|
|
cover_user = np.nanmin([coverHZ, coverKL]) # %
|
|||
|
|
|
|||
|
|
if (RESULT):
|
|||
|
|
printresults(case, P_HZ_user, P_KL_user, E_HZ_user, E_KL_user, cover_user, Pel_heatpump_user, Pel_chiller_user,
|
|||
|
|
Eel_heatpump_user, Eel_chiller_user, nBoreholes, H, BB_L3, BB_mean, 0., COP, EER, SCOP, SEER)
|
|||
|
|
if (RESULT):
|
|||
|
|
printratio(Efactor_user, Pfactor)
|
|||
|
|
user_defined_results = [case, P_HZ_user, P_KL_user, E_HZ_user, E_KL_user, cover_user, Pel_heatpump_user, Pel_chiller_user, Eel_heatpump_user,
|
|||
|
|
Eel_chiller_user, nBoreholes, H, BB_L3, BB_mean, 0., COP, EER, SCOP, SEER, Efactor_user, image_hash, image_hash_borefield, GTcalc]
|
|||
|
|
|
|||
|
|
# return automatic balanced results:
|
|||
|
|
P_HZ_bal = P_KL_bal = E_HZ_bal = E_KL_bal = cover_bal = 0
|
|||
|
|
if (BALANCED > 0):
|
|||
|
|
P_HZ_bal = PHZ_L3_bal*divider/1000.
|
|||
|
|
P_KL_bal = PKL_L3_bal*divider/1000.
|
|||
|
|
E_HZ_bal = P_HZ_bal * BS_HZ_bal
|
|||
|
|
E_KL_bal = P_KL_bal * BS_KL_bal
|
|||
|
|
Efactor_bal = Efactor_bal
|
|||
|
|
if (COP > 0):
|
|||
|
|
# electric power for heating kW
|
|||
|
|
Pel_heatpump_bal = P_HZ_bal / (COP - 1)
|
|||
|
|
if (EER > 0):
|
|||
|
|
# electric power for cooling kW
|
|||
|
|
Pel_chiller_bal = P_KL_bal / (EER + 1)
|
|||
|
|
if (SCOP_bal > 0):
|
|||
|
|
# electric power for heating kW
|
|||
|
|
Eel_heatpump_bal = E_HZ_bal / (SCOP_bal - 1)
|
|||
|
|
if (SEER_bal > 0):
|
|||
|
|
# electric power for heating kW
|
|||
|
|
Eel_chiller_bal = E_KL_bal / (SEER_bal + 1)
|
|||
|
|
if (P_HZ > 0):
|
|||
|
|
coverHZ = 100.0 / (P_HZ*BS_HZ) * E_HZ_bal # %
|
|||
|
|
if (P_KL < 0):
|
|||
|
|
coverKL = 100.0 / (P_KL*BS_KL) * E_KL_bal # %
|
|||
|
|
if (np.isfinite(coverHZ) or np.isfinite(coverKL)):
|
|||
|
|
cover_bal = np.nanmin([coverHZ, coverKL]) # %
|
|||
|
|
if (Efactor_user < 1): # more cooling
|
|||
|
|
cover_rise = (P_KL_bal-P_KL_user) / P_KL_user * 100.
|
|||
|
|
else: # more cooling
|
|||
|
|
cover_rise = (P_HZ_bal-P_HZ_user) / P_HZ_user * 100.
|
|||
|
|
|
|||
|
|
if (RESULT):
|
|||
|
|
printresults("AUTO-BALANCED", P_HZ_bal, P_KL_bal, E_HZ_bal, E_KL_bal, cover_bal, Pel_heatpump_bal, Pel_chiller_bal,
|
|||
|
|
Eel_heatpump_bal, Eel_chiller_bal, nBoreholes, H, BB_L3, BB_mean, cover_rise, COP, EER, SCOP_bal, SEER_bal)
|
|||
|
|
if (RESULT):
|
|||
|
|
printratio(Efactor_bal, Pfactor_bal)
|
|||
|
|
|
|||
|
|
automatic_results = [BALANCED, P_HZ_bal, P_KL_bal, E_HZ_bal, E_KL_bal, cover_bal, Pel_heatpump_bal, Pel_chiller_bal, Eel_heatpump_bal, Eel_chiller_bal,
|
|||
|
|
nBoreholes, H, BB_L3, BB_mean, cover_rise, COP, EER, SCOP_bal, SEER_bal, Efactor_bal, image_hash_bal, BS_HZ_bal, BS_KL_bal, T_radiator]
|
|||
|
|
line = list(map(str, user_defined_results)) + \
|
|||
|
|
list(map(str, automatic_results))
|
|||
|
|
print(line)
|
|||
|
|
|
|||
|
|
# results for post script:
|
|||
|
|
|
|||
|
|
# GTcalc
|
|||
|
|
# COP
|
|||
|
|
# EER
|
|||
|
|
|
|||
|
|
# P_HZ_user
|
|||
|
|
# P_KL_user
|
|||
|
|
# E_HZ_user
|
|||
|
|
# E_KL_user
|
|||
|
|
# cover_user
|
|||
|
|
# Pel_heatpump_user
|
|||
|
|
# Pel_chiller_user
|
|||
|
|
# Eel_heatpump_user
|
|||
|
|
# Eel_chiller_user
|
|||
|
|
# Efactor_user
|
|||
|
|
# SCOP
|
|||
|
|
# SEER
|
|||
|
|
|
|||
|
|
# P_HZ_bal
|
|||
|
|
# P_KL_bal
|
|||
|
|
# BS_HZ_bal
|
|||
|
|
# BS_KL_bal
|
|||
|
|
# E_HZ_bal
|
|||
|
|
# E_KL_bal
|
|||
|
|
# cover_bal
|
|||
|
|
# Pel_heatpump_bal
|
|||
|
|
# Pel_chiller_bal
|
|||
|
|
# Eel_heatpump_bal
|
|||
|
|
# Eel_chiller_bal
|
|||
|
|
# SCOP_bal
|
|||
|
|
# SEER_bal
|
|||
|
|
# cover_rise
|
|||
|
|
|
|||
|
|
# BB_mean
|
|||
|
|
# nBoreholes
|
|||
|
|
|
|||
|
|
# line = sys.argv[1:12] + list(map(str, [PHZ_L3, PKL_L3,...
|
|||
|
|
# log(line)
|
|||
|
|
|
|||
|
|
|
|||
|
|
def printdemand(text, PHZ, PKL, BS_HZ, BS_KL, COP, EER, T_radiator, GTcalc, lamda):
|
|||
|
|
print(f"\n########################### Ergebnisse #######################:")
|
|||
|
|
print(f"\nGebäudeleistung {text}:")
|
|||
|
|
if (COP > 0):
|
|||
|
|
print(f" PHZ = {BS_HZ:.1f} h x {PHZ/(1-1/COP):.1f} kW")
|
|||
|
|
else:
|
|||
|
|
print(f" PHZ = {BS_HZ:.1f} h x {PHZ:.1f} kW")
|
|||
|
|
if (EER > 0):
|
|||
|
|
print(f" PKL = {BS_KL:.1f} h x {PKL/(1+1/EER):.1f} kW")
|
|||
|
|
else:
|
|||
|
|
print(f" PKL = {BS_KL:.1f} h x {PKL:.1f} kW")
|
|||
|
|
print(f" T_Vorlauf_Heizung: {T_radiator:.1f} °C")
|
|||
|
|
print(f"Vorgabe für Erdsonden:")
|
|||
|
|
if (COP > 0):
|
|||
|
|
print(f" COP = {COP:.1f}")
|
|||
|
|
if (EER > 0):
|
|||
|
|
print(f" EER = {EER:.1f}")
|
|||
|
|
print(f" PHZ = {BS_HZ:.1f} h x {PHZ:.1f} kW")
|
|||
|
|
print(f" PKL = {BS_KL:.1f} h x {PKL:.1f} kW")
|
|||
|
|
print(f"Erdreichparameter:")
|
|||
|
|
print(f" mean ground temperature = {GTcalc:.1f}")
|
|||
|
|
print(f" mean thermal conductivity= {lamda:.1f}")
|
|||
|
|
|
|||
|
|
|
|||
|
|
def printratio(Efactor, Pfactor):
|
|||
|
|
if (Efactor > 1/1000. and Efactor < 1000.):
|
|||
|
|
print(f" energy ratio HZ/KL = {Efactor:.3f}")
|
|||
|
|
if (Pfactor > 1/1000. and Pfactor < 1000.):
|
|||
|
|
print(f" power ratio HZ/KL = {Pfactor:.3f}")
|
|||
|
|
|
|||
|
|
|
|||
|
|
def printresults(text, PHZ, PKL, EHZ, EKL, cover, P_heatpump, P_chiller, E_heatpump, E_chiller, nBoreholes, H, BB_L3, BB_mean, cover_rise, COP, EER, SCOP, SEER):
|
|||
|
|
if (np.isfinite(cover)):
|
|||
|
|
covert = (f"{cover:.1f} %")
|
|||
|
|
else:
|
|||
|
|
covert = ("not defined")
|
|||
|
|
|
|||
|
|
print(
|
|||
|
|
f"\nRESULTS {text} for {nBoreholes} BHEs with {H} m depth and average spacing {BB_mean} {BB_L3}:")
|
|||
|
|
print(
|
|||
|
|
f" heat extraction power from BHE = {PHZ:.1f} kW")
|
|||
|
|
print(
|
|||
|
|
f" + electr. power of heat pump (COP = {COP:.1f}) = {P_heatpump:.1f} kW")
|
|||
|
|
print(
|
|||
|
|
f" = heating power for building = {PHZ+P_heatpump:.1f} kW")
|
|||
|
|
print(
|
|||
|
|
f" yearly heat extraction from BHE = {EHZ/1000.:.1f} MWh/a")
|
|||
|
|
print(
|
|||
|
|
f" + electr. energy of heat pump (SCOP = {SCOP:.1f}) = {E_heatpump/1000.:.1f} MWh/a")
|
|||
|
|
print(
|
|||
|
|
f" = yearly heating energy for building = {(EHZ+E_heatpump)/1000.:.1f} MWh/a")
|
|||
|
|
print("\n")
|
|||
|
|
print(
|
|||
|
|
f" heat injection power to BHE = {PKL:.1f} kW")
|
|||
|
|
print(
|
|||
|
|
f" - electr. power of heat pump (EER = {EER:.1f}) = {P_chiller:.1f} kW")
|
|||
|
|
print(
|
|||
|
|
f" = cooling power for building = {PKL-P_chiller:.1f} kW")
|
|||
|
|
print(
|
|||
|
|
f" yearly heat injection to BHE = {EKL/1000.:.1f} MWh/a")
|
|||
|
|
print(
|
|||
|
|
f" - electr. energy of heat pump (SEER = {SEER:.1f}) = {E_chiller/1000.:.1f} MWh/a")
|
|||
|
|
print(
|
|||
|
|
f" = yearly cooling energy for building = {(EKL-E_chiller)/1000.:.1f} MWh/a")
|
|||
|
|
print("\n")
|
|||
|
|
print(f" Coverage from BHEs = {covert}")
|
|||
|
|
if (cover_rise > 0):
|
|||
|
|
print(f" {cover_rise:.1f} % more energy use @ balanced mode")
|
|||
|
|
|
|||
|
|
|
|||
|
|
def calculateL3(T_g, lamda, BS_HZ, BS_KL, Tmin, Tmax, years, D, H, r_b, cv, Pfactor, m_flow_borehole, boreField, PNGPath, pngname2, T_radiator, DEBUG, RESULT):
|
|||
|
|
def log(s):
|
|||
|
|
if DEBUG:
|
|||
|
|
print(s)
|
|||
|
|
|
|||
|
|
if (BS_HZ == 0):
|
|||
|
|
Efactor = 0.
|
|||
|
|
elif (BS_KL > 0 and Pfactor < 1000.):
|
|||
|
|
Efactor = Pfactor * BS_HZ / BS_KL
|
|||
|
|
else:
|
|||
|
|
Efactor, BS_KL = 1000., 0.
|
|||
|
|
|
|||
|
|
# Simulation parameters
|
|||
|
|
dt = 1.*3600. # Time step (s) for hourly simulation
|
|||
|
|
|
|||
|
|
# Pipe dimensions
|
|||
|
|
r_out = 0.016 # Pipe outer radius (m)
|
|||
|
|
r_in = 0.014 # Pipe inner radius (m)
|
|||
|
|
D_s = 0.04 # Shank spacing (m)
|
|||
|
|
epsilon = 1.0e-6 # Pipe roughness (m)
|
|||
|
|
|
|||
|
|
# Pipe positions
|
|||
|
|
# Double U-tube [(x_in1, y_in1), (x_in2, y_in2),
|
|||
|
|
# (x_out1, y_out1), (x_out2, y_out2)]
|
|||
|
|
pos = [(-D_s, 0.), (0., -D_s), (D_s, 0.), (0., D_s)]
|
|||
|
|
|
|||
|
|
# Ground properties
|
|||
|
|
cv = 2.2 # vol heat capacity of earth (MJ/m³/K)
|
|||
|
|
alpha = lamda / cv * 1e-6 # Ground thermal diffusivity (m2/s)
|
|||
|
|
|
|||
|
|
# Grout properties
|
|||
|
|
k_g = 2.0 # Grout thermal conductivity (W/m.K)
|
|||
|
|
|
|||
|
|
# Pipe properties
|
|||
|
|
k_p = 0.35 # Pipe thermal conductivity (W/m.K)
|
|||
|
|
|
|||
|
|
# The fluid is propylene-glycol (20 %) at 20 degC
|
|||
|
|
|
|||
|
|
# ‘Water’ - Complete water solution
|
|||
|
|
# ‘MEG’ - Ethylene glycol mixed with water
|
|||
|
|
# ‘MPG’ - Propylene glycol mixed with water
|
|||
|
|
# ‘MEA’ - Ethanol mixed with water
|
|||
|
|
# ‘MMA’ - Methanol mixed with water
|
|||
|
|
|
|||
|
|
# fluid = gt.media.Fluid('MEA', 30.)
|
|||
|
|
mix = 'MEA'
|
|||
|
|
percent = 12
|
|||
|
|
log(f" fluid type: {mix}")
|
|||
|
|
|
|||
|
|
fluid = gt.media.Fluid(mix, percent)
|
|||
|
|
cp_f = fluid.cp # Fluid specific isobaric heat capacity (J/kg.K)
|
|||
|
|
rho_f = fluid.rho # Fluid density (kg/m3)
|
|||
|
|
mu_f = fluid.mu # Fluid dynamic viscosity (kg/m.s)
|
|||
|
|
k_f = fluid.k # Fluid thermal conductivity (W/m.K)
|
|||
|
|
log(fluid)
|
|||
|
|
|
|||
|
|
# g-Function calculation options
|
|||
|
|
options = {'nSegments': 8, 'disp': DEBUG}
|
|||
|
|
|
|||
|
|
# -------------------------------------------------------------------------
|
|||
|
|
# Initialize bore field and pipe models
|
|||
|
|
# -------------------------------------------------------------------------
|
|||
|
|
|
|||
|
|
nBoreholes = len(boreField)
|
|||
|
|
divider = nBoreholes * H
|
|||
|
|
|
|||
|
|
# Total fluid mass flow rate (kg/s)
|
|||
|
|
m_flow_network = m_flow_borehole*nBoreholes
|
|||
|
|
|
|||
|
|
# Pipe thermal resistance
|
|||
|
|
log(f"calculate pipe resistance")
|
|||
|
|
R_p = gt.pipes.conduction_thermal_resistance_circular_pipe(
|
|||
|
|
r_in, r_out, k_p)
|
|||
|
|
|
|||
|
|
# Fluid to inner pipe wall thermal resistance (Double U-tube in parallel)
|
|||
|
|
m_flow_pipe = m_flow_borehole/2
|
|||
|
|
h_f = gt.pipes.convective_heat_transfer_coefficient_circular_pipe(
|
|||
|
|
m_flow_pipe, r_in, mu_f, rho_f, k_f, cp_f, epsilon)
|
|||
|
|
R_f = 1.0/(h_f*2*pi*r_in)
|
|||
|
|
R_f_p = R_f+R_p
|
|||
|
|
log(f" Rf {R_f:.3f} + Rp {R_p:.3f} = {R_f+R_p:.3f} m.K/W")
|
|||
|
|
|
|||
|
|
# Double U-tube (parallel), same for all boreholes in the bore field
|
|||
|
|
log(f"define U-tubes")
|
|||
|
|
UTubes = []
|
|||
|
|
for borehole in boreField:
|
|||
|
|
UTube = gt.pipes.MultipleUTube(
|
|||
|
|
pos, r_in, r_out, borehole, lamda, k_g, R_f_p, nPipes=2, config='parallel')
|
|||
|
|
UTubes.append(UTube)
|
|||
|
|
|
|||
|
|
image_hash_borefield = ""
|
|||
|
|
if (len(pngname2) > 0):
|
|||
|
|
image_hash_borefield = plotborefield(boreField, pngname2, RESULT)
|
|||
|
|
# Build a network object from the list of UTubes
|
|||
|
|
log(f"building pipe network")
|
|||
|
|
network = gt.networks.Network(
|
|||
|
|
boreField, UTubes, m_flow_network=m_flow_network, cp_f=cp_f)
|
|||
|
|
|
|||
|
|
# -------------------------------------------------------------------------
|
|||
|
|
# Calculate g-function for hourly Simulation
|
|||
|
|
# -------------------------------------------------------------------------
|
|||
|
|
|
|||
|
|
if (Efactor >= 1): # more heating than cooling
|
|||
|
|
BS_first = BS_HZ
|
|||
|
|
BS_second = BS_KL
|
|||
|
|
else: # more cooling than heating
|
|||
|
|
BS_first = BS_KL
|
|||
|
|
BS_second = BS_HZ
|
|||
|
|
|
|||
|
|
tmax = ((years)*8760. + BS_first/2) * 3600. # Maximum time (s)
|
|||
|
|
x1 = BS_first/2
|
|||
|
|
x5 = 8760/2-BS_second/2
|
|||
|
|
x6 = 8760/2+BS_second/2
|
|||
|
|
x10 = 8760-BS_first/2
|
|||
|
|
x11 = 8760.
|
|||
|
|
|
|||
|
|
xout = int(3*8760/2+BS_second/2)-1
|
|||
|
|
|
|||
|
|
Nt = int(np.floor(tmax/dt)) # Number of time steps
|
|||
|
|
time = dt * np.arange(1, Nt+1)
|
|||
|
|
|
|||
|
|
# Load aggregation scheme
|
|||
|
|
log("Load Aggregation Scheme")
|
|||
|
|
LoadAgg = gt.load_aggregation.ClaessonJaved(dt, tmax)
|
|||
|
|
# Get time values needed for g-function evaluation
|
|||
|
|
time_req = LoadAgg.get_times_for_simulation()
|
|||
|
|
# Calculate g-function
|
|||
|
|
log(f"calculate g-function")
|
|||
|
|
gFunc = gt.gfunction.gFunction(
|
|||
|
|
network, alpha, time_req, boundary_condition='MIFT', options=options)
|
|||
|
|
|
|||
|
|
# Initialize load aggregation scheme
|
|||
|
|
LoadAgg.initialize(gFunc.gFunc/(2*pi*lamda))
|
|||
|
|
|
|||
|
|
# -------------------------------------------------------------------------
|
|||
|
|
# Potential
|
|||
|
|
# -------------------------------------------------------------------------
|
|||
|
|
|
|||
|
|
# Evaluate the effective bore field thermal resistance (m.K/W)
|
|||
|
|
R_bmf = gt.networks.network_thermal_resistance(
|
|||
|
|
network, m_flow_network, cp_f)
|
|||
|
|
log(f" Rbmf = {R_bmf:5.3f} m.K/W (effective internal resistance)")
|
|||
|
|
|
|||
|
|
# calculate maximum deltaT
|
|||
|
|
deltaT_HZ = T_g - Tmin
|
|||
|
|
deltaT_KL = T_g - Tmax # negative
|
|||
|
|
|
|||
|
|
if (round(Efactor, 2) < 1): # more cooling than heating: Beginn with cooling period
|
|||
|
|
log(" cooling dominant")
|
|||
|
|
# calculate Resistance of Earth with g-functions after cooling period, year 2
|
|||
|
|
time_superpos_sort, time_superpos_ind_rev, Pfakt = prepare_times_for_gfunction_2J(
|
|||
|
|
Pfactor, Efactor, BS_HZ, BS_KL, years)
|
|||
|
|
# gfun2 = gt.gfunction.gFunction(network, alpha, time=time_superpos_sort*3600., boundary_condition='MIFT',options=options)
|
|||
|
|
gfun3 = np.interp(time_superpos_sort*3600., time_req, gFunc.gFunc)
|
|||
|
|
# gfun3 = np.array(gfun2.gFunc)
|
|||
|
|
gfun3_unsort = gfun3[time_superpos_ind_rev]
|
|||
|
|
product = np.multiply(gfun3_unsort, Pfakt)
|
|||
|
|
R_g_2J = abs(sum(product)/(2*pi*lamda))
|
|||
|
|
log(f" Rg = {R_g_2J:5.3f} m.K/W")
|
|||
|
|
|
|||
|
|
# calculate Resistance of Earth with g-functions after first heating period, year 21
|
|||
|
|
time_superpos_sort, time_superpos_ind_rev, Pfakt, time_sort = prepare_times_for_gfunction_20J(
|
|||
|
|
Pfactor, Efactor, BS_HZ, BS_KL, years)
|
|||
|
|
# gfun2 = gt.gfunction.gFunction(network, alpha, time=time_superpos_sort*3600., boundary_condition='MIFT',options=options)
|
|||
|
|
gfun3 = np.interp(time_superpos_sort*3600., time_req, gFunc.gFunc)
|
|||
|
|
# gfun3 = np.array(gfun2.gFunc)
|
|||
|
|
gfun3_unsort = gfun3[time_superpos_ind_rev]
|
|||
|
|
product = np.multiply(gfun3_unsort, Pfakt)
|
|||
|
|
R_g = abs(sum(product)/(2*pi*lamda))
|
|||
|
|
log(f" Rg = {R_g:5.3f} m.K/W")
|
|||
|
|
|
|||
|
|
# calculate maximum power from Heating and Cooling in respect to given Pfactor
|
|||
|
|
PKL = deltaT_KL / (R_g+R_bmf)
|
|||
|
|
PHZ = -PKL * Pfactor
|
|||
|
|
PHZ2 = deltaT_HZ / (R_g_2J+R_bmf)
|
|||
|
|
|
|||
|
|
if (abs(PHZ2) < abs(PHZ)):
|
|||
|
|
log(
|
|||
|
|
f" Heizlimit im Jahr 2 überschritten, reduziere Heizleistung um Faktor: {PHZ2/PHZ:.3f} !")
|
|||
|
|
log(f" Pot PHZ precalc = {PHZ:6.2f} W/m")
|
|||
|
|
PHZ = PHZ2
|
|||
|
|
PKL = -PHZ2/Pfactor
|
|||
|
|
else:
|
|||
|
|
log(f" Pot PHZ2 precalc = {PHZ2:6.2f} W/m")
|
|||
|
|
Tf_mean_max = T_g - (R_g+R_bmf)*PKL
|
|||
|
|
Tf_mean_min = T_g - (R_g_2J+R_bmf)*PHZ
|
|||
|
|
|
|||
|
|
Tb_pot_HZ = R_bmf * PHZ + Tf_mean_min
|
|||
|
|
Tb_pot_KL = R_bmf * PKL + Tf_mean_max
|
|||
|
|
|
|||
|
|
else: # more heating than cooling: Beginn with heating period
|
|||
|
|
log(" heating dominant")
|
|||
|
|
# calculate Resistance of Earth with g-functions after cooling period, year 2
|
|||
|
|
time_superpos_sort, time_superpos_ind_rev, Pfakt = prepare_times_for_gfunction_2J(
|
|||
|
|
Pfactor, Efactor, BS_HZ, BS_KL, years)
|
|||
|
|
# gfun2 = gt.gfunction.gFunction(network, alpha, time=time_superpos_sort*3600., boundary_condition='MIFT',options=options)
|
|||
|
|
gfun3 = np.interp(time_superpos_sort*3600., time_req, gFunc.gFunc)
|
|||
|
|
# gfun3 = np.array(gfun2.gFunc)
|
|||
|
|
gfun3_unsort = gfun3[time_superpos_ind_rev]
|
|||
|
|
product = np.multiply(gfun3_unsort, Pfakt)
|
|||
|
|
R_g_2J = abs(sum(product)/(2*pi*lamda))
|
|||
|
|
log(f" Rg_2J = {R_g_2J:5.3f} m.K/W")
|
|||
|
|
|
|||
|
|
# calculate Resistance of Earth with g-functions after first heating period, year 21
|
|||
|
|
time_superpos_sort, time_superpos_ind_rev, Pfakt, time_sort = prepare_times_for_gfunction_20J(
|
|||
|
|
Pfactor, Efactor, BS_HZ, BS_KL, years)
|
|||
|
|
# gfun2 = gt.gfunction.gFunction(network, alpha, time=time_superpos_sort*3600., boundary_condition='MIFT',options=options)
|
|||
|
|
# gfun3 = np.array(gfun2.gFunc)
|
|||
|
|
gfun3 = np.interp(time_superpos_sort*3600., time_req, gFunc.gFunc)
|
|||
|
|
gfun3_unsort = gfun3[time_superpos_ind_rev]
|
|||
|
|
product = np.multiply(gfun3_unsort, Pfakt)
|
|||
|
|
R_g = abs(sum(product)/(2*pi*lamda))
|
|||
|
|
log(f" Rg_20J = {R_g:5.3f} m.K/W")
|
|||
|
|
|
|||
|
|
# calculate maximum power from Heating and Cooling in respect to given Pfactor
|
|||
|
|
PHZ = deltaT_HZ / (R_g+R_bmf)
|
|||
|
|
PKL = -PHZ / Pfactor
|
|||
|
|
PKL2 = deltaT_KL / (R_g_2J+R_bmf)
|
|||
|
|
|
|||
|
|
if (abs(PKL2) < abs(PKL)):
|
|||
|
|
log(
|
|||
|
|
f" Kühllimit im Jahr 2 überschritten, reduziere Heizleistung um Faktor: {PKL2/PKL:.3f} !")
|
|||
|
|
log(f" Pot PHZ precalc = {PHZ:6.2f} W/m")
|
|||
|
|
PKL = PKL2
|
|||
|
|
PHZ = -PKL2*Pfactor
|
|||
|
|
|
|||
|
|
Tf_mean_max = T_g - (R_g_2J+R_bmf)*PKL
|
|||
|
|
Tf_mean_min = T_g - (R_g+R_bmf)*PHZ
|
|||
|
|
|
|||
|
|
Tb_pot_HZ = R_bmf * PHZ + Tf_mean_min
|
|||
|
|
Tb_pot_KL = R_bmf * PKL + Tf_mean_max
|
|||
|
|
|
|||
|
|
log(f" deltaT_HZ = {deltaT_HZ:6.2f} K")
|
|||
|
|
log(f" deltaT_KL = {deltaT_KL:6.2f} K")
|
|||
|
|
log(f" Efaktor = {Efactor:6.4f} ")
|
|||
|
|
log(f" Pfaktor = {Pfactor:6.4f} ")
|
|||
|
|
|
|||
|
|
log("\nresults from potential calculations:")
|
|||
|
|
log(f" Pot PHZ = {PHZ:6.2f} W/m")
|
|||
|
|
log(f" Pot PKL = {PKL:6.2f} W/m")
|
|||
|
|
log(f" Pot Tbmin = {Tb_pot_HZ:6.2f} °C")
|
|||
|
|
log(f" Pot Tbmax = {Tb_pot_KL:6.2f} °C")
|
|||
|
|
log(f" Pot Tf_mean_min = {Tf_mean_min:6.2f} °C")
|
|||
|
|
log(f" Pot Tf_mean_max = {Tf_mean_max:6.2f} °C")
|
|||
|
|
|
|||
|
|
# -------------------------------------------------------------------------
|
|||
|
|
# hourly Simulation
|
|||
|
|
# -------------------------------------------------------------------------
|
|||
|
|
|
|||
|
|
# Evaluate heat extraction rate
|
|||
|
|
P_tot = np.zeros(Nt)
|
|||
|
|
|
|||
|
|
T_b = np.zeros(Nt)
|
|||
|
|
# Tf_in = np.zeros(Nt)
|
|||
|
|
# Tf_out = np.zeros(Nt)
|
|||
|
|
Tf_m_sim = np.zeros(Nt)
|
|||
|
|
# Tf_m_sim2= np.zeros(Nt)
|
|||
|
|
COP_h = np.tile(np.nan, Nt)
|
|||
|
|
EER_h = np.tile(np.nan, Nt)
|
|||
|
|
SEER, SCOP = 0., 0.
|
|||
|
|
|
|||
|
|
for i, (t) in enumerate(time):
|
|||
|
|
# Increment time step by (1)
|
|||
|
|
# log(t/3600)
|
|||
|
|
hh = t/3600.
|
|||
|
|
xx = hh-np.floor(hh/8760.)*8760 # Stunde im Jahr
|
|||
|
|
|
|||
|
|
LoadAgg.next_time_step(t)
|
|||
|
|
|
|||
|
|
if xx <= x1:
|
|||
|
|
if (Efactor < 1):
|
|||
|
|
P_tot[i] = PKL
|
|||
|
|
else:
|
|||
|
|
P_tot[i] = PHZ
|
|||
|
|
elif xx <= x5:
|
|||
|
|
P_tot[i] = 0.
|
|||
|
|
elif xx <= x6:
|
|||
|
|
if (Efactor < 1):
|
|||
|
|
P_tot[i] = PHZ
|
|||
|
|
else:
|
|||
|
|
P_tot[i] = PKL
|
|||
|
|
elif xx <= x10:
|
|||
|
|
P_tot[i] = 0.
|
|||
|
|
elif xx <= x11:
|
|||
|
|
if (Efactor < 1):
|
|||
|
|
P_tot[i] = PKL
|
|||
|
|
else:
|
|||
|
|
P_tot[i] = PHZ
|
|||
|
|
else:
|
|||
|
|
log(f"WARNING at time {t/3600.}")
|
|||
|
|
P_tot[i] = 0.0
|
|||
|
|
|
|||
|
|
# Apply current load (in watts per meter of borehole)
|
|||
|
|
LoadAgg.set_current_load(P_tot[i])
|
|||
|
|
|
|||
|
|
# calculate borehole and fluid temperatures
|
|||
|
|
# Evaluate borehole wall temperature
|
|||
|
|
deltaT_b = LoadAgg.temporal_superposition()
|
|||
|
|
T_b[i] = T_g - deltaT_b
|
|||
|
|
|
|||
|
|
Tf_m_sim[i] = T_b[i]-P_tot[i]*R_bmf
|
|||
|
|
|
|||
|
|
if (PHZ > 0. and P_tot[i] == PHZ):
|
|||
|
|
Qc, Qe, Pel = get_COP_potential(T_radiator, Tf_m_sim[i]-1.5, 1.)
|
|||
|
|
COP_h[i] = Qc/Pel
|
|||
|
|
if (PKL < 0. and P_tot[i] == PKL):
|
|||
|
|
if (Tf_m_sim[i] > 19.5):
|
|||
|
|
# Qc, Qe, Pel = get_COP_potential(30., 19.5 , 1.)
|
|||
|
|
# EER_h[i] = Qe/Pel
|
|||
|
|
EER_h[i] = 6.9
|
|||
|
|
else:
|
|||
|
|
EER_h[i] = 20.
|
|||
|
|
|
|||
|
|
# Evaluate inlet fluid temperature (all boreholes are the same)
|
|||
|
|
# Tf_in[i] = network.get_network_inlet_temperature(P_tot[i]*divider, T_b[i], m_flow_network, cp_f, nSegments=1)
|
|||
|
|
|
|||
|
|
# Evaluate outlet fluid temperature
|
|||
|
|
# Tf_out[i] = network.get_network_outlet_temperature(Tf_in[i], T_b[i], m_flow_network, cp_f, nSegments=1)
|
|||
|
|
# Tf_m_sim2[i]=(Tf_in[i]+Tf_out[i])/2
|
|||
|
|
|
|||
|
|
last = T_b.size-1
|
|||
|
|
|
|||
|
|
Tm1 = min(Tf_m_sim[last], Tf_m_sim[xout])
|
|||
|
|
Tm2 = max(Tf_m_sim[last], Tf_m_sim[xout])
|
|||
|
|
Tb1 = min(T_b[last], T_b[xout])
|
|||
|
|
Tb2 = max(T_b[last], T_b[xout])
|
|||
|
|
if (np.any(np.isfinite(COP_h))):
|
|||
|
|
SCOP = np.nanmean(COP_h)
|
|||
|
|
if (np.any(np.isfinite(EER_h))):
|
|||
|
|
SEER = np.nanmean(EER_h)
|
|||
|
|
# np.savetxt(PNGPath+".txt",COP_h)
|
|||
|
|
|
|||
|
|
log("\nresults from hourly simulation run:")
|
|||
|
|
log(f" Sim Tbmin = {Tb1:6.2f} °C")
|
|||
|
|
log(f" Sim Tbmax = {Tb2:6.2f} °C")
|
|||
|
|
log(f" Sim Tf_mean_min = {Tm1:6.2f} °C")
|
|||
|
|
log(f" Sim Tf_mean_max = {Tm2:6.2f} °C")
|
|||
|
|
log(f" Sim SCOP = {SCOP:6.1f}")
|
|||
|
|
log(f" Sim SEER = {SEER:6.1f}")
|
|||
|
|
|
|||
|
|
# -------------------------------------------------------------------------
|
|||
|
|
# Plot hourly heat extraction rates and temperatures
|
|||
|
|
# -------------------------------------------------------------------------
|
|||
|
|
# plotgraf2(PNGPath,time/3600./24./365.,P_tot,"Betriebszeit [Jahre]","spezifische Sondenleistung [W/lm]",time/3600./24./365.,Tf_m_sim,"Betriebszeit [Jahre]","Temperatur [°C]", T_g, "Leistung pro Bohrmeter (Heizen = positiv, Kühlen/Regeneration = negativ)","Entwicklung der mittleren Fluidtemperatur der Sonden (Mittelwert Vor- und Rücklauf)", BS_HZ, BS_KL, RESULT)
|
|||
|
|
image_hash = plotgraf2(PNGPath, time/3600./24./365., P_tot, "Betriebszeit [Jahre]", "spezifische Sondenleistung [W/lm]", time/3600./24./365., Tf_m_sim, "Betriebszeit [Jahre]", "Temperatur [°C]",
|
|||
|
|
T_g, "Leistung pro Bohrmeter (Heizen = positiv, Kühlen/Regeneration = negativ)", "Entwicklung der mittleren Fluidtemperatur der Sonden (Mittelwert Vor- und Rücklauf)", BS_HZ, BS_KL, RESULT)
|
|||
|
|
|
|||
|
|
return PHZ, PKL, Efactor, SCOP, SEER, image_hash, image_hash_borefield
|
|||
|
|
|
|||
|
|
|
|||
|
|
def prepare_times_for_gfunction_20J(Pfactor, Efactor, BS_HZ, BS_KL, years):
|
|||
|
|
|
|||
|
|
# -------------------------------------------------------------------------
|
|||
|
|
# Potential
|
|||
|
|
# -------------------------------------------------------------------------
|
|||
|
|
# create array with Pfactors to multiply with gfunction values
|
|||
|
|
if Efactor >= 1:
|
|||
|
|
|
|||
|
|
Pfakt = np.array([1., -1., -1/Pfactor, 1/Pfactor])
|
|||
|
|
Pfakt = np.tile(Pfakt, int(years))
|
|||
|
|
# add last heating season
|
|||
|
|
Pfakt = np.concatenate((Pfakt, np.array([1.])), axis=0)
|
|||
|
|
|
|||
|
|
BS_first = BS_HZ
|
|||
|
|
BS_second = BS_KL
|
|||
|
|
else:
|
|||
|
|
Pfakt = np.array([-1., 1., Pfactor, -Pfactor])
|
|||
|
|
Pfakt = np.tile(Pfakt, int(years))
|
|||
|
|
# add last cooling season
|
|||
|
|
Pfakt = np.concatenate((Pfakt, np.array([-1.])), axis=0)
|
|||
|
|
|
|||
|
|
BS_first = BS_KL
|
|||
|
|
BS_second = BS_HZ
|
|||
|
|
|
|||
|
|
# create heating cooling time array
|
|||
|
|
t_HZ = np.arange(BS_first/2, 8760.0*(years+1)+BS_first/2, 8760)
|
|||
|
|
t_HZS = np.arange(8760/2-BS_second/2, 8760.0*(years), 8760)
|
|||
|
|
t_KL = np.arange(8760/2+BS_second/2, 8760.0*(years), 8760)
|
|||
|
|
t_KLS = np.arange(8760-BS_first/2, 8760.0*(years), 8760)
|
|||
|
|
|
|||
|
|
time2 = np.concatenate((t_HZ, t_HZS, t_KL, t_KLS), axis=0)
|
|||
|
|
# index array for sorting in increasing time order
|
|||
|
|
time_sort_ind = time2.argsort()
|
|||
|
|
# sort in increasing time order
|
|||
|
|
time_sort = time2[time_sort_ind]
|
|||
|
|
|
|||
|
|
# add zero value at beginning
|
|||
|
|
time_superpos = np.concatenate([[0.], time_sort])
|
|||
|
|
|
|||
|
|
Len = time2.size # length of arrays
|
|||
|
|
# get highest value = time point to evaluate power potential
|
|||
|
|
Last = time_sort[Len-1]
|
|||
|
|
|
|||
|
|
# prepare time points for g-function
|
|||
|
|
time_superpos = Last-time_superpos
|
|||
|
|
# remove last zero entry
|
|||
|
|
time_superpos = time_superpos[0:time_superpos.size-1]
|
|||
|
|
|
|||
|
|
# index array for sorting in increasing time order
|
|||
|
|
time_superpos_ind = time_superpos.argsort()
|
|||
|
|
# index array for sorting in decreasing time order
|
|||
|
|
time_superpos_ind_rev = time_superpos_ind.argsort()
|
|||
|
|
time_superpos_sort = time_superpos[time_superpos_ind]
|
|||
|
|
|
|||
|
|
# log(time_sort)
|
|||
|
|
# log(time_superpos)
|
|||
|
|
# log(Pfakt)
|
|||
|
|
# input()
|
|||
|
|
|
|||
|
|
return time_superpos_sort, time_superpos_ind_rev, Pfakt, time_sort
|
|||
|
|
|
|||
|
|
|
|||
|
|
def prepare_times_for_gfunction_2J(Pfactor, Efactor, BS_HZ, BS_KL, years):
|
|||
|
|
|
|||
|
|
# create array with Pfactors to multiply with gfunction values
|
|||
|
|
# plusminus = 1 for dominant heating and -1 for dominant cooling
|
|||
|
|
if Efactor >= 1:
|
|||
|
|
Pfakt = np.array([Pfactor, -Pfactor, -1., 1.,
|
|||
|
|
Pfactor, -Pfactor, -1.])
|
|||
|
|
BS_first = BS_HZ
|
|||
|
|
BS_second = BS_KL
|
|||
|
|
|
|||
|
|
else:
|
|||
|
|
Pfakt = np.array([-1/Pfactor, 1/Pfactor, 1., -1.,
|
|||
|
|
-1/Pfactor, 1/Pfactor, 1.])
|
|||
|
|
BS_first = BS_KL
|
|||
|
|
BS_second = BS_HZ
|
|||
|
|
# create same arrays to evaluate after end of heating phase in second year
|
|||
|
|
time2 = np.array([BS_first/2, 8760/2-BS_second/2, 8760/2+BS_second/2, 8760-BS_first/2,
|
|||
|
|
8760+BS_first/2, 8760+8760/2-BS_second/2, 8760+8760/2+BS_second/2])
|
|||
|
|
|
|||
|
|
# index array for sorting in increasing time order
|
|||
|
|
time_sort_ind = time2.argsort()
|
|||
|
|
# sort in increasing time order
|
|||
|
|
time_sort = time2[time_sort_ind]
|
|||
|
|
|
|||
|
|
time_superpos = np.concatenate([[0.], time_sort])
|
|||
|
|
Len = time_sort.size # length of arrays
|
|||
|
|
Last = time_sort[Len-1]
|
|||
|
|
time_superpos = Last - time_superpos
|
|||
|
|
time_superpos = time_superpos[0:time_superpos.size-1]
|
|||
|
|
|
|||
|
|
time_superpos_ind = time_superpos.argsort()
|
|||
|
|
time_superpos_ind_rev = time_superpos_ind.argsort()
|
|||
|
|
time_superpos_sort = time_superpos[time_superpos_ind]
|
|||
|
|
|
|||
|
|
# log(time_sort)
|
|||
|
|
# log(time_superpos)
|
|||
|
|
# log(Pfakt)
|
|||
|
|
# input()
|
|||
|
|
|
|||
|
|
return time_superpos_sort, time_superpos_ind_rev, Pfakt
|
|||
|
|
|
|||
|
|
|
|||
|
|
def plotborefield(boreField, pngname, RESULT):
|
|||
|
|
# Configure figure and axes
|
|||
|
|
if (RESULT):
|
|||
|
|
print(f"create plot {pngname}")
|
|||
|
|
fig = gt.boreholes.visualize_field(
|
|||
|
|
boreField, viewTop=True, view3D=False, labels=True) # , showTilt=True)
|
|||
|
|
|
|||
|
|
# save figure to binary stream and convert to base 64 string
|
|||
|
|
binary_stream = io.BytesIO()
|
|||
|
|
fig.savefig(binary_stream, format='png', dpi=120)
|
|||
|
|
binary_stream.seek(0)
|
|||
|
|
img_hash = base64.b64encode(binary_stream.read()).decode('utf-8')
|
|||
|
|
return img_hash
|
|||
|
|
|
|||
|
|
|
|||
|
|
def plotgraf2(pngname, x1, y1, x1_label, y1_label, x2, y2, x2_label, y2_label, Tg, title1, title2, BS_HZ, BS_KL, RESULT):
|
|||
|
|
# -------------------------------------------------------------------------
|
|||
|
|
# Plot hourly heat extraction rates and temperatures
|
|||
|
|
# -------------------------------------------------------------------------
|
|||
|
|
if (RESULT):
|
|||
|
|
print(f"create plot {pngname}")
|
|||
|
|
xmin = np.floor(min(x1))
|
|||
|
|
xmax = np.ceil(max(x1))
|
|||
|
|
xticks_major = np.arange(xmin+1, xmax)
|
|||
|
|
xticks_minor = np.arange(xmin+1, xmax, 0.5)
|
|||
|
|
|
|||
|
|
# Configure figure and axes
|
|||
|
|
plt.rc('figure')
|
|||
|
|
fig = gt.utilities._initialize_figure()
|
|||
|
|
# fig.suptitle(title, fontsize=8)
|
|||
|
|
ax1 = fig.add_subplot(211)
|
|||
|
|
# Axis labels
|
|||
|
|
ax1.set_xlabel(x1_label, fontsize=8)
|
|||
|
|
ax1.set_ylabel(y1_label, fontsize=8)
|
|||
|
|
ax1.set_title(title1, fontsize=9)
|
|||
|
|
ax1.set_xlim([xmin, xmax])
|
|||
|
|
# Specify tick label size
|
|||
|
|
ax1.set_xticks(xticks_minor, minor=True)
|
|||
|
|
ax1.set_xticks(xticks_major)
|
|||
|
|
|
|||
|
|
# gt.utilities._format_axes(ax1)
|
|||
|
|
#
|
|||
|
|
# Plot heat extraction rates
|
|||
|
|
lab = (
|
|||
|
|
f"Heizen {max(y1):.1f} W/m / Kühlen {min(y1):.1f} W/m \nHeizen {BS_HZ} h/a / Kühlen {BS_KL} h/a")
|
|||
|
|
ax1.plot(x1, y1, label=lab, color='teal')
|
|||
|
|
ax1.legend(fontsize=8, loc='lower center')
|
|||
|
|
|
|||
|
|
ax2 = fig.add_subplot(212)
|
|||
|
|
# Axis labels
|
|||
|
|
ax2.set_xlabel(x2_label, fontsize=8)
|
|||
|
|
ax2.set_ylabel(y2_label, fontsize=8)
|
|||
|
|
ax2.set_title(title2, fontsize=9)
|
|||
|
|
ax2.set_xlim([xmin, xmax])
|
|||
|
|
ax2.set_xticks(xticks_major)
|
|||
|
|
ax2.set_xticks(xticks_minor, minor=True)
|
|||
|
|
# gt.utilities._format_axes(ax2)
|
|||
|
|
|
|||
|
|
# Plot temperatures
|
|||
|
|
lab = (
|
|||
|
|
f"mittlere Sondenfluidtemperatur {min(y2[x2>1]):.1f} °C / {max(y2[x2>1]):.1f} °C")
|
|||
|
|
ax2.plot(x2, y2, label=lab, color='sienna')
|
|||
|
|
x3 = np.array([min(x2), max(x2)])
|
|||
|
|
y3 = np.array([Tg, Tg])
|
|||
|
|
ax2.plot(x3, y3, label="ungestörte Untergrundtemperatur",
|
|||
|
|
color='slategrey', linestyle='dotted')
|
|||
|
|
ax2.legend(fontsize=8, loc='upper center')
|
|||
|
|
plt.grid(True)
|
|||
|
|
|
|||
|
|
# Adjust to plot window
|
|||
|
|
plt.tight_layout()
|
|||
|
|
# fig.savefig(pngname)
|
|||
|
|
|
|||
|
|
# save figure to binary stream and encode it to base64 string
|
|||
|
|
binary_stream = io.BytesIO()
|
|||
|
|
fig.savefig(binary_stream, format='png', dpi=120)
|
|||
|
|
binary_stream.seek(0)
|
|||
|
|
img_hash = base64.b64encode(binary_stream.read()).decode('utf-8')
|
|||
|
|
return img_hash
|
|||
|
|
|
|||
|
|
|
|||
|
|
def get_COP_potential(T_sink_high, T_Source_low, P_SOW35):
|
|||
|
|
# polynomial fit from measurement data project Geofit and ZWEIFELDSPEICHER
|
|||
|
|
# Autor: Michael Lauermann (AIT 2015, 2022)
|
|||
|
|
# Qc = Q_heat condenser[kW]
|
|||
|
|
# Qe = Q_cool evaporator[kW]
|
|||
|
|
# Pel= electric power [kW]
|
|||
|
|
# a, b, c = linear coefficients
|
|||
|
|
# for delta T of heat sink = 5K, 10K and 15K
|
|||
|
|
# general formula: ( a + b * T_source_low + c * T_sink_high ) * P_SO-W35 / 44.14
|
|||
|
|
|
|||
|
|
def hpfunc(a, b, c, T_source_low, T_sink_high):
|
|||
|
|
result = (a + b * T_source_low + c * T_sink_high)
|
|||
|
|
return result
|
|||
|
|
|
|||
|
|
Qc_a_10K = 49.761136
|
|||
|
|
Qc_b_10K = 1.265924
|
|||
|
|
Qc_c_10K = -0.058782
|
|||
|
|
Qe_a_10K = 52.002178
|
|||
|
|
Qe_b_10K = 1.298574
|
|||
|
|
Qe_c_10K = -0.358732
|
|||
|
|
Pel_a_10K = -1.016948
|
|||
|
|
Pel_b_10K = 0.078546
|
|||
|
|
Pel_c_10K = 0.300213
|
|||
|
|
|
|||
|
|
Qc = hpfunc(Qc_a_10K, Qc_b_10K, Qc_c_10K, T_Source_low,
|
|||
|
|
T_sink_high) * P_SOW35 / 44.14
|
|||
|
|
Qe = hpfunc(Qe_a_10K, Qe_b_10K, Qe_c_10K, T_Source_low,
|
|||
|
|
T_sink_high) * P_SOW35 / 44.14
|
|||
|
|
Pel = hpfunc(Pel_a_10K, Pel_b_10K, Pel_c_10K,
|
|||
|
|
T_Source_low, T_sink_high) * P_SOW35 / 44.14
|
|||
|
|
# COP = Qc / Pel
|
|||
|
|
# EER = Qe / Pel
|
|||
|
|
|
|||
|
|
return Qc, Qe, Pel
|
|||
|
|
|
|||
|
|
|
|||
|
|
# Main function
|
|||
|
|
if __name__ == '__main__':
|
|||
|
|
freeze_support()
|
|||
|
|
import time as t
|
|||
|
|
seconds = t.time()
|
|||
|
|
|
|||
|
|
main()
|
|||
|
|
seconds = t.time()-seconds
|
|||
|
|
minutes = int(seconds/60)
|
|||
|
|
hours = int(minutes/60)
|
|||
|
|
sec = seconds-minutes*60-hours*3600
|
|||
|
|
# print(f"Total Running Time hhhh:mm:ss: {hours:04d}:{minutes:02d}:{sec:2.2f}")
|