We have contributing standards detailed in CONTRIBUTING.md.

Perhaps we should also request that a standard be followed for commit messages?

When attempting to run a whole script which doesn't have an if __name__ == '__main__' call, the pygfunction.gfunction.uniform_temperature function gets stuck in an infinite loop while undergoing the "Identifying similarities" step. This occurs on Windows with Python 2.7, both with a "normal" Python installation (Anaconda2) and with IPython. However, running line-by-line doesn't pose any problems.

If an if __name__ == '__main__' call is indeed used, then the problem disappears. Therefore, the following .py script will fail:

from__future__ import division, print_function, absolute_import
import pygfunction as gt
   
D = 4.
H = 175.
r_b = 0.075
B = 7.5
alpha = 1.0e-6
    
dt = 100*3600.
tmax = 3000. * 8760. * 3600.
Nt = 50
time = gt.utilities.time_geometric(dt, tmax, Nt)
N_1 = 1
N_2 = 1
boreField1 = gt.boreholes.rectangle_field(N_1, N_2, B, B, H, D, r_b)
gfunc = gt.gfunction.uniform_temperature(boreField1, time, alpha, disp=True)

while the following will work properly:

from __future__ import division, print_function, absolute_import
import pygfunction as gt

def main():    
    D = 4.
    H = 175.
    r_b = 0.075
    B = 7.5
    alpha = 1.0e-6
    
    dt = 100*3600.
    tmax = 3000. * 8760. * 3600.
    Nt = 50
    time = gt.utilities.time_geometric(dt, tmax, Nt)
    N_1 = 1
    N_2 = 1
    boreField1 = gt.boreholes.rectangle_field(N_1, N_2, B, B, H, D, r_b)
    gfunc = gt.gfunction.uniform_temperature(boreField1, time, alpha, disp=True)

if __name__ == '__main__':
    main()

This issue concerns the implementation of the infinite line source and cylindrical heat source solutions to correct the short-term response of geothermal bore fields.

Add a function to evaluate g-functions considering the connections between boreholes.

• Add load aggregation algorithms
  • MLAA
  • Liu
  • Claesson and Javed
• Write documentation
  • MLAA
  • Liu
  • Claesson and Javed
• Implement examples
  • docstring examples
  • 1 scheme
  • compare schemes

The input parameter to any load aggregation class method LoadAgg.initialize(g_d) is a 3d array of size (nSources, nSources, Nt).

In cases where nSources = 1, this class method should accept a 1D array of size (Nt,).

Thanks so much for your work. I’m having an issue as my borefield has 13 different circuits (A-M), I’m not sure how to go about creating the bore_connectivity variable as it would not identify the right inlets and outlets.

A workaround I was thinking of using was to create 13 individual bore_connectivity variables and then subsequently create network variables for each circuit. However, I cannot later merge the network variables as they are classes, hence the g-function computation would not be accurate.

How would you work around this problem and have you encountered it before? Thanks a lot for your help.

Right now, pipe models recalculate the coefficients necessary to calculate outputs (temperatures and heat extraction rates) at every call of the class methods. These coefficients should be stored in memory and recalculated only when conditions (e.g. flow rate) change.

Here is the log:

Running Sphinx v1.8.5

Traceback (most recent call last):
  File "/home/docs/checkouts/readthedocs.org/user_builds/pygfunction/envs/latest/lib/python2.7/site-packages/sphinx/cmd/build.py", line 303, in build_main
    args.tags, args.verbosity, args.jobs, args.keep_going)
  File "/home/docs/checkouts/readthedocs.org/user_builds/pygfunction/envs/latest/lib/python2.7/site-packages/sphinx/application.py", line 201, in __init__
    self.config = Config.read(self.confdir, confoverrides or {}, self.tags)
  File "/home/docs/checkouts/readthedocs.org/user_builds/pygfunction/envs/latest/lib/python2.7/site-packages/sphinx/config.py", line 203, in read
    namespace = eval_config_file(filename, tags)
  File "/home/docs/checkouts/readthedocs.org/user_builds/pygfunction/envs/latest/lib/python2.7/site-packages/sphinx/config.py", line 373, in eval_config_file
    raise ConfigError(msg % err)
