Simulation of a borehole using load aggregation

This example demonstrates the use of the load agggregation module to predict the borehole wall temperature of a single temperature with known heat extraction rates.

The g-function of a single borehole is first calculated. Then, the borehole wall temperature variations are calculated using the load aggregation scheme of Claesson and Javed 1. The time-variation of heat extraction rates is given by the synthetic load profile of Bernier et al. 2.

The following script validates the load aggregation scheme with the exact solution obtained from convolution in the Fourier domain (see ref. 3).

The script is located in: pygfunction/examples/load_aggregation.py

  1# -*- coding: utf-8 -*-
  2""" Example of simulation of a geothermal system.
  3
  4    The g-function of a single borehole is calculated for boundary condition of
  5    uniform borehole wall temperature along the borehole. Then, the borehole
  6    wall temperature variations resulting from a time-varying load profile
  7    are simulated using the aggregation method of Claesson and Javed (2012).
  8
  9"""
 10import matplotlib.pyplot as plt
 11import numpy as np
 12from scipy.constants import pi
 13from scipy.interpolate import interp1d
 14from scipy.signal import fftconvolve
 15
 16import pygfunction as gt
 17
 18
 19def main():
 20    # -------------------------------------------------------------------------
 21    # Simulation parameters
 22    # -------------------------------------------------------------------------
 23
 24    # Borehole dimensions
 25    D = 4.0             # Borehole buried depth (m)
 26    H = 150.0           # Borehole length (m)
 27    r_b = 0.075         # Borehole radius (m)
 28
 29    # Ground properties
 30    alpha = 1.0e-6      # Ground thermal diffusivity (m2/s)
 31    k_s = 2.0           # Ground thermal conductivity (W/m.K)
 32    T_g = 10.0          # Undisturbed ground temperature (degC)
 33
 34    # g-Function calculation options
 35    options = {'nSegments': 8,
 36               'disp': True}
 37
 38    # Simulation parameters
 39    dt = 3600.                  # Time step (s)
 40    tmax = 20.*8760. * 3600.    # Maximum time (s)
 41    Nt = int(np.ceil(tmax/dt))  # Number of time steps
 42    time = dt * np.arange(1, Nt+1)
 43
 44    # Evaluate heat extraction rate
 45    Q_b = synthetic_load(time/3600.)
 46
 47    # Load aggregation scheme
 48    LoadAgg = gt.load_aggregation.ClaessonJaved(dt, tmax)
 49
 50    # -------------------------------------------------------------------------
 51    # Calculate g-function
 52    # -------------------------------------------------------------------------
 53
 54    # The field contains only one borehole
 55    boreField = [gt.boreholes.Borehole(H, D, r_b, x=0., y=0.)]
 56    # Get time values needed for g-function evaluation
 57    time_req = LoadAgg.get_times_for_simulation()
 58    # Calculate g-function
 59    gFunc = gt.gfunction.gFunction(
 60        boreField, alpha, time=time_req, options=options)
 61    # Initialize load aggregation scheme
 62    LoadAgg.initialize(gFunc.gFunc/(2*pi*k_s))
 63
 64    # -------------------------------------------------------------------------
 65    # Simulation
 66    # -------------------------------------------------------------------------
 67
 68    T_b = np.zeros(Nt)
 69    for i, (t, Q_b_i) in enumerate(zip(time, Q_b)):
 70        # Increment time step by (1)
 71        LoadAgg.next_time_step(t)
 72
 73        # Apply current load
 74        LoadAgg.set_current_load(Q_b_i/H)
 75
 76        # Evaluate borehole wall temeprature
 77        deltaT_b = LoadAgg.temporal_superposition()
 78        T_b[i] = T_g - deltaT_b
 79
 80    # -------------------------------------------------------------------------
 81    # Calculate exact solution from convolution in the Fourier domain
 82    # -------------------------------------------------------------------------
 83
 84    # Heat extraction rate increment
 85    dQ = np.zeros(Nt)
 86    dQ[0] = Q_b[0]
 87    dQ[1:] = Q_b[1:] - Q_b[:-1]
 88    # Interpolated g-function
 89    g = interp1d(time_req, gFunc.gFunc)(time)
 90
 91    # Convolution in Fourier domain
 92    T_b_exact = T_g - fftconvolve(dQ, g/(2.0*pi*k_s*H), mode='full')[0:Nt]
 93
 94    # -------------------------------------------------------------------------
 95    # plot results
 96    # -------------------------------------------------------------------------
 97
 98    # Configure figure and axes
 99    fig = gt.utilities._initialize_figure()
100
101    ax1 = fig.add_subplot(311)
102    # Axis labels
103    ax1.set_xlabel(r'$t$ [hours]')
104    ax1.set_ylabel(r'$Q_b$ [W]')
105    gt.utilities._format_axes(ax1)
106
107    hours = np.arange(1, Nt+1) * dt / 3600.
108    ax1.plot(hours, Q_b)
109
110    ax2 = fig.add_subplot(312)
111    # Axis labels
112    ax2.set_xlabel(r'$t$ [hours]')
113    ax2.set_ylabel(r'$T_b$ [degC]')
114    gt.utilities._format_axes(ax2)
115
116    ax2.plot(hours, T_b)
117    ax2.plot(hours, T_b_exact, 'k.')
118
119    ax3 = fig.add_subplot(313)
120    # Axis labels
121    ax3.set_xlabel(r'$t$ [hours]')
122    ax3.set_ylabel(r'Error [degC]')
123    gt.utilities._format_axes(ax3)
124
125    ax3.plot(hours, T_b - T_b_exact)
126
127    # Adjust to plot window
128    plt.tight_layout()
129
130    return
131
132
133def synthetic_load(x):
134    """
135    Synthetic load profile of Bernier et al. (2004).
136
137    Returns load y (in watts) at time x (in hours).
138    """
139    A = 2000.0
140    B = 2190.0
141    C = 80.0
142    D = 2.0
143    E = 0.01
144    F = 0.0
145    G = 0.95
146
147    func = (168.0-C)/168.0
148    for i in [1,2,3]:
149        func += 1.0/(i*pi)*(np.cos(C*pi*i/84.0)-1.0) \
150                          *(np.sin(pi*i/84.0*(x-B)))
151    func = func*A*np.sin(pi/12.0*(x-B)) \
152           *np.sin(pi/4380.0*(x-B))
153
154    y = func + (-1.0)**np.floor(D/8760.0*(x-B))*abs(func) \
155      + E*(-1.0)**np.floor(D/8760.0*(x-B))/np.sign(np.cos(D*pi/4380.0*(x-F))+G)
156    return -y
157
158
159# Main function
160if __name__ == '__main__':
161    main()

References

1

Claesson, J., & Javed, S. (2011). A load-aggregation method to calculate extraction temperatures of borehole heat exchangers. ASHRAE Transactions, 118 (1): 530–539.

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.

3

Marcotte, D., & Pasquier, P. (2008). Fast fluid and ground temperature computation for geothermal ground-loop heat exchanger systems. Geothermics, 37 (6) : 651-665.