g-Functions for variable fluid mass flow rate and reversible flow direction

[fuma01]

I have more a question, which may result in an enhancment: Is it already possible to simulate a borefield with given heat extraction rate and different network flow values per time step?

my script is based on the example: https://pygfunction.readthedocs.io/en/stable/examples/fluid_temperature_multiple_boreholes.html#simulation-of-fluid-temperatures-in-a-field-of-multiple-boreholes

until now I used constant network flow, and the algorithm works really well eg: _m_flow_network_w = 0.35 N1_w N2w

(milesteps of the script): _LoadAgg_w = gt.load_aggregation.ClaessonJaved(dt, tmax) boreField_w = gt.boreholes.rectangle_field(N1_w, N2_w, B, B, H_w, D, r_b) network_w = gt.networks.Network(boreField_w, UTubes_w, m_flow_network=m_flow_network_w, cp_f=cp_fluid) gFunc_w = gt.gfunction.gFunction(network_w, alpha, time=time_req, boundary_condition='MIFT', options=options) LoadAgg_w.initialize(gFunc_w.gFunc/(2pik_s)) for i, (t) in enumerate(time): LoadAgg_w.next_time_step(t) Tf_in_w[i] = network_w.get_network_inlet_temperature(P_tot_w[i], Tb_w[i], m_flow_network_w, cpfluid, nSegments=1)

Now, I want to implement 3 different mass flow rates (e.g. 0.1, 0.25 and 0.5 kg/s per borehole), depending on actual power value (P_tot_w). The output should be inlet and outlet temperature of the field. as the gfunctions itself and hence the initialization of the Load Aggregation is dependent on the network flow, I think it do not work properly only changing the flow value in "_get_network_inlettemperature"... is there already a possibility to use several network flow steps? thank you

[MassimoCimmino]

Hello @fuma01.

We found in an earlier study that g-functions are not strongly dependent on mass flow rate (as long as it is above a threshold value) [1]. For parallel-connected boreholes, the variation is often negligible [2]. This means that the same g-function (at nominal mass flow rate) can be used for a simulation with variable mass flow rate.

One way to verify is to evaluate and compare g-functions at your mininum an maximum mass flow rates up to a time equal to the maximum simulation time. If there is no significant difference, then all you need is to use the variable mass flow rate when calling network.get_<> class methods.

It is also possible to consider different g-functions with different mass flow rates using a non-stationary combination [3]. This is not supported by pygfunction and would need to be adapted to time-marching simulations before it can be implemented.

---

[1] Cimmino, M. (2019). Semi-Analytical Method for g-Function Calculation of bore fields with series-and parallel-connected boreholes. Science and Technology for the Built Environment, 25(8), 1007-1022. [2] Cimmino, M. (2015). The effects of borehole thermal resistances and fluid flow rate on the g-functions of geothermal bore fields. International Journal of Heat and Mass Transfer, 91, 1119-1127. [3] Beaudry, G., Pasquier, P., & Marcotte, D. (2021). A fast convolution-based method to simulate time-varying flow rates in closed-loop and standing column well ground heat exchangers. Renewable Energy, 174, 55-72.

[fuma01]

Dear Massimo,

Thank you for the quick response! Indeed the g-function is not very sensitiv to flow values, at least if the field stays in turbulent mode (< 0.6 % in my simulation) but a bit more in low flow mode (up to 10 % difference) I testet it with three flow rates (0.1, 0.3 and 0.5 kg/s), which represents my target range and plottet the g-function and the difference to 0.5 kg/s, related to Eskilson time: I use hourly time steps up to 20 years. rectangular borefield 7x9 with 150 m length and 5 m spacing

now I am thinking about, if the error will be significant when using variable mass flow rate just at"network.get_" method, as suggested? also, I am wondering, that the computation time needs 5 times longer with switching flow rates @ network.get.

regards from Vienna, Martin

[MassimoCimmino]

Accuracy

The accuracy of your simulation will depend on your simulation objectives and how often the flow rate is 0.1 kg/s.

A difference in g-functions at low ln(t/ts) translates to an error in borehole wall temperature variations in the short-term (i.e. how quickly the temperature varies following a heat pulse). The variation of fluid temperatures is a combination of two factors: (1) the temperature variation at the borehole wall, and (2) the effective borehole thermal resistance. The second also varies with mass flow and is captured by the Network class.

If you are simulating for heat pump energy use, then an error in temperature affects the performance of the heat pump. If you are simulating to size the field in peak conditions then I figure the peak coincides with the maximum mass flow rate. The magnitude of these impacts strongly depend on your simulation case.

Here, I would suggest running and comparing simulations using g-functions evaluated with the 0.1 and 0.5 kg/s mass flow rates. This would give an estimation of the expected magnitude of the errors caused by using only one g-function.

Factor in that g-functions do not capture short-term effects (i.e. the heat capacity of the grout and the travel time of the fluid). This is a limitation of steady-state resistance models such as the one implemented in the Network class.

Calculation time

When the mass flow rate changes, the Network class re-evaluates all borehole coefficients. the class only remembers coefficients for the latest mass flow rate.

