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

La PET corrigée en téléchargement

Voici le code de la PET modifiée, ainsi que décrit dans l’article suivant . Les corrections mentionnées sont incluses.

##################################################################
#_______  ______    _______  _______ 
#|   _   ||    _ |  |       ||       |
#|  |_|  ||   | ||  |    ___||    _  |
#|       ||   |_||_ |   |___ |   |_| |
#|       ||    __  ||    ___||    ___|
#|   _   ||   |  | ||   |___ |   |    
#|__| |__||___|  |_||_______||___|    
# AREP - 16, av. d'Ivry, 75013 Paris, FRANCE
##################################################################
# based on: Peter Hoeppe PET fortran code, from:
# "Urban climatic map and standards for wind environment - Feasibility study, Technical Input Report No.1",
# Edouard Walther and Quentin Goestchel
# Most of the changes in the implementaion are explained in the resolution function comments                                  #
##################################################################
# Last modification : 10/04/2018                                 #
##################################################################
# Copyright (C) 2018 Édouard 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]com                                
##################################################################


import os
import numpy as np
import math as math
import scipy.optimize as optimize

## Implementation of the skin and core temperatures set values #######
tc_set=36.6 # 36.8
tsk_set=34 # 33.7
tbody_set=0.1*tsk_set+0.9*tc_set # Calculation of the body temperature through a weighted average


## Skin blood flow calculation function: #######
def vasoC(tcore,tsk):
    #Set value signals to consider every cases:
    sig_skin = tsk_set - tsk
    sig_core = tcore - tc_set
    if sig_core<0:
        # In this case, Tcore<Tc_set: the body needs to keep the heat --> the blood flow is reduced
        sig_core=0.
    if sig_skin<0:
        # In this case, Tsk>Tsk_set: the body needs to loose heat --> the blood flow is increased
        sig_skin=0.
    # 6.3 L/m^2/h is the set value of the blood flow
    qmblood = (6.3 + 75. * sig_core) / (1. + 0.5 * sig_skin)
    # 90 L/m^2/h is the blood flow upper limit, not sustainable for a human being
    if qmblood>90:
        qmblood=90.
    # Alpha can be used to calculate tbody in ameliorated models
    #alfa = 0.04177 + 0.74518 / (qmblood + 0.585417)
    alfa = 0.1
    return (qmblood,alfa)


## Sweating flow calculation function: #######
def Suda(tbody,tsk):
    sig_body = tbody - tbody_set
    sig_skin = tsk - tsk_set
    if sig_body<0:
        #In this case, Tbody<Tbody_set: the body needs to keep the heat --> The sweat flow is 0
        sig_body=0.
    if sig_skin<0:
        # In this case, Tsk<Tsk_set: the body needs to keep the heat --> the sweat flow is reduced
        sig_skin=0.
    #qmsw = 170 * sig_body * math.exp((sig_skin) / 10.7)  # [g/m2/h] Expression from Gagge Model
    qmsw = 304.94*10**(-3) * sig_body
    # 90 L/m^2/h is the blood flow upper limit, not sustainable for a human being
    if qmsw > 500:
        qmsw = 500
    return (qmsw)


