Negative g-function values

[tblanke]

Hi,

I have the issue that, under some circumstances, the g-function values are getting negative. Is this a model-based issue which should be captured, or is this issue solvable? Here an example:

import numpy as np
from pygfunction.boreholes import rectangle_field
from pygfunction.gfunction import gFunction
from pygfunction.load_aggregation import ClaessonJaved
def main():
    field = rectangle_field(10,7,2,2,150,2,0.2)
    time = ClaessonJaved(3600, 3600*8760*20).get_times_for_simulation()
    g_func = gFunction(field, 1/5000/1000, time).gFunc
    assert np.all(g_func > 0)
# Main function
if __name__ == '__main__':
    main()

Kind regards

[MassimoCimmino]

This seems to stem from a numerical issue related to your "large" radius r_b=0.2. In #231, we introduced a threshold time for the evaluation of g-functions to avoid numerical errors (linear_threshold = r_b_max**2 / (25 * alpha)). This is not sufficient in your case.

I was able to evaluate your g-function by increasing the threshold:

import numpy as np
from pygfunction.boreholes import rectangle_field
from pygfunction.gfunction import gFunction
from pygfunction.load_aggregation import ClaessonJaved
def main():
    field = rectangle_field(10,7,2,2,150,2,0.2)
    time = ClaessonJaved(3600, 3600*8760*20).get_times_for_simulation()
    options = {'linear_threshold': 24 * 3600.}
    g_func = gFunction(field, 1/5000/1000, time, options=options).gFunc
    assert np.all(g_func > 0)
# Main function
if __name__ == '__main__':
    main()

This issue could be resolved by increasing the default threshold. There is also the possibilty that #44 (which is overdue) fixes the issue for all cases.

[tblanke]

Thanks, that works perfectly.

[wouterpeere]

Hi @MassimoCimmino Hi @tblanke

I have a question when implementing the fix above in https://github.com/wouterpeere/GHEtool/issues/187. As far as I understood it correctly, when it was implemented in #231, the linear threshold does interpolate all the gfunction values for times below it, leaving the other gfunction values untouched.

Now, when I run the code below, I see that all gfunctions are different, also the ones above the threshold. Is this correct behaviour?

from pygfunction.boreholes import rectangle_field
from pygfunction.gfunction import gFunction
from pygfunction.load_aggregation import ClaessonJaved
import matplotlib.pyplot as plt

def main():
    field = rectangle_field(10,7,2,2,150,2,0.2)
    time = ClaessonJaved(3600, 3600*8760*20).get_times_for_simulation()
    g_func_without_linear_threshold = gFunction(field, 1/5000/1000, time).gFunc
    options = {'linear_threshold': 24 * 3600.}
    g_func_with_linear_threshold = gFunction(field, 1 / 5000 / 1000, time, options=options).gFunc

    plt.plot(time, g_func_with_linear_threshold - g_func_without_linear_threshold, label="Difference")
    plt.xlabel('Time in seconds')
    plt.ylabel('Gfunction')
    plt.legend()
    plt.show()
# Main function
if __name__ == '__main__':
    main()

[MassimoCimmino]

In you example above, the g-function without the threshold returns negative values, which causes the differences (e.g. r_b=0.75 shows negligible difference above the threshold). pygfunction is not presently well suited for "geothermal pile" geometries that have a large borehole radius. To remedy the issue, we need to implement the cylindrical correction in #44.

[wouterpeere]