ConfigError: There is a syntax error in your configuration file: invalid syntax (boreholes.py, line 140)

Configuration error:
There is a syntax error in your configuration file: invalid syntax (boreholes.py, line 140)
Command time: 0s Return: 2

The reference is to:

It may also fail on the next few lines:

The function that computes the effective borehole thermal resistance does not have a citation. Does this come from the literature somewhere?

Add the option to define a bore field from an external text file, having the following structure:

# x   y     H     D     rb
0.    0.    100.  2.5   0.075
5.    0.    100.  2.5   0.075
0.    5.    100.  2.5   0.075
0.    10.   100.  2.5   0.075
0.    20.   100.  2.5   0.075

TODO:

• Implement function boreholes.field_from_file()
• Implement example
• Add documentation

Hello Mr. Cimmino,

I am using the pygfunction module for a simulation of a larger geothermal field for the moment. I like many of the implemented features and I was able to already get decent results when modeling the entire field as a network. In my case however, it turned out that the pipes connecting the heat exchangers to the heat pump can largly influence the fluid temperatures. Especially when heating in summer or cooling in winter. Some of these pipes have a length of up to 60 meters and are buried at about 1m below surface. I found no pipe class for connection pipes in the repository, so my question would be:

Is it possible to include the thermal response of those single pipes connected to the inlet and outlet of the BHE as a series connection into the gfunction?

Is there an approach, you would recommend me, if I needed to add this function myself to the code?

Instead of Example of calculation of g-functions using uniform heat extraction rates, it should indicate using uniform temperatures.

A method for a uniformly distributed field (filled circular or hexagonal, in the same way as the DST model) is missing.

A very raw implementation done by one of my students is shown in the method below

def circularField(N_b,B,H,D,r_b):

  coo = np.array([[0,0]])
  a = 6
  B2 = B
  counter = 0
  shift=0
  for i in range(N_b-1):
      i = i+1
      counter = counter + 1
      if counter > a:
          a = a+6
          B2 = B2 + B
          #shift = shift + 2*np.pi/3.
          counter = 0
      coo = np.append(coo, np.array([[B2*np.sin(2*np.pi*counter/a+shift),B2*np.cos(2*np.pi*counter/a+shift)]]),0)

  field = [gt.boreholes.Borehole(H, D, r_b, x, y) for (x, y) in coo]
    
  return field

However, such implementation seems to result in a singular system of equations.

See this interactive example design script. to reproduce the problem.

Going forward, unit tests will be needed to validate future modifications and implementations.

@MassimoCimmino Would you be interested in adding a Fluid object? Here's something I previously wrote up.

This object takes in a fluid temperature, pressure, mixing fluid and percent mixture. The object directly accesses https://github.com/CoolProp/CoolProp to get the fluid properties.

