g-Function becomes negative for coaxial pipe depending on fluid

[OskarRaftegard]

I've was playing around with a coaxial pipe flows/dimensions, when the g-function suddenly became negative.

G-function is positive for 20 degC water (I think it is a laminar flow in the inner pipe as well as the outer annulus) It is negative for 60 degC water (I think it is laminar flow in the inner pipe and turbulent in the outer annulus)

I based my code on the examples in the documentation, but I could of course have made mistakes.

Below is my code:

# -*- coding: utf-8 -*-

import numpy as np
import pygfunction as gt
from scipy.constants import pi
import matplotlib.pyplot as plt
from matplotlib.ticker import AutoMinorLocator

# -------------------------------------------------------------------------
#  EXAMPLE OF SETTINGS - 
# -------------------------------------------------------------------------

B = 4               # Distance between boreholes (m)
D = 4.0             # Borehole buried depth (m)
r_b = 0.1143 / 2    # Borehole radius (m)
H = 200.            # Borehole length (m)

# ---- DETAILED SETTINGS -----

# g-Function calculation options
nSegments = 8
method = "similarities"
boundary_condition= "MIFT"

Jmultipoles = 4 # Used by pipes module

# Simulation parameters
dt = 3600.                  # Time step (s)
tmax = 100.*8760. * 3600.    # Maximum time (s)

# Load aggregation scheme
cells_per_level = 10 # 5 is defualt value
nSources=1

# Borefield geometry
N_1 = 1           # Number of boreholes in the x-direction (columns)
N_2 = 1             # Number of boreholes in the y-direction (rows)

# Ground properties
T_g = 10.0          # Undisturbed ground temperature (degC)
k_s = 3.1           # Ground thermal conductivity (W/m.K)
Cp_vol = 2160000    # Volumetric heat capacity (J/m3.K)
alpha = k_s/Cp_vol      # Ground thermal diffusivity (m2/s)

# Pipe dimensions and properties:
# pos, r_in, r_out, borehole, k_s, k_g, R_ff, R_fp, J

# Pipe properties
epsilon = 1.5e-6    # Pipe roughness (m)
k_p = 0.42           # Pipe thermal conductivity (W/m.K)

# Pipe positions
# Coaxial pipe (x, y)
pos = (0., 0.)

# Pipe dimensions (coaxial)
# Not: Outer pipe is 0.5 mm smaller than borehole and 0.001 mm thick
r_in_in = 90/2/1000    # Inside pipe inner radius (m)
r_in_out = r_in_in + 6 / 1000    # Inside pipe outer radius (m)
r_out_out = r_b - 0.5 / 1000   # Outer pipe outside radius (m)
r_out_in = r_out_out - 0.001 /1000    # Outer pipe inside radius (m)

# Grout properties
k_g = 1.5           # Grout thermal conductivity (W/m.K)

# Fluid properties
m_flow_borehole = 0.14      # Total fluid mass flow rate per borehole (kg/s)
media = "Water"
Tmedia = 60. # Media average tmperature (degC)

nBoreholes = N_1 * N_2

# g-function 
segment_ratios=gt.utilities.segment_ratios(nSegments, end_length_ratio=0.02)
options = {'nSegments': nSegments,
           'segment_ratios' : segment_ratios,
           'disp': True}

# Fluid properties
m_flow_network = m_flow_borehole*N_1*N_2    # Total fluid mass flow rate (kg/s)

fluid = gt.media.Fluid(media, 0, Tmedia)
cp_f = fluid.cp     # Fluid specific isobaric heat capacity (J/kg.K)
rho_f = fluid.rho   # Fluid density (kg/m3)
mu_f = fluid.mu     # Fluid dynamic viscosity (kg/m.s)
k_f = fluid.k       # Fluid thermal conductivity (W/m.K)

# -------------------------------------------------------------------------
# Define a rectangular borehole field
# -------------------------------------------------------------------------

# Rectangular field
boreField = gt.boreholes.rectangle_field(N_1, N_2, B, B, H, D, r_b)
nBoreholes=len(boreField)

