Skip to content

Update PET_fsolve_code.py #2

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
kastnerp wants to merge 1 commit into
base: master
Choose a base branch
from
Open

Update PET_fsolve_code.py #2

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
248 changes: 137 additions & 111 deletions PET_fsolve_code.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12,107 +12,118 @@
import scipy.optimize as optimize

# 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
tc_set = 36.6 # 36.8
tsk_set = 34 # 33.7
# Calculation of the body temperature through a weighted average
tbody_set = 0.1*tsk_set+0.9*tc_set


# Skin blood flow calculation function:
def vasoC(tcore,tsk):
#Set value signals


def vasoC(tcore, tsk):
# Set value signals
sig_skin = tsk_set - tsk
sig_core = tcore - tc_set
if sig_core<0:
if sig_core < 0:
# In this case, Tcore<Tc_set --> the blood flow is reduced
sig_core=0.
if sig_skin<0:
sig_core = 0.
if sig_skin < 0:
# In this case, Tsk>Tsk_set --> the blood flow is increased
sig_skin=0.
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
if qmblood>90:
qmblood=90.
if qmblood > 90:
qmblood = 90.
# in the transient model, alpha is used to update tbody
#alpha = 0.04177 + 0.74518 / (qmblood + 0.585417)
alpha = 0.1
return (qmblood,alpha)
return (qmblood, alpha)


# Sweating calculation function
def Suda(tbody,tsk):
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 sweat flow is 0
sig_body=0.
if sig_skin<0:
if sig_body < 0:
# In this case, Tbody<Tbody_set --> The sweat flow is 0
sig_body = 0.
if sig_skin < 0:
# In this case, Tsk<Tsk_set --> the sweat flow is reduced
sig_skin=0.
#qmsw = 170 * sig_body * math.exp((sig_skin) / 10.7) # [g/m2/h] is the expression from Gagge's model
sig_skin = 0.
# qmsw = 170 * sig_body * math.exp((sig_skin) / 10.7) # [g/m2/h] is the expression from Gagge's model
qmsw = 304.94*10**(-3) * sig_body
# 500 g/m^2/h is the upper sweat rate limit
if qmsw > 500:
qmsw = 500
return (qmsw)

# Vectorial MEMI balance calculation function
def Syst(T, Ta, Tmrt, HR, v, age, sex, ht, mbody, pos, M, icl,mode):


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)) #required for the vectorial expression of the balance
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
# required for the vectorial expression of the balance
enbal_vec = np.zeros((3, 1))

# Area parameters of the body:
Adu = 0.203 * mbody ** 0.425 * ht ** 0.725
feff=0.725
feff = 0.725
if pos == 1 or pos == 3:
feff = 0.725
if pos == 2:
feff = 0.696
# Calculation of the Burton surface increase coefficient, k = 0.31 for Hoeppe:
fcl = 1 + (0.31 * icl) # Increase heat exchange surface depending on clothing level
facl = (173.51 * icl - 2.36 - 100.76 * icl * icl + 19.28 * icl ** 3.0) / 100
# Increase heat exchange surface depending on clothing level
fcl = 1 + (0.31 * icl)
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: # mode=True is the calculation of the actual environment
vpa = HR / 100.0 * 6.105 * math.exp(17.27 * Ta / (237.7 + Ta )) #[hPa]
else: # mode=False means we are calculating the PET
vpa= 12 # [hPa] vapour pressure of the standard environment

# Partial pressure of water in the air depending on relative humidity and air temperature:
if mode: # mode=True is the calculation of the actual environment
vpa = HR / 100.0 * 6.105 * math.exp(17.27 * Ta / (237.7 + Ta)) # [hPa]
else: # mode=False means 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)
hc = 2.67 + (6.5 * v**0.67)
if pos == 2:
hc = 2.26 + (7.42 *v**0.67)
hc = 2.26 + (7.42 * v**0.67)
if pos == 3:
hc = 8.6 * (v ** 0.513)
# modification of hc with the total pressure
hc = hc * (p / po) ** 0.55
hc = hc * (p / po) ** 0.55

# Base 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))
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
if mode==True: # = actual environment
metab = (M + metab_male)/Adu
fec = (M + metab_female)/Adu
else:# False=reference environment
if mode == True: # = actual environment
metab = (M + metab_male)/Adu
fec = (M + metab_female)/Adu
else: # False=reference environment
metab = (80 + metab_male)/Adu
fec = (80 + metab_female)/Adu

he = 0.0
# Attribution of internal energy depending on the sex of the subject
if sex == 1:
he = metab
elif sex == 2:
he = fec
h = he *(1.0 - eta) # [W/m2]
h = he * (1.0 - eta) # [W/m2]

