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