# -----------
# Calculate Coaxial BHE  (based on  example in PyG-function):
# ---------------

# Vectors of inner and outer pipe radii
# Note : The dimensions of the inlet pipe are the first elements of
#        the vectors. In this example, the inlet pipe is the inside pipe.
r_inner = np.array([r_in_in, r_out_in])     # Inner pipe radii (m)
r_outer = np.array([r_in_out, r_out_out])   # Outer pip radii (m)

# Pipe thermal resistance
# (the two pipes have the same thermal conductivity, k_p)
# Inner pipe
R_p_in = gt.pipes.conduction_thermal_resistance_circular_pipe(
    r_in_in, r_in_out, k_p)
# Outer pipe
R_p_out = gt.pipes.conduction_thermal_resistance_circular_pipe(
    r_out_in, r_out_out, k_p)

# Fluid-to-fluid thermal resistance
# Inner pipe
h_f_in = gt.pipes.convective_heat_transfer_coefficient_circular_pipe(
    m_flow_borehole, r_in_in, mu_f, rho_f, k_f, cp_f, epsilon)
R_f_in = 1.0 / (h_f_in * 2 * pi * r_in_in)
# Outer pipe
h_f_a_in, h_f_a_out = \
    gt.pipes.convective_heat_transfer_coefficient_concentric_annulus(
        m_flow_borehole, r_in_out, r_out_in, mu_f, rho_f, k_f, cp_f,
        epsilon)

R_f_out_in = 1.0 / (h_f_a_in * 2 * pi * r_in_out)
R_ff = R_f_in + R_p_in + R_f_out_in

# Coaxial GHE in borehole
R_f_out_out = 1.0 / (h_f_a_out * 2 * pi * r_out_in)
R_fp = R_p_out + R_f_out_out

# All borehole collectors in the borefield (Coaxial borehole collector)
CoaxBHCs = []
for borehole in boreField:
    CoaxBHC = gt.pipes.Coaxial(pos, 
                               r_inner, r_outer, 
                               borehole, 
                               k_s, k_g, 
                               R_ff, R_fp, 
                               J=Jmultipoles)
    CoaxBHCs.append(CoaxBHC)

# Evaluate and print the effective borehole thermal resistance for the first borehole in borefield
R_b = CoaxBHCs[0].effective_borehole_thermal_resistance(
    m_flow_borehole, fluid.cp)

# Visualize the borehole geometry (and save the figure)
CoaxBHCs[0].visualize_pipes()

# Build a network object from the list of UTubes
network = gt.networks.Network(boreField,CoaxBHCs,
                              m_flow_network=m_flow_network, 
                              cp_f=cp_f,
                              nSegments=nSegments,
                              segment_ratios=segment_ratios)

# Simulation parameters
Nt = int(np.ceil(tmax/dt))  # Number of time steps
time = dt * np.arange(1, Nt+1)

# Load aggregation scheme
LoadAgg = gt.load_aggregation.ClaessonJaved(dt, tmax, 
                                            nSources=nSources,
                                            cells_per_level=cells_per_level)

# Get time values needed for g-function evaluation
time_req = LoadAgg.get_times_for_simulation()

# -------------------------------------------------------------------------
# Calculate g-function
# -------------------------------------------------------------------------

# Calculate g-function
gFunc = gt.gfunction.gFunction(network, alpha,time=time_req, 
                               method=method, boundary_condition=boundary_condition,
                               options=options)

gFunc.visualize_g_function()

print(f'Coaxial Borehole thermal resistance: {R_b:.4f} m.K/W')

[MassimoCimmino]

Thank you @OskarRaftegard for reporting this issue.

The issue seems to be related to the pipe model. I tested the evaluation of the fluid temperature profiles for the two temperatures, continued from your code. There is a non-physical increase of the fluid temperature at the bottom of the borehole when Tmedia=60 and the flow is turbulent. In that case, both the fluid to fluid (inner to outer) and the fluid to ground (outer to wall) are small and of the same order of magnitude (a few 0.01 m-K/W).

See the code and figures below.

