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))
