Simulation of fluid temperatures in a field of series-connected boreholes and reversible flow direction

This example demonstrates the use of the networks module to predict the fluid temperature variations in a bore field with series-connected boreholes and reversible flow direction.

The variable fluid mass flow rates g-functions of a bore field are first calculated using the mixed inlet fluid temperature boundary condition [1]. Then, the effective borehole wall temperature variations are calculated using the load aggregation scheme of Claesson and Javed [2]. g-Functions used in the temporal superposition are interpolated with regards to the fluid mass flow rate at the moment of heat extraction.

The script is located in: pygfunction/examples/fluid_temperature_reversible_flow_direction.py

  1# -*- coding: utf-8 -*-
  2""" Example of simulation of a geothermal system with multiple series-
  3    connected boreholes and reversible flow direction.
  4
  5    The g-function of a bore field is calculated for boundary condition of
  6    mixed inlet fluid temperature into the boreholes and a combination of two
  7    fluid mass flow rates in two opposing flow direction. Then, the borehole
  8    wall temperature variations resulting from a time-varying load profile
  9    are simulated using the aggregation method of Claesson and Javed (2012).
 10    Predicted outlet fluid temperatures of double U-tube borehole are
 11    evaluated.
 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    # Bore field geometry (rectangular array)
 31    N_1 = 1             # Number of boreholes in the x-direction (columns)
 32    N_2 = 6             # Number of boreholes in the y-direction (rows)
 33    B = 7.5             # Borehole spacing, in both directions (m)
 34
 35    # Pipe dimensions
 36    r_out = 33.6e-3/2   # Pipe outer radius (m)
 37    r_in = 27.4e-3/2    # Pipe inner radius (m)
 38    D_s = 0.11/2        # Shank spacing (m)
 39    epsilon = 1.0e-6    # Pipe roughness (m)
 40
 41    # Pipe positions
 42    # Single U-tube [(x_in, y_in), (x_out, y_out)]
 43    pos_pipes = [(-D_s, 0.), (D_s, 0.)]
 44
 45    # Ground properties
 46    alpha = 2.5/2.2e6   # Ground thermal diffusivity (m2/s)
 47    k_s = 2.5           # Ground thermal conductivity (W/m.K)
 48    T_g = 10.           # Undisturbed ground temperature (degC)
 49
 50    # Grout properties
 51    k_g = 1.5           # Grout thermal conductivity (W/m.K)
 52
 53    # Pipe properties
 54    k_p = 0.42          # Pipe thermal conductivity (W/m.K)
 55
 56    # Fluid properties
 57    m_flow_borehole_min = 0.5   # Minimum fluid mass flow rate per borehole (kg/s)
 58    m_flow_borehole_max = 1.2   # Maximum fluid mass flow rate per borehole (kg/s)
 59    m_flow_borehole = np.array([
 60        -m_flow_borehole_max,
 61        -m_flow_borehole_min,
 62        m_flow_borehole_min,
 63        m_flow_borehole_max
 64        ])
 65    nModes = len(m_flow_borehole)
 66    # Total fluid mass flow rate (kg/s)
 67    m_flow_network = m_flow_borehole
 68    # The fluid is propylene-glycol (20 %) at 20 degC
 69    fluid = gt.media.Fluid('MPG', 20.)
 70    cp_f = fluid.cp     # Fluid specific isobaric heat capacity (J/kg.K)
 71    rho_f = fluid.rho   # Fluid density (kg/m3)
 72    mu_f = fluid.mu     # Fluid dynamic viscosity (kg/m.s)
 73    k_f = fluid.k       # Fluid thermal conductivity (W/m.K)
 74
 75    # g-Function calculation options
 76    options = {'approximate_FLS': True,
 77               'nSegments': 8,
 78               'disp': True}
 79
 80    # Simulation parameters
 81    dt = 600.                   # Time step (s)
 82    tmax = 48 * 3600.           # Maximum time (s)
 83    Nt = int(np.ceil(tmax/dt))  # Number of time steps
 84    time = dt * np.arange(1, Nt+1)
 85    time_h = time / 3600.
 86
 87    # Load aggregation scheme
 88    LoadAgg = gt.load_aggregation.ClaessonJaved(dt, tmax, nSources=nModes)
 89
 90    # -------------------------------------------------------------------------
 91    # Initialize bore field and pipe models
 92    # -------------------------------------------------------------------------
 93
 94    # The field is a rectangular array
 95    borefield = gt.borefield.Borefield.rectangle_field(
 96        N_1, N_2, B, B, H, D, r_b)
 97    nBoreholes = len(borefield)
 98    H_tot = np.sum(borefield.H)
 99
100    # Boreholes are connected in series
101    bore_connectivity = [i-1 for i in range(nBoreholes)]
102
103    # Pipe thermal resistance
104    R_p = gt.pipes.conduction_thermal_resistance_circular_pipe(
105            r_in, r_out, k_p)
106
107    # Fluid to inner pipe wall thermal resistance (Double U-tube in parallel)
108    m_flow_pipe = np.max(np.abs(m_flow_borehole))
109    h_f = gt.pipes.convective_heat_transfer_coefficient_circular_pipe(
110            m_flow_pipe, r_in, mu_f, rho_f, k_f, cp_f, epsilon)
111    R_f = 1.0 / (h_f * 2 * np.pi * r_in)
112
113    # Double U-tube (parallel), same for all boreholes in the bore field
114    UTubes = []
115    for borehole in borefield:
116        UTube = gt.pipes.SingleUTube(
117            pos_pipes, r_in, r_out, borehole, k_s, k_g, R_f + R_p)
118        UTubes.append(UTube)
119    # Build a network object from the list of UTubes
120    network = gt.networks.Network(
121        borefield, UTubes, bore_connectivity=bore_connectivity)
122
123    # -------------------------------------------------------------------------
124    # Calculate g-function
125    # -------------------------------------------------------------------------
126
127    # Get time values needed for g-function evaluation
128    time_req = LoadAgg.get_times_for_simulation()
129    # Calculate g-function
130    gFunc = gt.gfunction.gFunction(
131        network, alpha, time=time_req, boundary_condition='MIFT',
132        m_flow_network=m_flow_network, cp_f=cp_f, options=options)
133    # Initialize load aggregation scheme
134    LoadAgg.initialize(gFunc.gFunc / (2 * np.pi * k_s))
135
136    # -------------------------------------------------------------------------
137    # Simulation
138    # -------------------------------------------------------------------------
139
140    # Sinusoidal variation between -80 and 80 W/m with minimum of 10 W/m
141    Q_min = -80. * H * nBoreholes
142    Q_max = 80. * H * nBoreholes
143    Q_small = 10. * H * nBoreholes
144    Q_avg = 0.5 * (Q_min + Q_max)
145    Q_amp = 0.5 * (Q_max - Q_min)
146    Q_tot = -Q_avg - Q_amp * np.sin(2 * np.pi * time_h / 24)
147    # Mass flow rate is proportional to the absolute load
148    m_flow_min = np.min(np.abs(m_flow_network))
149    m_flow_max = np.max(np.abs(m_flow_network))
150    m_flow = np.sign(Q_tot) * np.maximum(
151        m_flow_min,
152        np.minimum(m_flow_max, m_flow_max * np.abs(Q_tot) / Q_max))
153    # The minimum heat load is 10 W/m and the mass flow rate is 0 when Q_tot=0
154    Q_tot[np.abs(Q_tot) < Q_small] = 0.
155    m_flow[np.abs(Q_tot) < Q_small] = 0.
156
157    T_b = np.zeros(Nt)
158    T_f_in = np.zeros(Nt)
159    T_f_out = np.zeros(Nt)
160    for i, (t, Q_i) in enumerate(zip(time, Q_tot)):
161        # Increment time step by (1)
162        LoadAgg.next_time_step(t)
163
164        # Solve for temperatures and heat extraction rates (variable mass flow rate)
165        if np.abs(Q_tot[i]) >= Q_small:
166            xi = np.zeros(nModes)
167            iMode = np.argmax(m_flow[i] <= m_flow_network)
168            if iMode > 0:
169                xi[iMode-1] = (m_flow_network[iMode] - m_flow[i]) / (m_flow_network[iMode] - m_flow_network[iMode-1])
170                xi[iMode] = (m_flow[i] - m_flow_network[iMode-1]) / (m_flow_network[iMode] - m_flow_network[iMode-1])
171            else:
172                xi[0] = 1.
173            LoadAgg.set_current_load(Q_tot[i] / H_tot * xi)
174            deltaT_b = LoadAgg.temporal_superposition()
175            T_b[i] = T_g - deltaT_b @ xi
176            T_f_in[i] = network.get_network_inlet_temperature(
177                Q_tot[i], T_b[i], m_flow[i], cp_f, nSegments=1, segment_ratios=None)
178            T_f_out[i] = network.get_network_outlet_temperature(
179                T_f_in[i], T_b[i], m_flow[i], cp_f, nSegments=1, segment_ratios=None)
180        else:
181            T_b[i] = np.nan
182            T_f_in[i] = np.nan
183            T_f_out[i] = np.nan
184
185    # -------------------------------------------------------------------------
186    # Plot hourly heat extraction rates and temperatures
187    # -------------------------------------------------------------------------
188
189    # Configure figure and axes
190    fig = gt.utilities._initialize_figure()
191
192    ax1 = fig.add_subplot(221)
193    # Axis labels
194    ax1.set_xlabel(r'Time [hours]')
195    ax1.set_ylabel(r'Total heat extraction rate [W/m]')
196    gt.utilities._format_axes(ax1)
197
198    # Plot heat extraction rates
199    hours = np.arange(1, Nt+1) * dt / 3600.
200    ax1.plot(hours, Q_tot / H_tot)
201
202    ax2 = fig.add_subplot(222)
203    # Axis labels
204    ax2.set_xlabel(r'Time [hours]')
205    ax2.set_ylabel(r'Fluid mass flow rate [kg/s]')
206    gt.utilities._format_axes(ax2)
207
208    # Plot temperatures
209    ax2.plot(hours, m_flow)
210
211    ax3 = fig.add_subplot(223)
212    # Axis labels
213    ax3.set_xlabel(r'Time [hours]')
214    ax3.set_ylabel(r'Temperature [degC]')
215    gt.utilities._format_axes(ax3)
216
217    # Plot temperatures
218    ax3.plot(hours, T_b, label='Borehole wall')
219    ax3.plot(hours, T_f_out, '-.',
220             label='Outlet')
221    ax3.legend()
222
223    # Adjust to plot window
224    plt.tight_layout()
225
226    return
227
228
229# Main function
230if __name__ == '__main__':
231    main()

References