Une première explication exhaustive du modèle de la P.E.T. et une proposition de correction des erreurs qu’on y trouve (lien vers l’article).
La version corrigée a été soumise aux développeurs de LadyBug : affaire à suivre !
Une première explication exhaustive du modèle de la P.E.T. et une proposition de correction des erreurs qu’on y trouve (lien vers l’article).
La version corrigée a été soumise aux développeurs de LadyBug : affaire à suivre !
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))