Heat transfer & exchangers¶

ε-NTU heat exchanger calculations, the air-side fan + heat-exchanger model used by ASHP/ASHPB, and the borehole g-function and coupling abstractions used by GSHP/GSHPB.

ε-NTU heat exchanger¶

Heat transfer and fluid dynamics calculations.

tmhp.heat_transfer.TRIDIAG_MATRIX_ALGORITHM(a_M, a_P, a_E, a_W, b_P)[source]¶

Solve tridiagonal matrix system using Thomas algorithm.

Parameters:

a_M (list[float]) – Main diagonal (a_P in standard notation).
a_P (list[float]) – Not used directly, maintained for signature compatibility.
a_E (list[float]) – Upper diagonal.
a_W (list[float]) – Lower diagonal.
b_P (list[float]) – RHS vector.

Returns:

Solution vector.

Return type:

list[float]

tmhp.heat_transfer.calc_LMTD_counter_flow(Th_in, Th_out, Tc_in, Tc_out)[source]¶

Calculate Log-Mean Temperature Difference for counter-flow heat exchanger.

Parameters:

Th_in (float) – Hot stream inlet temp [K].
Th_out (float) – Hot stream outlet temp [K].
Tc_in (float) – Cold stream inlet temp [K].
Tc_out (float) – Cold stream outlet temp [K].

Returns:

LMTD [K].

Return type:

float

tmhp.heat_transfer.calc_LMTD_parallel_flow(Th_in, Th_out, Tc_in, Tc_out)[source]¶

Calculate Log-Mean Temperature Difference for parallel-flow heat exchanger.

Parameters:

Th_in (float) – Hot stream inlet temp [K].
Th_out (float) – Hot stream outlet temp [K].
Tc_in (float) – Cold stream inlet temp [K].
Tc_out (float) – Cold stream outlet temp [K].

Returns:

LMTD [K].

Return type:

float

tmhp.heat_transfer.calc_UA_tank_arr(arr_D_in, arr_D_out, arr_L, arr_k, h_in, h_out)[source]¶

Calculate thermal conductance (UA) of a multi-layer cylindrical tank.

Parameters:

arr_D_in (list[float]) – Inner diameters of layers [m].
arr_D_out (list[float]) – Outer diameters of layers [m].
arr_L (list[float]) – Lengths of layers [m].
arr_k (list[float]) – Thermal conductivities of layers [W/mK].
h_in (float) – Inner convection coefficient [W/m2K].
h_out (float) – Outer convection coefficient [W/m2K].

Returns:

Thermal conductance [W/K].

Return type:

float

tmhp.heat_transfer.calc_h_vertical_plate(T_s, T_inf, L, fluid='Air', is_active=True)[source]¶

Calculate natural convection heat transfer coefficient for a vertical plate.

Parameters:

T_s (float) – Surface temperature [K].
T_inf (float) – Fluid temperature [K].
L (float) – Characteristic length [m].
fluid (str) – Fluid name. Default is ‘Air’.
is_active (bool) – If False, returns np.nan.

Returns:

Heat transfer coefficient [W/m2K].

Return type:

float

tmhp.heat_transfer.calc_simple_tank_UA(r0=0.2, H=0.8, x_shell=0.01, x_ins=0.1, k_shell=25.0, k_ins=0.03, h_o=10.0)[source]¶

Calculate simple tank UA value.

Parameters:

r0 (float) – Tank radius [m]
H (float) – Tank height [m]
x_shell (float) – Shell thickness [m]
x_ins (float) – Insulation thickness [m]
k_shell (float) – Shell thermal conductivity [W/mK]
k_ins (float) – Insulation thermal conductivity [W/mK]
h_o (float) – External convective heat transfer coefficient [W/m²K]

Returns:

Tank UA value [W/K]

Return type:

float

tmhp.heat_transfer.darcy_friction_factor(Re, e, d, is_active=True)[source]¶

Calculate the Darcy friction factor.

Uses Haaland equation.

Parameters:

