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.