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