## Vectorial MEMI balance calculation function for the 3 node model: #######
def Syst(T, Ta, Tmrt, HR, v, age, sex, ht, mbody, pos, M, icl,mode):
    ## Conversion of T vector in an array
    arr = np.ones((3,1))
    arr[0,0]=T[0] #Corresponds to T_core
    arr[1,0]=T[1] #Corresponds to T_skin
    arr[2,0]=T[2] #Corresponds to T_clothes
    T=arr
    enbal_vec = np.zeros((3,1)) #Useful for the vectorial expression of the balance

    ## Area parameters of the body: #######
    Adu = 0.203 * mbody ** 0.425 * ht ** 0.725
    feff=0.725
    if pos == 1 or pos == 3:
        feff = 0.725
    if pos == 2:
        feff = 0.696
    # Calculation of the Burton coefficient, k = 0.31 for Hoeppe:
    fcl = 1 + (0.31 * icl) # Increase of the heat exchange surface  depending on clothing level
    facl = (173.51 * icl - 2.36 - 100.76 * icl * icl + 19.28 * icl ** 3.0) / 100
    Aclo = Adu * facl + Adu * (fcl - 1.0)
    Aeffr = Adu * feff  # Effective radiative area depending on the position of the subject
    # Partial pressure of water in the air depending on relative humidity and air temperature:
    if mode: # actual environment
        vpa = HR / 100.0 * 6.105 * math.exp(17.27 * Ta / (237.7 + Ta )) #[hPa]
    else:# mode==False implies we are calculating the PET
        vpa= 12 # [hPa] vapour pressure of the standard environment

    # Convection coefficient depending on wind velocity and subject position
    hc = 0
    if pos == 1:
        hc = 2.67 + (6.5 *v**0.67)
    if pos == 2:
        hc = 2.26 + (7.42 *v**0.67)
    if pos == 3:
        hc = 8.6 * (v ** 0.513)
    hc = hc * (p / po) ** 0.55

    # Basic metabolism for men and women in [W] #######
    metab_female = 3.19 * mbody**0.75 * (1.0 + 0.004 * (30.0 - age) + 0.018 * (ht * 100.0 / mbody**(1.0/3.0) - 42.1))
    metab_male = 3.45 * mbody**0.75 * (1.0 + 0.004 * (30.0 - age) + 0.01 * (ht * 100.0 / mbody**(1.0/3.0) - 43.4))
    # Source term : metabolic activity
    eswpot = (M + metab_male)/Adu
    fec = (M + metab_female)/Adu
    he = 0.0
    # Attribution of internal energy depending on the sex of the subject
    if sex == 1:
        he = eswpot
    elif sex == 2:
        he = fec
    h = he *(1.0 - eta)  # [W/m2]

    # Respiratory energy losses
    # Expired air temperature calculation:
    texp = 0.47 * Ta + 21.0  # [degC]

    # Pulmonary flow rate
    dventpulm = he * 1.44 * 10.0**(-6.0)

    # Sensible heat energy loss:
    eres = cair * (Ta - texp) * dventpulm  # [W/m2]

    # Latent heat energy loss:
    vpexp = 6.11 * 10.0**(7.45 * texp / (235.0 + texp))
    erel = 0.623 * Lvap / p * (vpa-vpexp) * dventpulm  # [W/m2]
    ere = eres + erel  # [W/m2]

    # Clothed fraction of the body approximation
    rcl = icl / 6.45  # Conversion in m2.K/W
    y = 0
    if facl > 1.0:
        facl = 1.0
    if icl >= 2.0:
        y = 1.0
    if icl > 0.6 and icl < 2.0:
        y = (ht - 0.2)/ht
    if icl <= 0.6 and icl > 0.3:
        y = 0.5
    if icl <= 0.3 and icl > 0.0:
        y = 0.1
    # calculation of the closing radius depending on the clothing level (6.28 = 2* pi !)
    r2 = Adu * (fcl - 1.0 + facl) / (6.28 * ht * y)  # External radius
    r1 = facl * Adu /(6.28 * ht * y)  # Internal radius
    di = r2 - r1

    # Calculation of the thermal resistance of the body
    alpha = vasoC(T[0,0],T[1,0])[1]
    tbody = alpha * T[1,0] + (1 - alpha) * T[0,0]
    htcl = (6.28 * ht * y * di)/(rcl * math.log(r2 / r1)*Aclo)  # [W/(m2.K)]

    # Calculation of sweat losses

    qmsw = Suda(tbody,T[1,0])
    # Lvap Latent heat of evaporation : 2400 [J/g] divided by 3600 for [g/m2/h] to [g/m2/s]
    esw = 2400 * qmsw / 3600  # [W/m2]
    # Saturation vapor pressure at temperature Tsk and for 100% HR
    Pvsk = 6.105 * math.exp((17.27 * (T[1,0]+273.15) - 4717.03)/ (237.7 + T[1,0])) # [hPa]

    rscl=0.155*icl
    Lw = 16.7 * 10 ** (-1)  # [K/hPa] Lewis factor
    he_diff = hc * Lw
    fecl=1/(1+0.92*hc*rscl)
    emax = he_diff * fecl * (Pvsk - vpa)
    w = esw / emax  # Dermal wetness
    if w > 1:
        w=1
        delta = esw-emax
        if delta < 0:
            esw=emax
    if esw < 0:
        esw=0
    i_m=0.3
    R_ecl=(1/(fcl*hc) + rscl)/(Lw*i_m)
    #R_ecl=0.79*1e7 #version Hoeppe
    ediff = (1 - w)*(Pvsk - vpa)/R_ecl  # Basal perspiration
    evap = -(ediff + esw)  # [W/m2]

    # Radiation losses
    # For bare skin area:
    rbare = Aeffr*(1.0 - facl) * emsk * sigm * ((Tmrt + 273.15)**(4.0) - (T[1,0] + 273.15)**(4.0))/Adu
    # For dressed area:
    rclo = feff * Aclo * emcl * sigm * ((Tmrt + 273.15)**(4.0) - (T[2,0] + 273.15)**(4.0))/Adu
    rsum = rclo+rbare


    ## Convection losses #######
    cbare = hc * (Ta - T[1,0]) * Adu * (1.0 - facl)/Adu  # [w/m^2]
    cclo = hc * (Ta - T[2,0]) * Aclo/Adu  # [W/m^2]
    csum = cclo+cbare

    ## Balance equations of the 3-nodes model
    enbal_vec[0,0] = h + ere - (vasoC(T[0,0],T[1,0])[0]/3600*cb+5.28)*(T[0,0]-T[1,0]) # Core balance [W/m^2]
    enbal_vec[1,0] = rbare + cbare + evap + (vasoC(T[0,0],T[1,0])[0]/3600*cb+5.28)*(T[0,0]-T[1,0]) - htcl*(T[1,0]-T[2,0])  # Skin balance [W/m^2]
    enbal_vec[2,0] = cclo + rclo + htcl*(T[1,0]-T[2,0]) # Clothes balance [W/m^2]
    enbal_scal = h + ere + rsum + csum +evap

    if mode:
        return [enbal_vec[0,0],enbal_vec[1,0],enbal_vec[2,0]] #List of the balances values
    else:
        return enbal_scal #Scalar balance used in the PET calculation

