Simulation of a borehole using load aggregation

This example demonstrates the use of the load aggregation 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.interpolate import interp1d
 13from scipy.signal import fftconvolve
 14
 15import pygfunction as gt
 16
 17
 18def main():
 19    # -------------------------------------------------------------------------
 20    # Simulation parameters
 21    # -------------------------------------------------------------------------
 22
 23    # Borehole dimensions
 24    D = 4.0             # Borehole buried depth (m)
 25    H = 150.0           # Borehole length (m)
 26    r_b = 0.075         # Borehole radius (m)
 27
 28    # Ground properties
 29    alpha = 1.0e-6      # Ground thermal diffusivity (m2/s)
 30    k_s = 2.0           # Ground thermal conductivity (W/m.K)
 31    T_g = 10.0          # Undisturbed ground temperature (degC)
 32
 33    # g-Function calculation options
 34    options = {'nSegments': 8,
 35               'disp': True}
 36
 37    # Simulation parameters
 38    dt = 3600.                  # Time step (s)
 39    tmax = 20.*8760. * 3600.    # Maximum time (s)
 40    Nt = int(np.ceil(tmax/dt))  # Number of time steps
 41    time = dt * np.arange(1, Nt+1)
 42
 43    # Evaluate heat extraction rate
 44    Q_b = synthetic_load(time/3600.)
 45
 46    # Load aggregation scheme
 47    LoadAgg = gt.load_aggregation.ClaessonJaved(dt, tmax)
 48
 49    # -------------------------------------------------------------------------
 50    # Calculate g-function
 51    # -------------------------------------------------------------------------
 52
 53    # The field contains only one borehole
 54    borehole = gt.boreholes.Borehole(H, D, r_b, x=0., y=0.)
 55    # Get time values needed for g-function evaluation
 56    time_req = LoadAgg.get_times_for_simulation()
 57    # Calculate g-function
 58    gFunc = gt.gfunction.gFunction(
 59        borehole, alpha, time=time_req, options=options)
 60    # Initialize load aggregation scheme
 61    LoadAgg.initialize(gFunc.gFunc / (2 * np.pi * k_s))
 62
 63    # -------------------------------------------------------------------------
 64    # Simulation
 65    # -------------------------------------------------------------------------
 66
 67    T_b = np.zeros(Nt)
 68    for i, (t, Q_b_i) in enumerate(zip(time, Q_b)):
 69        # Increment time step by (1)
 70        LoadAgg.next_time_step(t)
 71
 72        # Apply current load
 73        LoadAgg.set_current_load(Q_b_i/H)
 74
 75        # Evaluate borehole wall temperature
 76        deltaT_b = LoadAgg.temporal_superposition()
 77        T_b[i] = T_g - deltaT_b
 78
 79    # -------------------------------------------------------------------------
 80    # Calculate exact solution from convolution in the Fourier domain
 81    # -------------------------------------------------------------------------
 82
 83    # Heat extraction rate increment
 84    dQ = np.zeros(Nt)
 85    dQ[0] = Q_b[0]
 86    dQ[1:] = Q_b[1:] - Q_b[:-1]
 87    # Interpolated g-function
 88    g = interp1d(time_req, gFunc.gFunc)(time)
 89
 90    # Convolution in Fourier domain
 91    T_b_exact = T_g - fftconvolve(
 92        dQ, g / (2.0 * np.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*np.pi)*(np.cos(C*np.pi*i/84.0)-1.0) \
150                          *(np.sin(np.pi*i/84.0*(x-B)))
151    func = func*A*np.sin(np.pi/12.0*(x-B)) \
152           *np.sin(np.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*np.pi/4380.0*(x-F))+G)
156    return -y
157
158
159# Main function
160if __name__ == '__main__':
161    main()

References