Re (float) – Reynolds number.
e (float) – Surface roughness [m].
d (float) – Diameter [m].
is_active (bool) – If False, returns np.nan.

Returns:

Friction factor.

Return type:

float

Air-side fan & heat exchanger¶

Air-side fan UA and power utility functions.

Airflow-based correlations for the outdoor-coil heat exchanger: a velocity-dependent UA scaling law (calc_UA_from_dV_fan) and a fan power curve (calc_fan_power_from_dV_fan). HX performance solving lives in enex_functions and HP schedule checking in dynamic_context.

tmhp.hx_fan.is_generic_fan_curve(vsd_coeffs)[source]¶

Identify the unmodified ASHRAE Appendix G Method 2 correlation.

Parameters:: vsd_coeffs (dict | None)
Return type:: bool

tmhp.hx_fan.calc_ashrae_fan_power_ratio(flow_ratio)[source]¶

Raw Appendix G Method 2 fit for table parity, not a flow controller.

Table G3.1.3.15 Method 1 spans 0–1. Its x=0.1 row can be compared mathematically, but generic operating control separately enforces x>=0.15. Source: https://www.ashrae.org/file%20library/technical%20resources/standards%20and%20guidelines/standards%20addenda/90.1-2016/90_1_2016_be_bm_bn_bo_bp_br_bs_bu_bv_cf_cl_cm_cq_ct_cu_cv_cw_cy_20210324.pdf

Parameters:: flow_ratio (float)
Return type:: float

tmhp.hx_fan.calc_fan_operating_point(dV_fan, fan_params, vsd_coeffs, is_active=True)[source]¶

Clamp fan flow to its model/control limits, then compute power.

The default is the single-zone VAV surrogate with an independent 10% electrical floor. Its 15% airflow bound is a TMHP control assumption. The following fixed-SP curve remains an explicit alternative.

Source: ANSI/ASHRAE/IES Standard 90.1-2016, Appendix G, Table G3.1.3.15, Method 2: P* = 0.0013 + 0.1470*x + 0.9506*x**2 - 0.0998*x**3. x = actual airflow / fixed reference airflow. The generic correlation is restricted to 0.15 <= x <= 1.0. ASHRAE 2025 Fundamentals Chapter 19 cautions against extrapolation below a minimum ratio (example 0.15). 90.1-2022 Addendum u is supporting multizone-VAV turndown evidence; its informative foreword discusses 16% power at 15% airflow, not a normative 16% requirement or a universal electrical floor. The amended body 6.5.3.2.1(b) removes the old 30% power sentence without adding 16%. Its airflow turndown provision includes the design minimum outdoor air. Evidence: references/fan_model/01_ashrae_90_1_fan_curve.png through 04_ashrae_addendum_u_normative.png, captured in Windows Chrome via CDP. Generic 0.15–1.0 is a TMHP assumption, not a heat-pump requirement. Nondefault user coefficients define a custom curve with explicit limits. A raw demand is clamped here; HX solvers must separately close heat duty. An optional power_min_ratio implements a part-flow electrical power floor independently of those airflow bounds: P/P_ref=max(polynomial, floor). It does not invert the power curve to impose a minimum airflow. The SINGLE_ZONE_VAV_COEFFICIENTS alternative comes from PNNL-26917, PRM Reference Manual, Eq. (11) / Table 50 (printed pp. 3.153–3.154). That row is a single-zone/static-pressure-reset modeling surrogate, not a measured fan curve or a requirement for every heat pump.

Parameters:

dV_fan (float) – Current flow rate [m³/s].
fan_params (dict) – Must contain fan_design_flow_rate and fan_design_power.
vsd_coeffs (dict | None) – VSD Curve coefficients (c1 through c5).
is_active (bool) – If False, reports zero actual flow and np.nan power.

Returns:

Actual flow, ratio, power [W] and min/max limit flags.

Return type:

dict

tmhp.hx_fan.calc_fan_power_from_dV_fan(dV_fan, fan_params, vsd_coeffs, is_active=True)[source]¶

