g-functions for energy piles

[wouterpeere]

Hi Massimo,

I received a question w.r.t. energy piles and pygfunction. I looked at issue #269 and branch issue44b_shortTermCorrection, but I think there is still some numerical instability there.

Below is the code (which has to be run in issue44b_shortTermCorrection) I used, based on the code you mentioned in the last update on issue #269 and adopted for this specific energy pile with a diameter 120cm.
The results I get are also shown below.

Is this a bug, since I thought the cylindrical correction should have solved this?

Best regards,
Wouter

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

import pygfunction as gt

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

    # Borehole dimensions
    D = 0.8             # Borehole buried depth (m)
    H = 25              # Borehole length (m)
    r_b = 0.6           # 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_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_withCorrection = gt.gfunction.gFunction(
        field, alpha, time=time, options=options_withCorrection,
        method='equivalent')

    # Draw g-function
    ax = gfunc_withCorrection.visualize_g_function().axes[0]
    ax.set_title(f'Field of {len(field)} boreholes')

    plt.tight_layout()

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