Le calcul de la SET en Python

Toujours dans les modèles de confort, on donne ci-dessous le code de la SET*. L’indicateur est calculé selon les modèles les plus à jour décrits dans l’article « Modélisation du confort » de ce blog.

###################################################################
#_______  ______    _______  _______ 
#|   _   ||    _ |  |       ||       |
#|  |_|  ||   | ||  |    ___||    _  |
#|       ||   |_||_ |   |___ |   |_| |
#|       ||    __  ||    ___||    ___|
#|   _   ||   |  | ||   |___ |   |    
#|__| |__||___|  |_||_______||___|    
#
#  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