Simulation of fluid temperatures in a borehole¶
This example demonstrates the use of the pipes module to predict the fluid temperature variations in a borehole 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. Three pipe configurations are compared: (1) a single U-tube, using the model of Eskilson and Claesson 3, (2) a double U-tube in parallel, using the model of Cimmino 4, and (3) a double U-tube in series, using the model of Cimmino 4.
The script is located in: pygfunction/examples/fluid_temperature.py
1# -*- coding: utf-8 -*-
2""" Example of simulation of a geothermal system and comparison between
3 single and double U-tube pipe configurations.
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 method of Claesson and Javed (2012).
9 Predicted outlet fluid temperatures of three pipe configurations
10 (single U-tube, double U-tube in series, double U-tube in parallel) are
11 compared.
12
13"""
14import matplotlib.pyplot as plt
15import numpy as np
16from scipy.constants import pi
17
18import pygfunction as gt
19
20
21def main():
22 # -------------------------------------------------------------------------
23 # Simulation parameters
24 # -------------------------------------------------------------------------
25
26 # Borehole dimensions
27 D = 4.0 # Borehole buried depth (m)
28 H = 150.0 # Borehole length (m)
29 r_b = 0.075 # Borehole radius (m)
30
31 # Pipe dimensions
32 rp_out = 0.0211 # Pipe outer radius (m)
33 rp_in = 0.0147 # Pipe inner radius (m)
34 D_s = 0.052 # Shank spacing (m)
35 epsilon = 1.0e-6 # Pipe roughness (m)
36
37 # Pipe positions
38 # Single U-tube [(x_in, y_in), (x_out, y_out)]
39 pos_single = [(-D_s, 0.), (D_s, 0.)]
40 # Double U-tube [(x_in1, y_in1), (x_in2, y_in2),
41 # (x_out1, y_out1), (x_in2, y_in2)]
42 # Note: in series configuration, fluid enters pipe (in,1), exits (out,1),
43 # then enters (in,2) and finally exits (out,2)
44 pos_double = [(-D_s, 0.), (0., -D_s), (D_s, 0.), (0., D_s)]
45
46 # Ground properties
47 alpha = 1.0e-6 # Ground thermal diffusivity (m2/s)
48 k_s = 2.0 # Ground thermal conductivity (W/m.K)
49 T_g = 10.0 # Undisturbed ground temperature (degC)
50
51 # Grout properties
52 k_g = 1.0 # Grout thermal conductivity (W/m.K)
53
54 # Pipe properties
55 k_p = 0.4 # Pipe thermal conductivity (W/m.K)
56
57 # Fluid properties
58 m_flow = 0.25 # Total fluid mass flow rate (kg/s)
59 # The fluid is propylene-glycol (20 %) at 20 degC
60 fluid = gt.media.Fluid('MPG', 20.)
61 cp_f = fluid.cp # Fluid specific isobaric heat capacity (J/kg.K)
62 den_f = fluid.rho # Fluid density (kg/m3)
63 visc_f = fluid.mu # Fluid dynamic viscosity (kg/m.s)
64 k_f = fluid.k # Fluid thermal conductivity (W/m.K)
65
66 # g-Function calculation options
67 options = {'nSegments': 8,
68 'disp': True}
69
70 # Simulation parameters
71 dt = 3600. # Time step (s)
72 tmax = 1.*8760. * 3600. # Maximum time (s)
73 Nt = int(np.ceil(tmax/dt)) # Number of time steps
74 time = dt * np.arange(1, Nt+1)
75
76 # Evaluate heat extraction rate
77 Q = synthetic_load(time/3600.)
78
79 # Load aggregation scheme
80 LoadAgg = gt.load_aggregation.ClaessonJaved(dt, tmax)
81
82 # -------------------------------------------------------------------------
83 # Calculate g-function
84 # -------------------------------------------------------------------------
85
86 # The field contains only one borehole
87 borehole = gt.boreholes.Borehole(H, D, r_b, x=0., y=0.)
88 boreField = [borehole]
89 # Get time values needed for g-function evaluation
90 time_req = LoadAgg.get_times_for_simulation()
91 # Calculate g-function
92 gFunc = gt.gfunction.gFunction(
93 boreField, alpha, time=time_req, options=options)
94 # gt.gfunction.uniform_temperature(boreField, time_req, alpha,
95 # nSegments=nSegments)
96 # Initialize load aggregation scheme
97 LoadAgg.initialize(gFunc.gFunc/(2*pi*k_s))
98
99 # -------------------------------------------------------------------------
100 # Initialize pipe models
101 # -------------------------------------------------------------------------
102
103 # Pipe thermal resistance
104 R_p = gt.pipes.conduction_thermal_resistance_circular_pipe(
105 rp_in, rp_out, k_p)
106 # Fluid to inner pipe wall thermal resistance (Single U-tube and double
107 # U-tube in series)
108 h_f = gt.pipes.convective_heat_transfer_coefficient_circular_pipe(
109 m_flow, rp_in, visc_f, den_f, k_f, cp_f, epsilon)
110 R_f_ser = 1.0/(h_f*2*pi*rp_in)
111 # Fluid to inner pipe wall thermal resistance (Double U-tube in parallel)
112 h_f = gt.pipes.convective_heat_transfer_coefficient_circular_pipe(
113 m_flow/2, rp_in, visc_f, den_f, k_f, cp_f, epsilon)
114 R_f_par = 1.0/(h_f*2*pi*rp_in)
115
116 # Single U-tube
117 SingleUTube = gt.pipes.SingleUTube(pos_single, rp_in, rp_out,
118 borehole, k_s, k_g, R_f_ser + R_p)
119 # Double U-tube (parallel)
120 DoubleUTube_par = gt.pipes.MultipleUTube(pos_double, rp_in, rp_out,
121 borehole, k_s, k_g, R_f_par + R_p,
122 nPipes=2, config='parallel')
123 # Double U-tube (series)
124 DoubleUTube_ser = gt.pipes.MultipleUTube(pos_double, rp_in, rp_out,
125 borehole, k_s, k_g, R_f_ser + R_p,
126 nPipes=2, config='series')
127
128 # -------------------------------------------------------------------------
129 # Simulation
130 # -------------------------------------------------------------------------
131
132 T_b = np.zeros(Nt)
133 T_f_in_single = np.zeros(Nt)
134 T_f_in_double_par = np.zeros(Nt)
135 T_f_in_double_ser = np.zeros(Nt)
136 T_f_out_single = np.zeros(Nt)
137 T_f_out_double_par = np.zeros(Nt)
138 T_f_out_double_ser = np.zeros(Nt)
139 for i, (t, Q_b_i) in enumerate(zip(time, Q)):
140 # Increment time step by (1)
141 LoadAgg.next_time_step(t)
142
143 # Apply current load
144 LoadAgg.set_current_load(Q_b_i/H)
145
146 # Evaluate borehole wall temperature
147 deltaT_b = LoadAgg.temporal_superposition()
148 T_b[i] = T_g - deltaT_b
149
150 # Evaluate inlet fluid temperature
151 T_f_in_single[i] = SingleUTube.get_inlet_temperature(
152 Q[i], T_b[i], m_flow, cp_f)
153 T_f_in_double_par[i] = DoubleUTube_par.get_inlet_temperature(
154 Q[i], T_b[i], m_flow, cp_f)
155 T_f_in_double_ser[i] = DoubleUTube_ser.get_inlet_temperature(
156 Q[i], T_b[i], m_flow, cp_f)
157
158 # Evaluate outlet fluid temperature
159 T_f_out_single[i] = SingleUTube.get_outlet_temperature(
160 T_f_in_single[i], T_b[i], m_flow, cp_f)
161 T_f_out_double_par[i] = DoubleUTube_par.get_outlet_temperature(
162 T_f_in_double_par[i], T_b[i], m_flow, cp_f)
163 T_f_out_double_ser[i] = DoubleUTube_ser.get_outlet_temperature(
164 T_f_in_double_ser[i], T_b[i], m_flow, cp_f)
165
166 # -------------------------------------------------------------------------
167 # Plot hourly heat extraction rates and temperatures
168 # -------------------------------------------------------------------------
169
170 # Configure figure and axes
171 fig = gt.utilities._initialize_figure()
172
173 ax1 = fig.add_subplot(211)
174 # Axis labels
175 ax1.set_xlabel(r'Time [hours]')
176 ax1.set_ylabel(r'Total heat extraction rate [W]')
177 gt.utilities._format_axes(ax1)
178
179 # Plot heat extraction rates
180 hours = np.arange(1, Nt+1) * dt / 3600.
181 ax1.plot(hours, Q)
182
183 ax2 = fig.add_subplot(212)
184 # Axis labels
185 ax2.set_xlabel(r'Time [hours]')
186 ax2.set_ylabel(r'Temperature [degC]')
187 gt.utilities._format_axes(ax2)
188
189 # Plot temperatures
190 ax2.plot(hours, T_b, 'k-', lw=1.5, label='Borehole wall')
191 ax2.plot(hours, T_f_out_single, '--',
192 label='Outlet, single U-tube')
193 ax2.plot(hours, T_f_out_double_par, '-.',
194 label='Outlet, double U-tube (parallel)')
195 ax2.plot(hours, T_f_out_double_ser, ':',
196 label='Outlet, double U-tube (series)')
197 ax2.legend()
198
199 # Adjust to plot window
200 plt.tight_layout()
201
202 # -------------------------------------------------------------------------
203 # Plot fluid temperature profiles
204 # -------------------------------------------------------------------------
205
206 # Evaluate temperatures at nz evenly spaced depths along the borehole
207 # at the (it+1)-th time step
208 nz = 20
209 it = 8724
210 z = np.linspace(0., H, num=nz)
211 T_f_single = SingleUTube.get_temperature(z,
212 T_f_in_single[it],
213 T_b[it],
214 m_flow,
215 cp_f)
216 T_f_double_par = DoubleUTube_par.get_temperature(z,
217 T_f_in_double_par[it],
218 T_b[it],
219 m_flow,
220 cp_f)
221 T_f_double_ser = DoubleUTube_ser.get_temperature(z,
222 T_f_in_double_ser[it],
223 T_b[it],
224 m_flow,
225 cp_f)
226
227 # Configure figure and axes
228 fig = gt.utilities._initialize_figure()
229
230 ax3 = fig.add_subplot(131)
231 # Axis labels
232 ax3.set_xlabel(r'Temperature [degC]')
233 ax3.set_ylabel(r'Depth from borehole head [m]')
234 gt.utilities._format_axes(ax3)
235
236 # Plot temperatures
237 ax3.plot(np.array([T_b[it], T_b[it]]), np.array([0., H]), 'k--')
238 ax3.plot(T_f_single, z, 'b-')
239 ax3.legend(['Borehole wall', 'Fluid'])
240
241 ax4 = fig.add_subplot(132)
242 # Axis labels
243 ax4.set_xlabel(r'Temperature [degC]')
244 ax4.set_ylabel(r'Depth from borehole head [m]')
245 gt.utilities._format_axes(ax4)
246
247 # Plot temperatures
248 ax4.plot(T_f_double_par, z, 'b-')
249 ax4.plot(np.array([T_b[it], T_b[it]]), np.array([0., H]), 'k--')
250
251 ax5 = fig.add_subplot(133)
252 # Axis labels
253 ax5.set_xlabel(r'Temperature [degC]')
254 ax5.set_ylabel(r'Depth from borehole head [m]')
255 gt.utilities._format_axes(ax5)
256
257 # Plot temperatures
258 ax5.plot(T_f_double_ser, z, 'b-')
259 ax5.plot(np.array([T_b[it], T_b[it]]), np.array([0., H]), 'k--')
260
261 # Reverse y-axes
262 ax3.set_ylim(ax3.get_ylim()[::-1])
263 ax4.set_ylim(ax4.get_ylim()[::-1])
264 ax5.set_ylim(ax5.get_ylim()[::-1])
265 # Adjust to plot window
266 plt.tight_layout()
267
268 return
269
270
271def synthetic_load(x):
272 """
273 Synthetic load profile of Bernier et al. (2004).
274
275 Returns load y (in watts) at time x (in hours).
276 """
277 A = 2000.0
278 B = 2190.0
279 C = 80.0
280 D = 2.0
281 E = 0.01
282 F = 0.0
283 G = 0.95
284
285 func = (168.0-C)/168.0
286 for i in [1, 2, 3]:
287 func += 1.0/(i*pi)*(np.cos(C*pi*i/84.0)-1.0) \
288 *(np.sin(pi*i/84.0*(x-B)))
289 func = func*A*np.sin(pi/12.0*(x-B)) \
290 *np.sin(pi/4380.0*(x-B))
291
292 y = func + (-1.0)**np.floor(D/8760.0*(x-B))*abs(func) \
293 + E*(-1.0)**np.floor(D/8760.0*(x-B))/np.sign(np.cos(D*pi/4380.0*(x-F))+G)
294 return -y
295
296
297# Main function
298if __name__ == '__main__':
299 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
Eskilson, P., & Claesson, J. (1988). Simulation model for thermally interacting heat extraction boreholes. Numerical Heat Transfer 13 : 149-165.
- 4(1,2)
Cimmino, M. (2016). Fluid and borehole wall temperature profiles in vertical geothermal boreholes with multiple U-tubes. Renewable Energy 96 : 137-147.