At the moment, this is done to all boreholes individually even if all borehole are geometrically identical.

[fuma01]

in my current investigation, I have 5 heat pumps parallel to one (or two) BHE fields. Input is the heating/cooling power demand of the building and solar module. Then the field size is adapted, that the fluid temperatures are in specified range. The output is the electric power of the heat pumps (calculated with characteristic curve) in proportion to solar and BHE power and relation to field size. In dependence of the requested power load the needed number of heat pumps are activated. In current simulation the mass flow rate is set constant, hence the delta T is varying. The idea was to set a variably flow rate and get more or less constant delta T, which will also affect input temperature of the heat pump and hence the COP (or the needed electric power).

now, I understand, that the mass flow rate is implementet in the g-function, which is calculated in advance. By switching laminar/turbulent flow mode, this will probably not be neclectable. As a workaround, I will make two simulations at minimum and maximum constant flow rate, to see the efffects. The increase of calculation time by changing flow rate later in network class, is too high for me...

I think, it would be exciting if an implementation of "different g-functions with different mass flow rates using a non-stationary combination [3]" in pygfunction would be implemented in future releases!

[MassimoCimmino]

I developed g-functions for the simulation of bore fields with variable mass flow rate and reversible flow direction. This makes it possible to simulate borehole thermal energy storage systems using g-functions. The solution is presented in a paper now available online [1].

The are several tasks to complete for the implementation into pygfunction:

