from __future__ import annotations
import math
import numpy as np
from numba import float64, jit, vectorize
from pythermalcomfort.classes_input import GaggeTwoNodesInputs
from pythermalcomfort.classes_return import SET, GaggeTwoNodes
from pythermalcomfort.utilities import Postures, met_to_w_m2, p_sat_torr
[docs]
def two_nodes_gagge(
tdb: float | list[float],
tr: float | list[float],
v: float | list[float],
rh: float | list[float],
met: float | list[float],
clo: float | list[float],
wme: float | list[float] = 0,
body_surface_area: float | list[float] = 1.8258,
p_atm: float | list[float] = 101325,
position: str = Postures.standing.value,
max_skin_blood_flow: float | list[float] = 90,
round_output: bool = True,
max_sweating: float | list[float] = 500,
w_max: float | list[float] = False,
calculate_ce: bool = False,
) -> SET | GaggeTwoNodes:
"""Gagge Two-node model of human temperature regulation Gagge et al (1986)
[Gagge1986]_.
[Gagge1986]_ This model can be used to calculate a variety of indices, including:
* Gagge's version of Fanger's Predicted Mean Vote (PMV). This function uses the Fanger's PMV equations but it replaces the heat loss and gain terms with those calculated by the two-node model developed by Gagge et al. (1986) [Gagge1986]_.
* PMV SET and the predicted thermal sensation based on SET [Gagge1986]_. This function is similar in all aspects to the :py:meth:`pythermalcomfort.models.pmv_gagge` however, it uses the :py:meth:`pythermalcomfort.models.set` equation to calculate the dry heat loss by convection.
* Thermal discomfort (DISC) as the relative thermoregulatory strain necessary to restore a state of comfort and thermal equilibrium by sweating [Gagge1986]_. DISC is described numerically as: comfortable and pleasant (0), slightly uncomfortable but acceptable (1), uncomfortable and unpleasant (2), very uncomfortable (3), limited tolerance (4), and intolerable (5). The range of each category is ± 0.5 numerically. In the cold, the classical negative category descriptions used for Fanger's PMV apply [Gagge1986]_.
* Heat gains and losses via convection, radiation, and conduction.
* The Standard Effective Temperature (SET).
* The New Effective Temperature (ET).
* The Predicted Thermal Sensation (TSENS).
* The Predicted Percent Dissatisfied Due to Draft (PD).
* Predicted Percent Satisfied With the Level of Air Movement (PS).
Parameters
----------
tdb : float or list of floats
Dry bulb air temperature, [°C].
tr : float or list of floats
Mean radiant temperature, [°C].
v : float or list of floats
Air speed, [m/s].
rh : float or list of floats
Relative humidity, [%].
met : float or list of floats
Metabolic rate, [met].
clo : float or list of floats
Clothing insulation, [clo].
wme : float or list of floats, optional
External work, [met]. Defaults to 0.
body_surface_area : float or list of floats, optional
Body surface area, default value 1.8258 [m2]. Defaults to 1.8258.
.. note::
The body surface area can be calculated using the function
:py:meth:`pythermalcomfort.utilities.body_surface_area`.
p_atm : float or list of floats, optional
Atmospheric pressure, default value 101325 [Pa]. Defaults to 101325.
position : str, optional
Select either "sitting" or "standing". Defaults to "standing".
max_skin_blood_flow : float or list of floats, optional
Maximum blood flow from the core to the skin, [kg/h/m2]. Defaults to 90.
round_output : bool, optional
If True, rounds output value. If False, it does not round it. Defaults to True.
max_sweating : float or list of floats, optional
Maximum rate at which regulatory sweat is generated, [kg/h/m2]. Defaults to 500.
w_max : float or list of floats, optional
Maximum skin wettedness (w) adimensional. Ranges from 0 and 1. Defaults to False.
Returns
-------
GaggeTwoNodes
A dataclass containing the results of the two-node model of human temperature regulation.
See :py:class:`~pythermalcomfort.classes_return.GaggeTwoNodes` for more details.
To access the results, use the corresponding attributes of the returned `GaggeTwoNodes` instance, e.g., `result.e_skin`.
Examples
--------
.. code-block:: python
from pythermalcomfort.models import two_nodes_gagge
result = two_nodes_gagge(tdb=25, tr=25, v=0.1, rh=50, clo=0.5, met=1.2)
print(result.w) # 100.0
result = two_nodes_gagge(tdb=[25, 25], tr=25, v=0.3, rh=50, met=1.2, clo=0.5)
print(result.e_skin) # [100.0, 100.0]
"""
# Validate inputs using the TwoNodesInputs class
GaggeTwoNodesInputs(
tdb=tdb,
tr=tr,
v=v,
rh=rh,
met=met,
clo=clo,
wme=wme,
body_surface_area=body_surface_area,
p_atm=p_atm,
position=position,
max_skin_blood_flow=max_skin_blood_flow,
round_output=round_output,
max_sweating=max_sweating,
w_max=w_max,
)
tdb = np.asarray(tdb)
tr = np.asarray(tr)
v = np.asarray(v)
rh = np.asarray(rh)
met = np.asarray(met)
clo = np.asarray(clo)
wme = np.asarray(wme)
position = np.asarray(position)
vapor_pressure = rh * p_sat_torr(tdb) / 100
if calculate_ce:
result = _gagge_two_nodes_optimized_return_set(
tdb,
tr,
v,
met,
clo,
vapor_pressure,
wme,
body_surface_area,
p_atm,
1,
)
return SET(set=result)
(
_set,
e_skin,
e_rsw,
e_max,
q_sensible,
q_skin,
q_res,
t_core,
t_skin,
m_bl,
m_rsw,
w,
w_max,
et,
pmv_gagge,
pmv_set,
disc,
t_sens,
) = np.vectorize(_gagge_two_nodes_optimized, cache=True)(
tdb=tdb,
tr=tr,
v=v,
met=met,
clo=clo,
vapor_pressure=vapor_pressure,
wme=wme,
body_surface_area=body_surface_area,
p_atm=p_atm,
position=position,
calculate_ce=calculate_ce,
max_skin_blood_flow=max_skin_blood_flow,
max_sweating=max_sweating,
w_max=w_max,
)
output = {
"e_skin": e_skin,
"e_rsw": e_rsw,
"e_max": e_max,
"q_sensible": q_sensible,
"q_skin": q_skin,
"q_res": q_res,
"t_core": t_core,
"t_skin": t_skin,
"m_bl": m_bl,
"m_rsw": m_rsw,
"w": w,
"w_max": w_max,
"set": _set,
"et": et,
"pmv_gagge": pmv_gagge,
"pmv_set": pmv_set,
"disc": disc,
"t_sens": t_sens,
}
if round_output:
for key in output:
output[key] = np.around(output[key], 2)
return GaggeTwoNodes(**output)
@jit(nopython=True, cache=True)
def _gagge_two_nodes_optimized(
tdb,
tr,
v,
met,
clo,
vapor_pressure,
wme,
body_surface_area,
p_atm,
position,
calculate_ce=False,
max_skin_blood_flow=90,
max_sweating=500,
w_max=None,
) -> tuple[
float,
float,
float,
float,
float,
float,
float,
float,
float,
float,
float,
float,
float,
float,
float,
float,
float,
float,
]:
# Initial variables as defined in the ASHRAE 55-2020
air_speed = max(v, 0.1)
k_clo = 0.25
body_weight = 70 # body weight in kg
met_factor = 58.2 # met conversion factor
sbc = 0.000000056697 # Stefan-Boltzmann constant (W/m2K4)
c_sw = 170 # driving coefficient for regulatory sweating
c_dil = 120 # driving coefficient for vasodilation ashrae says 50 see page 9.19
c_str = 0.5 # driving coefficient for vasoconstriction
temp_skin_neutral = 33.7
temp_core_neutral = 36.8
alfa = 0.1
temp_body_neutral = alfa * temp_skin_neutral + (1 - alfa) * temp_core_neutral
skin_blood_flow_neutral = 6.3
t_skin = temp_skin_neutral
t_core = temp_core_neutral
m_bl = skin_blood_flow_neutral
# initialize some variables
e_skin = 0.1 * met # total evaporative heat loss, W
q_sensible = 0 # total sensible heat loss, W
w = 0 # skin wettedness
_set = 0 # standard effective temperature
e_rsw = 0 # heat lost by vaporization sweat
e_diff = 0 # vapor diffusion through skin
e_max = 0 # maximum evaporative capacity
m_rsw = 0 # regulatory sweating
q_res = 0 # heat loss due to respiration
et = 0 # effective temperature
e_req = 0 # evaporative heat loss required for tmp regulation
r_ea = 0
r_ecl = 0
c_res = 0 # convective heat loss respiration
pressure_in_atmospheres = p_atm / 101325
length_time_simulation = 60 # length time simulation
n_simulation = 1
r_clo = 0.155 * clo # thermal resistance of clothing, C M^2 /W
f_a_cl = 1.0 + 0.15 * clo # increase in body surface area due to clothing
lr = 2.2 / pressure_in_atmospheres # Lewis ratio
rm = (met - wme) * met_factor # metabolic rate
m = met * met_factor # metabolic rate
e_comfort = 0.42 * (rm - met_factor) # evaporative heat loss during comfort
e_comfort = max(e_comfort, 0)
i_cl = 1.0 # permeation efficiency of water vapour naked skin
if clo > 0:
i_cl = 0.45 # permeation efficiency of water vapour through the clothing layer
if not w_max: # if the user did not pass a value of w_max
w_max = 0.38 * pow(air_speed, -0.29) # critical skin wettedness naked
if clo > 0:
w_max = 0.59 * pow(air_speed, -0.08) # critical skin wettedness clothed
# h_cc corrected convective heat transfer coefficient
h_cc = 3.0 * pow(pressure_in_atmospheres, 0.53)
# h_fc forced convective heat transfer coefficient, W/(m2 °C)
h_fc = 8.600001 * pow((air_speed * pressure_in_atmospheres), 0.53)
h_cc = max(h_cc, h_fc)
if not calculate_ce and met > 0.85:
h_c_met = 5.66 * (met - 0.85) ** 0.39
h_cc = max(h_cc, h_c_met)
h_r = 4.7 # linearized radiative heat transfer coefficient
h_t = h_r + h_cc # sum of convective and radiant heat transfer coefficient W/(m2*K)
r_a = 1.0 / (f_a_cl * h_t) # resistance of air layer to dry heat
t_op = (h_r * tr + h_cc * tdb) / h_t # operative temperature
t_body = alfa * t_skin + (1 - alfa) * t_core # mean body temperature, °C
# respiration
q_res = 0.0023 * m * (44.0 - vapor_pressure) # latent heat loss due to respiration
c_res = 0.0014 * m * (34.0 - tdb) # sensible convective heat loss respiration
while n_simulation < length_time_simulation:
n_simulation += 1
iteration_limit = 150 # for following while loop
# t_cl temperature of the outer surface of clothing
t_cl = (r_a * t_skin + r_clo * t_op) / (r_a + r_clo) # initial guess
n_iterations = 0
tc_converged = False
while not tc_converged:
# 0.95 is the clothing emissivity from ASHRAE fundamentals Ch. 9.7 Eq. 35
if position == Postures.sitting.value:
# 0.7 ratio between radiation area of the body and the body area
h_r = 4.0 * 0.95 * sbc * ((t_cl + tr) / 2.0 + 273.15) ** 3.0 * 0.7
else: # if standing
# 0.73 ratio between radiation area of the body and the body area
h_r = 4.0 * 0.95 * sbc * ((t_cl + tr) / 2.0 + 273.15) ** 3.0 * 0.73
h_t = h_r + h_cc
r_a = 1.0 / (f_a_cl * h_t)
t_op = (h_r * tr + h_cc * tdb) / h_t
t_cl_new = (r_a * t_skin + r_clo * t_op) / (r_a + r_clo)
if abs(t_cl_new - t_cl) <= 0.01:
tc_converged = True
t_cl = t_cl_new
n_iterations += 1
if n_iterations > iteration_limit:
max_iteration_exceeded = "Max iterations exceeded"
raise StopIteration(max_iteration_exceeded)
q_sensible = (t_skin - t_op) / (r_a + r_clo) # total sensible heat loss, W
# hf_cs rate of energy transport between core and skin, W
# 5.28 is the average body tissue conductance in W/(m2 C)
# 1.163 is the thermal capacity of blood in Wh/(L C)
hf_cs = (t_core - t_skin) * (5.28 + 1.163 * m_bl)
s_core = m - hf_cs - q_res - c_res - wme # rate of energy storage in the core
s_skin = hf_cs - q_sensible - e_skin # rate of energy storage in the skin
tc_sk = 0.97 * alfa * body_weight # thermal capacity skin
tc_cr = 0.97 * (1 - alfa) * body_weight # thermal capacity core
d_t_sk = (s_skin * body_surface_area) / (
tc_sk * 60.0
) # rate of change skin temperature °C per minute
d_t_cr = (
s_core * body_surface_area / (tc_cr * 60.0)
) # rate of change core temperature °C per minute
t_skin = t_skin + d_t_sk
t_core = t_core + d_t_cr
t_body = alfa * t_skin + (1 - alfa) * t_core
# sk_sig thermoregulatory control signal from the skin
sk_sig = t_skin - temp_skin_neutral
warm_sk = (sk_sig > 0) * sk_sig # vasodilation signal
colds = ((-1.0 * sk_sig) > 0) * (-1.0 * sk_sig) # vasoconstriction signal
# c_reg_sig thermoregulatory control signal from the skin, °C
c_reg_sig = t_core - temp_core_neutral
# c_warm vasodilation signal
c_warm = (c_reg_sig > 0) * c_reg_sig
# c_cold vasoconstriction signal
c_cold = ((-1.0 * c_reg_sig) > 0) * (-1.0 * c_reg_sig)
# bd_sig thermoregulatory control signal from the body
bd_sig = t_body - temp_body_neutral
warm_b = (bd_sig > 0) * bd_sig
m_bl = (skin_blood_flow_neutral + c_dil * c_warm) / (1 + c_str * colds)
m_bl = min(m_bl, max_skin_blood_flow)
m_bl = max(m_bl, 0.5)
m_rsw = c_sw * warm_b * math.exp(warm_sk / 10.7) # regulatory sweating
m_rsw = min(m_rsw, max_sweating)
e_rsw = 0.68 * m_rsw # heat lost by vaporization sweat
r_ea = 1.0 / (lr * f_a_cl * h_cc) # evaporative resistance air layer
r_ecl = r_clo / (lr * i_cl)
e_req = (
rm - q_res - c_res - q_sensible
) # evaporative heat loss required for tmp regulation
e_max = (math.exp(18.6686 - 4030.183 / (t_skin + 235.0)) - vapor_pressure) / (
r_ea + r_ecl
)
if e_max == 0: # added this otherwise e_rsw / e_max cannot be calculated
e_max = 0.001
p_rsw = e_rsw / e_max # ratio heat loss sweating to max heat loss sweating
w = 0.06 + 0.94 * p_rsw # skin wetness
e_diff = w * e_max - e_rsw # vapor diffusion through skin
if w > w_max:
w = w_max
p_rsw = w_max / 0.94
e_rsw = p_rsw * e_max
e_diff = 0.06 * (1.0 - p_rsw) * e_max
if e_max < 0:
e_diff = 0
e_rsw = 0
w = w_max
e_skin = (
e_rsw + e_diff
) # total evaporative heat loss sweating and vapor diffusion
m_rsw = (
e_rsw / 0.68
) # back calculating the mass of regulatory sweating as a function of e_rsw
met_shivering = 19.4 * colds * c_cold # met shivering W/m2
m = rm + met_shivering
alfa = 0.0417737 + 0.7451833 / (m_bl + 0.585417)
q_skin = q_sensible + e_skin # total heat loss from skin, W
# p_s_sk saturation vapour pressure of water of the skin
p_s_sk = math.exp(18.6686 - 4030.183 / (t_skin + 235.0))
# standard environment - where _s at end of the variable names stands for standard
h_r_s = h_r # standard environment radiative heat transfer coefficient
h_c_s = 3.0 * pow(pressure_in_atmospheres, 0.53)
if not calculate_ce and met > 0.85:
h_c_met = 5.66 * (met - 0.85) ** 0.39
h_c_s = max(h_c_s, h_c_met)
h_c_s = max(h_c_s, 3.0)
h_t_s = (
h_c_s + h_r_s
) # sum of convective and radiant heat transfer coefficient W/(m2*K)
r_clo_s = (
1.52 / ((met - wme / met_factor) + 0.6944) - 0.1835
) # thermal resistance of clothing, °C M^2 /W
r_cl_s = 0.155 * r_clo_s # thermal insulation of the clothing in M2K/W
f_a_cl_s = 1.0 + k_clo * r_clo_s # increase in body surface area due to clothing
f_cl_s = 1.0 / (
1.0 + 0.155 * f_a_cl_s * h_t_s * r_clo_s
) # ratio of surface clothed body over nude body
i_m_s = 0.45 # permeation efficiency of water vapour through the clothing layer
i_cl_s = (
i_m_s * h_c_s / h_t_s * (1 - f_cl_s) / (h_c_s / h_t_s - f_cl_s * i_m_s)
) # clothing vapor permeation efficiency
r_a_s = 1.0 / (f_a_cl_s * h_t_s) # resistance of air layer to dry heat
r_ea_s = 1.0 / (lr * f_a_cl_s * h_c_s)
r_ecl_s = r_cl_s / (lr * i_cl_s)
h_d_s = 1.0 / (r_a_s + r_cl_s)
h_e_s = 1.0 / (r_ea_s + r_ecl_s)
# calculate Standard Effective Temperature (SET)
delta = 0.0001
dx = 100.0
set_old = round(t_skin - q_skin / h_d_s, 2)
while abs(dx) > 0.01:
err_1 = (
q_skin
- h_d_s * (t_skin - set_old)
- w
* h_e_s
* (p_s_sk - 0.5 * (math.exp(18.6686 - 4030.183 / (set_old + 235.0))))
)
err_2 = (
q_skin
- h_d_s * (t_skin - (set_old + delta))
- w
* h_e_s
* (
p_s_sk
- 0.5 * (math.exp(18.6686 - 4030.183 / (set_old + delta + 235.0)))
)
)
_set = set_old - delta * err_1 / (err_2 - err_1)
dx = _set - set_old
set_old = _set
# calculate Effective Temperature (ET)
h_d = 1 / (r_a + r_clo)
h_e = 1 / (r_ea + r_ecl)
et_old = t_skin - q_skin / h_d
delta = 0.0001
dx = 100.0
while abs(dx) > 0.01:
err_1 = (
q_skin
- h_d * (t_skin - et_old)
- w
* h_e
* (p_s_sk - 0.5 * (math.exp(18.6686 - 4030.183 / (et_old + 235.0))))
)
err_2 = (
q_skin
- h_d * (t_skin - (et_old + delta))
- w
* h_e
* (p_s_sk - 0.5 * (math.exp(18.6686 - 4030.183 / (et_old + delta + 235.0))))
)
et = et_old - delta * err_1 / (err_2 - err_1)
dx = et - et_old
et_old = et
tbm_l = (
0.194 / met_to_w_m2
) * rm + 36.301 # lower limit for evaporative regulation
tbm_h = (
0.347 / met_to_w_m2
) * rm + 36.669 # upper limit for evaporative regulation
t_sens = 0.4685 * (t_body - tbm_l) # predicted thermal sensation
if (t_body >= tbm_l) & (t_body < tbm_h):
t_sens = w_max * 4.7 * (t_body - tbm_l) / (tbm_h - tbm_l)
elif t_body >= tbm_h:
t_sens = w_max * 4.7 + 0.4685 * (t_body - tbm_h)
if t_sens > 0 > e_max * w_max - e_comfort - e_diff:
disc = 6.0
else:
disc = (
4.7 * (e_rsw - e_comfort) / (e_max * w_max - e_comfort - e_diff)
) # predicted thermal discomfort
if disc <= 0:
disc = t_sens
if disc > 6:
disc = 6.0
# PMV Gagge
pmv_gagge = (0.303 * math.exp(-0.036 * m) + 0.028) * (e_req - e_comfort - e_diff)
# PMV SET
dry_set = h_d_s * (t_skin - _set)
e_req_set = rm - c_res - q_res - dry_set
pmv_set = (0.303 * math.exp(-0.036 * m) + 0.028) * (e_req_set - e_comfort - e_diff)
# Predicted Percent Satisfied With the Level of Air Movement
# ps = 100 * (1.13 * (t_op**0.5) - 0.24 * t_op + 2.7 * (v**0.5) - 0.99 * v)
return (
_set,
e_skin,
e_rsw,
e_max,
q_sensible,
q_skin,
q_res,
t_core,
t_skin,
m_bl,
m_rsw,
w,
w_max,
et,
pmv_gagge,
pmv_set,
disc,
t_sens,
)
@vectorize(
[ # signature: output(float64)
float64(
float64,
float64,
float64,
float64,
float64,
float64,
float64,
float64,
float64,
float64,
),
],
cache=True,
)
def _gagge_two_nodes_optimized_return_set(
tdb,
tr,
v,
met,
clo,
vapor_pressure,
wme,
body_surface_area,
p_atm,
position,
):
return _gagge_two_nodes_optimized(
tdb=tdb,
tr=tr,
v=v,
met=met,
clo=clo,
vapor_pressure=vapor_pressure,
wme=wme,
body_surface_area=body_surface_area,
p_atm=p_atm,
position=position,
calculate_ce=True,
)[0]