S.E.T. calculation in Python
Stil in the comfort models, the SET* comfort index code is given below. The index is calculated according to the most up-to-date models described in the article Comfort Modeling on this website.
Show code
###################################################################
#_______ ______ _______ _______
#| _ || _ | | || |
#| |_| || | || | ___|| _ |
#| || |_||_ | |___ | |_| |
#| || __ || ___|| ___|
#| _ || | | || |___ | |
#|__| |__||___| |_||_______||___|
#
# AREP, 16 avenue d'Ivry,75013 Paris, FRANCE
##################################################################
#Calcul de la SET en Python
##################################################################
# Last modification : 12/03/2018 #
##################################################################
# Copyright (C) 2018 Edouard Walther #
# This program is free software; you can redistribute it and/or #
# modify it under the terms of the GNU General Public License #
# as published by the Free Software Foundation; either version #
# of the License, or (at your option) any later version. #
##################################################################
# Contact : edouard[dot]walther[at]arep[dot]fr
##################################################################
import math
import random
def pv_sat(T):
if T >= 0:
pv_sat2 = pow(10, (2.7877 + (7.625*T)/(241 + T)))
else :
pv_sat2 = pow(10,(2.7877 + (9.756*T)/(272.7 + T)))
return pv_sat2
def pv_calc(T, RH):
pv_calc2 = (RH * pv_sat(T))/100
return pv_calc2
def w(T, RH, p):
w2 = 0.622 * (pv_calc(T, RH)/(p - pv_calc(T, RH)))
return w2
def Cp_ah(T, RH, p):
cpa = 1006
cpv = 1830
water = w(T, RH, p)
Cp_ah2 = (cpa + water * cpv)/(1 + water)
return Cp_ah2
def v_spe(T, RH, p):
v_spe2 = (461.24 * (T + 273.15) * (0.622 + w(T, RH, p)))/p
return v_spe2
def h_radiation(T_rad, T_surf):
h_radiation2 = 0.72 * 5.67 * 0.00000001 * ((T_rad + T_surf) + 2 * 273.15) * (pow((T_rad + 273.15), 2) + pow((T_surf + 273.15), 2)) * 0.97 #<--- emissivite = 0,97
return h_radiation2
###############################################################################
############################ FONCTION SET* ####################################
###############################################################################
def modele_metabolique_SET(RadTempMtx, WindSpeedMtx, T_ambient, phi_ambient, p_ambient, hauteur, masse, fat, Cst_dilat,Cst_sweat, Cst_constr, T_core_set, T_skin_set, SkinBloodFlow_set,U_muscle_fat_skin, C_shiv):
T_skin = T_skin_set
T_core = T_core_set
dT = 60
#p_ambient = 101325
#met = 1.1
#clo = 1
#i_m = 0.45
# metabolisme masculin en W
age=30
surface = 0.203*pow(hauteur,0.725)*pow(masse,0.425) #Surface exterieure du sujet [m2]
#genre=random.random()
genre = 0.1 # homme
if genre<0.5:
metabolisme_W = 3.45 * math.pow(masse, 0.75) * (1.0 + 0.004 * (30.0 - age) + 0.01 * (hauteur * 100.0 / math.pow(hauteur, 1.0 / 3.0) - 43.4))
else:
metabolisme_W = 3.19 * math.pow(masse, 0.75) * (1.0 + 0.004 * (30.0 - age) + 0.018* (hauteur * 100.0 / math.pow(hauteur, 1.0 / 3.0) - 42.1))
met = metabolisme_W/surface/58.2
#print surface, met
SkinBloodFlow = 6.3
minutes_metab = 60
minutes = minutes_metab - 1
duree=(minutes + 1) * 60
temps=0
while temps < duree:
#variables fonctionnelles
compteur = 0 #variable de boucle
temps = temps + dT
i_m = 0.45
i_m_static = i_m
clo = 1 #0.155 m2.K/W
clo_static = clo
#evolution metabolique dynamique
#if temps < duree/2:
# met = 2.2
# v_walk = 4/3.6
#else:
# met = 1.2
# v_walk = 0
#WindSpeedMtx = 0.2
#clo dynamique
#if v_walk < 0.7:
# v_walk = 0.0052 * (met * 58.2 - 58)
#corr_T = math.exp(0.042 - 0.398 * WindSpeedMtx + 0.066 * WindSpeedMtx**2 - 0.378 * v_walk + 0.094 * v_walk**2)
#if WindSpeedMtx > 3.5:
# corr_T = 0.582
# ATTENTION : ICI CHANGER "WindSpeedMtx" car sinon on le change a chaque pas de temps !
#WindSpeedMtx = math.sqrt(WindSpeedMtx**2 + v_walk**2)
#clo = clo_static * corr_T
#i_m = i_m_static * (4.9 - 6.5 * corr_T + 2.6 * corr_T**2)
#Constantes du corps humain
KCLO = 0.25 #coefficient augmentation surface d'echange
# masse = vect[1] #masse moyenne sur une population
# R_muscle_fat_skin = 5.28
#Cp_body = 0.97*3600 #capacite calorifique du corps humain [J/(kg.K)]
body_mass = masse
# fat = 15 # pourcentage masse graisseuse
fat_mass = fat/100*body_mass
Cp_body = fat_mass/body_mass*2510 + (body_mass - fat_mass)/body_mass*3650 # Modele HAVENITH
#Constantes de regulation de l'organisme
SBC = 0.0000000567 #Constante de Stefann-Boltzmann
# Cst_sweat = 170
# Cst_dilat = 200
# Cst_constr = 0.5
#Valeurs de consignes de la regulation du corps humain
# T_skin_set = 33.7 #temperature de peau
# T_core_set = 36.8 #temperature interne
T_body_set = 0.1*T_skin_set + 0.9*T_core_set #temperature corporelle
# SkinBloodFlow_set = 6.3 #debit sanguin [L/m2.h]
#Conversion d'unites
p_atmosphere = p_ambient/1000 #conversion en kPa
p_atmosphere = p_atmosphere*0.009869 #conversion en atm
R_clo = 0.155*clo
#correction de la veture ##########
if clo < 0.5:
f_surf_clo = 1 + 0.2*clo
else:
f_surf_clo = 1 + 0.15*clo
#calcul du nombre de Lewis
Lewis = 2434 * v_spe(T_ambient, phi_ambient*100, p_ambient)/(Cp_ah(T_ambient, phi_ambient*100, p_ambient) * 1.04 * pow(0.83, (2/3))) * (18/8.32/(T_ambient + 273.15))
#Calculs initiaux du metabolisme
RM = met * 58.2
Metab = met * 58.2
w_crit = 0.59 * pow(WindSpeedMtx, (-0.08))
#Calcul des coefficients d'echange convectif
h_c = 3 * pow(p_atmosphere, 0.53)
h_c_vent = 8.600001 * pow(WindSpeedMtx * p_atmosphere, 0.53)
h_c = max(h_c, h_c_vent)
#Coefficient d'echange radiatif
h_r = 4.7
#Coefficient d'echange global
h_g = h_r + h_c
#Resistance thermique convective+radiative
R_air = 1/(f_surf_clo * h_g)
#Temperature operative
T_op = (h_r * RadTempMtx + h_c * T_ambient)/h_g #Sous forme de matrice
#Temperature superficielle de veture
T_clo = T_op + (T_skin - T_op)/(h_g * (R_air + R_clo))
T_clo_OLD = T_clo + 0.5
############### Boucle calcul T_clo ###########################################
while abs(T_clo - T_clo_OLD) > 0.001:
T_clo_OLD = T_clo
h_r = h_radiation(RadTempMtx, T_clo)# Sous forme de matrice
h_g = h_r + h_c
R_air = 1/(f_surf_clo * h_g)
T_op = (h_r * RadTempMtx + h_c * T_ambient)/h_g# Sous forme de matrice
T_clo = (R_air * T_skin + R_clo * T_op)/(R_air + R_clo)
compteur = compteur + 1
if compteur > 20:
break
###############################################################################
#SkinBloodFlow = SkinBloodFlow_set
#Temperature corporelle
#alpha = 0.0417737 + 0.7451833/(SkinBloodFlow + 0.585417)
#T_body = alpha * T_skin + (1 - alpha) * T_core
################ Modele de regulation du corps humain #######################
#Skin signal
signal_skin = T_skin - T_skin_set
if signal_skin > 0:
warm_skin = signal_skin
cold_skin = 0
else :
warm_skin = 0
cold_skin = -signal_skin
# Core signal
signal_core = T_core - T_core_set
if signal_core > 0:
warm_core = signal_core
cold_core = 0
else:
warm_core = 0
cold_core = -signal_core
#Debit sanguin
SkinBloodFlow = (SkinBloodFlow_set + Cst_dilat * warm_core)/(1 + Cst_constr * cold_skin)
if SkinBloodFlow > 90 :
SkinBloodFlow = 90
if SkinBloodFlow < 0.5 :
SkinBloodFlow = 0.5
#Temperature corporelle
alpha = 0.0417737 + 0.7451833/(SkinBloodFlow + 0.585417)
T_body = alpha * T_skin + (1 - alpha) * T_core
#Corps/Body
signal_body = T_body - T_body_set
if signal_body > 0 :
warm_body = signal_body
cold_body = 0
else :
warm_body = 0
cold_body = -signal_body
#Debit sudation
qm_sweat = Cst_sweat * warm_body * math.exp(warm_skin/10.7)
if qm_sweat > 500:
qm_sweat = 500
#Chaleur latente maximale echangee
R_vap_tot = (R_clo + R_air)/(Lewis * i_m)
E_max = (pv_sat(T_skin) - phi_ambient * pv_sat(T_ambient))/R_vap_tot
h_e = (2.2 * h_c)/(1 + 0.928 * R_clo * h_c)/133.322
#E_max = h_e * (pv_sat(T_skin) - phi_ambient * pv_sat(T_ambient))
#Chaleur latente sudation
E_sweat = 0.68 * qm_sweat
#Ratios
pcent_sweat = E_sweat/E_max #Part d'energie echangee due a la sudation adimensionnel
pcent_wet = 0.06 + 0.94 * pcent_sweat #Part de la surface du corps mouille
#Chaleur latente echangee par diffusion de l'eau a travers la couche cutanee
E_diff = pcent_wet * E_max - E_sweat
#Chaleur latente totale echangee par la peau
E_skin = E_sweat + E_diff
#Sudation supercritique
if pcent_wet > w_crit :
pcent_wet = w_crit
pcent_sweat = w_crit/0.94
E_sweat = pcent_sweat * E_max
E_diff = 0.06 * (1 - pcent_sweat) * E_max
E_skin = E_sweat + E_diff
drip_cond_nope = 1
#Condensation (la pression de vapeur saturante a la temperature de la peau est inferieure a la pression de vapeur de l'air ambiant
elif E_max < 0 :
E_diff = 0
E_sweat = 0
pcent_wet = w_crit
pcent_sweat = w_crit
E_skin = E_max
drip_cond_nope = -2
else:
drip_cond_nope = 0
w_skin = pcent_wet
#Frisson/shivering
M_shiv = C_shiv * cold_skin * cold_core
Metab = RM + M_shiv
w_sweat_global = E_sweat/E_skin #Part d'humidite due a la sudation
#Flux echange entre l'interieur du corps (noyau) et la peau
Flx_core_skin = (T_core - T_skin) * (U_muscle_fat_skin + 1.163 * SkinBloodFlow) # 1.163 = 4 200/3 600 * 1000
#Chaleur sensible echangee entre la peau et l'exterieur
DRY = (T_skin - T_op)/(R_air + R_clo)
#Chaleur sensible echangee par la respiration
T_expir = 32.5 + 0.066 * T_ambient + 1.99 * 0.000001 * phi_ambient * pv_sat(T_ambient)
C_resp = 0.0014 * Metab * (T_expir - T_ambient) #
#Chaleur latente echangee par la respiration
E_resp = 0.000017251 * Metab * (pv_sat(35.5) - phi_ambient * pv_sat(T_ambient))
#Accumulation d'energie par la peau
SSK = Flx_core_skin - DRY - E_skin
Accumulation_skin = SSK
#Accumulation d'energie par le corps
SCR = Metab - Flx_core_skin - E_resp - C_resp
Accumulation_core = SCR
#Modification des temperatures par l'effet de l'accumulation
#Methode 1
dT_skin = SSK * surface * dT /(alpha * masse * Cp_body)
dT_core = SCR * surface * dT / ((1-alpha) * masse * Cp_body)
#Methode 2
TCSK = Cp_body/3600 * alpha * masse
TCCR = Cp_body/3600 * (1 - alpha) * masse
DTSK = (SSK * surface)/(TCSK * dT)
DTCR = (SCR * surface)/(TCCR * dT)
T_skin = T_skin + DTSK
T_core = T_core + DTCR
T_body = alpha * T_skin + (1 - alpha) * T_core
#Energie totale echangee par la peau
H_skin = DRY + E_skin
#Metabolisme net
RN = Metab
#Chaleur latente echangee par sudation en etat de confort
E_comf = 0.42 * (RN - (1*58.2))
if E_comf < 0:
E_comf = 0
E_max = E_max * w_crit
#Chaleur latente evaporative requise pour la thermoregulation
E_req = RN - E_resp - C_resp - DRY
E_sweat_global = E_sweat
#Coefficient d'echange chaleur sensible
HD = 1/(R_air + R_clo)
#Coefficient d'echange evaporatif
HE = 1/R_vap_tot
#Pression de vapeur saturante a la temperature cutanee
PSSK = pv_sat(T_skin)
#Coefficients d'echange
h_r_SET = h_r
if met < 0.85:
h_c_SET = 3
else:
h_c_SET = 5.66 * pow((met - 0.85), 0.39)
if h_c_SET < 3:
h_c_SET = 3
h_g_SET = h_c_SET + h_r_SET
#CLO metabolique
RCLOS = 1.52/(met + 0.6944) - 0.1835
#Resistance thermique de la veture
RCLS = 0.155 * RCLOS
#Correction de la surface d'echange due a la veture
f_surf_clo_SET = 1 + KCLO * RCLOS
#Facteur d'efficacite de BURTON
F_clo_SET = 1/(1 + 0.155 * f_surf_clo_SET * h_g_SET * RCLOS)
#Index de permeabilite de la veture
i_m_SET = 0.38
#Indice de permeation de la veture
i_clo_SET = i_m_SET * h_c_SET/h_g_SET * (1 - F_clo_SET)/(h_c_SET/h_g_SET - F_clo_SET * i_m_SET)
#Resistance convective + radiative corrigee
R_air_SET = 1/(f_surf_clo_SET * h_g_SET)
#Resistances evaporatives
REAS = 1/(Lewis * f_surf_clo_SET * h_c_SET)
RECLS = RCLS/(Lewis * i_clo_SET)
#Resistance totale au transfert de chaleur sensible
HD_S = 1/(R_air_SET + RCLS)
HE_S = 1/(REAS + RECLS)
#Variables de resolution de SET/ET
Delta = 0.0001
dx = 100
X_OLD = T_skin - H_skin/HD_S
while abs(dx) > 0.001:
#compteurr=compteurr+1
ERR1 = H_skin - HD_S * (T_skin - X_OLD) - w_skin * HE_S * (PSSK - 0.5 * pv_sat(X_OLD))
ERR2 = H_skin - HD_S * (T_skin - (X_OLD + Delta)) - w_skin * HE_S * (PSSK - 0.5 * pv_sat(X_OLD + Delta))
#if ERR2==ERR1:
# print("attention sortie brutale de boucle...")
# break
x = X_OLD - Delta * ERR1/(ERR2 - ERR1)
dx = x - X_OLD
X_OLD = x
ET_global = x
return ET_global