# Respiratory energy losses
# Expired air temperature calculation:
Expand Down Expand Up @@ -141,119 +152,134 @@ def Syst(T, Ta, Tmrt, HR, v, age, sex, ht, mbody, pos, M, icl,mode):
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
r1 = facl * Adu / (6.28 * ht * y) # Internal radius
di = r2 - r1
# Calculation of the equivalent thermal resistance of body tissues
alpha = vasoC(T[0,0],T[1,0])[1]
tbody = alpha * T[1,0] + (1 - alpha) * T[0,0]
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])
qmsw = Suda(tbody, T[1, 0])
# Lvap/1000 = 2400 000[J/kg] divided by 1000 = [J/g] // qwsw/3600 for [g/m2/h] to [g/m2/s]
esw = Lvap/1000* qmsw/3600 # [W/m2]
esw = Lvap/1000 * qmsw/3600 # [W/m2]
# Saturation vapor pressure at temperature Tsk
Pvsk = 6.105*math.exp((17.27 * (T[1,0]+273.15) - 4717.03)/ (237.7 + T[1,0])) # [hPa]
# Calculation of vapour transfer
Lw = 16.7*10 **(-1) # [K/hPa] Lewis factor
he_diff = hc * Lw # diffusion coefficient of air layer
fecl=1/(1+0.92*hc*rcl) # Burton efficiency factor
emax = he_diff * fecl * (Pvsk - vpa) # maximum diffusion at skin surface
Pvsk = 6.105 * \
math.exp((17.27 * (T[1, 0]+273.15) - 4717.03) /
(237.7 + T[1, 0])) # [hPa]
# Calculation of vapour transfer
Lw = 16.7*10 ** (-1) # [K/hPa] Lewis factor
he_diff = hc * Lw # diffusion coefficient of air layer
fecl = 1/(1+0.92*hc*rcl) # Burton efficiency factor
emax = he_diff * fecl * (Pvsk - vpa) # maximum diffusion at skin surface
w = esw / emax # skin wettedness
if w > 1:
w=1
w = 1
delta = esw-emax
if delta < 0:
esw=emax
esw = emax
if esw < 0:
esw=0
i_m=0.38 # Woodcock's ratio
R_ecl=(1/(fcl*hc) + rcl)/(Lw*i_m) # clothing vapour transfer resistance after Woodcock's method
#R_ecl=0.79*1e7 # Hoeppe's method for E_diff
esw = 0
i_m = 0.38 # Woodcock's ratio
# clothing vapour transfer resistance after Woodcock's method
R_ecl = (1/(fcl*hc) + rcl)/(Lw*i_m)
# R_ecl=0.79*1e7 # Hoeppe's method for E_diff
ediff = (1 - w)*(Pvsk - vpa)/R_ecl # diffusion heat transfer
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
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
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]
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
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

#returning either the calculated core,skin,clo temperatures or the PET
# returning either the calculated core,skin,clo temperatures or the PET
if mode:
# if we solve for the system we need to return 3 temperatures
return [enbal_vec[0,0],enbal_vec[1,0],enbal_vec[2,0]]
# if we solve for the system we need to return 3 temperatures
return [enbal_vec[0, 0], enbal_vec[1, 0], enbal_vec[2, 0]]
else:
# solving for the PET requires the scalar balance only
# solving for the PET requires the scalar balance only
return enbal_scal

# Solving the 3 equation non-linear 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))
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):
# 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 = Tmin # Start of the search interval
Tf = Tmax # End of the search interval
return Syst(Tstable, Tx, Tx, 50, 0.1, age, sex, ht, mbody, pos, M, 0.9, False)
Ti = Tmin # Start of the search interval
Tf = Tmax # End of the search interval
pet = 0
while Tf-Ti>eps: # Dichotomy loop
if f(Ti)*f(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
# definition of 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 # Clothing 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
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 # Clothing 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**(-6) # numerical tolerance
T = [38, 40, 40]
eps = 10**(-6) # numerical tolerance
# Dichotomy search interval (a=min / b=max)
Tmin = -40
Tmax = 60

# Input data for the PET
Ta=30 # Air temperature in [oC]
Tmrt=60 # Mean radiant temperature in [oC]
HR=50 # Air relative humidity [%]
v=1 # Wind velocity [m/s]
# Input data for the PET
Ta = 30 # Air temperature in [oC]
Tmrt = 60 # Mean radiant temperature in [oC]
HR = 50 # Air relative humidity [%]
v = 1 # Wind velocity [m/s]
age = 35
sex = 1 # 1 for men and 2 for women
sex = 1 # 1 for men and 2 for women
pos = 1
mbody = 75 #[kg]
ht = 1.80 #[m]
p = 1013.25 #[hPa]
M = 80 # [W] Metabolic activity level
icl = 0.5 # [clo] Clothing level

# Results
Tstable = resolution(Ta,Tmrt,HR,v,age,sex,ht,mbody,pos,M,icl,T)[0]
print("Nodes temperature [T_core, T_skin, T_clo]",Tstable)
print('Thermal Balance', Syst(Tstable, Ta, Tmrt, HR, v, age, sex, ht, mbody, pos, M, icl,True)[0])
print('PET:', round(PET(age, sex, ht, mbody, pos, M, icl, Tstable, Tmin, Tmax, eps),2))
mbody = 75 # [kg]
ht = 1.80 # [m]
p = 1013.25 # [hPa]
M = 80 # [W] Metabolic activity level
icl = 0.5 # [clo] Clothing level

# Results
Tstable = resolution(Ta, Tmrt, HR, v, age, sex, ht, mbody, pos, M, icl, T)[0]
print("Nodes temperature [T_core, T_skin, T_clo]", Tstable)
print('Thermal Balance', Syst(Tstable, Ta, Tmrt, HR,
v, age, sex, ht, mbody, pos, M, icl, True)[0])
print('PET:', round(PET(age, sex, ht, mbody, pos, M, icl, Tstable, Tmin, Tmax, eps), 2))