Thermal interference between two borefields

[wouterpeere]

Dear @MassimoCimmino

I have two somewhat related questions about what's (currently) possible with pygfunction.

1. Given a certain borefield, you can calculate the g-function of this borefield with pygfunction. However, I am wondering if it is possible to calculate the g-function for a point (x, y) outside of the borefield. Namely, image you have a constant heat extraction out of a certain borefield and you want to plot the temperature evolution over time from a location (x, y) outside of the borefield, can you do this with pygfunction?
2. Image you have two neighbouring houses each with 2 boreholes drilled in their garden. We can calculate the temperature response of each pair of boreholes, but since these two houses are close to each other, there is a thermal interference between both. Is it possible to determine these interference (namely the difference between the thermal response of one pair of boreholes without its neighbour and the thermal response of the boreholes given the neighbouring pair).

Best regards,
@wouterpeere

[mitchute]

Re: 2 - I will point out a paper that we published a few years back that proposes some methods for doing this. @MassimoCimmino may have a more elegant solution, though.

Mitchell_et_al_2020_WGC_Cross_g_functions.pdf

[MassimoCimmino]

No.1

There is nothing currently implemented. However, the heat extraction rates of all segments are saved when using the option 'profiles': True. The weigthed sum of the FLS of all segments onto a location would give the thermal response factor related to the temperature at that location. I have recently received requests for such a feature by email. This is a topic for another issue.

No.2

The method of @mitchute is what I was going to recommend. Thank you @mitchute.

It should be noted that the cross-g-functions calculated this way would be calculated using the 'UHTR' boundary condition. The 'UBWT' condition could be calculated using the heat extraction rates of the segments calculated when solving the system of equations (which are saved with the option 'profiles': True).

The following snippet reproduces fig.2 of the paper with pygfunction and the 'equivalent' solver:

from matplotlib import pyplot as plt
import numpy as np

import pygfunction as gt


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

    # Borehole dimensions
    D = 1.0             # Borehole buried depth (m)
    H = 100.0           # Borehole length (m)
    r_b = 0.05          # Borehole radius (m)
    B = 5.0             # Borehole spacing (m)

    # Thermal properties
    alpha = 1.0e-6      # Ground thermal diffusivity (m2/s)

    # g-Function calculation options
    options = {'nSegments': 1,
               'segment_ratios': None,
               'disp': True,
               'kClusters': 0}

    # Geometrically expanding time vector.
    dt = 100*3600.                  # Time step
    tmax = 3000. * 8760. * 3600.    # Maximum time
    Nt = 25                         # Number of time steps
    ts = H**2/(9.*alpha)            # Bore field characteristic time
    time = gt.utilities.time_geometric(dt, tmax, Nt)
    lntts = np.log(time/ts)

    # -------------------------------------------------------------------------
    # Borehole fields
    # -------------------------------------------------------------------------

    # Borefield A
    N_1 = 2
    N_2 = 1
    field_A = gt.boreholes.rectangle_field(N_1, N_2, B, B, H, D, r_b)
    # Move field 5m north
    for b in field_A: b.y = b.y + B

    # Borefield B
    N_1 = 2
    N_2 = 1
    field_B = gt.boreholes.rectangle_field(N_1, N_2, B, B, H, D, r_b)
    # Slightly change the radius of the boreholes. This tricks pygfunction
    # into classifying the boreholes of fields A and B into different groups
    for b in field_B: b.r_b = b.r_b * 1.0001

    # Borefield C
    N_1 = 2
    N_2 = 2
    field_C = gt.boreholes.rectangle_field(N_1, N_2, B, B, H, D, r_b)

    # -------------------------------------------------------------------------
    # Evaluate g-functions for field C
    # -------------------------------------------------------------------------
    gfunc_C = gt.gfunction.gFunction(
        field_C, alpha, time=time, options=options, method='equivalent', boundary_condition='UHTR')

    # -------------------------------------------------------------------------
    # Evaluate self- and cross- g-functions for fields A and B
    # -------------------------------------------------------------------------
    field_AB = field_A + field_B
    gfunc_AB = gt.gfunction.gFunction(
        field_AB, alpha, options=options, method='equivalent')
    # Number of equivalent boreholes in fields A and B
    nEqBoreholes_A = np.max(gfunc_AB.solver.clusters[:len(field_A)]) + 1
    nEqBoreholes_B = gfunc_AB.solver.nEqBoreholes - nEqBoreholes_A
    # Number of boreholes represented by each equivalent borehole in each group
    nBoreholes_A = np.array([b.nBoreholes for b in gfunc_AB.solver.boreholes[:nEqBoreholes_A]])
    nBoreholes_B = np.array([b.nBoreholes for b in gfunc_AB.solver.boreholes[nEqBoreholes_A:nEqBoreholes_A+nEqBoreholes_B]])
    # Matrix of thermal response factors (eq. borehole to eq. borehole)
    h_AB = gfunc_AB.solver.thermal_response_factors(time, alpha).y[:,:,1:]
    # Weighted sum of the borehole-to-borehole thermal response factors
    g_AA = np.zeros(Nt)
    g_AB = np.zeros(Nt)
    g_BA = np.zeros(Nt)
    g_BB = np.zeros(Nt)
    for k in range(Nt):
        g_AA[k] = np.sum(nBoreholes_A @ h_AB[:nEqBoreholes_A,:nEqBoreholes_A,k]) / np.sum(nBoreholes_A)
        g_AB[k] = np.sum(nBoreholes_A @ h_AB[:nEqBoreholes_A,nEqBoreholes_A:nEqBoreholes_A+nEqBoreholes_B,k]) / np.sum(nBoreholes_A)
        g_BA[k] = np.sum(nBoreholes_B @ h_AB[nEqBoreholes_A:nEqBoreholes_A+nEqBoreholes_B,:nEqBoreholes_A,k]) / np.sum(nBoreholes_B)
        g_BB[k] = np.sum(nBoreholes_B @ h_AB[nEqBoreholes_A:nEqBoreholes_A+nEqBoreholes_B,nEqBoreholes_A:nEqBoreholes_A+nEqBoreholes_B,k]) / np.sum(nBoreholes_B)

    # -------------------------------------------------------------------------
    # Plot g-functions
    # -------------------------------------------------------------------------
    ax = gfunc_C.visualize_g_function().axes[0]
    ax.plot(lntts, g_AA)
    ax.plot(lntts, g_AB)
    ax.plot(lntts, g_AA + g_AB, 'ks')
    ax.legend(['C->C (Ref.)', 'A->A', 'B->A', 'A->A + B->A'])
    plt.tight_layout()