class Fluid:
    # 4/21/20
    def __init__(self, T: float = 293., P: float = 239220., mixer: str = 'MEG', percent: float = 0):
        """

        :param T: temperature of the water (Kelvin)
        :param P: pressure of the water (Pascal)
        :param mixer: the substance being mixed with water
        :param percent: the percentage of the mixture which is not water

        :example:
        >>> import OSU_ORNL_Tool as szgtl
        >>> T = 20 + 273  # Temp at 20 C for now to match GLHEPRO
        >>> gage_P = 20  # PsiG
        >>> atm_P = 14.69595
        >>> P = (gage_P + atm_P) * 6894.75728  # Pressure in Pa

        >>> # pure water
        >>> fluid = szgtl.ghe.Fluid(T=T, P=P)
        >>> print(fluid.__repr__())

        >>> # 20 % propylene glycol mixed with water
        >>> mix = 'MPG'
        >>> percent = 20
        >>> fluid = szgtl.ghe.Fluid(T=T, P=P, mixer=mix, percent=percent)

        >>> # 60% ethylene glycol mixed with water
        >>> mix = 'MEG'
        >>> percent = 60
        >>> fluid = szgtl.ghe.Fluid(T=T, P=P, mixer=mix, percent=percent)
        >>> print(fluid.__repr__() + '\n')

        >>> # 5% methanol mixed with water water
        >>> mix = 'MMA'
        >>> percent = 5
        >>> fluid = szgtl.ghe.Fluid(T=T, P=P, mixer=mix, percent=percent)
        >>> print(fluid.__repr__() + '\n')

        >>> # ethanol / water
        >>> mix = 'MEA'
        >>> percent = 10
        >>> fluid = szgtl.ghe.Fluid(T=T, P=P, mixer=mix, percent=percent)
        >>> print(fluid.__repr__() + '\n')
        """
        self.fluid_mix = 'INCOMP::' + mixer + '-' + str(percent) + '%'
        if mixer == 'MEG':
            pass
        elif mixer == 'MPG':
            pass
        elif mixer == 'MMA':
            pass
        elif mixer == 'MEA':
            pass
        # add in else if's
        else:
            warnings.warn('It is unknown whether or not cool props has the mixing fluid requested, '
                          'proceed with caution.')
        self.T = T                              # temperature of the fluid {K}
        self.P = P                              # pressure of the fluid in {Pa}
        self.rho = self.density()               # Density, {kg/m^3}
        self.mu = self.dynamic_viscosity()      # Dynamic Viscosity, {Pa s} or {N s/ m^2}
        self.nu = self.kinematic_viscosity()    # Kinematic Viscosity, {m^2/s}
        self.cp = self.mass_heat_capacity()     # Mass specific Cp Specific Heat, {J/kg/K}
        self.rhoCp = self.volumetric_heat_capacity()    # Volumetric heat capacity, {kJ/m3/K}
        self.k = self.thermal_conductivity()    # Thermal conductivity, {W/m/K}
        self.Pr = self.Prandlt_number()         # Prandlt number

    def __repr__(self):
        return str(self.__class__) + '\n' + '\n'.join(
            (str(item) + ' = ' + '{}'.format(self.__dict__[item]) for item in sorted(self.__dict__)))

    def append_to_dict(self, dnary):
        if len(list(dnary.keys())) == 0:
            for item in sorted(self.__dict__):
                dnary[item] = []
        for item in sorted(self.__dict__):
            dnary[item].append('{:.5E}'.format(self.__dict__[item]))

    def density(self):
        return PropsSI('D', 'T', self.T, 'P', self.P, self.fluid_mix)

    def dynamic_viscosity(self):
        return PropsSI('V', 'T', self.T, 'P', self.P, self.fluid_mix)

    def kinematic_viscosity(self):
        return self.mu / self.rho

    def mass_heat_capacity(self):
        return PropsSI('C', 'T', self.T, 'P', self.P, self.fluid_mix)

    def volumetric_heat_capacity(self):
        return self.rho * self.cp / 1000

    def thermal_conductivity(self):
        return PropsSI('L', 'T', self.T, 'P', self.P, self.fluid_mix)

    def Prandlt_number(self):
        return PropsSI('PRANDTL', 'T', self.T, 'P', self.P, self.fluid_mix)

The borehole object does not have an instance for drawing a borehole. Visually representing a borehole upon creation will provide a quick and easy check to make sure the borehole physically possible.

The set_stored_coefficients() class method of the Network class is never called to store matrix coefficients. Like in the pipe classes, this is required to lower the calculation time. Repeated calls using the same mass flow rate should return the previous values instead of recomputing.

This is necesseray for all coefficients_[...] class methods.

Add a function to visualize bore field:

• Implement boreholes.visualize_field()
• Update examples

The function thermal_response_factors() is currently in the gfunction module. It should properly be placed in the heat_transfer module.

With the addition of g-function of non-parallel bore fields in #32, it will be useful to have an interface similar to pipes for the calculations involving networks of boreholes.

A new network module will be developed to introduce a network class that takes as parameters lists of pipes (forming a bore field) and the connectivity map of the network. Class methods similar to pipe classes can then be defined, based on the work presented here.

Following this, gfunction.equal_inlet_temperature and mixed_inlet_temperature can be simplified.

If two boreholes are defined on top of each other and the g-function is calculated, the results are erroneous.

check_duplicates(field) -> duplicate_pairs: finds all boreholes that are less than borehole radius away, places into a list of tuples
remove_duplicates(field, duplicate_pairs) -> new_field: removes any of the duplicates of the input field found from the check_duplicates function

g-function.uniform_temperature() and g-function.equal_inlet_temperature() return an error when called with a single time value (both with np.array() and float). This seems to be due to the use of scipy.interpolate.interp1d().