Electrical power [W] at the bounded fan operating point.

ASHRAE Appendix G Method 2 uses its unchanged empirical coefficients only over 0.15–1.0 of reference airflow. See calc_fan_operating_point for actual flow and limit flags. The default single-zone surrogate has a 10% electrical floor; custom settings may explicitly supply an independent power_min_ratio.

Parameters:

dV_fan (float)
fan_params (dict)
vsd_coeffs (dict | None)
is_active (bool)

Return type:

float

Borehole g-function¶

Borehole g-function and air property helpers.

Provides: - Finite Line Source (FLS) g-function for borehole heat exchangers - Air dynamic viscosity (Sutherland’s formula) and Prandtl number

tmhp.g_function.f(x)[source]¶

Helper function for G-function calculation.

Parameters:: x (float) – Input value
Returns:: f(x) = x*erf(x) - (1-exp(-x²))/√π
Return type:: float

tmhp.g_function.chi(s, rb, H, z0=0)[source]¶

Helper function for G-function calculation.

Parameters:

s (float) – Integration variable
rb (float) – Borehole radius [m]
H (float) – Borehole height [m]
z0 (float, optional) – Reference depth [m] (default: 0)

Returns:

chi function value

Return type:

float

tmhp.g_function.G_FLS(t, ks, as_, rb, H)[source]¶

Calculate the g-function for finite line source (FLS) model.

This function calculates the g-function used in ground source heat pump analysis. Results are cached for performance.

Parameters:

t (float | ndarray) – Time [s]
ks (float) – Ground thermal conductivity [W/mK]
as (float) – Ground thermal diffusivity [m²/s]
rb (float) – Borehole radius [m]
H (float) – Borehole height [m]
as_ (float)

Returns:

g-function value [mK/W]. Returns scalar for single time value, array for multiple time values.

Return type:

float | ndarray

tmhp.g_function.precompute_gfunction(N_1, N_2, B, H_b, D_b, r_b, alpha_s, k_s, t_max_s, dt_s, boundary_condition='UBWT')[source]¶

Precompute g-function using pygfunction and return an interpolator.

Creates a rectangular borehole field and computes the g-function for log-spaced time steps up to t_max_s (plus an extra margin). Returns a callable interp1d object predicting the g-function [mK/W].

The g-function is built from the finite line source (FLS) solution of Claesson and Javed (2011), extended to boreholes at different vertical positions by Cimmino and Bernier (2014). Because a borehole has finite length, the wall temperature is not uniform along its depth, so an axial boundary condition must be assumed:

'UHTR' (uniform heat transfer rate): every depth exchanges the same heat per unit length; the wall temperature then varies with depth and the g-function is its depth-average.
'UBWT' (uniform borehole wall temperature): the wall is at a single temperature over the whole depth and the heat flux varies instead. This is Eskilson’s original definition and is the more realistic of the two for a circulated borehole.

Parameters:

N_1 (int) – Number of boreholes in x-direction.
N_2 (int) – Number of boreholes in y-direction.
B (float) – Borehole spacing [m].
H_b (float) – Borehole depth/length [m].
D_b (float) – Buried depth [m].
r_b (float) – Borehole radius [m].
alpha_s (float) – Ground thermal diffusivity [m²/s].
k_s (float) – Ground thermal conductivity [W/mK].
t_max_s (float) – Maximum simulation time [s].
dt_s (float) – Simulation timestep [s].
boundary_condition (str) – Axial boundary condition along the borehole. Either 'UBWT' (uniform borehole wall temperature) or 'UHTR' (uniform heat transfer rate). Default is 'UBWT'.

Returns:

Interpolator function mapping time [s] to g-function [mK/W].

Return type:

interp1d

References

Claesson, J., & Javed, S. (2011). An analytical method to calculate borehole fluid temperatures for time scales from minutes to decades. ASHRAE Transactions, 117(2), 279-288.