if __name__ == '__main__':
    main()

[wouterpeere]

@mitchute Thank you for this interesting paper. I didn't saw it before, but it is indeed exactly what I was looking for!

@MassimoCimmino : I do not get what you mean by:

It should be noted that the cross-g-functions calculated this way would be calculated using the 'UHTR' boundary condition. The 'UBWT' condition could be calculated using the heat extraction rates of the segments calculated when solving the system of equations (which are saved with the option 'profiles': True).'

Is there a way to use the cross-g-functions concept in combination with the UBWT boundary condition after all?

[MassimoCimmino]

The method assumes, like the 'UHTR' condition, that all boreholes extract heat at the same rate, uniformly along their lengths. When computing g-functions with the 'UBWT' condition, the heat transfer rate is variable along the boreholes, from borehole to borehole, and varies with time. This could be considered in the calculation of the cross-g-functions, since we already evaluate the heat tranfer rates of all segments.

One thing to note is that the impact of the choice of boundary condition is most prominent inside the borehole field. If the two fields are not too close, the choice of boundary condition shouldn't be too significant (intuitively, this remains to be shown). Another thing to keep in mind is that such a 'UBWT' condition is not necessarily more accurate than the 'UHTR' condition for cross-g-functions. They would need to both be compared to a more detailed model for validation.

I am working on the g-function module for #215 and this issue could be incorporated.

[smo999]

Dear @MassimoCimmino ,
thank you very much for the impressive work on the g-functions! I like very much how you compute them.

I wondered whether you had some updates on this question. I would like to simulate two circuits to use some boreholes for heat extraction and other for regeneration (heat injection).

Could you maybe tell me how I could do ?

That would be awesome.

Kind regards,
Sébastien

[wouterpeere]

I just wanted to mention another paper on this topic (Letizia Fascì et al., 2021). It is available here: https://www.sciencedirect.com/science/article/pii/S1359431121006785

[MassimoCimmino]

Thanks @wouterpeere. This is useful.

@smo999 - Somehow I had missed a notification for your comment. This is still far from being implemented but I will look into it after #215. The reason for the delay is that we try to fully integrate to new features to the existing ones in pygfunction. In the case of the present issue, it means that we need to differentiate between boundary conditions of uniform heat extraction rates ('UHTR'), uniform borehole wall temperature ('UBWT'), and mixed inlet fluid conditions ('MIFT'). The first two are at least partially solved in Mitchell et al. (2020) and Fascì et al. (2021). The last one could be obtained by an extension of the solution to #215.

This, like some other "feature" issues, is as much of a research challenge than an implementation challenge and is the reason why we have a fairly long development time.

In the meantime, the solution would depend on the specific geometry of your system. My script above could be use to calculate cross-g-functions based on the work of Mitchell et al. (2020). There is also the possibility of using the geothermal model that should be soon merged into the Modelica Buildings library : lbl-srg/modelica-buildings#3812