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.