Example:

from pygfunction.gfunction import uniform_temperature
from pygfunction.boreholes import rectangle_field
H = 150.           # Borehole length [m]
D = 4.             # Borehole buried depth [m]
r_b = 0.075        # Borehole radius [m]
B = 7.5            # Borehole spacing [m]
alpha = 1.0e-6     # Ground thermal diffusivity [m2/s]
N_1 = 3
N_2 = 2
boreField = rectangle_field(N_1, N_2, B, B, H, D, r_b)
time = 33554478*3600.
g = uniform_temperature(boreField, time, alpha)

Provide an object which has access to the dimensions of the HDPE pipe dimensions. There will be an option for one of the two types of pipes: SDR-11 or Schedule-40. There are nominal pipe sizes associated with each of these, the nominal size will be input; the inside and outside diameter will be returned.

Following the pipe selection and the nominal pipe size, the u-tube size will be known. The shank spacing needs to be supplied. Three shank spacing options will be provided by an A, B or C configuration from Paul's (1996) MS Thesis:

• A configuration - the pipes are touching at the center of the borehole
• B configuration - the pipes are evenly spaced in the borehole
• C configuration - the pipes are touching the outside wall of the borehole

Paul, N. D. (1996). The effect of grout thermal conductivity on vertical geothermal heat exchanger design and performance.

The function pipes.convective_heat_transfer_coefficient_circular_pipe() does not consider the case where Re < 2300.

The multipole method would give more accurate evaluations of the internal borehole resistances than the currently implemented line source approximation.

Add in borehole resistance calculation using Claesson and Javed (2019) for double u-tubes.

@MassimoCimmino Do you have any comments on this or how it differs from what you already have in the package? I see that you mention Cimmino (2016) in the MultipleUTube class.

Due date: February 28, 2021

Johan Claesson & Saqib Javed (2019) Explicit multipole formulas and thermal network models for calculating thermal resistances of double U-pipe borehole heat exchangers, Science and Technology for the Built Environment, 25:8, 980-992, DOI: https://doi.org/10.1080/23744731.2019.1620565

Cimmino, M. (2016). Fluid and borehole wall temperature profiles in vertical geothermal boreholes with multiple U-tubes. Renewable Energy, 96, 137–147. https://doi.org/10.1016/j.renene.2016.04.067

While github shows the list of contributors to the project, the documentation only mentions myself by name.

A section should be added to the documentation (and may be in a CONTRIBUTOR.md file) to credit contributors to the project.

It would be nice with an example with multiple boreholes using load aggregation with variable load over time where inlet and outlet temperatures are calculated for the boreholes. All examples using load aggregation currently use a single borehole, and it is not completely clear from the documentation how to proceed when using multiple boreholes. What I am looking for is an example e.g. with the following features:

Initialization:

• Set up a borehole field with e.g. 5 x 5 boreholes (boreholes.rectangle_field)
• Calculate g-functions (gfunction.uniform_temperature)
• Initialize load aggregation (LoadAgg.initialize, e.g. based on ClaessonJaved)
• Calculate pipe thermal resistance and heat transfer coefficients
• Loop over boreholes using e.g. pipes.SingleUTube

Loop over time:

• Set current load for load aggregation for each timestep (LoadAgg.set_current_load).
• Calculate inlet temperature (get_inlet_temperature) for each timestep, and correspondingly for outlet temperature

The uniform borehole wall temperature (UBHWT) documentation points to Cimmino and Bernier 2014. Specifically the uniform_temperature function is an implementation of Boundary Condition 3 (BC-III). Could it be placed in the documentation that it specifically applies to BC-III?

Could the documentation of the uniform_heat_extraction function point the boundary condition to which is applies from Cimmino and Bernier (2014)?

Hello everyone, I'm just getting started looking at this excellent piece of work. It looks promising!

I found some minor issues while playing around with the example scripts. I am using anaconda as the package manager for python libraries, and I installed this library through pip install. However there seems to be a problem with the media.py library, it is missing (see screenshot below).

Trying to solve it, I added manually the missing file to the folder by taking it out directly from the repo, but then I get an error when trying to import the library. This is the error line:

