From 8bc2ca9c85a824dc302c349167c345fdb157c813 Mon Sep 17 00:00:00 2001 From: Stephan Dueber <47361541+aflu@users.noreply.github.com> Date: Fri, 4 Feb 2022 11:13:26 +0100 Subject: [PATCH] Add files via upload --- Examples/Calc_FFT_Coeffs.py | 86 +++++++++++++++ Examples/Calc_SA_Coeffs.py | 87 +++++++++++++++ Examples/Calc_opt_periods.py | 23 ++++ Examples/example_error_coupling.py | 160 ++++++++++++++++++++++++++++ Examples/multiple_BHE_backward.py | 161 ++++++++++++++++++++++++++++ Examples/single_BHE_backward.py | 165 +++++++++++++++++++++++++++++ Examples/single_BHE_forward.py | 160 ++++++++++++++++++++++++++++ 7 files changed, 842 insertions(+) create mode 100644 Examples/Calc_FFT_Coeffs.py create mode 100644 Examples/Calc_SA_Coeffs.py create mode 100644 Examples/Calc_opt_periods.py create mode 100644 Examples/example_error_coupling.py create mode 100644 Examples/multiple_BHE_backward.py create mode 100644 Examples/single_BHE_backward.py create mode 100644 Examples/single_BHE_forward.py diff --git a/Examples/Calc_FFT_Coeffs.py b/Examples/Calc_FFT_Coeffs.py new file mode 100644 index 0000000..e7841b0 --- /dev/null +++ b/Examples/Calc_FFT_Coeffs.py @@ -0,0 +1,86 @@ +from __future__ import division + +import numpy as np +import matplotlib.pyplot as plt +import os +import sys +import timeit + +from scipy.optimize import curve_fit +from scipy.interpolate import interp1d as interp1d + +sourcepath = os.path.dirname(os.path.dirname(os.path.abspath(__file__))) +sys.path.append(sourcepath) + +import utilities + + +def calc_FFT_coeff(nt,gfunc,tlog,lin_gm,nBhe,nt_future,loads,dt): + + nsteps_total = loads.size + nt_future + Time = np.linspace(dt,nsteps_total*dt,nsteps_total) + + # convert loads + FFT_Loads = [] + for i in range(0,nBhe): + gload = np.concatenate([loads,np.zeros(nt_future)]) + gload[1:] -= np.roll(gload,1)[1:] + + FFT_Loads.append(utilities.FourierSci(gload)) + + # Temperatures at each borehole + T_borehole = np.ones([nBhe,Time.size])*10 + for i in range(0,nBhe): + for j in range(0,nBhe): + + # FFT Gfunc + fG = interp1d(tlog,gfunc) + FFT_Gfunc = utilities.FourierSci(fG(Time)) + + # Calc Tborehole + T_borehole[j,:] -= np.real(utilities.invFourierSci(FFT_Loads[i]*FFT_Gfunc)) + + return True + + +def main(): + + dt = 30 + nmax = 250000 + stepsize = 7500 + + steps = np.arange(stepsize,nmax,stepsize) + gfunc = np.random.rand(200)*10 + times_arr = np.zeros(steps.size) + + + for i in range(0,steps.size): + print(steps[i]) + nt = steps[i] + + tlog = np.geomspace(dt,3*nt*dt,200) + lin_g = interp1d(tlog,gfunc) + nBhe = 1 + loads = np.random.rand(nt) + nt_future = loads.size + + times_arr[i] = np.min(timeit.repeat(lambda :calc_FFT_coeff(nt,gfunc,tlog,lin_g,nBhe,nt_future,loads,dt),repeat = 20, number = 1,globals=globals()))/1 + + + print(steps) + popt_l, pcov = curve_fit(utilities.func_lin0,steps, times_arr) + + print('m: '+ str(popt_l[0])) + + plt.plot(steps,times_arr,label = 'measured',color = 'b') + plt.plot(steps,popt_l[0]*steps,label = 'lin fit', color = 'r') + + plt.ylabel('comp time') + plt.xlabel('steps') + plt.legend() + plt.show() + + +# Main function +if __name__ == '__main__': + main() \ No newline at end of file diff --git a/Examples/Calc_SA_Coeffs.py b/Examples/Calc_SA_Coeffs.py new file mode 100644 index 0000000..fc8f851 --- /dev/null +++ b/Examples/Calc_SA_Coeffs.py @@ -0,0 +1,87 @@ +from __future__ import division + +import numpy as np +import math + +import matplotlib.pyplot as plt +import timeit +from scipy.optimize import curve_fit +import os +import sys + +sourcepath = os.path.dirname(os.path.dirname(os.path.abspath(__file__))) +sys.path.append(sourcepath) + +import utilities + + + +def nsec(nt,psi_FFT,psi_sa): + return np.sqrt(psi_sa/psi_FFT*nt-1) + + +def calc_sa_coeff(T_borehole,loads,gMatrix,i,nSteps_soil,nBhe): + for j in range(0,nBhe): + for k in range(0,nBhe): + T_borehole[k][i:nSteps_soil] -= (loads[j][i]-loads[j][i-1]) * gMatrix[j,k,0:nSteps_soil-i] + +def main(): + + nt = 250000 + stepsize = 1000 + nBhe = 1 + nb = 130000 + + # random values for testing + T_borehole = np.random.rand(nBhe,nt) + gMatrix = np.random.rand(nBhe,nBhe,nt) + loads = np.random.rand(nBhe,nt) + times_arr = np.zeros(int(nt/stepsize)-1) + + for j in range(stepsize,nt,stepsize): + # measure computational times + times_vec = timeit.repeat(lambda: calc_sa_coeff(T_borehole,loads,gMatrix,j,nt,nBhe),repeat =50, number = 10,globals=globals()) + times_arr[int(j/stepsize)-1] = np.min(times_vec)/10 + + times_arr = np.flip(times_arr) + steps = np.arange(stepsize,nt,stepsize) + + # curve fitting + # fit nrnb + popt_l, pcov = curve_fit(utilities.func_lin,steps[steps>nb], times_arr[steps>nb]) + m_high = popt_l[0] + n_high = popt_l[1] + print('m_high:, n_high %d',m_high,n_high) + + # - - - - - - - - - - - - - - - - - - - - - - - + # plot results + # - - - - - - - - - - - - - - - - - - - - - - - + + fig = plt.figure(figsize=utilities.cm2inch(10, 6)) + ax1 = fig.add_subplot(111) + plt.subplots_adjust(bottom=0.22, top=0.90, left=0.18, right = 0.94) + + ax1.set_ylabel(r'$\mathrm{Computation \: time \: [} 10^{-3} \: \mathrm{s]}$') + ax1.set_xlabel(r'$\mathrm{remaining \: Steps \: [-]}$') + + ax1.plot(steps,1000*times_arr,label = r'$\mathrm{measured}$',color = 'k',linewidth = 0.7) + + ax1.plot(steps[stepsnb],1000*(steps[steps>nb]*m_high + n_high),linestyle = '--',color = 'grey',linewidth = 0.7,marker = 'o',markersize = 4,markevery=2,markerfacecolor='none') + + ax1.legend(loc=0,frameon=False) + plt.show() + + + + + +# Main function +if __name__ == '__main__': + main() \ No newline at end of file diff --git a/Examples/Calc_opt_periods.py b/Examples/Calc_opt_periods.py new file mode 100644 index 0000000..7642997 --- /dev/null +++ b/Examples/Calc_opt_periods.py @@ -0,0 +1,23 @@ +from __future__ import division +import os +import sys + +sourcepath = os.path.dirname(os.path.dirname(os.path.abspath(__file__))) +sys.path.append(sourcepath) + +import utilities + + +nb = 130000 # breakpoint calc time sa model +nt = 2109906 # total number of timesteps + +# computational time coefficients derived from measurements +psi_FFT = 1.58e-6 # coeff for FFT model +psi_sa_l = 5.8e-10 # coeff for sa model below breakpoint +psi_sa_h = 3.7e-9 # coeff for sa model above breakpoint +c_sa_l = 5.4e-6 # coeff for sa model below breakpoint +c_sa_h = -9.1e-5 # coeff for sa model above breakpoint + + +print(utilities.opt_nper(nt,psi_sa_l,psi_FFT,c_sa_l)) + diff --git a/Examples/example_error_coupling.py b/Examples/example_error_coupling.py new file mode 100644 index 0000000..4bf1ada --- /dev/null +++ b/Examples/example_error_coupling.py @@ -0,0 +1,160 @@ +from __future__ import division + +import numpy as np +import matplotlib.pyplot as plt +from scipy import interpolate +import os +import sys +import time as time +from multiprocessing import Pool + +sourcepath = os.path.dirname(os.path.dirname(os.path.abspath(__file__))) +sys.path.append(sourcepath) + +import bhe_models +import gfunctions +import utilities +import hybrid_model + + + +def main(): + + + BheData = { 'id' : 'BHE01', + 'type' : '2-U', + 'length': 100., # [m] length + 'diamB': 0.152, # [m] borehole diameter + 'lmG': 2.0, # [W/m/K] thermal conductivity grout + 'capG': 1000000, # [MJ/m³/K] volumetric heat capacity grout + 'lmP': 0.3, # [W/m/K] thermal conductivity pipe + 'odiamP': 0.032, # [m] outer diameter pipe + 'thickP': 0.0029, # [m] wall thickness pipe + 'distP': 0.04, # [m] pipe distance + 'densF': 1045, # [kg/m³] density fluid + 'dynviscF': 0.0035, # [m] dynamic viscosity fluid + 'lmF': 0.43, # [m] thermal conductivity fluid + 'capF': 3800000., # [m] volumetric heat capacity fluid + 'covH': 1., # [m] cover height + 'linkedHP': 'HP01', # [m] linked heat pump + 'xCoord' : 30000, # [m] x coordinate + 'yCoord' : 40000, # [m] y coordinate + 'Qf' : 1.92/3600., # [m³/s] flow rate + 'Tundist': 12.0, # [°C] undisturbed ground temperature + 'lm' : 2.3, # [W/m/K] thermal conductivity ground + 'Cm' : 2300000., # [MJ/m³/K] volumetric heat capacity ground + } + + + sim_setup = { 'nz': 5, # bhe cells vertical + 'dt_bhe': 30, # bhe time step + 'dt_soil': 3600, # soil time step + 'nt_part': 24, # soil time steps per period + 'n_parts': 1, # number of periods + 'error': 0.001, + 'type' : 'forward' # euler scheme + } + + + + #create random Tin and flow rate + M_Tin = np.ones(sim_setup['nt_part']*sim_setup['n_parts'])*50 + M_qf = np.ones(sim_setup['nt_part']*sim_setup['n_parts'])*BheData['Qf'] + + + x_pos = np.array([1]) + y_pos = np.array([1]) + nBhe = x_pos.size + + # ------------------------------------------------------------------------- + # calc gfunction and create time array tlog + # ------------------------------------------------------------------------- + + tstart = sim_setup['dt_soil'] + tstop = sim_setup['dt_soil']*sim_setup['nt_part']*sim_setup['n_parts']+1 + ntgfunc = 200 + + tlog = np.geomspace(tstart, tstop, ntgfunc) + gfuncs = gfunctions.g_Matrix(x_pos,y_pos,tlog,BheData) + + + gfunc_data = { 'gfuncs' : gfuncs, + 'gtime' : tlog + } + + # gfunc for semi-analytical model + tlin = np.linspace(sim_setup['dt_soil'],sim_setup['dt_soil']*sim_setup['nt_part'],sim_setup['nt_part']) # time array soil + gMatrix = np.zeros([nBhe,nBhe,sim_setup['nt_part']]) # matrix with gfuncs + for i in range(0,nBhe): + for j in range(0,nBhe): + fG = interpolate.interp1d(gfunc_data['gtime'],gfunc_data['gfuncs'][i,j]) + gMatrix[i,j,:] = fG(tlin) # gfunc matrix interpolated to linear time grid + + # ------------------------------------------------------------------------- + # setup for simulation + # ------------------------------------------------------------------------- + + res_loads = np.zeros([nBhe,sim_setup['nt_part']*sim_setup['n_parts']]) + res_Tins = np.zeros([nBhe,sim_setup['nt_part']*sim_setup['n_parts']]) + res_Touts = np.zeros([nBhe,sim_setup['nt_part']*sim_setup['n_parts']]) + + M_Tin = np.array([M_Tin]) + M_qf = np.array([M_qf]) + + T_borehole = np.ones([nBhe,sim_setup['nt_part']])*BheData['Tundist'] + T_borehole_ini = np.ones(nBhe)*BheData['Tundist'] + + BHEs = hybrid_model.init_BHEs(nBhe ,BheData,sim_setup['dt_bhe'],sim_setup['nz'], sim_setup['type']) + + # num first part + start = time.time() + (res_Tins[:,0:sim_setup['nt_part']], + res_Touts[:,0:sim_setup['nt_part']], + res_loads[:,0:sim_setup['nt_part']]) = hybrid_model.calc_sa_sec_U_forw(sim_setup, + gMatrix, + M_Tin[:,0:sim_setup['nt_part']], + M_qf[:,0:sim_setup['nt_part']], + T_borehole,BHEs,T_borehole_ini) + + for p in range(1,sim_setup['n_parts']): + + # analytical firt part + T_borehole = hybrid_model.calc_FFT_sec(gfunc_data['gfuncs'], + gfunc_data['gtime'], + res_loads[:,0:p*sim_setup['nt_part']], + sim_setup['nt_part'], + sim_setup['dt_soil'], + nBhe,BheData['Tundist']) + + T_borehole_ini = [T_borehole[k,p*sim_setup['nt_part']-1] for k in range(0,nBhe)] + + # num second part + (res_Tins[:,p*sim_setup['nt_part']:(p+1)*sim_setup['nt_part']], + res_Touts[:,p*sim_setup['nt_part']:(p+1)*sim_setup['nt_part']], + res_loads[:,p*sim_setup['nt_part']:(p+1)*sim_setup['nt_part']]) = hybrid_model.calc_sa_sec_U_forw(sim_setup, + gMatrix, + M_Tin[:,p*sim_setup['nt_part']:(p+1)*sim_setup['nt_part']], + M_qf[:,p*sim_setup['nt_part']:(p+1)*sim_setup['nt_part']], + T_borehole[:,p*sim_setup['nt_part']:],BHEs,T_borehole_ini) + + print('time: ',time.time() - start) + + + + # ------------------------------------------------------------------------- + # plot results + # ------------------------------------------------------------------------- + + plt.plot(M_Tin[0],label = 'Tin') + plt.plot(res_Touts[0],label = 'Tout') + plt.legend() + plt.show() + + + + + + +# Main function +if __name__ == '__main__': + main() diff --git a/Examples/multiple_BHE_backward.py b/Examples/multiple_BHE_backward.py new file mode 100644 index 0000000..7cac95d --- /dev/null +++ b/Examples/multiple_BHE_backward.py @@ -0,0 +1,161 @@ +from __future__ import division + +import numpy as np +import matplotlib.pyplot as plt +from scipy import interpolate +import os +import sys +import time as time +from multiprocessing import Pool + +sourcepath = os.path.dirname(os.path.dirname(os.path.abspath(__file__))) +sys.path.append(sourcepath) + +import bhe_models +import gfunctions +import utilities +import hybrid_model + + + +def main(): + + + BheData = { 'id' : 'BHE01', + 'type' : '2-U', + 'length': 100., # [m] length + 'diamB': 0.152, # [m] borehole diameter + 'lmG': 2.0, # [W/m/K] thermal conductivity grout + 'capG': 1000000, # [MJ/m³/K] volumetric heat capacity grout + 'lmP': 0.3, # [W/m/K] thermal conductivity pipe + 'odiamP': 0.032, # [m] outer diameter pipe + 'thickP': 0.0029, # [m] wall thickness pipe + 'distP': 0.04, # [m] pipe distance + 'densF': 1045, # [kg/m³] density fluid + 'dynviscF': 0.0035, # [m] dynamic viscosity fluid + 'lmF': 0.43, # [m] thermal conductivity fluid + 'capF': 3800000., # [m] volumetric heat capacity fluid + 'covH': 1., # [m] cover height + 'linkedHP': 'HP01', # [m] linked heat pump + 'xCoord' : 30000, # [m] x coordinate + 'yCoord' : 40000, # [m] y coordinate + 'Qf' : 1.92/3600., # [m³/s] flow rate + 'Tundist': 12.0, # [°C] undisturbed ground temperature + 'lm' : 2.3, # [W/m/K] thermal conductivity ground + 'Cm' : 2300000., # [MJ/m³/K] volumetric heat capacity ground + } + + + sim_setup = { 'nz': 5, # bhe cells vertical + 'dt_bhe': 30, # bhe time step + 'dt_soil': 30, # soil time step + 'nt_part': 14880, # soil time steps per period + 'n_parts': 6, # number of periods + 'error': 0.001, + 'type' : 'backward' # euler scheme + } + + + + + x_pos = np.array([1,2,3,4]) + y_pos = np.array([1,1,1,1]) + nBhe = x_pos.size + + #create random Tin and flow rate + Tin = np.ones(sim_setup['nt_part']*sim_setup['n_parts']) + M_Tin = [Tin*20,Tin*20,Tin*20,Tin*20,] + qf = np.ones(sim_setup['nt_part']*sim_setup['n_parts'])*BheData['Qf'] + M_qf = [qf,qf,qf,qf] + + # ------------------------------------------------------------------------- + # calc gfunction and create time array tlog + # ------------------------------------------------------------------------- + + tstart = sim_setup['dt_soil'] + tstop = sim_setup['dt_soil']*sim_setup['nt_part']*sim_setup['n_parts']+1 + ntgfunc = 200 + + tlog = np.geomspace(tstart, tstop, ntgfunc) + gfuncs = gfunctions.g_Matrix(x_pos,y_pos,tlog,BheData) + + + gfunc_data = { 'gfuncs' : gfuncs, + 'gtime' : tlog + } + + # gfunc for semi-analytical model + tlin = np.linspace(sim_setup['dt_soil'],sim_setup['dt_soil']*sim_setup['nt_part'],sim_setup['nt_part']) # time array soil + gMatrix = np.zeros([nBhe,nBhe,sim_setup['nt_part']]) # matrix with gfuncs + for i in range(0,nBhe): + for j in range(0,nBhe): + fG = interpolate.interp1d(gfunc_data['gtime'],gfunc_data['gfuncs'][i,j]) + gMatrix[i,j,:] = fG(tlin) # gfunc matrix interpolated to linear time grid + + # ------------------------------------------------------------------------- + # setup for simulation + # ------------------------------------------------------------------------- + + res_loads = np.zeros([nBhe,sim_setup['nt_part']*sim_setup['n_parts']]) + res_Tins = np.zeros([nBhe,sim_setup['nt_part']*sim_setup['n_parts']]) + res_Touts = np.zeros([nBhe,sim_setup['nt_part']*sim_setup['n_parts']]) + + + T_borehole = np.ones([nBhe,sim_setup['nt_part']])*BheData['Tundist'] + T_borehole_ini = np.ones(nBhe)*BheData['Tundist'] + + BHEs = hybrid_model.init_BHEs(nBhe ,BheData,sim_setup['dt_bhe'],sim_setup['nz'], sim_setup['type']) + + # num first part + start = time.time() + (res_Tins[:,0:sim_setup['nt_part']], + res_Touts[:,0:sim_setup['nt_part']], + res_loads[:,0:sim_setup['nt_part']]) = hybrid_model.calc_sa_sec_U_backw(sim_setup, + gMatrix, + [M_Tin[k][0:sim_setup['nt_part']] for k in range(0,nBhe)], + [M_qf[k][0:sim_setup['nt_part']] for k in range(0,nBhe)], + T_borehole,BHEs,T_borehole_ini) + + for p in range(1,sim_setup['n_parts']): + + # analytical firt part + T_borehole = hybrid_model.calc_FFT_sec(gfunc_data['gfuncs'], + gfunc_data['gtime'], + res_loads[:,0:p*sim_setup['nt_part']], + sim_setup['nt_part'], + sim_setup['dt_soil'], + nBhe,BheData['Tundist']) + + T_borehole_ini = [T_borehole[k,p*sim_setup['nt_part']-1] for k in range(0,nBhe)] + + # num second part + (res_Tins[:,p*sim_setup['nt_part']:(p+1)*sim_setup['nt_part']], + res_Touts[:,p*sim_setup['nt_part']:(p+1)*sim_setup['nt_part']], + res_loads[:,p*sim_setup['nt_part']:(p+1)*sim_setup['nt_part']]) = hybrid_model.calc_sa_sec_U_backw(sim_setup, + gMatrix, + [M_Tin[k][p*sim_setup['nt_part']:(p+1)*sim_setup['nt_part']] for k in range(0,nBhe)], + [M_qf[k][p*sim_setup['nt_part']:(p+1)*sim_setup['nt_part']] for k in range(0,nBhe)], + T_borehole[:,p*sim_setup['nt_part']:],BHEs,T_borehole_ini) + + print('time: ',time.time() - start) + + + + # ------------------------------------------------------------------------- + # plot results + # ------------------------------------------------------------------------- + + for i in range(0,nBhe): + plt.plot(M_Tin[i],label = 'Tin_'+str(i)) + plt.plot(res_Touts[i],label = 'Tout_'+str(i)) + plt.legend() + plt.show() + + + + + + +# Main function +if __name__ == '__main__': + main() diff --git a/Examples/single_BHE_backward.py b/Examples/single_BHE_backward.py new file mode 100644 index 0000000..0ab3860 --- /dev/null +++ b/Examples/single_BHE_backward.py @@ -0,0 +1,165 @@ +from __future__ import division + +import numpy as np +import matplotlib.pyplot as plt +from scipy import interpolate +import os +import sys +import time as time +from multiprocessing import Pool + +sourcepath = os.path.dirname(os.path.dirname(os.path.abspath(__file__))) +sys.path.append(sourcepath) + +import bhe_models +import gfunctions +import utilities +import hybrid_model + + + +def main(): + + + BheData = { 'id' : 'BHE01', + 'type' : '1-U', + 'length': 100., # [m] length + 'diamB': 0.152, # [m] borehole diameter + 'lmG': 2.0, # [W/m/K] thermal conductivity grout + 'capG': 1000000, # [J/m³/K] volumetric heat capacity grout + 'lmP': 0.3, # [W/m/K] thermal conductivity pipe + 'odiamP': 0.032, # [m] outer diameter pipe + 'thickP': 0.0029, # [m] wall thickness pipe + 'distP': 0.04, # [m] pipe distance + 'densF': 1045, # [kg/m³] density fluid + 'dynviscF': 0.0035, # [m] dynamic viscosity fluid + 'lmF': 0.43, # [m] thermal conductivity fluid + 'capF': 3800000., # [m] volumetric heat capacity fluid + 'covH': 1., # [m] cover height + 'linkedHP': 'HP01', # [m] linked heat pump + 'xCoord' : 3 , # [m] x coordinate + 'yCoord' : 4 , # [m] y coordinate + 'Qf' : 1.92/3600., # [m³/s] flow rate + 'Tundist': 12.0, # [°C] undisturbed ground temperature + 'lm' : 2.3, # [W/m/K] thermal conductivity ground + 'Cm' : 2300000., # [J/m³/K] volumetric heat capacity ground + } + + + sim_setup = { 'nz': 100, # bhe cells vertical + 'dt_bhe': 5, # bhe time step + 'dt_soil': 5, # soil time step + 'nt_part': 17280, # soil time steps per period + 'n_parts': 1, # number of periods + 'error': 0.001, + 'type' : 'backward' # euler scheme + } + + + + #create random Tin and flow rate + M_Tin = np.ones(sim_setup['nt_part']*sim_setup['n_parts'])*30 + M_qf = np.ones(sim_setup['nt_part']*sim_setup['n_parts'])*BheData['Qf'] + + + x_pos = np.array([1]) + y_pos = np.array([1]) + nBhe = x_pos.size + + # ------------------------------------------------------------------------- + # calc gfunction and create time array tlog + # ------------------------------------------------------------------------- + + tstart = 1 + tstop = sim_setup['dt_soil']*sim_setup['nt_part']*sim_setup['n_parts']+1 + ntgfunc = 200 + + tlog = np.geomspace(tstart, tstop, ntgfunc) + gfuncs = gfunctions.g_Matrix(x_pos,y_pos,tlog,BheData) + + + gfunc_data = { 'gfuncs' : gfuncs, + 'gtime' : tlog + } + + # gfunc for semi-analytical model + tlin = np.linspace(sim_setup['dt_soil'],sim_setup['dt_soil']*sim_setup['nt_part'],sim_setup['nt_part']) # time array soil + gMatrix = np.zeros([nBhe,nBhe,sim_setup['nt_part']]) # matrix with gfuncs + for i in range(0,nBhe): + for j in range(0,nBhe): + fG = interpolate.interp1d(gfunc_data['gtime'],gfunc_data['gfuncs'][i,j]) + gMatrix[i,j,:] = fG(tlin) # gfunc matrix interpolated to linear time grid + + # ------------------------------------------------------------------------- + # setup for simulation + # ------------------------------------------------------------------------- + + res_loads = np.zeros([nBhe,sim_setup['nt_part']*sim_setup['n_parts']]) + res_Tins = np.zeros([nBhe,sim_setup['nt_part']*sim_setup['n_parts']]) + res_Touts = np.zeros([nBhe,sim_setup['nt_part']*sim_setup['n_parts']]) + + M_Tin = np.array([M_Tin]) + M_qf = np.array([M_qf]) + + T_borehole = np.ones([nBhe,sim_setup['nt_part']])*BheData['Tundist'] + T_borehole_ini = np.ones(nBhe)*BheData['Tundist'] + + BHEs = hybrid_model.init_BHEs(nBhe ,BheData,sim_setup['dt_bhe'],sim_setup['nz'], sim_setup['type']) + + # num first part + start = time.time() + (res_Tins[:,0:sim_setup['nt_part']], + res_Touts[:,0:sim_setup['nt_part']], + res_loads[:,0:sim_setup['nt_part']]) = hybrid_model.calc_sa_sec_U_backw(sim_setup, + gMatrix, + M_Tin[:,0:sim_setup['nt_part']], + M_qf[:,0:sim_setup['nt_part']], + T_borehole,BHEs,T_borehole_ini) + + for p in range(1,sim_setup['n_parts']): + + # analytical firt part + T_borehole = hybrid_model.calc_FFT_sec(gfunc_data['gfuncs'], + gfunc_data['gtime'], + res_loads[:,0:p*sim_setup['nt_part']], + sim_setup['nt_part'], + sim_setup['dt_soil'], + nBhe,BheData['Tundist']) + + T_borehole_ini = [T_borehole[k,p*sim_setup['nt_part']-1] for k in range(0,nBhe)] + + # num second part + (res_Tins[:,p*sim_setup['nt_part']:(p+1)*sim_setup['nt_part']], + res_Touts[:,p*sim_setup['nt_part']:(p+1)*sim_setup['nt_part']], + res_loads[:,p*sim_setup['nt_part']:(p+1)*sim_setup['nt_part']]) = hybrid_model.calc_sa_sec_U_backw(sim_setup, + gMatrix, + M_Tin[:,p*sim_setup['nt_part']:(p+1)*sim_setup['nt_part']], + M_qf[:,p*sim_setup['nt_part']:(p+1)*sim_setup['nt_part']], + T_borehole[:,p*sim_setup['nt_part']:],BHEs,T_borehole_ini) + + print('time: ',time.time() - start) + + + + # ------------------------------------------------------------------------- + # plot results + # ------------------------------------------------------------------------- + + plt.plot(M_Tin[0],label = 'Tin') + plt.plot(res_Touts[0],label = 'Tout') + plt.legend() + plt.show() + + plt.plot(BHEs[0].Tf_in) + plt.plot(BHEs[0].Tf_out) + plt.show() + + + + + + + +# Main function +if __name__ == '__main__': + main() diff --git a/Examples/single_BHE_forward.py b/Examples/single_BHE_forward.py new file mode 100644 index 0000000..f95787e --- /dev/null +++ b/Examples/single_BHE_forward.py @@ -0,0 +1,160 @@ +from __future__ import division + +import numpy as np +import matplotlib.pyplot as plt +from scipy import interpolate +import os +import sys +import time as time +from multiprocessing import Pool + +sourcepath = os.path.dirname(os.path.dirname(os.path.abspath(__file__))) +sys.path.append(sourcepath) + +import bhe_models +import gfunctions +import utilities +import hybrid_model + + + +def main(): + + + BheData = { 'id' : 'BHE01', + 'type' : '2-U', + 'length': 100., # [m] length + 'diamB': 0.152, # [m] borehole diameter + 'lmG': 2.0, # [W/m/K] thermal conductivity grout + 'capG': 1000000, # [MJ/m³/K] volumetric heat capacity grout + 'lmP': 0.3, # [W/m/K] thermal conductivity pipe + 'odiamP': 0.032, # [m] outer diameter pipe + 'thickP': 0.0029, # [m] wall thickness pipe + 'distP': 0.04, # [m] pipe distance + 'densF': 1045, # [kg/m³] density fluid + 'dynviscF': 0.0035, # [m] dynamic viscosity fluid + 'lmF': 0.43, # [m] thermal conductivity fluid + 'capF': 3800000., # [m] volumetric heat capacity fluid + 'covH': 1., # [m] cover height + 'linkedHP': 'HP01', # [m] linked heat pump + 'xCoord' : 30000, # [m] x coordinate + 'yCoord' : 40000, # [m] y coordinate + 'Qf' : 1.92/3600., # [m³/s] flow rate + 'Tundist': 12.0, # [°C] undisturbed ground temperature + 'lm' : 2.3, # [W/m/K] thermal conductivity ground + 'Cm' : 2300000., # [MJ/m³/K] volumetric heat capacity ground + } + + + sim_setup = { 'nz': 5, # bhe cells vertical + 'dt_bhe': 30, # bhe time step + 'dt_soil': 30, # soil time step + 'nt_part': 14880, # soil time steps per period + 'n_parts': 6, # number of periods + 'error': 0.001, + 'type' : 'forward' # euler scheme + } + + + + #create random Tin and flow rate + M_Tin = np.ones(sim_setup['nt_part']*sim_setup['n_parts'])*50 + M_qf = np.ones(sim_setup['nt_part']*sim_setup['n_parts'])*BheData['Qf'] + + + x_pos = np.array([1]) + y_pos = np.array([1]) + nBhe = x_pos.size + + # ------------------------------------------------------------------------- + # calc gfunction and create time array tlog + # ------------------------------------------------------------------------- + + tstart = sim_setup['dt_soil'] + tstop = sim_setup['dt_soil']*sim_setup['nt_part']*sim_setup['n_parts']+1 + ntgfunc = 200 + + tlog = np.geomspace(tstart, tstop, ntgfunc) + gfuncs = gfunctions.g_Matrix(x_pos,y_pos,tlog,BheData) + + + gfunc_data = { 'gfuncs' : gfuncs, + 'gtime' : tlog + } + + # gfunc for semi-analytical model + tlin = np.linspace(sim_setup['dt_soil'],sim_setup['dt_soil']*sim_setup['nt_part'],sim_setup['nt_part']) # time array soil + gMatrix = np.zeros([nBhe,nBhe,sim_setup['nt_part']]) # matrix with gfuncs + for i in range(0,nBhe): + for j in range(0,nBhe): + fG = interpolate.interp1d(gfunc_data['gtime'],gfunc_data['gfuncs'][i,j]) + gMatrix[i,j,:] = fG(tlin) # gfunc matrix interpolated to linear time grid + + # ------------------------------------------------------------------------- + # setup for simulation + # ------------------------------------------------------------------------- + + res_loads = np.zeros([nBhe,sim_setup['nt_part']*sim_setup['n_parts']]) + res_Tins = np.zeros([nBhe,sim_setup['nt_part']*sim_setup['n_parts']]) + res_Touts = np.zeros([nBhe,sim_setup['nt_part']*sim_setup['n_parts']]) + + M_Tin = np.array([M_Tin]) + M_qf = np.array([M_qf]) + + T_borehole = np.ones([nBhe,sim_setup['nt_part']])*BheData['Tundist'] + T_borehole_ini = np.ones(nBhe)*BheData['Tundist'] + + BHEs = hybrid_model.init_BHEs(nBhe ,BheData,sim_setup['dt_bhe'],sim_setup['nz'], sim_setup['type']) + + # num first part + start = time.time() + (res_Tins[:,0:sim_setup['nt_part']], + res_Touts[:,0:sim_setup['nt_part']], + res_loads[:,0:sim_setup['nt_part']]) = hybrid_model.calc_sa_sec_U_forw(sim_setup, + gMatrix, + M_Tin[:,0:sim_setup['nt_part']], + M_qf[:,0:sim_setup['nt_part']], + T_borehole,BHEs,T_borehole_ini) + + for p in range(1,sim_setup['n_parts']): + + # analytical firt part + T_borehole = hybrid_model.calc_FFT_sec(gfunc_data['gfuncs'], + gfunc_data['gtime'], + res_loads[:,0:p*sim_setup['nt_part']], + sim_setup['nt_part'], + sim_setup['dt_soil'], + nBhe,BheData['Tundist']) + + T_borehole_ini = [T_borehole[k,p*sim_setup['nt_part']-1] for k in range(0,nBhe)] + + # num second part + (res_Tins[:,p*sim_setup['nt_part']:(p+1)*sim_setup['nt_part']], + res_Touts[:,p*sim_setup['nt_part']:(p+1)*sim_setup['nt_part']], + res_loads[:,p*sim_setup['nt_part']:(p+1)*sim_setup['nt_part']]) = hybrid_model.calc_sa_sec_U_forw(sim_setup, + gMatrix, + M_Tin[:,p*sim_setup['nt_part']:(p+1)*sim_setup['nt_part']], + M_qf[:,p*sim_setup['nt_part']:(p+1)*sim_setup['nt_part']], + T_borehole[:,p*sim_setup['nt_part']:],BHEs,T_borehole_ini) + + print('time: ',time.time() - start) + + + + # ------------------------------------------------------------------------- + # plot results + # ------------------------------------------------------------------------- + + plt.plot(M_Tin[0],label = 'Tin') + plt.plot(res_Touts[0],label = 'Tout') + plt.legend() + plt.show() + + + + + + +# Main function +if __name__ == '__main__': + main()