I will delve deeper into the coaxial pipe model.

# -------------------------------------------------------------------------
# Plot fluid temperature profiles
# -------------------------------------------------------------------------

# Evaluate temperatures at nz evenly spaced depths along the borehole
# at the (it+1)-th time step
nz = 20
z = np.linspace(0., H, num=nz)
T_b = 0.
T_f_in = 1.

T_f = CoaxBHCs[0].get_temperature(
    z, T_f_in, T_b, m_flow_borehole, cp_f)

# Configure figure and axes
fig = gt.utilities._initialize_figure()

ax = fig.add_subplot(111)
# Axis labels
ax.set_xlabel(r'Temperature [degC]')
ax.set_ylabel(r'Depth from borehole head [m]')
ax.set_title(f'Tmedia = {Tmedia:.0f} [degC]')
gt.utilities._format_axes(ax)

# Plot temperatures
ax.plot(np.array([T_b, T_b]), np.array([0., H]), 'k--')
ax.plot(T_f, z, 'b-')
ax.legend(['Borehole wall', 'Fluid'])

# Reverse y-axes
ax.set_ylim(ax.get_ylim()[::-1])
# Adjust to plot window
plt.tight_layout()

[MassimoCimmino]

Here is a higher resolution figure for Tmedia=60. There appears to be numerical noise, most likely due to a round-off error when subtracting large values calculated by pipe._general_solution() (which involves calls to exp(beta*z), sinh(beta*z) and cosh(beta*z)).

[OskarRaftegard]

Sounds reasonable with a numeric error. I have continued to play around, but with the double-U collector. I get negative temperatures with the double-U as well, but it do not look like the same kind of numerical error. T_f slowly drops below T_b with a few decimals. The fig is plotted with nz=1000. More or less same "settings" as above.

# Pipe dimensions and properties
r_out = 0.016      # Pipe outer radius (m)
r_in = r_out - 0.003       # Pipe inner radius (m)
D_s = 0.074 / 2         # Shank spacing radius (m)
epsilon = 1.5e-6    # Pipe roughness (m)
k_p = 0.42           # Pipe thermal conductivity (W/m.K)

# Pipe positions
# Double U-tube [(x_in1, y_in1), (x_in2, y_in2),
#                (x_out1, y_out1), (x_out2, y_out2)]
pos = [(-D_s, 0.), (0., -D_s), (D_s, 0.), (0., D_s)]

# Grout properties
k_g = 1.5           # Grout thermal conductivity (W/m.K)

# Fluid properties
m_flow_borehole = 0.14      # Total fluid mass flow rate per borehole (kg/s)
media = "Water"
Tmedia = 60. # Media average tmperature (degC) 

[MassimoCimmino]

I seem to a found a promising approximation for when gamma * z is large in pipes.SingleUTube.coefficients_temperature.

For now, it is only done for a_in:

# -------------------------------------------------------------------------
# Plot fluid temperature profiles
# -------------------------------------------------------------------------

# Evaluate temperatures at nz evenly spaced depths along the borehole
# at the (it+1)-th time step
nz = 200
z = np.linspace(0., H, num=nz)
T_b = 0.
T_f_in = 1.

T_f = CoaxBHCs[0].get_temperature(
    z, T_f_in, T_b, m_flow_borehole, cp_f)
a_in, a_b = CoaxBHCs[0].coefficients_temperature(
    z, m_flow_borehole, cp_f, nSegments)

delta = CoaxBHCs[0]._delta
gamma = CoaxBHCs[0]._gamma
beta1 = CoaxBHCs[0]._beta1
beta2 = CoaxBHCs[0]._beta2
beta = CoaxBHCs[0]._beta
beta12 = CoaxBHCs[0]._beta12
A = (1 - 0.5 * (beta1 + beta2) / gamma * np.tanh(gamma*H)) / (1 + 0.5 * (beta1 + beta2) / gamma * np.tanh(gamma*H))
z_crit = -0.5 / gamma * np.log(1e-8)
a_in[z > z_crit, 0, 0] = 0.5 * beta12 / gamma * np.exp((beta + gamma) * z[z > z_crit] - 2 * gamma * H) + 0.5 * (1 - A * beta12/gamma + delta) * np.exp((beta - gamma) * z[z > z_crit])
a_in[z > z_crit, 1, 0] = 0.5 * (1 + delta) * np.exp((beta + gamma) * z[z > z_crit] - 2 * gamma * H) + 0.5 * (A * (1 - delta) + beta12/gamma) * np.exp((beta - gamma) * z[z > z_crit])
T_f2 = a_in @ np.atleast_1d(T_f_in) + a_b @ np.full(nSegments, T_b)

