diff --git a/CHANGELOG.md b/CHANGELOG.md index 56e2dd1a..d4610574 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,7 @@ * [Commit 2bd12bd](https://github.com/MassimoCimmino/pygfunction/commit/2bd12bd254928889431366c2ddd38539e246ef05) - Implemented `UTube.visualize_pipes()` class method. * [Issue 30](https://github.com/MassimoCimmino/pygfunction/issues/30) - Laminar regime is now considered for calculation of convection heat transfer coefficient. +* [Issue 32](https://github.com/MassimoCimmino/pygfunction/issues/32) - g-Functions for bore fields with mixed series and parallel connections between boreholes. ### Bug fixes diff --git a/doc/source/example_bore_field_thermal_resistance.rst b/doc/source/example_bore_field_thermal_resistance.rst new file mode 100644 index 00000000..596bd629 --- /dev/null +++ b/doc/source/example_bore_field_thermal_resistance.rst @@ -0,0 +1,26 @@ +.. examples: + +****************************************************** +Calculation of effective bore field thermal resistance +****************************************************** + +This example demonstrates the use of the +:py:func:`.pipes.field_thermal_resistance` function to evaluate the bore field +thermal resistance. The concept effective bore field thermal is detailed by +Cimmino [1]_ + +The following script evaluates the effective bore field thermal resistance +for fields of 1 to 5 series-connected boreholes with fluid flow rates ranging +from 0.01 kg/s ti 1.00 kg/s. + +The script is located in: +`pygfunction/examples/bore_field_thermal_resistance.py` + +.. literalinclude:: ../../pygfunction/examples/bore_field_thermal_resistance.py + :language: python + :linenos: + +.. rubric:: References +.. [1] Cimmino, M. (2018). g-Functions for bore fields with mixed parallel and + series connections considering the axial fluid temperature variations. + IGSHPA Research Track, Stockholm. In review. diff --git a/doc/source/example_index.rst b/doc/source/example_index.rst index 3cff004d..2437edc3 100644 --- a/doc/source/example_index.rst +++ b/doc/source/example_index.rst @@ -18,3 +18,5 @@ Examples example_equal_inlet_temperature example_multiple_independent_Utubes example_multipole_temperature + example_bore_field_thermal_resistance + example_mixed_inlet_conditions diff --git a/doc/source/example_mixed_inlet_conditions.rst b/doc/source/example_mixed_inlet_conditions.rst new file mode 100644 index 00000000..534decb4 --- /dev/null +++ b/doc/source/example_mixed_inlet_conditions.rst @@ -0,0 +1,31 @@ +.. examples: + +********************************************************************* +Calculation of g-functions with mixed parallel and series connections +********************************************************************* + +This example demonstrates the use of the :doc:`g-function ` module +and the :doc:`pipes ` module to calculate *g*-functions considering the +piping connections between the boreholes, based on the method of Cimmino [1]_. +For boreholes connected in series, it is considered the outlet fluid temperature +of the upstream borehole is equal to the inlet fluid temperature of the +downstream borehole. The total rate of heat extraction in the bore field is +constant. + +The following script generates the *g*-functions of a field of 5 equally spaced +borehole on a straight line and connected in series. The boreholes have +different lengths. The *g*-function considering piping conections is compared to +the *g*-function obtained using a boundary condition of uniform borehole wall +temperature. + +The script is located in: +`pygfunction/examples/mixed_inlet_conditions.py` + +.. literalinclude:: ../../pygfunction/examples/mixed_inlet_conditions.py + :language: python + :linenos: + +.. rubric:: References +.. [1] Cimmino, M. (2018). g-Functions for bore fields with mixed parallel and + series connections considering the axial fluid temperature variations. + IGSHPA Research Track, Stockholm. In review. diff --git a/pygfunction/boreholes.py b/pygfunction/boreholes.py index 94ba11b9..6dad9327 100644 --- a/pygfunction/boreholes.py +++ b/pygfunction/boreholes.py @@ -481,3 +481,76 @@ def visualize_field(borefield): plt.tight_layout(rect=[0, 0.0, 0.95, 1.0]) return fig + + +def _path_to_inlet(bore_connectivity, bore_index): + """ + Returns the path from a borehole to the bore field inlet. + + This function raises an error if the supplied borehole connectivity is + invalid. + + Parameters + ---------- + bore_connectivity : list + Index of fluid inlet into each borehole. -1 corresponds to a borehole + connected to the bore field inlet. + bore_index : int + Index of borehole to evaluate path. + + Returns + ------- + path : list + List of boreholes leading to the bore field inlet, starting from + borehole bore_index + + """ + # Initialize path + path = [bore_index] + # Index of borehole feeding into borehole (bore_index) + index_in = bore_connectivity[bore_index] + # Stop when bore field inlet is reached (index_in == -1) + while not index_in == -1: + # Add index of upstream borehole to path + path.append(index_in) + # Get index of next upstream borehole + index_in = bore_connectivity[index_in] + + return path + + +def _verify_bore_connectivity(bore_connectivity, nBoreholes): + """ + Verifies that borehole connectivity is valid. + + This function raises an error if the supplied borehole connectivity is + invalid. + + Parameters + ---------- + bore_connectivity : list + Index of fluid inlet into each borehole. -1 corresponds to a borehole + connected to the bore field inlet. + nBoreholes : int + Number of boreholes in the bore field. + + """ + if not len(bore_connectivity) == nBoreholes: + raise ValueError( + 'The length of the borehole connectivity list does not correspond ' + 'to the number of boreholes in the bore field.') + # Cycle through each borehole and verify that connections lead to -1 + # (-1 is the bore field inlet) + for i in range(nBoreholes): + n = 0 # Initialize step counter + # Index of borehole feeding into borehole i + index_in = bore_connectivity[i] + # Stop when bore field inlet is reached (index_in == -1) + while not index_in == -1: + index_in = bore_connectivity[index_in] + n += 1 # Increment step counter + # Raise error if n exceeds the number of boreholes + if n > nBoreholes: + raise ValueError( + 'The borehole connectivity list is invalid.') + return diff --git a/pygfunction/examples/bore_field_thermal_resistance.py b/pygfunction/examples/bore_field_thermal_resistance.py new file mode 100644 index 00000000..a2fdf04f --- /dev/null +++ b/pygfunction/examples/bore_field_thermal_resistance.py @@ -0,0 +1,136 @@ +# -*- coding: utf-8 -*- +""" Example of calculation of effective bore field thermal resistance. + + The effective bore field thermal resistance of fields of up to 5 boreholes + of equal lengths connected in series is calculated for various fluid flow + rates. + +""" +from __future__ import division, print_function, absolute_import + +import matplotlib.pyplot as plt +from matplotlib.ticker import AutoMinorLocator +import numpy as np +from scipy import pi + +import pygfunction as gt + + +def main(): + # ------------------------------------------------------------------------- + # Simulation parameters + # ------------------------------------------------------------------------- + + # Number of boreholes + nBoreholes = 5 + # Borehole dimensions + D = 4.0 # Borehole buried depth (m) + # Borehole length (m) + H = 150. + r_b = 0.075 # Borehole radius (m) + B = 7.5 # Borehole spacing (m) + + # Pipe dimensions + rp_out = 0.02 # Pipe outer radius (m) + rp_in = 0.015 # Pipe inner radius (m) + D_s = 0.05 # Shank spacing (m) + epsilon = 1.0e-6 # Pipe roughness (m) + + # Pipe positions + # Single U-tube [(x_in, y_in), (x_out, y_out)] + pos_pipes = [(-D_s, 0.), (D_s, 0.)] + + # Ground properties + k_s = 2.0 # Ground thermal conductivity (W/m.K) + + # Grout properties + k_g = 1.0 # Grout thermal conductivity (W/m.K) + + # Pipe properties + k_p = 0.4 # Pipe thermal conductivity (W/m.K) + + # Fluid properties + # Total fluid mass flow rate per borehole (kg/s), from 0.01 kg/s to 1 kg/s + m_flow_boreholes = 10**np.arange(-2, 0.001, 0.05) + cp_f = 4000. # Fluid specific isobaric heat capacity (J/kg.K) + den_f = 1015. # Fluid density (kg/m3) + visc_f = 0.002 # Fluid dynamic viscosity (kg/m.s) + k_f = 0.5 # Fluid thermal conductivity (W/m.K) + + # ------------------------------------------------------------------------- + # Borehole field + # ------------------------------------------------------------------------- + + boreField = [] + bore_connectivity = [] + for i in range(nBoreholes): + x = i*B + borehole = gt.boreholes.Borehole(H, D, r_b, x, 0.) + boreField.append(borehole) + # Boreholes are connected in series: The index of the upstream + # borehole is that of the previous borehole + bore_connectivity.append(i - 1) + + # ------------------------------------------------------------------------- + # Evaluate the effective bore field thermal resistance + # ------------------------------------------------------------------------- + + # Initialize result array + R = np.zeros((nBoreholes, len(m_flow_boreholes))) + for i in range(nBoreholes): + for j in range(len(m_flow_boreholes)): + nBoreholes = i + 1 + m_flow = m_flow_boreholes[j] + + # Pipe thermal resistance + R_p = gt.pipes.conduction_thermal_resistance_circular_pipe( + rp_in, rp_out, k_p) + # Fluid to inner pipe wall thermal resistance (Single U-tube) + h_f = gt.pipes.convective_heat_transfer_coefficient_circular_pipe( + m_flow, rp_in, visc_f, den_f, k_f, cp_f, epsilon) + R_f = 1.0/(h_f*2*pi*rp_in) + + # Single U-tube, same for all boreholes in the bore field + UTubes = [] + for borehole in boreField: + SingleUTube = gt.pipes.SingleUTube(pos_pipes, rp_in, rp_out, + borehole, k_s, k_g, R_f + R_p) + UTubes.append(SingleUTube) + + # Effective bore field thermal resistance + R_field = gt.pipes.field_thermal_resistance( + UTubes[:nBoreholes], bore_connectivity[:nBoreholes], + m_flow, cp_f) + # Add to result array + R[i,j] = R_field + + # ------------------------------------------------------------------------- + # Plot bore field thermal resistances + # ------------------------------------------------------------------------- + + plt.rc('figure') + fig = plt.figure() + ax1 = fig.add_subplot(111) + # Bore field thermal resistances + ax1.plot(m_flow_boreholes, R[0,:], 'k-', lw=1.5, label='1 borehole') + ax1.plot(m_flow_boreholes, R[2,:], 'r--', lw=1.5, label='3 boreholes') + ax1.plot(m_flow_boreholes, R[4,:], 'b-.', lw=1.5, label='5 boreholes') + ax1.legend() + # Axis labels + ax1.set_xlabel(r'$\dot{m}$ [kg/s]') + ax1.set_ylabel(r'$R^*_{field}$ [m.K/W]') + # Axis limits + ax1.set_xlim([0., 1.]) + ax1.set_ylim([0., 1.]) + # Show minor ticks + ax1.xaxis.set_minor_locator(AutoMinorLocator()) + ax1.yaxis.set_minor_locator(AutoMinorLocator()) + # Adjust to plot window + plt.tight_layout() + + return + + +# Main function +if __name__ == '__main__': + main() diff --git a/pygfunction/examples/mixed_inlet_conditions.py b/pygfunction/examples/mixed_inlet_conditions.py new file mode 100644 index 00000000..47c81b1e --- /dev/null +++ b/pygfunction/examples/mixed_inlet_conditions.py @@ -0,0 +1,157 @@ +# -*- coding: utf-8 -*- +""" Example of calculation of g-functions using mixed inlet temperatures. + + The g-functions of a field of 5 boreholes of different lengths connected + in series are calculated for 2 boundary conditions: (a) uniform borehole + wall temperature, and (b) series connections between boreholes. The + g-function for case (b) is based on the effective borehole wall + temperature, rather than the average borehole wall temperature. + +""" +from __future__ import division, print_function, absolute_import + +import matplotlib.pyplot as plt +from matplotlib.ticker import AutoMinorLocator +import numpy as np +from scipy import pi + +import pygfunction as gt + + +def main(): + # ------------------------------------------------------------------------- + # Simulation parameters + # ------------------------------------------------------------------------- + + # Borehole dimensions + D = 4.0 # Borehole buried depth (m) + # Borehole length (m) + H_boreholes = np.array([75.0, 100.0, 125.0, 150.0, 75.0]) + H_mean = np.mean(H_boreholes) + r_b = 0.075 # Borehole radius (m) + B = 7.5 # Borehole spacing (m) + + # Pipe dimensions + rp_out = 0.02 # Pipe outer radius (m) + rp_in = 0.015 # Pipe inner radius (m) + D_s = 0.05 # Shank spacing (m) + epsilon = 1.0e-6 # Pipe roughness (m) + + # Pipe positions + # Single U-tube [(x_in, y_in), (x_out, y_out)] + pos_pipes = [(-D_s, 0.), (D_s, 0.)] + + # Ground properties + alpha = 1.0e-6 # Ground thermal diffusivity (m2/s) + k_s = 2.0 # Ground thermal conductivity (W/m.K) + + # Grout properties + k_g = 1.0 # Grout thermal conductivity (W/m.K) + + # Pipe properties + k_p = 0.4 # Pipe thermal conductivity (W/m.K) + + # Fluid properties + m_flow = 0.25 # Total fluid mass flow rate per borehole (kg/s) + cp_f = 4000. # Fluid specific isobaric heat capacity (J/kg.K) + den_f = 1015. # Fluid density (kg/m3) + visc_f = 0.002 # Fluid dynamic viscosity (kg/m.s) + k_f = 0.5 # Fluid thermal conductivity (W/m.K) + + # Number of segments per borehole + nSegments = 12 + + # Geometrically expanding time vector. + dt = 100*3600. # Time step + tmax = 3000. * 8760. * 3600. # Maximum time + Nt = 50 # Number of time steps + ts = H_mean**2/(9.*alpha) # Bore field characteristic time + time = gt.utilities.time_geometric(dt, tmax, Nt) + + # ------------------------------------------------------------------------- + # Borehole field + # ------------------------------------------------------------------------- + + boreField = [] + bore_connectivity = [] + for i in range(len(H_boreholes)): + H = H_boreholes[i] + x = i*B + borehole = gt.boreholes.Borehole(H, D, r_b, x, 0.) + boreField.append(borehole) + # Boreholes are connected in series: The index of the upstream + # borehole is that of the previous borehole + bore_connectivity.append(i - 1) + + # ------------------------------------------------------------------------- + # Initialize pipe model + # ------------------------------------------------------------------------- + + # Pipe thermal resistance + R_p = gt.pipes.conduction_thermal_resistance_circular_pipe(rp_in, + rp_out, + k_p) + # Fluid to inner pipe wall thermal resistance (Single U-tube) + h_f = gt.pipes.convective_heat_transfer_coefficient_circular_pipe(m_flow, + rp_in, + visc_f, + den_f, + k_f, + cp_f, + epsilon) + R_f = 1.0/(h_f*2*pi*rp_in) + + # Single U-tube, same for all boreholes in the bore field + UTubes = [] + for borehole in boreField: + SingleUTube = gt.pipes.SingleUTube(pos_pipes, rp_in, rp_out, + borehole, k_s, k_g, R_f + R_p) + UTubes.append(SingleUTube) + + # ------------------------------------------------------------------------- + # Evaluate the g-functions for the borefield + # ------------------------------------------------------------------------- + + # Calculate the g-function for uniform temperature + gfunc_Tb = \ + gt.gfunction.uniform_temperature( + boreField, time, alpha, + nSegments=nSegments, disp=True) + + # Calculate the g-function for mixed inlet fluid conditions + gfunc_equal_Tf_mixed = \ + gt.gfunction.mixed_inlet_temperature( + boreField, UTubes, bore_connectivity, m_flow, cp_f, time, alpha, + nSegments=nSegments, disp=True) + + # ------------------------------------------------------------------------- + # Plot g-functions + # ------------------------------------------------------------------------- + + plt.rc('figure') + fig = plt.figure() + ax1 = fig.add_subplot(111) + # g-functions + ax1.plot(np.log(time/ts), gfunc_Tb, + 'k-', lw=1.5, label='Uniform temperature') + ax1.plot(np.log(time/ts), gfunc_equal_Tf_mixed, + 'r-.', lw=1.5, label='Mixed inlet temperature') + ax1.legend() + # Axis labels + ax1.set_xlabel(r'$ln(t/t_s)$') + ax1.set_ylabel(r'g-function') + # Axis limits + ax1.set_xlim([-10.0, 5.0]) + ax1.set_ylim([0., 12.]) + # Show minor ticks + ax1.xaxis.set_minor_locator(AutoMinorLocator()) + ax1.yaxis.set_minor_locator(AutoMinorLocator()) + # Adjust to plot window + plt.tight_layout() + + return + + +# Main function +if __name__ == '__main__': + main() diff --git a/pygfunction/gfunction.py b/pygfunction/gfunction.py index e775d06f..cf3ee5b7 100644 --- a/pygfunction/gfunction.py +++ b/pygfunction/gfunction.py @@ -7,8 +7,9 @@ from scipy.constants import pi import time as tim -from .boreholes import Borehole +from .boreholes import Borehole, _path_to_inlet, _verify_bore_connectivity from .heat_transfer import thermal_response_factors +from .pipes import field_thermal_resistance def uniform_heat_extraction(boreholes, time, alpha, use_similarities=True, @@ -466,6 +467,240 @@ def equal_inlet_temperature(boreholes, UTubes, m_flow, cp, time, alpha, return gFunction +def mixed_inlet_temperature(boreholes, UTubes, bore_connectivity, m_flow, cp, + time, alpha, method='linear', nSegments=12, + use_similarities=True, disTol=0.1, tol=1.0e-6, + processes=None, disp=False): + """ + Evaluate the g-function with mixed inlet fluid temperatures. + + This function superimposes the finite line source (FLS) solution to + estimate the g-function of a geothermal bore field. Each borehole is + modeled as a series of finite line source segments, as proposed in + [#Cimmino2018]_. The piping configurations between boreholes can be any + combination of series and parallel connections. + + Parameters + ---------- + boreholes : list of Borehole objects + List of boreholes included in the bore field. + UTubes : list of pipe objects + Model for pipes inside each borehole. + bore_connectivity : list + Index of fluid inlet into each borehole. -1 corresponds to a borehole + connected to the bore field inlet. + m_flow : array + Fluid mass flow rate in each borehole (in kg/s). + cp : float + Fluid specific isobaric heat capacity (in J/kg.K) + time : float or array + Values of time (in seconds) for which the g-function is evaluated. + alpha : float + Soil thermal diffusivity (in m2/s). + nSegments : int, optional + Number of line segments used per borehole. + Default is 12. + method : string, optional + Interpolation method used for segment-to-segment thermal response + factors. See documentation for scipy.interpolate.interp1d. + Default is 'linear'. + use_similarities : bool, optional + True if similarities are used to limit the number of FLS evaluations. + Default is True. + disTol : float, optional + Absolute tolerance (in meters) on radial distance. Two distances + (d1, d2) between two pairs of boreholes are considered equal if the + difference between the two distances (abs(d1-d2)) is below tolerance. + Default is 0.1. + tol : float, optional + Relative tolerance on length and depth. Two lenths H1, H2 + (or depths D1, D2) are considered equal if abs(H1 - H2)/H2 < tol. + Default is 1.0e-6. + processes : int, optional + Number of processors to use in calculations. If the value is set to + None, a number of processors equal to cpu_count() is used. + Default is None. + disp : bool, optional + Set to true to print progression messages. + Default is False. + + Returns + ------- + gFunction : float or array + Values of the g-function + + Examples + -------- + >>> b1 = gt.boreholes.Borehole(H=150., D=4., r_b=0.075, x=0., y=0.) + >>> b2 = gt.boreholes.Borehole(H=150., D=4., r_b=0.075, x=5., y=0.) + >>> Utube1 = gt.pipes.SingleUTube(pos=[(-0.05, 0), (0, -0.05)], + r_in=0.015, r_out=0.02, + borehole=b1,k_s=2, k_g=1, R_fp=0.1) + >>> Utube2 = gt.pipes.SingleUTube(pos=[(-0.05, 0), (0, -0.05)], + r_in=0.015, r_out=0.02, + borehole=b1,k_s=2, k_g=1, R_fp=0.1) + >>> bore_connectivity = [-1, 0] + >>> time = np.array([1.0*10**i for i in range(4, 12)]) + >>> m_flow = 0.25 + >>> cp = 4000. + >>> alpha = 1.0e-6 + >>> gt.gfunction.mixed_inlet_temperature([b1, b2], [Utube1, Utube2], + bore_connectivity, + m_flow, cp, time, alpha) + array([ 0.63783569, 1.63305912, 2.72193357, 4.04093857, 5.98242643, + 7.77218495, 8.66198231, 8.77569636]) + + References + ---------- + .. [#Cimmino2018] Cimmino, M. (2018). g-Functions for bore fields with + mixed parallel and series connections considering the axial fluid + temperature variations. IGSHPA Research Track, Stockholm. In review. + + """ + if disp: + print(60*'-') + print('Calculating g-function for mixed inlet fluid temperatures') + print(60*'-') + # Initialize chrono + tic = tim.time() + # Number of boreholes + nBoreholes = len(boreholes) + # Total number of line sources + nSources = nSegments*nBoreholes + # Number of time values + nt = len(np.atleast_1d(time)) + # Initialize g-function + gFunction = np.zeros_like(np.atleast_1d(time)) + # Initialize segment heat extraction rates + Q = np.zeros((nSources, nt)) + + # If m_flow is supplied as float, apply m_flow to all boreholes + if np.isscalar(m_flow): + m_flow = np.tile(m_flow, nBoreholes) + m_flow_tot = sum([m_flow[i] for i in range(nBoreholes) + if bore_connectivity[i] == -1]) + + # Verify that borehole connectivity is valid + _verify_bore_connectivity(bore_connectivity, nBoreholes) + + # Split boreholes into segments + boreSegments = _borehole_segments(boreholes, nSegments) + # Vector of time values + t = np.atleast_1d(time).flatten() + # Calculate segment to segment thermal response factors + h_ij = thermal_response_factors( + boreSegments, t, alpha, use_similarities=use_similarities, + splitRealAndImage=True, disTol=disTol, tol=tol, processes=processes, + disp=disp) + toc1 = tim.time() + + if disp: + print('Building and solving system of equations ...') + # ------------------------------------------------------------------------- + # g-function. [A] is a coefficient matrix, [X] = [Qb,Tb,Tf_in] is a state + # Build a system of equation [A]*[X] = [B] for the evaluation of the + # space vector of the borehole heat extraction rates, borehole wall + # temperatures and inlet fluid temperature (into the bore field), + # [B] is a coefficient vector. + # ------------------------------------------------------------------------- + + # Segment lengths + Hb = np.array([b.H for b in boreSegments]) + # Vector of time steps + dt = np.hstack((t[0], t[1:] - t[:-1])) + if not np.isscalar(time) and len(time) > 1: + # Spline object for thermal response factors + h_dt = interp1d(np.hstack((0., t)), + np.dstack((np.zeros((nSources,nSources)), h_ij)), + kind=method, axis=2) + # Thermal response factors evaluated at t=dt + h_dt = h_dt(dt) + else: + h_dt = h_ij + # Thermal response factor increments + dh_ij = np.concatenate((h_ij[:,:,0:1], h_ij[:,:,1:]-h_ij[:,:,:-1]), axis=2) + + # Energy balance on borehole segments: + # [Q_{b,i}] = [a_in]*[T_{f,in}] + [a_{b,i}]*[T_{b,i}] + A_eq2 = np.hstack((-np.eye(nSources), np.zeros((nSources, nSources + 1)))) + B_eq2 = np.zeros(nSources) + for i in range(nBoreholes): + # Segment length + Hi = boreholes[i].H / nSegments + # Rows of equation matrix + j1 = i*nSegments + j2 = (i + 1)*nSegments + # Coefficients for current borehole + a_in, a_b = UTubes[i].coefficients_borehole_heat_extraction_rate( + m_flow[i], cp, nSegments) + # [a_b] is the coefficient matrix for [T_{b,i}] + n1 = i*nSegments + nSources + n2 = (i + 1)*nSegments + nSources + A_eq2[j1:j2, n1:n2] = a_b / (-2.0*pi*UTubes[i].k_s*Hi) + + # Assemble matrix coefficient for [T_{f,in}] and all [T_b] + path = _path_to_inlet(bore_connectivity, i) + b_in = a_in + for j in path[1:]: + # Coefficients for borehole j + c_in, c_b = UTubes[j].coefficients_outlet_temperature( + m_flow[j], cp, nSegments) + # Assign the coefficient matrix for [T_{b,j}] + n2 = (j + 1)*nSegments + nSources + n1 = j*nSegments + nSources + A_eq2[j1:j2, n1:n2] = b_in.dot(c_b)/(-2.0*pi*UTubes[i].k_s*Hi) + # Keep on building coefficient for [T_{f,in}] + b_in = b_in.dot(c_in) + A_eq2[j1:j2, -1:] = b_in / (-2.0*pi*UTubes[i].k_s*Hi) + + # Energy conservation: sum([Qb*Hb]) = sum([Hb]) + A_eq3 = np.hstack((Hb, np.zeros(nSources + 1))) + B_eq3 = np.atleast_1d(np.sum(Hb)) + + # Build and solve the system of equations at all times + for p in range(nt): + # Current thermal response factor matrix + h_ij_dt = h_dt[:,:,p] + # Reconstructed load history + Q_reconstructed = load_history_reconstruction(t[0:p+1], Q[:,0:p+1]) + # Borehole wall temperature for zero heat extraction at current step + Tb_0 = _temporal_superposition(dh_ij, Q_reconstructed) + # Spatial superposition: [Tb] = [Tb0] + [h_ij_dt]*[Qb] + A_eq1 = np.hstack((h_ij_dt, + -np.eye(nSources), + np.zeros((nSources, 1)))) + B_eq1 = -Tb_0 + # Assemble equations + B = np.hstack((B_eq1, B_eq2, B_eq3)) + A = np.vstack((A_eq1, A_eq2, A_eq3)) + # Solve the system of equations + X = np.linalg.solve(A, B) + # Store calculated heat extraction rates + Q[:,p] = X[0:nSources] + # The gFunction is equal to the average borehole wall temperature + Tf_in = X[-1] + Tf_out = Tf_in - 2*pi*UTubes[0].k_s*np.sum(Hb)/(m_flow_tot*cp) + Tf = 0.5*(Tf_in + Tf_out) + Rfield = field_thermal_resistance( + UTubes, bore_connectivity, m_flow, cp) + Tb_eff = Tf - 2*pi*UTubes[0].k_s*Rfield + gFunction[p] = Tb_eff + + toc2 = tim.time() + + if disp: + print('{} sec'.format(toc2 - toc1)) + print('Total time for g-function evaluation: {} sec'.format( + toc2 - tic)) + print(60*'-') + + # Return float if time is a scalar + if np.isscalar(time): + gFunction = np.asscalar(gFunction) + + return gFunction + + def load_history_reconstruction(time, Q): """ Reconstructs the load history. diff --git a/pygfunction/pipes.py b/pygfunction/pipes.py index 8b9e6f8b..dffeeeca 100644 --- a/pygfunction/pipes.py +++ b/pygfunction/pipes.py @@ -4,6 +4,8 @@ from scipy.constants import pi from scipy.special import binom +from .boreholes import _path_to_inlet, _verify_bore_connectivity + class _BasePipe(object): """ @@ -1696,6 +1698,105 @@ def thermal_resistances(pos, r_out, r_b, k_s, k_g, Rfp, J=2): return R, Rd +def borehole_thermal_resistance(pipe, m_flow, cp): + """ + Evaluate the effective borehole thermal resistance. + + Parameters + ---------- + pipe : pipe object + Model for pipes inside the borehole. + m_flow : float + Fluid mass flow rate (in kg/s). + cp : float + Fluid specific isobaric heat capacity (in J/kg.K) + + Returns + ------- + Rb : float + Effective borehole thermal resistance (m.K/W). + + """ + # Coefficient for T_{f,out} = a_out*T_{f,in} + [b_out]*[T_b] + a_out = np.asscalar( + pipe.coefficients_outlet_temperature(m_flow, cp, nSegments=1)[0]) + # Coefficient for Q_b = [a_Q]*T{f,in} + [b_Q]*[T_b] + a_Q = np.asscalar(pipe.coefficients_borehole_heat_extraction_rate( + m_flow, cp, nSegments=1)[0]) + # Borehole length + H = pipe.b.H + # Effective borehole thermal resistance + Rb = -0.5*H*(1. + a_out)/a_Q + + return Rb + + +def field_thermal_resistance(pipes, bore_connectivity, m_flow, cp): + """ + Evaluate the effective bore field thermal resistance. + + As proposed in [#Cimmino2018]_. + + Parameters + ---------- + pipes : list of pipe objects + Models for pipes inside each borehole. + bore_connectivity : list + Index of fluid inlet into each borehole. -1 corresponds to a borehole + connected to the bore field inlet. + m_flow : float or array + Fluid mass flow rate in each borehole (in kg/s). + cp : float + Fluid specific isobaric heat capacity (in J/kg.K) + + Returns + ------- + Rfield : float + Effective bore field thermal resistance (m.K/W). + + References + ---------- + .. [#Cimmino2018] Cimmino, M. (2018). g-Functions for bore fields with + mixed parallel and series connections considering the axial fluid + temperature variations. IGSHPA Research Track, Stockholm. In review. + + """ + # Number of boreholes + nBoreholes = len(pipes) + # If m_flow is supplied as float, apply m_flow to all boreholes + if np.isscalar(m_flow): + m_flow = np.tile(m_flow, nBoreholes) + # Total mass flow rate in the bore field + m_flow_tot = sum([m_flow[i] for i in range(nBoreholes) + if bore_connectivity[i] == -1]) + + # Total borehole length + H_tot = sum([pipes[i].b.H for i in range(nBoreholes)]) + + # Verify that borehole connectivity is valid + _verify_bore_connectivity(bore_connectivity, nBoreholes) + + + # Coefficients for T_{f,out} = A_out*T_{f,in} + [B_out]*[T_b], and + # Q_b = [A_Q]*T{f,in} + [B_Q]*[T_b] + A_Q = np.array([np.asscalar( + pipes[i].coefficients_borehole_heat_extraction_rate( + m_flow[i], cp, nSegments=1)[0]) + for i in range(len(pipes))]) + for i in range(len(pipes)): + path = _path_to_inlet(bore_connectivity, i) + for j in path[1:]: + a_out = np.asscalar(pipes[j].coefficients_outlet_temperature( + m_flow[j], cp, nSegments=1)[0]) + A_Q[i] = A_Q[i]*a_out + A_out = 1. + np.sum(A_Q)/(m_flow_tot*cp) + + # Effective bore field thermal resistance + Rfield = -0.5*H_tot*(1. + A_out)/np.sum(A_Q) + + return Rfield + + def fluid_friction_factor_circular_pipe(m_flow, r_in, visc, den, epsilon, tol=1.0e-6): """