Simulation of fluid temperatures in a field of multiple boreholes¶
This example demonstrates the use of the networks module to predict the fluid temperature variations in a bore field with known heat extraction rates.
The g-function of a bore field is first calculated using the equal inlet fluid temperature boundary condition 1. Then, the borehole wall temperature variations are calculated using the load aggregation scheme of Claesson and Javed 2. The time-variation of heat extraction rates is given by the synthetic load profile of Bernier et al. 3. Predicted inlet and outlet fluid temperatures of double U-tube boreholes are calculated using the model of Cimmino 4.
The script is located in: pygfunction/examples/fluid_temperature_multiple_boreholes.py
1# -*- coding: utf-8 -*-
2""" Example of simulation of a geothermal system with multiple boreholes.
3
4 The g-function of a bore field is calculated for boundary condition of
5 mixed inlet fluid temperature into the boreholes. 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 Predicted outlet fluid temperatures of double U-tube borehole are
9 evaluated.
10
11"""
12import matplotlib.pyplot as plt
13import numpy as np
14from scipy.constants import pi
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 # Bore field geometry (rectangular array)
30 N_1 = 6 # Number of boreholes in the x-direction (columns)
31 N_2 = 4 # Number of boreholes in the y-direction (rows)
32 B = 7.5 # Borehole spacing, in both directions (m)
33
34 # Pipe dimensions
35 r_out = 0.0211 # Pipe outer radius (m)
36 r_in = 0.0147 # Pipe inner radius (m)
37 D_s = 0.052 # Shank spacing (m)
38 epsilon = 1.0e-6 # Pipe roughness (m)
39
40 # Pipe positions
41 # Double U-tube [(x_in1, y_in1), (x_in2, y_in2),
42 # (x_out1, y_out1), (x_out2, y_out2)]
43 pos = [(-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_borehole = 0.25 # Total fluid mass flow rate per borehole (kg/s)
58 m_flow_network = m_flow_borehole*N_1*N_2 # 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 rho_f = fluid.rho # Fluid density (kg/m3)
63 mu_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 # Load aggregation scheme
77 LoadAgg = gt.load_aggregation.ClaessonJaved(dt, tmax)
78
79 # -------------------------------------------------------------------------
80 # Initialize bore field and pipe models
81 # -------------------------------------------------------------------------
82
83 # The field is a retangular array
84 boreField = gt.boreholes.rectangle_field(N_1, N_2, B, B, H, D, r_b)
85 nBoreholes = len(boreField)
86
87 # Pipe thermal resistance
88 R_p = gt.pipes.conduction_thermal_resistance_circular_pipe(
89 r_in, r_out, k_p)
90
91 # Fluid to inner pipe wall thermal resistance (Double U-tube in parallel)
92 m_flow_pipe = m_flow_borehole/2
93 h_f = gt.pipes.convective_heat_transfer_coefficient_circular_pipe(
94 m_flow_pipe, r_in, mu_f, rho_f, k_f, cp_f, epsilon)
95 R_f = 1.0/(h_f*2*pi*r_in)
96
97 # Double U-tube (parallel), same for all boreholes in the bore field
98 UTubes = []
99 for borehole in boreField:
100 UTube = gt.pipes.MultipleUTube(
101 pos, r_in, r_out, borehole, k_s, k_g, R_f + R_p,
102 nPipes=2, config='parallel')
103 UTubes.append(UTube)
104 # Build a network object from the list of UTubes
105 network = gt.networks.Network(
106 boreField, UTubes, m_flow_network=m_flow_network, cp_f=cp_f)
107
108 # -------------------------------------------------------------------------
109 # Calculate g-function
110 # -------------------------------------------------------------------------
111
112 # Get time values needed for g-function evaluation
113 time_req = LoadAgg.get_times_for_simulation()
114 # Calculate g-function
115 gFunc = gt.gfunction.gFunction(
116 network, alpha, time=time_req, boundary_condition='MIFT',
117 options=options)
118 # Initialize load aggregation scheme
119 LoadAgg.initialize(gFunc.gFunc/(2*pi*k_s))
120
121 # -------------------------------------------------------------------------
122 # Simulation
123 # -------------------------------------------------------------------------
124
125 # Evaluate heat extraction rate
126 Q_tot = nBoreholes*synthetic_load(time/3600.)
127
128 T_b = np.zeros(Nt)
129 T_f_in = np.zeros(Nt)
130 T_f_out = np.zeros(Nt)
131 for i, (t, Q_i) in enumerate(zip(time, Q_tot)):
132 # Increment time step by (1)
133 LoadAgg.next_time_step(t)
134
135 # Apply current load (in watts per meter of borehole)
136 Q_b = Q_i/nBoreholes
137 LoadAgg.set_current_load(Q_b/H)
138
139 # Evaluate borehole wall temperature
140 deltaT_b = LoadAgg.temporal_superposition()
141 T_b[i] = T_g - deltaT_b
142
143 # Evaluate inlet fluid temperature (all boreholes are the same)
144 T_f_in[i] = network.get_network_inlet_temperature(
145 Q_tot[i], T_b[i], m_flow_network, cp_f, nSegments=1)
146
147 # Evaluate outlet fluid temperature
148 T_f_out[i] = network.get_network_outlet_temperature(
149 T_f_in[i], T_b[i], m_flow_network, cp_f, nSegments=1)
150
151 # -------------------------------------------------------------------------
152 # Plot hourly heat extraction rates and temperatures
153 # -------------------------------------------------------------------------
154
155 # Configure figure and axes
156 fig = gt.utilities._initialize_figure()
157
158 ax1 = fig.add_subplot(211)
159 # Axis labels
160 ax1.set_xlabel(r'Time [hours]')
161 ax1.set_ylabel(r'Total heat extraction rate [W]')
162 gt.utilities._format_axes(ax1)
163
164 # Plot heat extraction rates
165 hours = np.arange(1, Nt+1) * dt / 3600.
166 ax1.plot(hours, Q_tot)
167
168 ax2 = fig.add_subplot(212)
169 # Axis labels
170 ax2.set_xlabel(r'Time [hours]')
171 ax2.set_ylabel(r'Temperature [degC]')
172 gt.utilities._format_axes(ax2)
173
174 # Plot temperatures
175 ax2.plot(hours, T_b, label='Borehole wall')
176 ax2.plot(hours, T_f_out, '-.',
177 label='Outlet, double U-tube (parallel)')
178 ax2.legend()
179
180 # Adjust to plot window
181 plt.tight_layout()
182
183 # -------------------------------------------------------------------------
184 # Plot fluid temperature profiles
185 # -------------------------------------------------------------------------
186
187 # Evaluate temperatures at nz evenly spaced depths along the borehole
188 # at the (it+1)-th time step
189 nz = 20
190 it = 8724
191 z = np.linspace(0., H, num=nz)
192 T_f = UTubes[0].get_temperature(
193 z, T_f_in[it], T_b[it], m_flow_borehole, cp_f)
194
195
196 # Configure figure and axes
197 fig = gt.utilities._initialize_figure()
198
199 ax3 = fig.add_subplot(111)
200 # Axis labels
201 ax3.set_xlabel(r'Temperature [degC]')
202 ax3.set_ylabel(r'Depth from borehole head [m]')
203 gt.utilities._format_axes(ax3)
204
205 # Plot temperatures
206 pltFlu = ax3.plot(T_f, z, 'b-', label='Fluid')
207 pltWal = ax3.plot(np.array([T_b[it], T_b[it]]), np.array([0., H]),
208 'k--', label='Borehole wall')
209 ax3.legend(handles=[pltFlu[0]]+pltWal)
210
211 # Reverse y-axes
212 ax3.set_ylim(ax3.get_ylim()[::-1])
213 # Adjust to plot window
214 plt.tight_layout()
215
216 return
217
218
219def synthetic_load(x):
220 """
221 Synthetic load profile of Bernier et al. (2004).
222
223 Returns load y (in watts) at time x (in hours).
224 """
225 A = 2000.0
226 B = 2190.0
227 C = 80.0
228 D = 2.0
229 E = 0.01
230 F = 0.0
231 G = 0.95
232
233 func = (168.0-C)/168.0
234 for i in [1, 2, 3]:
235 func += 1.0/(i*pi)*(np.cos(C*pi*i/84.0)-1.0) \
236 *(np.sin(pi*i/84.0*(x-B)))
237 func = func*A*np.sin(pi/12.0*(x-B)) \
238 *np.sin(pi/4380.0*(x-B))
239
240 y = func + (-1.0)**np.floor(D/8760.0*(x-B))*abs(func) \
241 + E*(-1.0)**np.floor(D/8760.0*(x-B))/np.sign(np.cos(D*pi/4380.0*(x-F))+G)
242 return -y
243
244
245# Main function
246if __name__ == '__main__':
247 main()
References
- 1
Cimmino, M. (2015). The effects of borehole thermal resistances and fluid flow rate on the g-functions of geothermal bore fields. International Journal of Heat and Mass Transfer, 91, 1119-1127.
- 2
Claesson, J., & Javed, S. (2011). A load-aggregation method to calculate extraction temperatures of borehole heat exchangers. ASHRAE Transactions, 118 (1): 530–539.
- 3
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.
- 4
Cimmino, M. (2016). Fluid and borehole wall temperature profiles in vertical geothermal boreholes with multiple U-tubes. Renewable Energy 96 : 137-147.