# Configure figure and axes
fig = gt.utilities._initialize_figure()

ax = fig.add_subplot(111)
# Axis labels
ax.set_xlabel(r'Temperature [degC]')
ax.set_ylabel(r'Depth from borehole head [m]')
ax.set_title(f'Tmedia = {Tmedia:.0f} [degC]')
gt.utilities._format_axes(ax)

# Plot temperatures
ax.plot(T_f, z, 'b-')
ax.plot(T_f2, z, 'r-')
ax.plot(np.array([T_b, T_b]), np.array([0., H]), 'k--')
ax.plot(np.array([0., 1.]), np.array([z_crit, z_crit]), 'k:')

# Reverse y-axes
ax.set_ylim(ax.get_ylim()[::-1])
# Adjust to plot window
plt.tight_layout()

In the figure below, the blue line is the original method, red is the approximation which only applied below the dotted line, and the dashed line is the borehole wall.

It is still needed for the coefficient a_b, for example if we set T_f_in = 0. and T_b = 1:

[MassimoCimmino]

To do:

• [x] Approximation of a_in
• [x] Approximation of a_b
• [x] Validation (using a finite difference method?)

[MassimoCimmino]

The new model seems to work well for the fluid temperature profile (for the case in the top comment).

Here is a test script for uniform and non-uniform borehole wall temperatures:

# -*- coding: utf-8 -*-
import numpy as np
import pygfunction as gt
from scipy.constants import pi
from scipy.interpolate import interp1d
import matplotlib.pyplot as plt
from matplotlib.ticker import AutoMinorLocator

# -------------------------------------------------------------------------
#  EXAMPLE OF SETTINGS - 
# -------------------------------------------------------------------------

D = 4.0             # Borehole buried depth (m)
r_b = 0.1143 / 2    # Borehole radius (m)
H = 200.            # Borehole length (m)

# Ground properties
k_s = 3.1           # Ground thermal conductivity (W/m.K)

# Pipe dimensions and properties:
# pos, r_in, r_out, borehole, k_s, k_g, R_ff, R_fp, J

# Pipe properties
epsilon = 1.5e-6    # Pipe roughness (m)
k_p = 0.42           # Pipe thermal conductivity (W/m.K)

# Pipe positions
# Coaxial pipe (x, y)
pos = (0., 0.)

# Pipe dimensions (coaxial)
# Not: Outer pipe is 0.5 mm smaller than borehole and 0.001 mm thick
r_in_in = 90/2/1000    # Inside pipe inner radius (m)
r_in_out = r_in_in + 6 / 1000    # Inside pipe outer radius (m)
r_out_out = r_b - 0.5 / 1000   # Outer pipe outside radius (m)
r_out_in = r_out_out - 0.001 /1000    # Outer pipe inside radius (m)

# Grout properties
k_g = 1.5           # Grout thermal conductivity (W/m.K)

# Fluid properties
m_flow_borehole = 0.14      # Total fluid mass flow rate per borehole (kg/s)
media = "Water"
Tmedia = 60. # Media average tmperature (degC)

fluid = gt.media.Fluid(media, 0, Tmedia)
cp_f = fluid.cp     # Fluid specific isobaric heat capacity (J/kg.K)
rho_f = fluid.rho   # Fluid density (kg/m3)
mu_f = fluid.mu     # Fluid dynamic viscosity (kg/m.s)
k_f = fluid.k       # Fluid thermal conductivity (W/m.K)