from CoolProp.CoolProp import PropsSI

So I couldn't solve the issue and I can't get the media library up and running.

Add pipe model for coaxial pipes.

Add functionnality to import and export bore fields and pipes, possibly using json.

Bore fields should be exportable/importable in both json and simple text formats.

This could require changes in the interfaces of Borehole and Pipe classes.

The segment response matrix h_ij is of size nq x nq x nt, where nq is the number of sources in the field, and nt is the number of time steps.

Two variants of the segment response matrix h_dt and dh_ij are computed. The factor increment matrix dh_ij is used in calculating the borehole wall temperatures based on all previous responses, equation 18 from Cimmino (2018). The h_dt matrix sets up the linear interpolation for the current thermal response factor matrix h_ij_dt, and can be seen in equation 18 in Cimmino (2018).

With the current methodology here, the segment response matrix h_ij is necessary for computing both h_dt and dh_ij. Therefore, it is not directly possible to get around having to hold 3 (nq x nq x nt) segment response factors in memory. With each of the segment responses being double precision (at 16 bytes).

Would it be possible to directly compute h_ij_dt and dh_ij[:, :, t] from h_ij without computing dh_ij and h_dt without losing speed performance?

Massimo Cimmino (2018) Fast calculation of the g-functions of geothermal borehole fields using similarities in the evaluation of the finite line source solution, Journal of Building Performance Simulation, 11:6, 655-668.

gfunction.uniform_heat_extraction, gfunction.uniform_temperature and gfunction.equal_inlet_temperature include too much duplicate code. Code should be refactored:

• function to evaluate segment-to-segment response factors with or without similarities
• function for load history reconstruction
• use coefficient methods from pipe models in gfunction.equal_inlet_temperature

Implement setup script to facilitate the use of pygfunction.

The g-function is returned as int if the time vector is an array or a list of int.

An if clause should be inserted at the start of the g-function functions to convert time to a float.

Inside the UBHWT g-function calculation I have added a couple of lines to time the total time spent in thermal_response_factors.

   # Calculate segment to segment thermal response factors
    toca = tim.time()
    h_ij = thermal_response_factors(
        boreSegments, t, alpha, use_similarities=use_similarities,
        splitRealAndImage=True, disTol=disTol, tol=tol, processes=processes,
        disp=disp)
    tocb = tim.time()
    print('Total time spent in thermal response: {}'.format(tocb - toca))

The output here for a 5x5 field with regular spacing and 30 number of time steps is Total time spent in thermal response: 43.6.

Additionally I have slightly modified lines in pygfunction/heat_transfer.py for the timing of the calculation of h_ij to occur after the pool.close().

    # Close pool of workers
    pool.close()
    pool.join()

    if disp:
        print('{} sec'.format(toc2 - tic))

The output here is 10.633503675460815 sec.

So what exactly is going on here? Why does a function that completes the hefty portion of the calculation in 10 seconds taking 43 seconds to return. It is because Python is having to copy every element back into the scope of gfunction.py.

