Incorrect outlet fluid temperature in pipes.SingleUTube

[MassimoCimmino]

An error in pipes.SingleUTube._continuity_condition_base() causes an error in functions that depend on it. The error, however, only causes errors for nSegments > 1 and segment_ratios that are not symmetric in relation with the borehole mid-length. This limits the impact of the error since there is no known advantage in using a non-symmetric discretization and it would have had to be generated manually by the unfortunate user.

The erroneous code is located here:

pygfunction/pygfunction/pipes.py

Lines 1104 to 1109 in 25afb35

	z = self.b._segment_edges(nSegments, segment_ratios=segment_ratios)[::-1]
	F4 = self._F4(z)
	dF4 = F4[:-1] - F4[1:]
	F5 = self._F5(z)
	dF5 = F5[:-1] - F5[1:]
	a_b[0, :] = (dF4 + dF5) / (self._f3(self.b.H) - self._f2(self.b.H))

From Hellström (1991), the integrated functions f_4 and f_5 should be evaluated at (H-z).

This explains why tests are failing in #283. The script shows the error for one of the failing tests.

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.075         # Borehole radius (m)
H = 150.            # Borehole length (m)
borehole = gt.boreholes.Borehole(H, D, r_b, 0., 0.)

# Pipe positions [m]
D_s = 0.05
pos_pipes = [(-D_s, 0.), (D_s, 0.)]
k_s = 2.0               # Ground thermal conductivity [W/m.K]
k_g = 1.0               # Grout thermal conductivity [W/m.K]
k_p = 0.4               # Pipe thermal conductivity [W/m.K]
r_out = 0.02            # Pipe outer radius [m]
r_in = 0.015            # Pipe inner radius [m]
epsilon = 1.0e-06       # Pipe surface roughness [m]
m_flow_borehole = 0.05  # Nominal fluid mass flow rate [kg/s]
m_flow_pipe = m_flow_borehole
# Fluid is propylene-glycol (20 %) at 20 degC
fluid = gt.media.Fluid('MPG', 20.)
# Pipe thermal resistance [m.K/W]
R_p = gt.pipes.conduction_thermal_resistance_circular_pipe(
    r_in, r_out, k_p)
# Convection heat transfer coefficient [W/m2.K]
h_f = gt.pipes.convective_heat_transfer_coefficient_circular_pipe(
    m_flow_pipe, r_in, fluid.mu, fluid.rho, fluid.k, fluid.cp,
    epsilon)
# Film thermal resistance [m.K/W]
R_f = 1.0/(h_f*2*np.pi*r_in)
# Initialize pipe
singleUTube = gt.pipes.SingleUTube(
        pos_pipes, r_in, r_out, borehole, k_s, k_g, R_f + R_p)

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

# Thermal parameters
m_flow_borehole = 0.2  # Nominal fluid mass flow rate [kg/s]
singleUTube._update_model_variables(m_flow_borehole, fluid.cp, 1, None)
delta = singleUTube._delta
gamma = singleUTube._gamma
beta1 = singleUTube._beta1
beta2 = singleUTube._beta2
beta = singleUTube._beta
beta12 = singleUTube._beta12

# Parameters
nz = 800
nSegments = 4
segment_ratios = np.array([0.1, 0.35, 0.40, 0.15])
z = np.linspace(0., H, num=nz)
z_u = singleUTube.b._segment_edges(nSegments, segment_ratios=segment_ratios)
T_b = np.array([1., 2., 3., 1.])
T_f_in = 5.


# Initial model
# -------------
T_f_initial = singleUTube.get_temperature(
    z, T_f_in, T_b, m_flow_borehole, fluid.cp, segment_ratios=segment_ratios)
T_f_out_initital = singleUTube.get_outlet_temperature(
    T_f_in, T_b, m_flow_borehole, fluid.cp, 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_out_initital, 0, 'bo')
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()

With symmetric segment_ratios ( segment_ratios = np.array([0.1, 0.05, 0.30, 0.10, 0.30, 0.05, 0.1]) and T_b = np.array([1., 2., 2., 3., 3., 1., 1.])):