# -----------
# Calculate Coaxial BHE  (based on  example in PyG-function):
# ---------------

# Vectors of inner and outer pipe radii
# Note : The dimensions of the inlet pipe are the first elements of
#        the vectors. In this example, the inlet pipe is the inside pipe.
r_inner = np.array([r_in_in, r_out_in])     # Inner pipe radii (m)
r_outer = np.array([r_in_out, r_out_out])   # Outer pip radii (m)

# Pipe thermal resistance
# (the two pipes have the same thermal conductivity, k_p)
# Inner pipe
R_p_in = gt.pipes.conduction_thermal_resistance_circular_pipe(
    r_in_in, r_in_out, k_p)
# Outer pipe
R_p_out = gt.pipes.conduction_thermal_resistance_circular_pipe(
    r_out_in, r_out_out, k_p)

# Fluid-to-fluid thermal resistance
# Inner pipe
h_f_in = gt.pipes.convective_heat_transfer_coefficient_circular_pipe(
    m_flow_borehole, r_in_in, mu_f, rho_f, k_f, cp_f, epsilon)
R_f_in = 1.0 / (h_f_in * 2 * pi * r_in_in)
# Outer pipe
h_f_a_in, h_f_a_out = \
    gt.pipes.convective_heat_transfer_coefficient_concentric_annulus(
        m_flow_borehole, r_in_out, r_out_in, mu_f, rho_f, k_f, cp_f,
        epsilon)

R_f_out_in = 1.0 / (h_f_a_in * 2 * pi * r_in_out)
R_ff = R_f_in + R_p_in + R_f_out_in

# Coaxial GHE in borehole
R_f_out_out = 1.0 / (h_f_a_out * 2 * pi * r_out_in)
R_fp = R_p_out + R_f_out_out

# Coaxial borehole collector
borehole = gt.boreholes.Borehole(H, D, r_b, 0., 0.)
coaxial = gt.pipes.Coaxial(
    pos, r_inner, r_outer, borehole, k_s, k_g, R_ff, R_fp)

# -------------------------------------------------------------------------
# Plot fluid temperature profiles
# -------------------------------------------------------------------------

# Thermal parameters
coaxial._update_model_variables(m_flow_borehole, cp_f, 1, None)
delta = coaxial._delta
gamma = coaxial._gamma
beta1 = coaxial._beta1
beta2 = coaxial._beta2
beta = coaxial._beta
beta12 = coaxial._beta12

# Parameters
nz = 400
nSegments = 1
# segment_ratios = None
segment_ratios = gt.utilities.segment_ratios(nSegments)
z = np.linspace(0., H, num=nz)
z_u = coaxial.b._segment_edges(nSegments, segment_ratios=segment_ratios)
T_b = -2. - np.random.rand(nSegments)
T_f_in = 1.
tol = 1e-8
z_crit = -0.5 / gamma * np.log(tol)

# New model
# ---------
A = 1 - delta + beta12 / gamma
B = 1 + delta - beta12 / gamma
C = (A + B * np.exp(-2 * gamma * H)) / (B + A * np.exp(-2 * gamma * H))
a_in = np.zeros((nz, 2, 1))
a_in[:,0,0] = (
    (0.5 * (A * beta12 / gamma + B * (1 - delta)) / (B + A * np.exp(-2 * gamma * H)) * np.exp((beta + gamma) * z) if H < z_crit else 0.)
    + (0.5 * (B * beta12 / gamma + A * (1 - delta)) / (B + A * np.exp(-2 * gamma * H)) * np.exp((beta + gamma) * z - 2 * gamma * H))
    + 0.5 * (1 + delta - beta12 / gamma * C) * np.exp((beta - gamma) * z)
    )
a_in[:,1,0] = (
    (0.5 * (A * (1 + delta) - B * beta12 / gamma) / (B + A * np.exp(-2 * gamma * H)) * np.exp((beta + gamma) * z) if H < z_crit else 0.)
    + (0.5 * (B * (1 + delta) - A * beta12 / gamma) / (B + A * np.exp(-2 * gamma * H)) * np.exp((beta + gamma) * z - 2 * gamma * H))
    + 0.5 * (beta12 / gamma + (1 - delta) * C) * np.exp((beta - gamma) * z)
    )