Cimmino, M., & Bernier, M. (2014). A semi-analytical method to generate g-functions for geothermal bore fields. International Journal of Heat and Mass Transfer, 70, 641-650.

tmhp.g_function.chi_mfls(s, r, H, x_prime, U, alpha_s, z0=0)[source]¶

Helper function for MFLS (Moving Finite Line Source) G-function calculation.

Ref: Molina-Giraldo et al. (2011), “A moving finite line source model to simulate borehole heat exchangers with groundwater advection”

tmhp.g_function.G_MFLS_Field(times, boreholes, v_gw, theta_gw, rho_w, c_w, alpha_s, k_s, rho_s, c_s)[source]¶

Calculate the spatial superposition of the MFLS response for a bore field.

Parameters:

times (ndarray) – Array of time values [s]
boreholes (list) – List of pygfunction Borehole objects
v_gw (float) – Darcy velocity of groundwater [m/s]
theta_gw (float) – Direction of groundwater flow [rad]
rho_w (float) – Density of groundwater [kg/m³]
c_w (float) – Specific heat capacity of groundwater [J/kgK]
alpha_s (float) – Ground thermal diffusivity [m²/s]
k_s (float) – Ground thermal conductivity [W/mK]
rho_s (float) – Density of ground [kg/m³]
c_s (float) – Specific heat capacity of ground [J/kgK]

Returns:

Dimensional g-values for the entire field over time [mK/W]

Return type:

ndarray

tmhp.g_function.precompute_gfunction_mls(N_1, N_2, B, H_b, D_b, r_b, alpha_s, k_s, rho_s, c_s, v_gw, theta_gw, rho_w, c_w, t_max_s, dt_s)[source]¶

Precompute the MFLS g-function and return an interpolator.

Parameters:

N_1 (int)
N_2 (int)
B (float)
H_b (float)
D_b (float)
r_b (float)
alpha_s (float)
k_s (float)
rho_s (float)
c_s (float)
v_gw (float)
theta_gw (float)
rho_w (float)
c_w (float)
t_max_s (float)
dt_s (float)

Return type:

interp1d

tmhp.g_function.air_dynamic_viscosity(T_K)[source]¶

Calculate air dynamic viscosity using Sutherland’s formula.

Parameters:

T_K (float) – Temperature [K]

Returns:

float – Dynamic viscosity [Pa·s]
Reference (Sutherland’s formula for air)
mu = mu0 * (T/T0)^1.5 * (T0 + S) / (T + S)
where mu0 = 1.716e-5 Pa·s at T0 = 273.15 K, S = 110.4 K

tmhp.g_function.air_prandtl_number(T_K)[source]¶

Calculate air Prandtl number.

Parameters:

T_K (float) – Temperature [K]

Returns:

float – Prandtl number [-]
Note (Pr ≈ 0.71 for air at typical temperatures (20-50°C))
Temperature dependence is weak, so using constant value.

tmhp.g_function.calc_submerged_coil_thermal_resistance(r_out, r_in, D_s, k_p, m_flow_pipe, rho_f, mu_f, cp_f, k_f, v_river=0.5)[source]¶

Calculate the local thermal resistance [mK/W] of a submerged surface water heat exchanger coil.

Uses the Churchill-Bernstein correlation for cross-flow forced convection over a cylinder to estimate the external (river water) convective heat transfer coefficient. It tricks pygfunction’s SingleUTube model into capturing this pure pipe resistance without any ground thermal mass by assigning exceptionally high thermal conductivities to the grout and ground.

Parameters:

r_out (float) – Pipe outer radius [m]
r_in (float) – Pipe inner radius [m]
D_s (float) – Shank spacing (half distance between pipes) [m]
k_p (float) – Pipe thermal conductivity [W/mK]
m_flow_pipe (float) – Mass flow rate per pipe [kg/s]
rho_f (float) – Internal fluid density [kg/m³]
mu_f (float) – Internal fluid dynamic viscosity [Pa·s]
cp_f (float) – Internal fluid specific heat capacity [J/kgK]
k_f (float) – Internal fluid thermal conductivity [W/mK]
v_river (float) – Velocity of the river water cross-flow [m/s]