Even with the cylindrical correction (if I take the code from #44), there still appears to be negative values. Was there something wrong with the implementation below?

Note: in contrast with #44, I changed the "y0(p 2)" to "y0(p u)" for I think this was a typo.

from scipy.special import j0, j1, y0, y1, erf, exp1
from scipy.integrate import quad
import numpy as np
from numpy import pi
def _shortTermCorrection(time, gFunc, r_b, alpha):

    def _CHS(u, Fo, p):
        CHS_integrand = 1. / (u ** 2 * pi ** 2) * (np.exp(-u ** 2 * Fo) - 1.0) / (j1(u) ** 2 + y1(u) ** 2) * (
                j0(p * u) * y1(u) - j1(u) * y0(p * u))
        return CHS_integrand

    def _ILS(t, alpha, dis):
        ILS = exp1(dis ** 2 / (4 * alpha * t))
        return ILS

    for i in range(len(time)):
        ILS = _ILS(time[i], alpha, r_b)
        CHS, err = quad(_CHS, 1e-12, 100., args=(alpha * time[i] / r_b ** 2, 1.))

        gFunc[i] = gFunc[i] + 2 * pi * CHS - 0.5 * ILS

    return gFunc

from pygfunction.boreholes import rectangle_field
from pygfunction.gfunction import gFunction
from pygfunction.load_aggregation import ClaessonJaved
import matplotlib.pyplot as plt

def main():
    field = rectangle_field(10,7,2,2,150,2,0.2)
    time = ClaessonJaved(3600, 3600*8760*20).get_times_for_simulation()
    gFunc = gFunction(field, 1/5000/1000, time).gFunc
    gFunc = _shortTermCorrection(time, gFunc, field[0].r_b, 1/5000/1000)
    plt.plot(time, gFunc, label="gFunc")
    plt.xlabel('Time in seconds')
    plt.ylabel('Gfunction')
    plt.legend()
    plt.show()
# Main function
if __name__ == '__main__':
    main()

[MassimoCimmino]

I am sorry. I should have provided a more in depth explanation.

The cause of the negative values is the badly conditioned system of equations due to the near-zero segment-to-segment thermal response factors. When time is small enough or r_b is large enough, heat does not have time to travel from the line source (at the borehole axis) to the borehole wall, and the finite line source solution is close to zero (it may also actually be negative due to numerical error).

For the purpose of solving the bug, the correction cannot be applied after the calculation of the g-function. The diagonal terms of the thermal response factor matrix (the impact of the segments on themselves) need to be corrected before solving the system.

I remember this is why I had postponed issue #44 to be addressed after v2.0 and the refactored g-function interface, since it has non-trivial implications on the gFunction class.

I quickly implemented a fix for the 'equivalent' solver on the branch issue44b_shortTermCorrection (notice the 'b' in the branch name). You can test it out with the following code (adapted from yours):

# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt

import pygfunction as gt

def main():
    # -------------------------------------------------------------------------
    # Simulation parameters
    # -------------------------------------------------------------------------

    # Borehole dimensions
    D = 2.0             # Borehole buried depth (m)
    H = 150.0           # Borehole length (m)
    r_b = 0.2           # Borehole radius (m)
    B = 2.0             # Borehole spacing (m)

    # Field of 10x7 (n=70) boreholes
    N_1 = 10
    N_2 = 7
    field = gt.boreholes.rectangle_field(N_1, N_2, B, B, H, D, r_b)

    # Thermal properties
    alpha = 2e-7        # Ground thermal diffusivity (m2/s)

    options_noCorrection = {
        'cylindrical_correction': False,
        'linear_threshold': 24 * 3600.}
    options_withCorrection = {
        'cylindrical_correction': True,
        'linear_threshold': 0.}

    # Geometrically expanding time vector.
    ts = H**2/(9.*alpha)        # Bore field characteristic time
    # dt = 3600.                  # Time step
    # tmax = ts * np.exp(5)       # Maximum time
    # Nt = 50                     # Number of time steps
    # time = gt.utilities.time_geometric(dt, tmax, Nt)
    time = gt.load_aggregation.ClaessonJaved(3600, 3600*8760*20).get_times_for_simulation()
    lntts = np.log(time/ts)

    # -------------------------------------------------------------------------
    # Evaluate g-functions
    # -------------------------------------------------------------------------
    gfunc_noCorrection = gt.gfunction.gFunction(
        field, alpha, time=time, options=options_noCorrection,
        method='equivalent')
    gfunc_withCorrection = gt.gfunction.gFunction(
        field, alpha, time=time, options=options_withCorrection,
        method='equivalent')

    # Draw g-function
    ax = gfunc_noCorrection.visualize_g_function().axes[0]
    ax.plot(lntts, gfunc_withCorrection.gFunc)
    ax.legend(['Linear threshold',
               'Cylindrical correction',])
    ax.set_title(f'Field of {len(field)} boreholes')
    plt.tight_layout()

    # Draw difference

    # Configure figure and axes
    fig = gt.utilities._initialize_figure()
    ax = fig.add_subplot(111)
    # Axis labels
    ax.set_xlabel(r'ln(t/t_s)')
    ax.set_ylabel(r'difference')
    gt.utilities._format_axes(ax)
    # Plot difference
    ax.plot(lntts, gfunc_withCorrection.gFunc - gfunc_noCorrection.gFunc)
    ax.set_title('Difference between g-functions')
    # Adjust to plot window
    plt.tight_layout()

# Main function
if __name__ == '__main__':
    main()