A4 = (beta1 - delta * beta1 - beta2 * beta12 / gamma) / (beta + gamma)
B4 = (beta1 + delta * beta1 + beta2 * beta12 / gamma) / (beta - gamma)
A5 = (beta2 + delta * beta2 + beta1 * beta12 / gamma) / (beta + gamma)
B5 = (beta2 - delta * beta2 - beta1 * beta12 / gamma) / (beta - gamma)
D = 1 / (B + A * np.exp(-2 * gamma * H))
a_b = np.zeros((nz, 2, nSegments + 1))
dz = np.subtract.outer(z, z_u)
index = dz >= 0.
dz_p = dz[index]
z_p = np.tile(z[:, np.newaxis], (1, nSegments+1))[index]
z_u_p = np.tile(z_u[np.newaxis, :], (nz, 1))[index]
dz_m = dz[~index]
z_m = np.tile(z[:, np.newaxis], (1, nSegments+1))[~index]
z_u_m = np.tile(z_u[np.newaxis, :], (nz, 1))[~index]
a_b[:,0,:][index] = 0.5 * (
    ((beta12 / gamma * D * (A4 + A5) + A4) * np.exp((beta + gamma) * dz_p) if H < z_crit else A4 * A / (B + A * np.exp(-2 * gamma * H)) * np.exp((beta + gamma) * dz_p - 2 * gamma * H))
    + B4 * np.exp((beta - gamma) * dz_p)
    - beta12 / gamma * D * (A4 + A5) * np.exp((beta - gamma) * z_p - (beta + gamma) * z_u_p)
    + beta12 / gamma * D * (B4 + B5) * np.exp((beta + gamma) * z_p - (beta - gamma) * z_u_p - 2 * gamma * H)
    - beta12 / gamma * D * (B4 + B5) * np.exp((beta - gamma) * z_p - (beta - gamma) * z_u_p - 2 * gamma * H)
    )
a_b[:,0,:][~index] = 0.5 * (
    (beta12 / gamma * D * (A4 + A5) * np.exp((beta + gamma) * dz_m) + A4)
    + B4
    - beta12 / gamma * D * (A4 + A5) * np.exp((beta - gamma) * z_m - (beta + gamma) * z_u_m)
    + beta12 / gamma * D * (B4 + B5) * np.exp((beta + gamma) * z_m - (beta - gamma) * z_u_m - 2 * gamma * H)
    - beta12 / gamma * D * (B4 + B5) * np.exp((beta - gamma) * z_m - (beta - gamma) * z_u_m - 2 * gamma * H)
    )
a_b[:,1,:][index] = 0.5 * (
    (((1 + delta) * D * (A4 + A5) - A5) * np.exp((beta + gamma) * dz_p) if H < z_crit else -A5 * A / (B + A * np.exp(-2 * gamma * H)) * np.exp((beta + gamma) * dz_p - 2 * gamma * H))
    - B5 * np.exp((beta - gamma) * dz_p)
    + (1 - delta) * D * (A4 + A5) * np.exp((beta - gamma) * z_p - (beta + gamma) * z_u_p)
    + (1 + delta) * D * (B4 + B5) * np.exp((beta + gamma) * z_p - (beta - gamma) * z_u_p - 2 * gamma * H)
    + (1 - delta) * D * (B4 + B5) * np.exp((beta - gamma) * z_p - (beta - gamma) * z_u_p - 2 * gamma * H)
    )
a_b[:,1,:][~index] = 0.5 * (
    ((1 + delta) * D * (A4 + A5) * np.exp((beta + gamma) * dz_m) - A5)
    - B5
    + (1 - delta) * D * (A4 + A5) * np.exp((beta - gamma) * z_m - (beta + gamma) * z_u_m)
    + (1 + delta) * D * (B4 + B5) * np.exp((beta + gamma) * z_m - (beta - gamma) * z_u_m - 2 * gamma * H)
    + (1 - delta) * D * (B4 + B5) * np.exp((beta - gamma) * z_m - (beta - gamma) * z_u_m - 2 * gamma * H)
    )