Returns:

Thermal resistance of the submerged coil [mK/W].

Return type:

float

Ground coupling¶

Ground-coupling abstraction for borehole-heat-exchanger temporal superposition.

This module separates the ground thermal response (how the borehole-wall temperature responds to a history of ground loads) from the heat-pump cycle physics, so the response backend becomes swappable behind one small contract:

AggregateGFunctionCoupler (default) wraps a single field-average g-function interpolator g(tau) [m·K/W] and reproduces the legacy inline temporal superposition byte-for-byte.
External packages (e.g. geolink) can implement the GroundCoupler protocol with a resolved multi-borehole network response — replacing the single lumped g-function with full borehole-to-borehole superposition — without this package depending on them (dependency inversion: tmhp defines the abstraction it needs; the richer implementation lives elsewhere and is injected).

Contract¶

A coupler owns the ground-load pulse history. The host (GroundSourceHeatPumpBoiler) calls GroundCoupler.reset() once per simulation and GroundCoupler.wall_temperature_rise() once per timestep with the current per-length ground load q_unit [W/m]; the coupler returns the borehole-wall temperature rise magnitude dT [K].

Sign convention (matches the legacy GSHPB)¶

q_unit > 0 denotes heat extraction from the ground; the returned dT is the positive magnitude such that the disturbed wall temperature is T_wall = T_undisturbed - dT. (Note this is the opposite sign labelling from geolink’s q' > 0 = injection convention; a geolink coupler maps the sign at its boundary so the value returned here keeps the extraction-positive meaning.)

class tmhp.ground_coupling.GroundCoupler(*args, **kwargs)[source]¶

Bases: Protocol

Pluggable ground-response backend for BHE temporal superposition.

reset(n_steps, time_arr)[source]¶

Initialise pulse-history state for a simulation of n_steps steps.

Parameters:

n_steps (int) – Number of timesteps in the upcoming simulation.
time_arr (ndarray) – Absolute time of each step [s] (uniform grid). Backends that precompute a response matrix at fixed times use this.

Return type:

None

wall_temperature_rise(n, time_arr, q_unit)[source]¶

Borehole-wall temperature rise [K] at step n for load q_unit.

Parameters:

n (int) – Current step index (0-based).
time_arr (ndarray) – Absolute time of each step [s] (same array passed to reset()).
q_unit (float) – Per-length ground load at this step [W/m], extraction-positive.

Return type:

float

__init__(*args, **kwargs)¶

class tmhp.ground_coupling.AggregateGFunctionCoupler(g_interp, pulse_tol=1e-06)[source]¶

Bases: object

Single field-average g-function temporal superposition (legacy default).

Reproduces the legacy _compute_bhe_superposition inner loop exactly: a pulse is recorded whenever the per-length load changes by more than pulse_tol, and the pulse train is convolved with g(tau).

Parameters:

g_interp (Callable[[np.ndarray], np.ndarray]) – Field-average g-function interpolator mapping lag time tau [s] to the dimensional response g [m·K/W].
pulse_tol (float) – Minimum load change [W/m] that registers a new pulse (default 1e-6).

__init__(g_interp, pulse_tol=1e-06)[source]¶

Parameters:: pulse_tol (float)

reset(n_steps, time_arr)[source]¶

Parameters:

n_steps (int)
time_arr (ndarray)

Return type:

None

preview_wall_temperature_rise(n, time_arr, q_unit)[source]¶

Evaluate an uncommitted candidate without modifying load history.

Optional extension for coupled flow control; the original GroundCoupler protocol remains unchanged for legacy external backends.

Parameters:

n (int)
time_arr (ndarray)
q_unit (float)

Return type:

float

wall_temperature_rise(n, time_arr, q_unit)[source]¶

Parameters:

n (int)
time_arr (ndarray)
q_unit (float)

Return type:

float