MassimoCimmino / pygfunction

An open-source toolbox for the evaluation of thermal response factors (g-functions) of geothermal borehole fields.
BSD 3-Clause "New" or "Revised" License
48 stars 21 forks source link

The converging value of g-function is related with the discretization scheme #307

Open metab0t opened 4 weeks ago

metab0t commented 4 weeks ago

I have a script to test the value of g-function when using a different number of segments and discretization schemes (uniform and non-uniform).

import numpy as np
import matplotlib.pyplot as plt

import pygfunction as gt

def define_borefield(N1, N2, Bx, By):
    # Borehole dimensions
    D = 1.0  # Borehole buried depth (m)
    H = 80.0  # Borehole length (m)
    r_b = 0.10  # Borehole radius (m)

    # The field is a retangular array
    boreField = gt.boreholes.rectangle_field(N1, N2, Bx, By, H, D, r_b)

    return boreField

field = define_borefield(2, 1, 6.8, 5.7)

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

dt = 3600 * 24 * 365.0
Nt = 5
ts = [dt * Nt]

N_segs = [16, 100, 500, 1500]
gv_nonuniform = []
gv_uniform = []

for N_seg in N_segs:
    options = {
        "nSegments": N_seg,
        "disp": False,
        "approximate_FLS": True,
    }
    gf = gt.gfunction.gFunction(
        field, alpha, boundary_condition="UBWT", options=options, method="detailed"
    )
    gv = gf.evaluate_g_function(ts)
    gv_nonuniform.append(gv[-1])

    options["segment_ratios"] = np.full(N_seg, 1.0 / N_seg)
    gf = gt.gfunction.gFunction(
        field, alpha, boundary_condition="UBWT", options=options, method="detailed"
    )
    gv = gf.evaluate_g_function(ts)
    gv_uniform.append(gv[-1])

fig, ax = plt.subplots()
ax.plot(N_segs, gv_nonuniform, "o-", label="Non-uniform segments")
ax.plot(N_segs, gv_uniform, "o-", label="Uniform segments")
ax.set_xlabel("Number of segments")
ax.set_ylabel("g-function value")
ax.legend()

plt.show()

The output is: pic

The result is a little surprising because as the number of segments increase, g-function s given by non-uniform segmentation and uniform segmentation converge to different values. I wonder which value is more accurate, and should be regarded as the reference.

MassimoCimmino commented 3 weeks ago

The non-uniform discretization has two parameters : the number of segments and the end-length-ratio (i.e. the relative legnth of the segments at the edges). Both need to be refined when conducting a discretization-independence study. A good way of doing it is to subdivide a base discretization by splitting each segment into two at each discretization level.

You can then use Richardson extrapolation to compare where the two curves converge towards.

import numpy as np
import matplotlib.pyplot as plt

import pygfunction as gt

def define_borefield(N1, N2, Bx, By):
    # Borehole dimensions
    D = 1.0  # Borehole buried depth (m)
    H = 80.0  # Borehole length (m)
    r_b = 0.10  # Borehole radius (m)

    # The field is a retangular array
    boreField = gt.boreholes.rectangle_field(N1, N2, Bx, By, H, D, r_b)

    return boreField

field = define_borefield(2, 1, 6.8, 5.7)

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

dt = 3600 * 24 * 365.0
Nt = 5
ts = [dt * Nt]

N_segs = [2 * 8, 4 * 8, 8 * 8, 16 * 8, 32 * 8, 64 * 8] # Must double at each term
gv_nonuniform = []
gv_uniform = []

segment_ratios_base = gt.utilities.segment_ratios(8, 0.02)

for N_seg in N_segs:
    segment_ratios = np.repeat(segment_ratios_base, N_seg / 8) * 8 / N_seg
    options = {
        "nSegments": N_seg,
        "disp": False,
        "approximate_FLS": True,
        "segment_ratios": segment_ratios
    }
    gf = gt.gfunction.gFunction(
        field, alpha, boundary_condition="UBWT", options=options, method="detailed"
    )
    gv = gf.evaluate_g_function(ts)
    gv_nonuniform.append(gv[-1])

    segment_ratios = None
    options["segment_ratios"] = segment_ratios
    gf = gt.gfunction.gFunction(
        field, alpha, boundary_condition="UBWT", options=options, method="detailed"
    )
    gv = gf.evaluate_g_function(ts)
    gv_uniform.append(gv[-1])

q_nonuniform = np.log((gv_nonuniform[-3] - gv_nonuniform[-2]) / (gv_nonuniform[-2] - gv_nonuniform[-1])) / np.log(2)
g_extrap_nonuniform = gv_nonuniform[-1] + (gv_nonuniform[-1] - gv_nonuniform[-2]) / (2**q_nonuniform - 1)
q_uniform = np.log((gv_uniform[-3] - gv_uniform[-2]) / (gv_uniform[-2] - gv_uniform[-1])) / np.log(2)
g_extrap_uniform = gv_uniform[-1] + (gv_uniform[-1] - gv_uniform[-2]) / (2**q_uniform - 1)

fig, ax = plt.subplots()
ax.plot(N_segs, gv_nonuniform, "o-", label="Non-uniform segments")
ax.plot(N_segs, gv_uniform, "o-", label="Uniform segments")
ax.plot([N_segs[0], N_segs[-1]], [g_extrap_nonuniform, g_extrap_nonuniform], "--", label="Non-uniform extrap")
ax.plot([N_segs[0], N_segs[-1]], [g_extrap_uniform, g_extrap_uniform], "--", label="Uniform extrap")
ax.set_xlabel("Number of segments")
ax.set_ylabel("g-function value")
ax.legend()

plt.show()

Figure_1

metab0t commented 3 weeks ago

Thanks!

The extrapolation indicates that their converging value have slight difference. Does it mean that the "exact" value of g-function depends on which discretization we use?

MassimoCimmino commented 2 weeks ago

The exact value of the g-function should be unique.

Here, I evaluated the extrapolated values using only the last three points in the increasing discretization. We should expect closer values using trend lines on log-scale graphs instead.