a_b = -np.diff(a_b, axis=2)

T_f_new = a_in @ np.atleast_1d(T_f_in) + a_b @ T_b

# Initial model
# -------------
T_f_initial = coaxial.get_temperature(
    z, T_f_in, T_b, m_flow_borehole, cp_f, segment_ratios=segment_ratios)

# Finite diference model
# ----------------------
T_b_func = interp1d(z_u, np.concatenate((T_b, [T_b[-1]])), kind='zero')
A_fd = np.zeros((2*nz, 2*nz))
B_fd = np.zeros(2*nz)
for i in range(1, nz):
    # Downward pipe
    A_fd[i,i] = 1. + 0.5 * (z[i] - z[i-1]) * (beta1 + beta12)
    A_fd[i,i-1] = -1. + 0.5 * (z[i] - z[i-1]) * (beta1 + beta12)
    A_fd[i,-i] = -0.5 * (z[i] - z[i-1]) * beta12
    A_fd[i,-i-1] = -0.5 * (z[i] - z[i-1]) * beta12
    B_fd[i] = (z[i] - z[i-1]) * beta1 * T_b_func(0.5 * (z[i] + z[i-1]))
    # Upward pipe
    A_fd[-i,-i] = 1. + 0.5 * (z[-i] - z[-i-1]) * (beta2 + beta12)
    A_fd[-i,-i-1] = -1. + 0.5 * (z[-i] - z[-i-1]) * (beta2 + beta12)
    A_fd[-i,i] = -0.5 * (z[i] - z[i-1]) * beta12
    A_fd[-i,i-1] = -0.5 * (z[i] - z[i-1]) * beta12
    B_fd[-i] = (z[i] - z[i-1]) * beta2 * T_b_func(0.5 * (z[i] + z[i-1]))
# Boundary conditions
A_fd[0,0] = 1.
B_fd[0] = T_f_in
A_fd[nz,nz] = 1.
A_fd[nz,nz-1] = -1.

T_f_fd = np.linalg.solve(A_fd, B_fd)

# Configure figure and axes
fig = gt.utilities._initialize_figure()

ax = fig.add_subplot(111)
# Axis labels
ax.set_xlabel(r'Temperature [degC]')
ax.set_ylabel(r'Depth from borehole head [m]')
ax.set_title(f'Tmedia = {Tmedia:.0f} [degC]')
gt.utilities._format_axes(ax)

# Plot temperatures
ax.plot(T_f_initial[:,0], z, 'b-', label='Old model')
ax.plot(T_f_initial[:,1], z, 'b-')
ax.plot(T_f_new[:,0], z, 'r-', label='New model')
ax.plot(T_f_new[:,1], z, 'r-')
ax.plot(T_f_fd[:nz], z, 'k:', label='Finite difference')
ax.plot(T_f_fd[-nz:], z[::-1], 'k:')
ax.plot(np.repeat(T_b, 2), np.repeat(z_u, 2)[1:-1], 'g--', label='Borehole wall')
ax.legend()

# Reverse y-axes
ax.set_ylim(ax.get_ylim()[::-1])
# Adjust to plot window
plt.tight_layout()

Uniform borehole wall temperature (nSegments=1)

Non-uniform (random) borehole wall temperature (nSegments=10)

[MassimoCimmino]

To do :

• [x] Cleanup and implement in SingleUTube class
• [x] Test g-function generation
• [ ] Decide on tolerance parameter (Not needed)
• [x] Check if new expressions are also needed for inlet/outlet temperatures and effective thermal resistance
• [x] Add tests

[MassimoCimmino]

The initial implementation is on branch https://github.com/MassimoCimmino/pygfunction/tree/issue266_negativeGFunctionsCoaxial

It seems to solve the issue of negative g-functions. Something interesting is that nSegments=8 does not seem do be enough in the case of the initial post. Here is the g-function for Tmedia=60. using nSegments=24.