• [x] Implement negative mass flow rates in pipes classes and networks classes to model reversed flow (#282).
• [x] Implement the evaluation of variable mass flow rate g-functions in the gfunction module.
• [x] Extend the load aggregation methods to use variable mass flow rate g-functions (eq. 25 of the paper).

---

[1] Cimmino, M. g-Functions for fields of series- and parallel-connected boreholes with variable fluid mass flow rate and reversible flow direction. Renewable Energy. [In Press.] May 15, 2024. https://doi.org/10.1016/j.renene.2024.120661

[fuma01]

Dear Massimo, This is great work! Thanks for the implementation! I'm looking forward to studying the paper and using the new functions! Best regards

[MassimoCimmino]

The following scripts reproduce the results from the paper (after 566f3fe)

• Parallel boreholes

  -

  Code

  ``` import numpy as np import matplotlib.pyplot as plt from scipy.signal import square from time import perf_counter import pygfunction as gt def main(): # ------------------------------------------------------------------------- # Simulation parameters # ------------------------------------------------------------------------- # Borehole dimensions D = 4.0 # Borehole buried depth (m) H = 150.0 # Borehole length (m) r_b = 0.075 # Borehole radius (m) B = 5.0 # Borehole spacing (m) # Pipe dimensions r_out = 33.6e-3/2 # Pipe outer radius (m) r_in = 27.4e-3/2 # Pipe inner radius (m) D_s = 0.11/2 # Shank spacing (m) epsilon = 1.0e-6 # Pipe roughness (m) # Pipe positions # Single U-tube [(x_in, y_in), (x_out, y_out)] pos_pipes = [(-D_s, 0.), (D_s, 0.)] # Ground properties alpha = 2.5/2.2e6 # Ground thermal diffusivity (m2/s) k_s = 2.5 # Ground thermal conductivity (W/m.K) T_g = 10. # Undisturbed ground temperatue (degC) # Grout properties k_g = 1.5 # Grout thermal conductivity (W/m.K) # Pipe properties k_p = 0.42 # Pipe thermal conductivity (W/m.K) # Fluid properties # The fluid is propylene-glycol (20 %) at 20 degC fluid = gt.media.Fluid('MPG', 20.) 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) m_flow_borehole = np.array([0.25, 0.60]) # Total fluid mass flow rate per borehole (kg/s) # Calculation options for thermal response factors and g-functions nSegments = 8 segment_ratios = gt.utilities.segment_ratios(nSegments) options = {'approximate_FLS': True, 'nSegments': nSegments, 'segment_ratios': segment_ratios, 'disp': True, 'linear_threshold': 0.} # Calculation options for plotting g-functions dt_gFunc = 100*3600. # Time step tmax_gFunc = 3000. * 8760. * 3600. # Maximum time Nt_gFunc = 50 # Number of time steps # Geometrically expanding time vector time_gFunc = gt.utilities.time_geometric(dt_gFunc, tmax_gFunc, Nt_gFunc) # Calculation options for simulations dt = 600. # Time step tmax = 48 * 3600. # Maximum time cells_per_level = 5 # Cells per level for load aggregation time = np.arange(dt, tmax + dt, dt) time_h = time / 3600. Nt = len(time) # Number of time steps # ------------------------------------------------------------------------- # Borehole field #1 # ------------------------------------------------------------------------- # Field of 6x4 (n=24) boreholes N_1 = 6 N_2 = 4 boreField = gt.boreholes.rectangle_field(N_1, N_2, B, B, H, D, r_b) nBoreholes = len(boreField) H_tot = H * nBoreholes m_flow_network = m_flow_borehole * nBoreholes # ------------------------------------------------------------------------- # Initialize pipe model # ------------------------------------------------------------------------- # Always use the maximum fluid mass flow rate for thermal resistances m_flow_max = np.max(m_flow_borehole) # Pipe thermal resistance R_p = gt.pipes.conduction_thermal_resistance_circular_pipe( r_in, r_out, k_p) # Fluid to inner pipe wall thermal resistance (Single U-tube) m_flow_pipe = m_flow_max 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 * np.pi * r_in) # Single U-tube, same for all boreholes in the bore field UTubes = [] for borehole in boreField: SingleUTube = gt.pipes.SingleUTube(pos_pipes, r_in, r_out, borehole, k_s, k_g, R_f + R_p) UTubes.append(SingleUTube) network = gt.networks.Network(boreField, UTubes) # ------------------------------------------------------------------------- # Evaluate and plot variable mass flow rate g-functions # ------------------------------------------------------------------------- # Calculate the g-function for uniform heat extraction rate gfunc = gt.gfunction.gFunction( network, alpha, time=time_gFunc, m_flow_network=m_flow_network, cp_f=cp_f, options=options) # Plot g-functions plt.rc('figure', figsize=(90.0/25.4, 70.0/25.4)) fig = gfunc.visualize_g_function() plt.tight_layout() # ------------------------------------------------------------------------- # Evaluate the g-functions and thermal response factors for simulations # ------------------------------------------------------------------------- nSources = nSegments*nBoreholes # Number of borehole segments nModes = len(m_flow_borehole) # Number of operation states # Load aggregation object for simulations using variable mass flow rate g-functions # (The number of aggregated loads is equal to the number of operation states) loadAgg = gt.load_aggregation.ClaessonJaved( dt, tmax, nSources=nModes, cells_per_level=cells_per_level) # Times required for g-function and thermal resposne factors evaluation time_req = loadAgg.get_times_for_simulation() # Variable mass flow rate g-functions gfunc = gt.gfunction.gFunction( network, alpha, time=time_req, m_flow_network=m_flow_network, cp_f=cp_f, options=options) loadAgg.initialize(gfunc.gFunc / (2 * np.pi * k_s)) # ------------------------------------------------------------------------- # Flow rates and loads # ------------------------------------------------------------------------- # Linear variation from 0 to -80 W/m over 48 hours Q_tot = -80 * time_h / 48 * H * nBoreholes # Alternating between minimum and maximum mass flow rate every 4 hours m_flow_min = np.min(m_flow_borehole) m_flow_max = np.max(m_flow_borehole) m_flow_tot = m_flow_min + (m_flow_max - m_flow_min) * (0.5 + 0.5 * square(2*np.pi*(time_h - dt / 3600.) / 8, duty=0.5)) # ------------------------------------------------------------------------- # Simulation using variable mass flow rate g-functions # ------------------------------------------------------------------------- print(f'Simulation using variable mass flow rate g-functions: ', end='') # Initialize chrono tic = perf_counter() T_b = np.zeros(Nt) T_f_in = np.zeros(Nt) T_f_out = np.zeros(Nt) for i, t in enumerate(time): m_flow = m_flow_tot[i] * nBoreholes # Increment time step by (1) loadAgg.next_time_step(t) # Solve for temperatures and heat extraction rates (variable mass flow rate) xi = np.zeros(nModes) iMode = np.argmax(m_flow_tot[i] <= m_flow_borehole) if iMode > 0: xi[iMode-1] = (m_flow_borehole[iMode] - m_flow_tot[i]) / (m_flow_borehole[iMode] - m_flow_borehole[iMode-1]) xi[iMode] = (m_flow_tot[i] - m_flow_borehole[iMode-1]) / (m_flow_borehole[iMode] - m_flow_borehole[iMode-1]) else: xi[0] = 1. loadAgg.set_current_load(Q_tot[i]/H_tot * xi) deltaT_b = loadAgg.temporal_superposition() T_b[i] = T_g - deltaT_b @ xi T_f_in[i] = network.get_network_inlet_temperature( Q_tot[i], T_b[i], m_flow, cp_f, nSegments=1, segment_ratios=None) T_f_out[i] = network.get_network_outlet_temperature( T_f_in[i], T_b[i], m_flow, cp_f, nSegments=1, segment_ratios=None) toc = perf_counter() print(f' {toc - tic:.3f} sec') # ------------------------------------------------------------------------- # Plot simulation results # ------------------------------------------------------------------------- # Configure figure and axes plt.rc('figure', figsize=(140.0/25.4, 105.0/25.4)) fig = gt.utilities._initialize_figure() ax11 = fig.add_subplot(221) ax12 = fig.add_subplot(222) ax21 = fig.add_subplot(223) # Results - Flow rates # -------------------- # Axis labels ax11.set_xlabel(r'$t$ [h]') ax11.set_ylabel(r'$\dot{m}$ [kg/s]') gt.utilities._format_axes(ax11) # Plot temperatures t_plot = np.repeat(np.concatenate(([0.], time_h)), 2)[1:-1] m_flow_plot = np.repeat(m_flow_tot * nBoreholes, 2) # ax11.plot(time_h, m_flow_tot * nBoreholes, 'k-') ax11.plot(t_plot, m_flow_plot, 'k-') # Adjust to plot window plt.tight_layout() # Results - Heat injection rate # ----------------------------- # Axis labels ax12.set_xlabel(r'$t$ [h]') ax12.set_ylabel(r'$\bar{Q}^\prime$ [W/m]') gt.utilities._format_axes(ax12) # Plot temperatures t_plot = np.repeat(np.concatenate(([0.], time_h)), 2)[1:-1] Q_plot = np.repeat(Q_tot / (nBoreholes * H), 2) ax12.plot(t_plot, Q_plot, 'k-') # Adjust to plot window plt.tight_layout() # Results - Temperatures # ---------------------- # Axis labels ax21.set_xlabel(r'$t$ [h]') ax21.set_ylabel(r'$T_{f,out}$ [degC]') gt.utilities._format_axes(ax21) # Plot temperatures ax21.plot(time_h, T_f_out, 'k-') # Adjust to plot window plt.tight_layout() return # Main function if __name__ == '__main__': main() ```

  • Figures

  

• Series boreholes

  -

  Code

  ``` import numpy as np import matplotlib.pyplot as plt from scipy.signal import square from time import perf_counter import pygfunction as gt def main(): # ------------------------------------------------------------------------- # Simulation parameters # ------------------------------------------------------------------------- # Borehole dimensions D = 4.0 # Borehole buried depth (m) H = 150.0 # Borehole length (m) r_b = 0.075 # Borehole radius (m) B = 5.0 # Borehole spacing (m) # Pipe dimensions r_out = 33.6e-3/2 # Pipe outer radius (m) r_in = 27.4e-3/2 # Pipe inner radius (m) D_s = 0.11/2 # Shank spacing (m) epsilon = 1.0e-6 # Pipe roughness (m) # Pipe positions # Single U-tube [(x_in, y_in), (x_out, y_out)] pos_pipes = [(-D_s, 0.), (D_s, 0.)] # Ground properties alpha = 2.5/2.2e6 # Ground thermal diffusivity (m2/s) k_s = 2.5 # Ground thermal conductivity (W/m.K) T_g = 10. # Undisturbed ground temperatue (degC) # Grout properties k_g = 1.5 # Grout thermal conductivity (W/m.K) # Pipe properties k_p = 0.42 # Pipe thermal conductivity (W/m.K) # Fluid properties # The fluid is propylene-glycol (20 %) at 20 degC fluid = gt.media.Fluid('MPG', 20.) 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) m_flow_borehole = np.array([0.5, 1.2]) # Total fluid mass flow rate per borehole (kg/s) # Add reversed flows to list of mass flow rates m_flow_borehole = np.concatenate([-m_flow_borehole[::-1], m_flow_borehole]) # Calculation options for thermal response factors and g-functions nSegments = 8 segment_ratios = gt.utilities.segment_ratios(nSegments) options = {'approximate_FLS': True, 'nSegments': nSegments, 'segment_ratios': segment_ratios, 'disp': True, 'linear_threshold': 0.} # Calculation options for plotting g-functions dt_gFunc = 100*3600. # Time step tmax_gFunc = 3000. * 8760. * 3600. # Maximum time Nt_gFunc = 50 # Number of time steps # Geometrically expanding time vector time_gFunc = gt.utilities.time_geometric(dt_gFunc, tmax_gFunc, Nt_gFunc) # Calculation options for simulations dt = 600. # Time step tmax = 48 * 3600. # Maximum time cells_per_level = 5 # Cells per level for load aggregation time = np.arange(dt, tmax + dt, dt) time_h = time / 3600. Nt = len(time) # Number of time steps # ------------------------------------------------------------------------- # Borehole field #2 # ------------------------------------------------------------------------- # Field of 6x1 (n=6) boreholes N_1 = 6 N_2 = 1 boreField = gt.boreholes.rectangle_field(N_1, N_2, B, B, H, D, r_b) nBoreholes = len(boreField) H_tot = H * nBoreholes m_flow_network = m_flow_borehole # ------------------------------------------------------------------------- # Initialize pipe model # ------------------------------------------------------------------------- # Always use the maximum fluid mass flow rate for thermal resistances m_flow_max = np.max(np.abs(m_flow_borehole)) # Pipe thermal resistance R_p = gt.pipes.conduction_thermal_resistance_circular_pipe( r_in, r_out, k_p) # Fluid to inner pipe wall thermal resistance (Single U-tube) m_flow_pipe = m_flow_max 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 * np.pi * r_in) # Single U-tube, same for all boreholes in the bore field UTubes = [] for borehole in boreField: SingleUTube = gt.pipes.SingleUTube(pos_pipes, r_in, r_out, borehole, k_s, k_g, R_f + R_p) UTubes.append(SingleUTube) bore_connectivity=[i-1 for i in range(nBoreholes)] network = gt.networks.Network(boreField, UTubes, bore_connectivity=bore_connectivity) # ------------------------------------------------------------------------- # Evaluate and plot variable mass flow rate g-functions # ------------------------------------------------------------------------- # Calculate the g-function for uniform heat extraction rate gfunc = gt.gfunction.gFunction( network, alpha, time=time_gFunc, m_flow_network=m_flow_network, cp_f=cp_f, options=options) # Plot g-functions gfunc.visualize_g_function(which=[(i,i) for i in range(len(m_flow_network))]) gfunc.visualize_g_function(which=[(0,i) for i in range(len(m_flow_network))]) # ------------------------------------------------------------------------- # Evaluate the g-functions and thermal response factors for simulations # ------------------------------------------------------------------------- nSources = nSegments*nBoreholes # Number of borehole segments nModes = len(m_flow_borehole) # Number of operation states # Load aggregation object for simulations using variable mass flow rate g-functions # (The number of aggregated loads is equal to the number of operation states) loadAgg = gt.load_aggregation.ClaessonJaved( dt, tmax, nSources=nModes, cells_per_level=cells_per_level) # Times required for g-function and thermal resposne factors evaluation time_req = loadAgg.get_times_for_simulation() # Variable mass flow rate g-functions gfunc = gt.gfunction.gFunction( network, alpha, time=time_req, m_flow_network=m_flow_network, cp_f=cp_f, options=options).gFunc loadAgg.initialize(gfunc/ (2 * np.pi * k_s)) # ------------------------------------------------------------------------- # Flow rates and loads # ------------------------------------------------------------------------- # Sinusoidal variation between -80 and 80 W/m with minimum of 10 W/m Q_min = -80 * H * nBoreholes Q_max = 80 * H * nBoreholes Q_small = 10 * H * nBoreholes Q_avg = 0.5 * (Q_min + Q_max) Q_amp = 0.5 * (Q_max - Q_min) Q_tot = -Q_avg - Q_amp * np.sin(2 * np.pi * time_h / 24) # Mass flow rate is proportional to the absolute load m_flow_min = np.min(np.abs(m_flow_borehole)) m_flow_max = np.max(np.abs(m_flow_borehole)) m_flow_tot = np.sign(Q_tot) * np.maximum(m_flow_min, np.minimum(m_flow_max, m_flow_max * np.abs(Q_tot) / Q_max)) m_flow_tot[np.abs(Q_tot) < Q_small] = 0. Q_tot[np.abs(Q_tot) < Q_small] = 0. # ------------------------------------------------------------------------- # Simulation using variable mass flow rate g-functions # ------------------------------------------------------------------------- print(f'Simulation using variable mass flow rate g-functions: ', end='') # Initialize chrono tic = perf_counter() T_b = np.zeros(Nt) T_f_in = np.zeros(Nt) T_f_out = np.zeros(Nt) for i, t in enumerate(time): m_flow = m_flow_tot[i] # Increment time step by (1) loadAgg.next_time_step(t) # Solve for temperatures and heat extraction rates (variable mass flow rate) if np.abs(Q_tot[i]) >= Q_small: xi = np.zeros(nModes) iMode = np.argmax(m_flow_tot[i] <= m_flow_network) if iMode > 0: xi[iMode-1] = (m_flow_network[iMode] - m_flow_tot[i]) / (m_flow_network[iMode] - m_flow_network[iMode-1]) xi[iMode] = (m_flow_tot[i] - m_flow_network[iMode-1]) / (m_flow_network[iMode] - m_flow_network[iMode-1]) else: xi[0] = 1. loadAgg.set_current_load(Q_tot[i]/H_tot * xi) deltaT_b = loadAgg.temporal_superposition() T_b[i] = T_g - deltaT_b @ xi T_f_in[i] = network.get_network_inlet_temperature( Q_tot[i], T_b[i], m_flow, cp_f, nSegments=1, segment_ratios=None) T_f_out[i] = network.get_network_outlet_temperature( T_f_in[i], T_b[i], m_flow, cp_f, nSegments=1, segment_ratios=None) else: T_b[i] = np.nan T_f_in[i] = np.nan T_f_out[i] = np.nan toc = perf_counter() print(f' {toc - tic:.3f} sec') # ------------------------------------------------------------------------- # Plot simulation results # ------------------------------------------------------------------------- # Configure figure and axes plt.rc('figure', figsize=(140.0/25.4, 105.0/25.4)) fig = gt.utilities._initialize_figure() ax11 = fig.add_subplot(221) ax12 = fig.add_subplot(222) ax21 = fig.add_subplot(223) # Results - Flow rates # -------------------- # Axis labels ax11.set_xlabel(r'$t$ [h]') ax11.set_ylabel(r'$\dot{m}$ [kg/s]') gt.utilities._format_axes(ax11) # Plot temperatures t_plot = np.repeat(np.concatenate(([0.], time_h)), 2)[1:-1] m_flow_plot = np.repeat(m_flow_tot, 2) ax11.plot(t_plot, m_flow_plot, 'k-') # Adjust to plot window plt.tight_layout() # Results - Heat injection rate # ----------------------------- # Axis labels ax12.set_xlabel(r'$t$ [h]') ax12.set_ylabel(r'$\bar{Q}^\prime$ [W/m]') gt.utilities._format_axes(ax12) # Plot temperatures t_plot = np.repeat(np.concatenate(([0.], time_h)), 2)[1:-1] Q_plot = np.repeat(Q_tot / (nBoreholes * H), 2) ax12.plot(t_plot, Q_plot, 'k-') # Adjust to plot window plt.tight_layout() # Results - Temperatures # ---------------------- # Axis labels ax21.set_xlabel(r'$t$ [h]') ax21.set_ylabel(r'$T_{f,out}$ [degC]') gt.utilities._format_axes(ax21) # Plot temperatures ax21.plot(time_h, T_f_out, 'k-') # Adjust to plot window plt.tight_layout() return # Main function if __name__ == '__main__': main() ```

  • Figures

  

• Borehole thermal energy storage

  -

  Code

  ``` import numpy as np import matplotlib.pyplot as plt from scipy.signal import square from scipy.interpolate import interp1d from time import perf_counter import pygfunction as gt import pvlib as pv def main(): # ------------------------------------------------------------------------- # Simulation parameters # ------------------------------------------------------------------------- # Building properties T_eq = 18.0 # Base temperature of the buildings [degC] # eps weather file for Montreal based on CWEC2020 weather_file = './src/CAN_QC_MONTREAL-INTL-A_7025251_CWEC.epw' data, metadata = pv.iotools.read_epw(weather_file) time_data = np.arange(3600., 3600. * 8761, 3600.) T_amb = data['temp_air'].to_numpy() # Solar irradiance on a 45 degrees tilted surface facing south solar_position = pv.solarposition.get_solarposition( data.index, metadata['latitude'], metadata['longitude'], metadata['altitude']) G_t = pv.irradiance.get_total_irradiance( 45., 180., solar_position['zenith'], solar_position['azimuth'], data['dni'], data['ghi'], data['dhi'])['poa_global'].to_numpy() # return data, solar_position # Borehole dimensions D = 4.0 # Borehole buried depth (m) H = 150.0 # Borehole length (m) r_b = 0.075 # Borehole radius (m) B = 5.0 # Borehole spacing (m) # Pipe dimensions r_out = 33.6e-3/2 # Pipe outer radius (m) r_in = 27.4e-3/2 # Pipe inner radius (m) D_s = 0.11/2 # Shank spacing (m) epsilon = 1.0e-6 # Pipe roughness (m) # Pipe positions # Single U-tube [(x_in, y_in), (x_out, y_out)] pos_pipes = [(-D_s, 0.), (D_s, 0.)] # Ground properties alpha = 2.5/2.2e6 # Ground thermal diffusivity (m2/s) k_s = 2.5 # Ground thermal conductivity (W/m.K) T_g = 10. # Undisturbed ground temperatue (degC) # Grout properties k_g = 1.5 # Grout thermal conductivity (W/m.K) # Pipe properties k_p = 0.42 # Pipe thermal conductivity (W/m.K) # Fluid properties # The fluid is propylene-glycol (20 %) at 20 degC fluid = gt.media.Fluid('MPG', 20.) 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) m_flow_borehole = np.array([0.5, 1.2]) # Total fluid mass flow rate per borehole (kg/s) # Calculation options for thermal response factors and g-functions nSegments = 8 segment_ratios = gt.utilities.segment_ratios(nSegments) options = {'approximate_FLS': True, 'nSegments': nSegments, 'segment_ratios': segment_ratios, 'disp': True, 'linear_threshold': 0.} # Calculation options for plotting g-functions dt_gFunc = 100*3600. # Time step tmax_gFunc = 3000. * 8760. * 3600. # Maximum time Nt_gFunc = 50 # Number of time steps # Geometrically expanding time vector time_gFunc = gt.utilities.time_geometric(dt_gFunc, tmax_gFunc, Nt_gFunc) # Calculation options for simulations dt = 3600. # Time step tmax = 20 * 8760 * 3600. # Maximum time # tmax = 8760 * 3600. # Maximum time # tmax = 48 * 3600. # Maximum time cells_per_level = 5 # Cells per level for load aggregation time = np.arange(dt, tmax + dt, dt) time_h = time / 3600. Nt = len(time) # Number of time steps # ------------------------------------------------------------------------- # Borehole field #3 # ------------------------------------------------------------------------- # Field of 144 borehole in 24 parallel branches of 6 boreholes boreField = gt.boreholes.field_from_file('./src/BTES.txt') nBoreholes = len(boreField) nBranches = 24 nBorPerBranch = 6 H_tot = H * nBoreholes # ------------------------------------------------------------------------- # Initialize pipe model # ------------------------------------------------------------------------- bore_connectivity = [i+1 for i in range(nBoreholes)] for i in range(nBorPerBranch-1, nBoreholes, nBorPerBranch): bore_connectivity[i] = -1 # Add reversed flows to list of mass flow rates m_flow_borehole = np.concatenate([-m_flow_borehole[::-1], m_flow_borehole]) # Always use the maximum fluid mass flow rate for thermal resistances m_flow_max = np.max(m_flow_borehole) # Pipe thermal resistance R_p = gt.pipes.conduction_thermal_resistance_circular_pipe( r_in, r_out, k_p) # Fluid to inner pipe wall thermal resistance (Single U-tube) m_flow_pipe = m_flow_max 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 * np.pi * r_in) # Single U-tube, same for all boreholes in the bore field UTubes = [] for borehole in boreField: SingleUTube = gt.pipes.SingleUTube(pos_pipes, r_in, r_out, borehole, k_s, k_g, R_f + R_p) UTubes.append(SingleUTube) network = gt.networks.Network( boreField, UTubes, nSegments=nSegments, bore_connectivity=bore_connectivity) # ------------------------------------------------------------------------- # Evaluate and plot variable mass flow rate g-functions # ------------------------------------------------------------------------- # # Calculate the g-function for uniform heat extraction rate m_flow_network = m_flow_borehole * nBranches gfunc = gt.gfunction.gFunction( network, alpha, m_flow_network=m_flow_network, cp_f=cp_f, time=time_gFunc, options=options) # Plot g-functions gfunc.visualize_g_function(which=[(i,i) for i in range(len(m_flow_network))]) gfunc.visualize_g_function(which=[(0,i) for i in range(len(m_flow_network))]) plt.tight_layout() # ------------------------------------------------------------------------- # Evaluate the g-functions and thermal response factors for simulations # ------------------------------------------------------------------------- nSources = nSegments*nBoreholes # Number of borehole segments nModes = len(m_flow_borehole) # Number of operation states # Load aggregation object for simulations using variable mass flow rate g-functions # (The number of aggregated loads is equal to the number of operation states) loadAgg = gt.load_aggregation.ClaessonJaved( dt, tmax, nSources=nModes, cells_per_level=cells_per_level) # Variable mass flow rate g-functions time_req = loadAgg.get_times_for_simulation() gfunc = gt.gfunction.gFunction( network, alpha, time=time_req, m_flow_network=m_flow_network, cp_f=cp_f, options=options).gFunc loadAgg.initialize(gfunc / (2 * np.pi * k_s)) # ------------------------------------------------------------------------- # Flow rates and loads # ------------------------------------------------------------------------- # Sinusoidal variation between -80 and 80 W/m with minimum of 10 W/m Q_max = 80 * H * nBoreholes Q_small = 10 * H * nBoreholes Q_buildings = np.maximum( T_eq - np.interp(time, time_data, T_amb, period=8760*3600.), 0.) a = Q_max / np.max(Q_buildings) Q_buildings = Q_buildings * Q_max / np.max(Q_buildings) Q_solar = np.maximum( np.interp(np.ceil(time/3600.)*3600., time_data, G_t, period=8760*3600.), 0.) b = np.sum(Q_buildings) / np.sum(Q_solar) Q_solar = Q_solar * np.sum(Q_buildings) / np.sum(Q_solar) Q_tot = Q_buildings - Q_solar # Mass flow rate is proportional to the absolute load m_flow_min = np.min(np.abs(m_flow_borehole)) m_flow_max = np.max(np.abs(m_flow_borehole)) m_flow_tot = np.sign(Q_tot) * np.maximum(m_flow_min, np.minimum(m_flow_max, m_flow_max * np.abs(Q_tot) / Q_max)) m_flow_tot[np.abs(Q_tot) < Q_small] = 0. Q_tot[np.abs(Q_tot) < Q_small] = 0. # ------------------------------------------------------------------------- # Simulation using variable mass flow rate g-functions # ------------------------------------------------------------------------- print(f'Simulation using variable mass flow rate g-functions: ', end='') # Initialize chrono tic = perf_counter() T_b = np.zeros(Nt) T_f_in = np.zeros(Nt) T_f_out = np.zeros(Nt) iModes = np.zeros(Nt, dtype=int) xis = np.zeros((nModes, Nt)) # Initialize chrono tic = perf_counter() m_flow_all = np.concatenate((np.linspace(m_flow_network[0], m_flow_network[1]), np.linspace(m_flow_network[2], m_flow_network[3]))) R_field_all = [gt.networks.network_thermal_resistance(network, m_flow, cp_f) for m_flow in m_flow_all] f_R_field = interp1d(m_flow_all, R_field_all) for i, t in enumerate(time): m_flow = m_flow_tot[i] * nBranches # Increment time step by (1) loadAgg.next_time_step(t) # Solve for temperatures and heat extraction rates (variable mass flow rate) if np.abs(Q_tot[i]) >= Q_small: xi = np.zeros(nModes) iMode = np.argmax(m_flow_tot[i] <= m_flow_borehole) iModes[i] = iMode if iMode > 0: xi[iMode-1] = (m_flow_borehole[iMode] - m_flow_tot[i]) / (m_flow_borehole[iMode] - m_flow_borehole[iMode-1]) xi[iMode] = (m_flow_tot[i] - m_flow_borehole[iMode-1]) / (m_flow_borehole[iMode] - m_flow_borehole[iMode-1]) else: xi[0] = 1. xis[:,i] = xi loadAgg.set_current_load(Q_tot[i] / H_tot * xi) deltaT_b = loadAgg.temporal_superposition() T_b[i] = T_g - deltaT_b @ xi # T_f_in[i] = network.get_network_inlet_temperature( # Q_tot[i], T_b[i], m_flow, cp_f, nSegments=1, segment_ratios=None) # T_f_out[i] = network.get_network_outlet_temperature( # T_f_in[i], T_b[i], m_flow, cp_f, nSegments=1, segment_ratios=None) T_f = T_b[i] - Q_tot[i] / H_tot * f_R_field(m_flow) T_f_in[i] = T_f - 0.5 * Q_tot[i] / (np.abs(m_flow) * cp_f) T_f_out[i] = T_f + 0.5 * Q_tot[i] / (np.abs(m_flow) * cp_f) else: T_b[i] = np.nan T_f_in[i] = np.nan T_f_out[i] = np.nan toc = perf_counter() print(f' {toc - tic:.3f} sec') # ------------------------------------------------------------------------- # Plot simulation results # ------------------------------------------------------------------------- # Configure figure and axes plt.rc('figure', figsize=(140.0/25.4, 105.0/25.4)) fig = gt.utilities._initialize_figure() ax11 = fig.add_subplot(221) ax12 = fig.add_subplot(222) ax21 = fig.add_subplot(223) # Results - Flow rates # -------------------- # Axis labels ax11.set_xlabel(r'$t$ [h]') ax11.set_ylabel(r'$\dot{m}$ [kg/s]') gt.utilities._format_axes(ax11) # Plot mass flow rates t_plot = np.repeat(np.concatenate(([0.], time_h)), 2)[1:-1] m_flow_plot = np.repeat(m_flow_tot, 2) ax11.plot(t_plot[-2*8760:], m_flow_plot[-2*8760:] * nBranches, 'k-') # Adjust to plot window plt.tight_layout() # Results - Heat injection rate # ----------------------------- # Axis labels ax12.set_xlabel(r'$t$ [h]') ax12.set_ylabel(r'$\bar{Q}^\prime$ [W/m]') gt.utilities._format_axes(ax12) # Plot temperatures t_plot = np.repeat(np.concatenate(([0.], time_h)), 2)[1:-1] Q_plot = np.repeat(Q_tot / (nBoreholes * H), 2) ax12.plot(t_plot[-2*8760:], Q_plot[-2*8760:], 'k-') # Adjust to plot window plt.tight_layout() # Results - Temperatures # ---------------------- # Axis labels ax21.set_xlabel(r'$t$ [h]') ax21.set_ylabel(r'$T_{f,out}$ [degC]') gt.utilities._format_axes(ax21) # Plot temperatures ax21.plot(time_h[-8760:], T_f_out[-8760:]) ax21.legend() # Adjust to plot window plt.tight_layout() # ------------------------------------------------------------------------- # Plot simulation results # ------------------------------------------------------------------------- # Configure figure and axes plt.rc('figure', figsize=(140.0/25.4, 105.0/25.4)) fig = gt.utilities._initialize_figure() ax11 = fig.add_subplot(221) ax12 = fig.add_subplot(222) ax21 = fig.add_subplot(223) # Results - Flow rates # -------------------- # Axis labels ax11.set_xlabel(r'$t$ [h]') ax11.set_ylabel(r'$\dot{m}$ [kg/s]') gt.utilities._format_axes(ax11) # Plot mass flow rates t_plot = np.repeat(np.concatenate(([0.], time_h)), 2)[1:-1] m_flow_plot = np.repeat(m_flow_tot, 2) ax11.plot(t_plot[7200*24*2:7207*24*2], m_flow_plot[7200*24*2:7207*24*2] * nBranches, 'k-') # Adjust to plot window plt.tight_layout() # Results - Heat injection rate # ----------------------------- # Axis labels ax12.set_xlabel(r'$t$ [h]') ax12.set_ylabel(r'$\bar{Q}^\prime$ [W/m]') gt.utilities._format_axes(ax12) # Plot temperatures t_plot = np.repeat(np.concatenate(([0.], time_h)), 2)[1:-1] Q_plot = np.repeat(Q_tot / (nBoreholes * H), 2) ax12.plot(t_plot[7200*24*2:7207*24*2], Q_plot[7200*24*2:7207*24*2], 'k-') # Adjust to plot window plt.tight_layout() # Results - Temperatures # ---------------------- # Axis labels ax21.set_xlabel(r'$t$ [h]') ax21.set_ylabel(r'$T_{f,out}$ [degC]') gt.utilities._format_axes(ax21) # Plot temperatures ax21.plot(time_h[7200*24:7207*24], T_f_out[7200*24:7207*24]) ax21.legend() # Adjust to plot window plt.tight_layout() return # Main function if __name__ == '__main__': out = main() ```

  • Figures

  

  BTES.txt