Compare the accuracy and speed of different load aggregation algorithms¶
This example compares the simulation times and the accuracy of borehole wall temperature predictions of different load aggregation algorithms implemented into the load agggregation module.
The g-function of a single borehole is first calculated. Then, the borehole wall temperature variations are calculated using the load aggregation schemes of Bernier et al. 1, Liu 2, and Claesson and Javed 3,. The time-variation of heat extraction rates is given by the synthetic load profile of Bernier et al. 1.
The following script validates the load aggregation schemes with the exact solution obtained from convolution in the Fourier domain (see ref. 4).
The script is located in: pygfunction/examples/comparison_load_aggregation.py
1# -*- coding: utf-8 -*-
2""" Comparison of the accuracy and computational speed of different load
3 aggregation algorithms.
4
5 The g-function of a single borehole is calculated for boundary condition of
6 uniform borehole wall temperature along the borehole. Then, the borehole
7 wall temperature variations resulting from a time-varying load profile
8 are simulated using the aggregation methods of Bernier et al. (2004),
9 Liu (2005), and Claesson and Javed (2012). Results are compared to the
10 exact solution obtained by convolution in the Fourier domain.
11
12 Default parameters are used for each of the aggregation schemes.
13
14"""
15from time import perf_counter
16
17import matplotlib.pyplot as plt
18import numpy as np
19from scipy.constants import pi
20from scipy.interpolate import interp1d
21from scipy.signal import fftconvolve
22
23import pygfunction as gt
24
25
26def main():
27 # -------------------------------------------------------------------------
28 # Simulation parameters
29 # -------------------------------------------------------------------------
30
31 # Borehole dimensions
32 D = 4.0 # Borehole buried depth (m)
33 H = 150.0 # Borehole length (m)
34 r_b = 0.075 # Borehole radius (m)
35
36 # Ground properties
37 alpha = 1.0e-6 # Ground thermal diffusivity (m2/s)
38 k_s = 2.0 # Ground thermal conductivity (W/m.K)
39 T_g = 10.0 # Undisturbed ground temperature (degC)
40
41 # g-Function calculation options
42 options = {'nSegments': 8,
43 'disp': True}
44
45 # Simulation parameters
46 dt = 3600. # Time step (s)
47 tmax = 20.*8760. * 3600. # Maximum time (s)
48 Nt = int(np.ceil(tmax/dt)) # Number of time steps
49 time = dt * np.arange(1, Nt+1)
50
51 # Evaluate heat extraction rate
52 Q_b = synthetic_load(time/3600.)
53
54 # Load aggregation schemes
55 ClaessonJaved = gt.load_aggregation.ClaessonJaved(dt, tmax)
56 MLAA = gt.load_aggregation.MLAA(dt, tmax)
57 Liu = gt.load_aggregation.Liu(dt, tmax)
58 LoadAggSchemes = [ClaessonJaved, MLAA, Liu]
59 loadAgg_labels = ['Claesson and Javed', 'MLAA', 'Liu']
60 loadAgg_lines = ['b-', 'k--', 'r-.']
61
62 # -------------------------------------------------------------------------
63 # Calculate g-function
64 # -------------------------------------------------------------------------
65
66 # The field contains only one borehole
67 boreField = [gt.boreholes.Borehole(H, D, r_b, x=0., y=0.)]
68 # Evaluate the g-function on a geometrically expanding time grid
69 time_gFunc = gt.utilities.time_geometric(dt, tmax, 50)
70 # Calculate g-function
71 gFunc = gt.gfunction.gFunction(
72 boreField, alpha, time=time_gFunc, options=options)
73
74 # -------------------------------------------------------------------------
75 # Simulation
76 # -------------------------------------------------------------------------
77 nLoadAgg = len(LoadAggSchemes)
78 T_b = np.zeros((nLoadAgg, Nt))
79
80 t_calc = np.zeros(nLoadAgg)
81 for n, (LoadAgg, label) in enumerate(zip(LoadAggSchemes, loadAgg_labels)):
82 print(f'Simulation using {label} ...')
83 # Interpolate g-function at required times
84 time_req = LoadAgg.get_times_for_simulation()
85 gFunc_int = interp1d(np.hstack([0., time_gFunc]),
86 np.hstack([0., gFunc.gFunc]),
87 kind='cubic',
88 bounds_error=False,
89 fill_value=(0., gFunc.gFunc[-1]))(time_req)
90 # Initialize load aggregation scheme
91 LoadAgg.initialize(gFunc_int/(2*pi*k_s))
92
93 tic = perf_counter()
94 for i in range(Nt):
95 # Increment time step by (1)
96 LoadAgg.next_time_step(time[i])
97
98 # Apply current load
99 LoadAgg.set_current_load(Q_b[i]/H)
100
101 # Evaluate borehole wall temeprature
102 deltaT_b = LoadAgg.temporal_superposition()
103 T_b[n,i] = T_g - deltaT_b
104 toc = perf_counter()
105 t_calc[n] = toc - tic
106
107 # -------------------------------------------------------------------------
108 # Calculate exact solution from convolution in the Fourier domain
109 # -------------------------------------------------------------------------
110
111 # Heat extraction rate increment
112 dQ = np.zeros(Nt)
113 dQ[0] = Q_b[0]
114 dQ[1:] = Q_b[1:] - Q_b[:-1]
115 # Interpolated g-function
116 g = interp1d(time_gFunc, gFunc.gFunc)(time)
117
118 # Convolution in Fourier domain
119 T_b_exact = T_g - fftconvolve(dQ, g/(2.0*pi*k_s*H), mode='full')[0:Nt]
120
121 # -------------------------------------------------------------------------
122 # plot results
123 # -------------------------------------------------------------------------
124
125 # Configure figure and axes
126 fig = gt.utilities._initialize_figure()
127
128 ax1 = fig.add_subplot(311)
129 # Axis labels
130 ax1.set_xlabel(r'$t$ [hours]')
131 ax1.set_ylabel(r'$Q_b$ [W]')
132 gt.utilities._format_axes(ax1)
133 hours = np.array([(j+1)*dt/3600. for j in range(Nt)])
134 ax1.plot(hours, Q_b)
135
136 ax2 = fig.add_subplot(312)
137 # Axis labels
138 ax2.set_xlabel(r'$t$ [hours]')
139 ax2.set_ylabel(r'$T_b$ [degC]')
140 gt.utilities._format_axes(ax2)
141 for T_b_n, line, label in zip(T_b, loadAgg_lines, loadAgg_labels):
142 ax2.plot(hours, T_b_n, line, label=label)
143 ax2.plot(hours, T_b_exact, 'k.', label='exact')
144 ax2.legend()
145
146 ax3 = fig.add_subplot(313)
147 # Axis labels
148 ax3.set_xlabel(r'$t$ [hours]')
149 ax3.set_ylabel(r'Error [degC]')
150 gt.utilities._format_axes(ax3)
151 for T_b_n, line, label in zip(T_b, loadAgg_lines, loadAgg_labels):
152 ax3.plot(hours, T_b_n - T_b_exact, line, label=label)
153 # Adjust to plot window
154 plt.tight_layout()
155
156 # -------------------------------------------------------------------------
157 # Print performance metrics
158 # -------------------------------------------------------------------------
159
160 # Maximum errors in evaluation of borehole wall temperatures
161 maxError = np.array([np.max(np.abs(T_b_n-T_b_exact)) for T_b_n in T_b])
162 # Print results
163 print('Simulation results')
164 for label, maxError_n, t_calc_n in zip(loadAgg_labels, maxError, t_calc):
165 print()
166 print((f' {label} ').center(60, '-'))
167 print(f'Maximum absolute error : {maxError_n:.3f} degC')
168 print(f'Calculation time : {t_calc_n:.3f} sec')
169
170 return
171
172
173def synthetic_load(x):
174 """
175 Synthetic load profile of Bernier et al. (2004).
176
177 Returns load y (in watts) at time x (in hours).
178 """
179 A = 2000.0
180 B = 2190.0
181 C = 80.0
182 D = 2.0
183 E = 0.01
184 F = 0.0
185 G = 0.95
186
187 func = (168.0-C)/168.0
188 for i in [1, 2, 3]:
189 func += 1.0/(i*pi)*(np.cos(C*pi*i/84.0) - 1.0) \
190 *(np.sin(pi*i/84.0*(x - B)))
191 func = func*A*np.sin(pi/12.0*(x - B)) \
192 *np.sin(pi/4380.0*(x - B))
193
194 y = func + (-1.0)**np.floor(D/8760.0*(x - B))*abs(func) \
195 + E*(-1.0)**np.floor(D/8760.0*(x - B)) \
196 /np.sign(np.cos(D*pi/4380.0*(x - F)) + G)
197 return -y
198
199
200# Main function
201if __name__ == '__main__':
202 main()
References
- 1(1,2)
Bernier, M., Pinel, P., Labib, R. and Paillot, R. (2004). A multiple load aggregation algorithm for annual hourly simulations of GCHP systems. HVAC&R Research 10 (4): 471–487.
- 2
Liu, X. (2005). Development and experimental validation of simulation of hydronic snow melting systems for bridges. Ph.D. Thesis. Oklahoma State University.
- 3
Claesson, J., & Javed, S. (2011). A load-aggregation method to calculate extraction temperatures of borehole heat exchangers. ASHRAE Transactions, 118 (1): 530–539.
- 4
Marcotte, D., & Pasquier, P. (2008). Fast fluid and ground temperature computation for geothermal ground-loop heat exchanger systems. Geothermics, 37 (6) : 651-665.