# Solving the system
def resolution(Ta, Tmrt, HR, v, age, sex, ht, mbody, pos, M, icl, Tx):
    Tn = optimize.fsolve(Syst,Tx ,args=(Ta, Tmrt, HR, v, age, sex, ht, mbody, pos, M, icl,True))
    return (Tn, 1)


# PET calculation with dichotomy method 
def PET (age, sex, ht, mbody, pos, M, icl, Tstable,a,b,eps):
    # Definition of a function with the input variables of the PET reference situation
    def f(Tx):
        return Syst(Tstable, Tx, Tx, 50, 0.1, age, sex, ht, mbody, pos, M, 0.9,False)
    Ti = a # Start index of the browsing interval
    Tf = b # End index of the browsing interval
    pet = 0
    while Tf-Ti>eps: # Dichotomy loop
        if f(Ti)*f(pet)<0:
            Tf = pet
        else:
            Ti = pet
        pet = (Ti + Tf) / 2.0
    return pet

# Input data
# Environment constants
po = 1013.25 #[hPa]
rob = 1.06 # Blood density kg/L
cb = 3.64 * 1000. # Blood specific heat [j/kg/k]
cair = 1.01 * 1000. # Air specific heat  [J./kg/K-]
emsk = 0.99 # Skin emissivity
emcl = 0.95 # Clothes emissivity
Lvap = 2.42 * 10. ** 6. # Latent heat of evaporation [J/Kg]
sigm = 5.67 * 10. ** (-8.) # Stefan-Boltzmann constant [W/(m2*K^(-4))]
eta = 0. # Body efficiency

# Initialisation of Temperature vector with respectively: Tcore, Tskin, Tcl
T = [38,40,40]
eps = 10**(-3) #Numerical precision

# Dichotomy browsning parameters
a = -40
b = 60
# Input data of the PET 
Ta=50 # Air temperature in oC
Tmrt=50 # Mean radiant temperature in oC
HR=50 # Air relative humidity %
v=0.1 # Wind velocity m/s
age = 35
sex = 1 # 1 for men and 2 for women
pos = 1
mbody = 75 #[Kg]
ht = 1.80 #[m]
p = 1013.25 #[hPa]
M_activity = 80 #[W] Activity
icl = 0.9# [clo] Clothing level

#initialisation pour le premier calcul
Tstable = resolution(Ta,Tmrt,HR,v,age,sex,ht,mbody,pos,M_activity,icl,T)[0]
print(PET(age, sex, ht, mbody, pos, M_activity, icl, Tstable, -30, 90, eps))

IBPSA 2018

AREP a participé à la conférence IBPSA France les 15-16 mai 2018 @ Bordeaux. Trois posters différents ont été présentés, disponibles ci-dessous.

Indicateur de confort P.E.T. : Une revue critique (prix du meilleur poster & pitch)

Download the PDF file .

Calcul trigonométrique du flux solaire reçu par un individu.

Download the PDF file .

Simulation spatiale du confort aéraulique : influence de la rose des vents.

Download the PDF file .