
AREP et le CSTB étaient présent à la 3è édition de la conférence Foam’U, des utilisateurs francophones d’OpenFOAM, à l’Université de Valenciennes, les 23 & 24 mai 2018

Modélisation de l’effet piston en gare :
AREP et le CSTB étaient présent à la 3è édition de la conférence Foam’U, des utilisateurs francophones d’OpenFOAM, à l’Université de Valenciennes, les 23 & 24 mai 2018
Modélisation de l’effet piston en gare :
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 !
Article sur la réduction des incertitudes de modélisation de la ventilation naturelle en STD, dans REHVA Journal.
Le confort dans les espaces semi-ouverts, dans REHVA Journal.
Une extension de l’article précédent aux PM10 et PM2,5 concernant la modélisation de la qualité d’air en gares souterraines : modèle à deux classes de particules.
Un premier article sur le modélisation de la qualité d’air en gares souterraines : modèle à une classe de particules.
Et une routine de calcul du flux solaire incident sur un individu, représenté par un cylindre (méthode décrite dans cet article).
################################################################### #_______ ______ _______ _______ #| _ || _ | | || | #| |_| || | || | ___|| _ | #| || |_||_ | |___ | |_| | #| || __ || ___|| ___| #| _ || | | || |___ | | #|__| |__||___| |_||_______||___| # # AREP, 16 avenue d'Ivry,75013 Paris, FRANCE ################################################################## #Calcul du flux recu par un cylindre # calcul integral pour le flux direct # modelisation anisotrope pour le flux diffus, # d'apres la methode de (Perez et. al 1990) ################################################################## # Last modification : 01/08/2017 # ################################################################## # Copyright (C) 2018 El Mehdi Hamdani et 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 : elmehdi[dot]hamdani[at]arep[dot]fr # edouard[dot]walther[at]arep[dot]fr ################################################################## import math global pi pi = 3.141592653589793 # Cette fonction permet de retourner la valeur du flux diffus circumsolaire, a partir des valeurs du flux diffus horizontal, du flux direct normal, de la hauteur solaire et de l'angle d'incidence def flux_diffus_circumsolaire_perez(diffus_horizontal, direct_normal, hauteur_solaire, angle_incidence): # Valeur de la constante solaire (flux solaire moyen qui parvient sur la terre a la limite de l'atmosphere) en W/m² I0 = 1368 # angle solaire zenithal ksi_z = 90 - math.degrees(hauteur_solaire) # composante decrivant la transparence du ciel epsilon = (((diffus_horizontal + direct_normal)/diffus_horizontal) + 5.535 * 0.000001 * math.pow(ksi_z, 3))/(1 + 5.535 * 0.000001 * math.pow(ksi_z, 3)) # Masse atmospherique (Kasten and Young 1989) a = max(0, math.cos(angle_incidence)) b = max(math.cos(math.radians(85)), math.cos(math.radians(ksi_z))) ma = 1/(math.cos(math.radians(ksi_z)) + 0.50572 * ((96.07995 - ksi_z)**-1.6364)) # Delta represente l'eclairement du ciel delta = ma * diffus_horizontal / I0 # coefficients donnes par PEREZ et al (1990) if epsilon <= 1.065: f11 = 0.013 f12 = 0.764 f13 = -0.1 f21 = -0.058 f22 = 0.127 f23 = -0.023 elif epsilon <= 1.23: f11 = 0.095 f12 = 0.920 f13 = -0.152 f21 = 0 f22 = 0.051 f23 = -0.02 elif epsilon <= 1.5: f11 = 0.464 f12 = 0.421 f13 = -0.28 f21 = 0.064 f22 = -0.051 f23 = -0.002 elif epsilon <= 1.95: f11 = 0.759 f12 = -0.009 f13 = -0.373 f21 = 0.201 f22 = -0.382 f23 = 0.01 elif epsilon <= 2.8: f11 = 0.976 f12 = -0.4 f13 = -0.436 f21 = 0.271 f22 = -0.638 f23 = 0.051 elif epsilon <= 4.5: f11 = 1.176 f12 = -1.254 f13 = -0.462 f21 = 0.295 f22 = -0.975 f23 = 0.129 elif epsilon <= 6.2: f11 = 1.106 f12 = -1.563 f13 = -0.398 f21 = 0.301 f22 = -1.442 f23 = 0.212 else: f11 = 0.934 f12 = -1.501 f13 = -0.271 f21 = 0.42 f22 = -2.917 f23 = 0.249 # coefficient ponderant les densites du flux circumsolaire K1 = max(0, f11 + f12 * delta + pi * ksi_z * f13 / 180) # coefficient ponderant les densites du flux de l'horizon K2 = f21 + f22 * delta + pi * ksi_z * f23 / 180 # flux diffus provenant du rayonnement circumsolaire W/m² diffus_circumsolaire = diffus_horizontal * K1 * a/b #retourne la valeur du flux diffus circumsolaire return diffus_circumsolaire def flux_diffus_horizon_voute_perez(diffus_horizontal, direct_normal, hauteur_solaire, angle_incidence): # Cette fonction permet de retourner la valeur du flux diffus circumsolaire, a partir des valeurs du flux diffus horizontal, du flux direct normal, de la hauteur solaire et de l'angle d'incidence # Valeur de la constante solaire (flux solaire moyen qui parvient sur la terre a la limite de l'atmosphere) en W/m² I0 = 1368 # angle solaire zenithal ksi_z = 90 - math.degrees(hauteur_solaire) # composante decrivant la transparence du ciel epsilon = (((diffus_horizontal + direct_normal)/diffus_horizontal) + 5.535 * 0.000001 * math.pow(ksi_z, 3))/(1 + 5.535 * 0.000001 * math.pow(ksi_z, 3)) # Masse atmospherique (Kasten and Young 1989) a = max(0, math.cos(angle_incidence)) b = max(math.cos(math.radians(85)), math.cos(math.radians(ksi_z))) ma = 1/(math.cos(math.radians(ksi_z)) + 0.50572 * ((96.07995 - ksi_z)**-1.6364)) # Delta represente l'eclairement du ciel delta = ma * diffus_horizontal / I0 # coefficients donnes par PEREZ et al (1990) if epsilon <= 1.065: f11 = 0.013 f12 = 0.764 f13 = -0.1 f21 = -0.058 f22 = 0.127 f23 = -0.023 elif epsilon <= 1.23: f11 = 0.095 f12 = 0.920 f13 = -0.152 f21 = 0 f22 = 0.051 f23 = -0.02 elif epsilon <= 1.5: f11 = 0.464 f12 = 0.421 f13 = -0.28 f21 = 0.064 f22 = -0.051 f23 = -0.002 elif epsilon <= 1.95: f11 = 0.759 f12 = -0.009 f13 = -0.373 f21 = 0.201 f22 = -0.382 f23 = 0.01 elif epsilon <= 2.8: f11 = 0.976 f12 = -0.4 f13 = -0.436 f21 = 0.271 f22 = -0.638 f23 = 0.051 elif epsilon <= 4.5: f11 = 1.176 f12 = -1.254 f13 = -0.462 f21 = 0.295 f22 = -0.975 f23 = 0.129 elif epsilon <= 6.2: f11 = 1.106 f12 = -1.563 f13 = -0.398 f21 = 0.301 f22 = -1.442 f23 = 0.212 else: f11 = 0.934 f12 = -1.501 f13 = -0.271 f21 = 0.42 f22 = -2.917 f23 = 0.249 # coefficient ponderant les densites du flux circumsolaire K1 = max(0, f11 + f12 * delta + pi * ksi_z * f13 / 180) # coefficient ponderant les densites du flux de l'horizon K2 = f21 + f22 * delta + pi * ksi_z * f23 / 180 # flux diffus provenant de l'horizon et de la voute celeste W/m² diffus_horizon_voute = diffus_horizontal * ((1 - K1) * (0.5 * (1 + math.cos(pi/2))) + K2 * math.sin(pi/2)) # retourne la valeur du flux diffus de l'horizon et de la voute celeste return diffus_horizon_voute # Cette fonction permet de retourner la valeur du flux solaire absorbe par un individu, # a partir des valeurs des flux solaires diffus et direct horizontaux, de l'albedo du sol, du jour de l'annee, de l'heure, du numero du fuseau, de la latitude et de la longitude du site def flux_cylindre_anisotrope (diffus_horizontal, direct_horizontal, albedo, jour, heure, fuseau, latitude, longitude_est): # source : Wikipedia. L'equation du temps permet de calculer le temps solaire vrai. equation_temps = 7.678 * math.sin(1.374 + (2 * pi * (jour - 81)/365)) - 9.87 * math.sin(2 * (2 * pi * (jour - 81)/365)) # le temps solaire vrai permet de calculer la position du soleil temps_solaire_vrai = heure + (longitude_est/15) - fuseau - (equation_temps/60) # l'unite de l'angle horaire doit être en radians car dans le langage python l'angle des fonctions trigonometriques est considere comme etant en radians angle_horaire = math.radians((temps_solaire_vrai - 12) * 15) # declinaison en radians declinaison = math.asin(0.4 * math.sin(math.radians(0.986 * jour - 80))) # hauteur solaire en radians hauteur_solaire = math.asin(math.sin(math.radians(latitude)) * math.sin(declinaison) + math.cos(math.radians(latitude)) * math.cos(declinaison) * math.cos(angle_horaire)) #Afin d'eviter des erreurs dans la suite du calcul nous limitons certaines valeurs de hauteur solaire a une valeur nulle if hauteur_solaire < math.radians(2): hauteur_solaire = 0 #L'azimut est l'angle compte a partir du sud, positivement vers l'ouest, negativement vers l'est (temps_solaire_vrai < 12), entre le plan vertical du soleil a l'instant donne et le plan meridien local. if temps_solaire_vrai < 12: azimut = -1 * math.acos((math.sin(math.radians(latitude)) * math.cos(declinaison) * math.cos(angle_horaire) - math.cos(math.radians(latitude)) * math.sin(declinaison)) / math.cos(hauteur_solaire)) else: azimut = math.acos((math.sin(math.radians(latitude)) * math.cos(declinaison) * math.cos(angle_horaire) - math.cos(math.radians(latitude)) * math.sin(declinaison)) / math.cos(hauteur_solaire)) #Cette formule permet de calculer le cosinus de l'angle d'incidence # (angle entre la normale au plan recepteur et le rayon solaire incident) # Le math.cos(0) signifie que le plan suit la trajectoire du soleil. cosinus_angle_incidence = math.cos(hauteur_solaire)*math.sin(0.5*pi)*math.cos(0)+math.sin(hauteur_solaire)*math.cos(0.5*pi) #pour eviter une division par 0 dans la formule ligne 224, lorsque la hauteur solaire est egale a 0 on considere que le flux direct normal est nul if hauteur_solaire == 0: direct_normal = 0 else: direct_normal = direct_horizontal/math.sin(hauteur_solaire) if diffus_horizontal == 0: flux_diffus_circumsolaire = 0 flux_diffus_voute_horizon = 0 else: flux_diffus_circumsolaire = flux_diffus_circumsolaire_perez(diffus_horizontal, direct_normal, hauteur_solaire, math.acos(cosinus_angle_incidence)) # calcul de l'azimut de la face arriere de la paroi verticale en degre if math.degrees(azimut) < 0: azimut_arr = math.degrees(azimut) + 180 else: azimut_arr = math.degrees(azimut) - 180 azimut_arriere = math.radians(azimut_arr) #cosinus de l'angle d'incidence de la face arriere de la paroi verticale cosinus_angle_incidence_arriere = math.cos(hauteur_solaire) * math.sin(0.5 * pi) * math.cos(azimut - azimut_arriere) + math.sin(hauteur_solaire) * math.cos(0.5 * pi) flux_diffus_voute_horizon = flux_diffus_horizon_voute_perez(diffus_horizontal, direct_normal, hauteur_solaire, math.acos(cosinus_angle_incidence)) + flux_diffus_horizon_voute_perez(diffus_horizontal, direct_normal, hauteur_solaire, math.acos(cosinus_angle_incidence_arriere)) if flux_diffus_circumsolaire < 0: flux_diffus_circumsolaire = 0 if flux_diffus_voute_horizon < 0: flux_diffus_voute_horizon = 2 * diffus_horizontal * 0.5 if hauteur_solaire == 0: flux_direct = 0 else: flux_direct = direct_horizontal * cosinus_angle_incidence / math.sin(hauteur_solaire) flux_reflechi_direct = albedo * direct_horizontal * cosinus_angle_incidence * (1 - math.cos(pi/2))/2 flux_reflechi_diffus = albedo * diffus_horizontal * (1 - math.cos(pi/2))/2 # geometrie du cylindre rayon_cylindre = 0.17 hauteur_cylindre = 1.73 # calcul des flux recus par le cylindre # > direct = anisotrope (calcul integral) flux_direct_cylindre = 2 * flux_direct * rayon_cylindre * hauteur_cylindre # > diffus = anisotrope circumsolaire (calcul integral) + isotrope diffus et horizon flux_diffus_cylindre = (2 * flux_diffus_circumsolaire + pi * flux_diffus_voute_horizon) * rayon_cylindre * hauteur_cylindre # > reflechi = anisotrope (calcul integral) + isotrope flux_reflechi_cylindre = (2 * flux_reflechi_direct + 2 * pi * flux_reflechi_diffus) * rayon_cylindre * hauteur_cylindre #Valeurs en Watts return(flux_direct_cylindre, flux_diffus_cylindre, flux_reflechi_cylindre) def temperature_mrt_fluxsolaire(temperature_mrt_classique, diffus_horizontal, direct_horizontal, albedo, jour, heure, fuseau, latitude, longitude_est): sigma = 5.67 * 0.00000001 #constante de stefan boltzmann epsilon = 0.97 #emissivite du corps humain #temperature radiante en Kelvin temperature_radiante = math.pow(temperature_mrt_classique + 273.15, 4) f_eff = 0.75 #facteur de surface de rayonnement effectif pour un individu debout d'après JENDRITZKY and NUBLER alpha_clo = 0.8 #coefficient d'absorption des vêtements # geometrie cylindre rayon_cylindre = 0.17 hauteur_cylindre = 1.73 S_cylindre = 2*pi*rayon_cylindre*hauteur_cylindre #flux direct, diffus, reflechi avec diffus anisotrope en W/m² flux_direct, flux_diffus, flux_reflechi = flux_cylindre_anisotrope(diffus_horizontal, direct_horizontal, albedo, jour, heure, fuseau, latitude, longitude_est) / S_cylindre temperature_mrt_solaire = math.pow(temperature_radiante + (alpha_clo * flux_direct + f_eff * alpha_clo * flux_diffus + f_eff * alpha_clo * flux_reflechi)/(sigma * epsilon), 0.25) - 273.15 return temperature_mrt_solaire
Scilab intègre un algorithme génétique dans la fonction optim_GA.
Pour gagner du temps lors de l’évaluation de chaque individu d’une génération, il est possible de paralléliser les évaluations : celles-ci sont indépendantes.
Nous vous proposons ici une modification de la routine Scilab optim_ga
permettant d’évaluer en parallèle un nombre d’individus égal au nombre de cœurs de l’ordinateur (la réduction du temps de calcul est ainsi proportionnelle au nombre de cœurs de votre machine).
Ce fonctionnement est valide sous Linux uniquement, la version Windows de Scilab 6 ne supportant pas le parallélisme : « In this current version of Scilab, parallel_run
uses only one core on Windows platforms. » (source).
///////////////////////////////////////////////////////////////////// //_______ ______ _______ _______ //| _ || _ | | || | //| |_| || | || | ___|| _ | //| || |_||_ | |___ | |_| | //| || __ || ___|| ___| //| _ || | | || |___ | | //|__| |__||___| |_||_______||___| // AREP - 16, av. d'Ivry, 7501/3 Paris, FRANCE ///////////////////////////////////////////////////////////////////// // Scilab (www.scilab.org) - This file is part of Scilab // Copyright (C) 2008 - Yann COLLETTE <yann[dot]collette[at]renault[dot]com> // Copyright (C) 2014 - Michael Baudin <michael[dot]baudin[at]contrib[dot]scilab[dot]org> // // 26/6/2017 : added the "parallel_run" capability for the evaluation of individuals (Edouard Walther) // // This file must be used under the terms of the CeCILL. // This source file is licensed as described in the file COPYING, which // you should have received as part of this distribution. The terms // are also available at // http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt // //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// // Last modification : 26/06/2017 // //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// // Contact : edouard[dot]walther[at]arep[dot]com //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// function [pop_opt, fobj_pop_opt, pop_init, fobj_pop_init] = optim_GA_parallel(ga_f, pop_size, nb_generation, p_mut, p_cross, Log, param) [nargout, nargin] = argn(); if ~isdef("param", "local") then param = []; end [codage_func, err] = get_param(param, "codage_func", coding_ga_identity); [init_func, err] = get_param(param, "init_func", init_ga_default); [crossover_func, err] = get_param(param, "crossover_func", crossover_ga_default); [mutation_func, err] = get_param(param, "mutation_func", mutation_ga_default); [selection_func, err] = get_param(param, "selection_func", selection_ga_elitist); [nb_couples, err] = get_param(param, "nb_couples", 100); [pressure, err] = get_param(param, "pressure", 0.05); [output_func, err] = get_param(param, "output_func", output_ga_default); if ~isdef("ga_f", "local") then error(sprintf(gettext("%s: ga_f is mandatory"), "optim_ga")); else if typeof(ga_f) == "list" then deff("y = _ga_f(x)", "y = ga_f(1)(x, ga_f(2:$))"); else deff("y = _ga_f(x)", "y = ga_f(x)"); end end if ~isdef("pop_size", "local") then pop_size = 100; end if ~isdef("nb_generation", "local") then nb_generation = 10; end if ~isdef("p_mut", "local") then p_mut = 0.1; end if ~isdef("p_cross", "local") then p_cross = 0.7; end if ~isdef("Log", "local") then Log = %F; end // Initialization of the population Pop = list(); Pop = init_func(pop_size, param); if (nargout >= 3) then pop_init = Pop; end // Code the individuals Pop = codage_func(Pop, "code", param); // Getting the objective function for each individual //disp("Extraction list..."); [vec_param_pop]=Pop(1); vec_param_total=vec_param_pop; for i = 2:length(Pop) [vec_param_pop]=Pop(i); vec_param_total=cat(2,vec_param_total,vec_param_pop); end // First launch //disp("Execution for the first population..."); [FObj_Pop]=parallel_run(vec_param_total,_ga_f); FObj_Pop=FObj_Pop'; if (nargout == 4) then fobj_pop_init = FObj_Pop; end FObj_Pop_Max = max(FObj_Pop); FObj_Pop_Min = min(FObj_Pop); // Normalization of the efficiency Efficiency = (1 - pressure) * (FObj_Pop_Max - FObj_Pop) / max([FObj_Pop_Max - FObj_Pop_Min %eps]) + pressure; //disp("Starting optimisation with GA..."); // The genetic algorithm for i = 1:nb_generation // // Selection // Indiv1 = list(); Indiv2 = list(); Wheel = cumsum(Efficiency); for j = 1:nb_couples // Selection of the first individual in the couple Shoot = grand(1, 1, "unf", 0, Wheel($)); Index = find(Shoot <= Wheel, 1); Indiv1(j) = Pop(Index); FObj_Indiv1(j) = FObj_Pop(Index); // Selection of the second individual in the couple Shoot = grand(1, 1, "unf", 0, Wheel($)); Index = 1; Index = find(Shoot <= Wheel, 1); Indiv2(j) = Pop(Index); FObj_Indiv2(j) = FObj_Pop(Index); end // // Crossover // for j = 1:nb_couples if (p_cross>grand(1, 1, "def")) then [x1, x2] = crossover_func(Indiv1(j), Indiv2(j), param); Indiv1(j) = x1; Indiv2(j) = x2; ToCompute_I1(j) = %T; ToCompute_I2(j) = %T; else ToCompute_I1(j) = %F; ToCompute_I2(j) = %F; end end // // Mutation // for j = 1:nb_couples if (p_mut>grand(1, 1, "def")) then x1 = mutation_func(Indiv1(j), param); Indiv1(j) = x1; ToCompute_I1(j) = %T; end if (p_mut>grand(1, 1, "def")) then x2 = mutation_func(Indiv2(j), param); Indiv2(j) = x2; ToCompute_I2(j) = %T; end end // // Computation of the objective functions k=0;kk=0; // counters to iterate for j = 1:nb_couples // for all couples in the population if ToCompute_I1(j) then// if to be computed k=k+1; if k==1 then // create the first vector of parameters [vec_param_pop1]=Indiv1(j); else // concatenate for parallel_run [vec_param_indiv1]=Indiv1(j); indices_indiv1(k)=j; vec_param_pop1=cat(2,vec_param_pop1,vec_param_indiv1); end end if ToCompute_I2(j) then// if to be computed kk=kk+1; if kk==1 then [vec_param_pop2]=Indiv2(j); else [vec_param_indiv2]=Indiv2(j); indices_indiv2(kk)=j; vec_param_pop2=cat(2,vec_param_pop2,vec_param_indiv2); end end end // Parallel_run //disp("Parallel launch for Indiv1..."); [objectifs_Indiv1]=parallel_run(vec_param_pop1,_ga_f); objectifs_Indiv1=objectifs_Indiv1'; //disp("Parallel launch for Indiv2..."); [objectifs_Indiv2]=parallel_run(vec_param_pop2,_ga_f); objectifs_Indiv2=objectifs_Indiv2'; // Updating indexes //disp("Updating FObj1 ..."); for k=1:length(objectifs_Indiv1) if indices_indiv1(k)<> 0 then FObj_Indiv1(indices_indiv1(k))= objectifs_Indiv1(k); end; end for k=1:length(objectifs_Indiv2) if indices_indiv2(k)<> 0 then FObj_Indiv2(indices_indiv2(k))= objectifs_Indiv2(k); end; end // Reinit ToCompute lists ToCompute_I1 = ToCompute_I1 & %F; ToCompute_I2 = ToCompute_I2 & %F; // Recombination [Pop, FObj_Pop] = selection_func(Pop, Indiv1, Indiv2, FObj_Pop, FObj_Indiv1, FObj_Indiv2, [], [], [], param); // Callback for plotting / printing intermediate results or stopping the algorithm if (Log) then stop = output_func(i, nb_generation, Pop, FObj_Pop, param); if (stop) then break end end end pop_opt = Pop; pop_opt = codage_func(pop_opt, "decode", param); fobj_pop_opt = FObj_Pop; endfunction
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
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))