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