To fix this issue, h_ij must be declared in gfunction.py and passed into thermal_response_factors(...). The (very large) cube will be passed by reference, returning from the function will be seamless. (h_ij could be an optional argument, if h_ij == None is true, then go about it this way else h_ij pass by ref.

How to solve 2 extraction wells with a remote reinjection well of about 100m each with g-function

• Refactor pipe models
  • _BasePipe
  • SingleUTube
  • MultipleUTube
  • IndependentMultipleUTube
• Example for fluid temperature profiles (use to test implementation)
  • SingleUTube
  • MultipleUTube
  • IndependentMultipleUTube
• Add class method(s) to calculate inlet temperatures from Q and Tb
  • _BasePipe
• Examples should be implemented for their use in simulations.
  • SingleUTube
• The documentation for pipe model classes is not in line with the other modules.
• The equal_inlet_temperature g-function algorithm should be modified according to the implemented changes.
• Example calculation of g-functions with equal_inlet_temperature (compare with other boundary conditions)

Hi Massimo,

Correct me if I missed it, but I could not find any means to extract the kappa vector with the dimensional thermal response factors in pygfunction.
I did a small tweak on my branch of pygfunction to extract the kappa vector. In the load_aggregation.py file, within the class ClaessonJaved I added the following function

    def get_kappa(self):

        return self.dg

It is a small addition, but perhaps can be useful in the code. If you are interested I can add it to the classes MLAA and Liu and make a pull request.

Best,
Iago

With the multiplication of borehole boundary conditions, a centralize function would simplify the use of pygfunction.

Increment files to v1.1.2 and upload to Pypi.

Implement a utilities module that contains functions to generate time vectors with commonly encountered grid spacings. This way the g-function can be subsambled and interpolated at the required times.

• Claesson and Javed (2012) scheme
• Marcotte and Pasquier (2008) scheme
• Geometric expansion
• Update examples
• Documentation

Does the most current version account for inclined boreholes?

Given that the segment-to-segment response factor is calculated by equation 8 (Cimmino and Bernier 2014), how would the distance formula d_{ij} be implemented?

For example, currently calculation of d_{ij} is given by the equation 11.

Would equation 11 need to become d_{uv}? Given the implementation of inclined boreholes, could it now be that x_i == x_u is false, and x_j == x_v is false?

Additionally, the location of segment u, (x_u, y_u) could vary slightly along the segment. Would the method to calculate the distance d_{uv} result in the distance from top-to-top, center-to-center, bottom-to-bottom or some other approach?

I'm assuming this implementation would occur in _borehole_segments. There has been extensive validation of the analytical implementation of the uniform borehole wall temperature (UBWHT) boundary condition when compared to Eskilson's numerical method for vertical boreholes. Given implementation, inclined boreholes would need compared to some of Eskilson's inclined borehole g-functions.

Cimmino, Massimo, and Michel Bernier. "A Semi-analytical Method to Generate G-functions for Geothermal Bore Fields." International Journal of Heat and Mass Transfer 70 (2014): 641-50. Web.

The g-function methods use similarities to accelerate the calculation of thermal response factors (i.e. the finite line source solution) by applying the same values of thermal response factors to borehole pairs (of equal dimensions) that share similar distances, under some tolerance. At the moment the tolerance on radial distances is absolute with a default value of 10 cm. This means that distances are deemed equal if their values are within 10 cm of each other, no matter if the distance is 10 m or 100 m.

A relative tolerance (e.g. 1 %) may yield results that are as accurate much faster.

This issue is to replace the absolute tolerance on distances by a relative tolerance and validate that there is no loss in accuracy.

Consider using quadpy or other specialized integration packages to speed up the calculation of integrals. quadpy supports vectorization for the integral bounds and has many options for integration schemes (some of which might be better adapted to the evaluation of the FLS).

I just checked and scipy seems to have support vectorization of the integrand. It could be used to vectorize over distances, whereas quadpy can vectorize over time.

This issue is to implement the calculation of g-functions of bore fields with uniform and equal fluid temperature along the borehole lengths.

This boundary condition was proposed by Monzó et al., here, and requires fewer parameters than the equal inlet temperature condition.

The documentation does not seem to be clear with respect to using mixed series/parallel-coupled boreholes with load aggregation, i.e. combining the example "mixed_inlet_conditions.py" with the example "fluid_temperature_multiple_boreholes.py". Is this possible at all, or are there restrictions here? I could not find such restrictions in the documentation.

Example: As there are N inlets (i.e. N series of coupled boreholes), I would have suspected that one should define the load-aggregation with nSources=N. That works fine in itself, but then Loadagg.initialize expects a 54xN array (at least according to section 2.4 in the documentation; 54 is here due to the number of timesteps), while the gfunction for the network (via gt.gfunction.mixed_inlet_temperature) gives an 54x1 array as output instead of a 54xN array. I've also tried using an array input of mflow to gt.gfunction.mixed_inlet_temperature, but that gives an error message and it looks like the function is not prepared to output a 54xN array.

If one instead uses nSources=1, one gets only one deltaT_b (ground temperature), which can not be correct? It is also not clear how one should calculate the ground temperature in each borehole with the functions as they are now (e.g. how to evaluate the correct input to Utubes[i].get_temperature()).

pygfunction makes heavy use of numpy and scipy. Both project dropped support for Python 2.7.

The development of pygfunction is done using Python3 with recent versions of numpy and scipy. It makes sense to also drop support for Python 2.7.