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