Option Explicit


Public Function coefficient_j(ByVal vitesse As Single, ByVal diametre As Single, _
ByVal rho As Single, Optional epsilon As Single = 0.0001, Optional nu As Single = 0.000016) As Single 'definition des variables pour éviter les erreurs

Dim re As Single,compteur  As Integer, lambda0 As Single, lambda1 As Single, lambda2 As Single
'rugosités typiques en mm
'cuivre , plomb, laiton, inox - 0,001 à 0,002
'tube PVC - 0,0015
'Acier inox - 0,015
'tube acier du commerce - 0,045 à 0,09
'Acier étiré - 0,015
'Acier soudé - 0,045
'Acier galvanisé - 0,15
'Acier rouillé - 0,1 à 1
'fonte neuve - 0,25 à 0,8
'fonte usagée - 0,8 à 1,5
'fonte incrustée - 1,5 à 2,5
'tôle ou fonte asphaltée - 0,01 à 0,015
'ciment bien lissé - 0,3
'béton ordinaire - 1
'béton grossier - 5
'bois bien raboté - 5
'bois ordinaire - 1

re = vitesse * diametre * 0.001 / nu
lambda0 = 64/re
lambda1 = 2*lambda0
compteur = 0
Do While Abs(lambda1 - lambda0) / lambda0 > 0.02 Or compteur < 100
lambda2 = (-2 * Log(epsilon / (0.001 * 3.7 * diametre) + 2.51 / (re * lambda0 ^ 0.5)) / Log(10)) ^ -2
lambda1 = lambda0
lambda0 = lambda2
compteur = compteur+1
Loop
If compteur > 100 Then
coefficient_j = "pas de solution trouvée"
Else
coefficient_j = 1000 * lambda0 / diametre * rho * vitesse ^ 2 / 2
End If
End Function

Public Function diametre_comm(ByVal diametre As Single, ByVal seuil As Single)
Dim serie1 As Variant, serie2 As Variant, diametre_comm1 As Integer, diametre_comm2 As Integer, index1 As Integer, index2 As Integer
If diametre < 125 Then
diametre_comm = 125
ElseIf diametre > 710 Then
diametre_comm = -1
    Else
    diametre = diametre + 0.1 'pour éviter de tomber pile-poil sur un diamètre commercial existant (génère des erreurs)
    serie1 = Array(125, 160, 200, 250, 315, 355, 400, 450, 500, 560, 630, 710) 'série des diamètres commerciaux
    serie2 = Array(710, 630, 560, 500, 450, 400, 355, 315, 250, 200, 160, 125)
    index1 = Application.WorksheetFunction.Match(diametre, serie1, 1)  'utilisation des fonctions d'excel
    index2 = Application.WorksheetFunction.Match(diametre, serie2, -1)
    diametre_comm1 = serie1(index1)
    diametre_comm2 = serie2(index2)
    seuil = 1 - seuil
    If Abs(diametre_comm1 - diametre) / Abs(diametre_comm1 - diametre_comm2) <= seuil Then
        diametre_comm = diametre_comm1
    Else
        diametre_comm = diametre_comm2
    End If
End If
End Function

Public Function puiss_emetteur(ByVal DT As Single, ByVal KS As Single _
, ByVal n As Single, Optional Q As Single = 1, Optional m As Single = 0)
puiss_emetteur = KS * Q ^ m * DT ^ n
End Function

Public Function rho_air(ByVal temperature_celsius As Single, Optional pression As Single = 101325)
rho_air = pression * 0.029 / ((273.15 + temperature_celsius) * 8.314)
End Function

Public Function pression_vap_sat(ByVal temperature As Single)
If temperature < -20 Or temperature > 35 Then
pression_vap_sat = -1
ElseIf temperature >= 0 And temperature <= 35 Then
pression_vap_sat = 288.68 * (1.098 + (temperature / 100)) ^ 8.02
ElseIf temperature < 0 And temperature >= -20 Then
pression_vap_sat = 4.689 * (1.486 + temperature / 100) ^ 12.3
End If

End Function
Public Function pression_vap(ByVal temperature As Single, ByVal humidite_relative As Single)
pression_vap = humidite_relative * pression_vap_sat(temperature)
End Function
Public Function t_rosee(ByVal temperature As Single, ByVal humidite_relative As Single)
t_rosee = (Exp(Log(pression_vap(temperature, humidite_relative) / 288.68) / 8.02) - 1.098) * 100
End Function
Public Function hum_abs(ByVal temperature As Single, ByVal humidite_relative As Single, Optional pression As Single = 101325)
hum_abs = 622 * pression_vap(temperature, humidite_relative) / (pression - pression_vap(temperature, humidite_relative))
End Function
Public Function enthalpie(ByVal temperature As Single, ByVal humidite_relative As Single)
enthalpie = 1.006 * temperature + hum_abs(temperature, humidite_relative) / 1000 * (2450 + 1.86 * temperature)
End Function
Public Function vol_spe(ByVal temperature As Single, ByVal humidite_relative As Single)
vol_spe = 461.5 * (temperature + 273) * (hum_abs(temperature, humidite_relative) / 1000 + 0.622) / 101300
End Function
Public Function hum_rel(ByVal temperature As Single, ByVal humid_abs As Single, Optional pression As Single = 101325)
hum_rel = pression * humid_abs / (622 * pression_vap_sat(temperature))
End Function

Public Function pression_vap_sat2(ByVal humid_rel As Single, ByVal humid_abs As Single, Optional pression As Single = 101325)
pression_vap_sat2 = pression * humid_abs / (622 * humid_rel)
End Function

Public Function temper(ByVal press_vap_sat As Single)
'température de rosée
temper = ((press_vap_sat / 288.68) ^ (1 / 8.02) - 1.098) * 100
If temper > 35 Then
    temper = "erreur"
End If
If temper < 0 Then
temper = 100 * ((press_vap_sat / 4.689) ^ (1 / 12.3) - 1.486)
    If temper < -20 Then
        temper = "erreur"
    End If
End If
End Function

Public Function temper2(ByVal enth As Single, ByVal humid_abs As Single)
'température d'air
    temper2 = (enth - 2450 * humid_abs) / (1.006 + 1.86 * humid_abs)
End Function
Public Function temp_bulb_hum(ByVal T As Single, ByVal hum_rel As Single)
Dim ecart As Single, T2 As Single, enthalp As Single, enthalp2 As Single, compteur As Integer
enthalp = enthalpie(T, hum_rel)
compteur = 0
T2 = T
ecart = 1
Do While ecart >= 0.1
enthalp2 = enthalpie(T2, 1)
ecart = Abs(enthalp2 - enthalp)
T2 = T2 - 0.05
If compteur > 100 Then
ecart = 0
End If
compteur = compteur + 1
Loop
temp_bulb_hum = T2
End Function

Public Function rho_eau(ByVal T As Single)
'masse volumique de l'eau
If T > 90 Then
    rho_eau = -1
ElseIf T < 0 Then
    rho_eau = -1
Else
    rho_eau = -0.0038 * T ^ 2 - 0.0531 * T + 1000.3
End If
End Function
Public Function lambda_eau(ByVal T As Single)
'conductivité thermique de l'eau
If T > 90 Then
    lambda_eau = -1
ElseIf T < 0 Then
    lambda_eau = -1
Else
    lambda_eau = -0.0000101515 * T ^ 2 + 0.002242121 * T + 0.554836364
End If
End Function
Public Function Cp_eau(ByVal T As Single)
'chaleur massique de l'eau
If T > 90 Then
    Cp_eau = -1
ElseIf T < 0 Then
    Cp_eau = -1
Else
    Cp_eau = 0.0000043852 * T ^ 4 - 0.000998155 * T ^ 3 + 0.088200758 * T ^ 2 - 3.181915307 * T + 4216.402098
End If
End Function
Public Function nu_eau(ByVal T As Single)
'viscosité de l'eau
If T > 90 Then
    nu_eau = -1
ElseIf T < 0 Then
    nu_eau = -1
Else
    nu_eau = 3.99767E-14 * T ^ 4 - 1.04207E-11 * T ^ 3 + 0.0000000010636 * T ^ 2 - 0.0000000567069 * T + 0.00000178708
End If
End Function
Public Function beta_eau(ByVal T As Single)
If T > 90 Then
    beta_eau = -1
ElseIf T < 0 Then
    beta_eau = -1
Else
    beta_eau = -2.64423E-11 * T ^ 4 + 0.00000000534334 * T ^ 3 - 0.000000384923 * T ^ 2 + 0.0000198047 * T - 0.0000722692
End If
End Function

Public Function lambda_air(ByVal T As Single)
'conductivité thermique de l'air
If T > 90 Then
    lambda_air = -1
ElseIf T < 0 Then
    lambda_air = -1
Else
    lambda_air = 0.0000707143 * T + 0.024280952
End If
End Function
Public Function Cp_air(ByVal T As Single)
'chaleur massique de l'air
If T > 90 Then
    Cp_air = -1
ElseIf T < 0 Then
    Cp_air = -1
Else
    Cp_air = 1009
End If
End Function
Public Function nu_air(ByVal T As Single)
'viscosité de l'air
If T > 90 Then
    nu_air = -1
ElseIf T < 0 Then
    nu_air = -1
Else
    nu_air = 0.0000000975714 * T + 0.0000131714
End If
End Function
Public Function beta_air(ByVal T As Single)
If T > 90 Then
    beta_air = -1
ElseIf T < 0 Then
    beta_air = -1
Else
    beta_air = -0.00000992857 * T + 0.003631429

End If
End Function

Public Function hcv_nat(ByVal T_pl As Single, ByVal T_inf As Single, ByVal D As Single, ByVal fluide As Integer)
'fluide 1 = air et fluide 2 = eau
Dim pr As Single, gr As Single, nu As Single, lambda As Single
If fluide = 1 Then
pr = rho_air(T_pl) * nu_air(T_pl) * Cp_air(T_pl) / lambda_air(T_pl)
gr = 9.81 * beta_air(T_pl) * (T_pl - T_inf) * D ^ 3 / nu_air(T_pl) ^ 2
lambda = lambda_air(T_pl)
Else
pr = rho_eau(T_pl) * nu_eau(T_pl) * Cp_eau(T_pl) / lambda_eau(T_pl)
gr = 9.81 * beta_eau(T_pl) * (T_pl - T_inf) * D ^ 3 / nu_eau(T_pl) ^ 2
lambda = lambda_eau(T_pl)
End If
If gr * pr < 0.05 Then
nu = 1.18 * (pr * gr) ^ (1 / 8)
ElseIf gr * pr < 2 * 10 ^ 7 Then
nu = 0.54 * (pr * gr) ^ (1 / 4)
ElseIf gr * pr < 10 ^ 13 Then
nu = 0.13 * (pr * gr) ^ (1 / 3)
Else
nu = -1
End If
If nu = -1 Then
hcv_nat = -1
Else
hcv_nat = lambda * nu / D
End If
End Function

Public Function T_bulb_humid(ByVal enthalp As Single)
'Cette fonction calcule la température de bulbe humide suivant une corrélation (2,5% d'erreur max ; T de -15°C à 30°C et Patm=101325 Pa)
T_bulb_humid = -0.0023 * enthalp ^ 2 + 0.5955 * enthalp - 5.9859
End Function

'------------------------------------------------------------------------------------------------------
' Predicted Mean Vote (PMV) and Predicted Percentage of Dissatisfied (PPD) in accordance with ISO 7730
'---------   December 18 2002, editted by Takahiro SATO, Tanabe Lab., Waseda Univ. ---------------
'------------------------------------------------------------------------------------------------------
Function FNPS(T)
    FNPS = Exp(16.6536 - 4030.183 / (T + 235))         'Saturated Vapor Pressure, [kPa]
End Function

Function PMV(ByVal Ta As Single, ByVal Tr As Single, ByVal Vel As Single, ByVal RH As Single, ByVal CLO As Single, ByVal MET As Single, Optional EW As Single = 0)   'Definition of the Function "PMV" by 7 factors
           
'  Ta  : Air Temperature,                [deg.C]
'@Tr@: Mean Radiant Temperature, @    [deg.C]
'@Vel : Relative Air Velocity,          [m/s]
'  RH  : Relative Humidity,              [%]
'  CLO : Clothing,                       [clo]
'  MET : Metabolic Rate,                 [met]
 
'  EW : External Work,                   [met] (=normally around 0)
'  PA  : Water Vapor Pressure,          [Pa]
   PA = RH * 10 * FNPS(Ta)   '[Pa]=(RH/100)*1000*[kPa]

'---METABORIC RATE---
 m = MET * 58.15:    'Metabolic Rate,    [W/m2]
 W = EW * 58.15:     'External Work,     [W/m2]
 MW = m - W          'internal heat production in the human body
  
'---CLOTHING---
 ICL = 0.155 * CLO:  'thermal insulation of the Clothing, [m2K/W]
 If ICL < 0.078 Then FCL = 1 + 1.29 * ICL Else FCL = 1.05 + 0.645 * ICL  'clothing area factor

'---CONVECTION---
 HCF = 12.1 * Sqr(Vel):  'convective heat transfer coefficient by forced convection
 TaA = Ta + 273:         'Air Temperature in Kelvin [K]
 TrA = Tr + 273:         'Mean Radiant Temperature in Kelvin [K]
    
'CALCULATE SURFACE TEMPERATURE OF CLOTHING BY ITERATION
 TCLA = TaA + (35.5 - Ta) / (3.5 * (6.45 * ICL + 0.1))
'first guess for surface temperature of clothing
 P1 = ICL * FCL:  'calculation term
 P2 = P1 * 3.96:  'calculation term
 P3 = P1 * 100:   'calculation term
 P4 = P1 * TaA:   'calculation term
 P5 = 308.7 - 0.028 * MW + P2 * (TrA / 100) ^ 4 'calculation term
 XN = TCLA / 100
 XF = XN
 n = 0:           'N: number of iterations
 EPS = 0.00015:   'stop criteria in iteration

Do
XF = (XF + XN) / 2

'convective heat Transf. coeff. by natural convection
 HCN = 2.38 * Abs(100 * XF - TaA) ^ 0.25
 If HCF > HCN Then HC = HCF Else HC = HCN
 XN = (P5 + P4 * HC - P2 * XF ^ 4) / (100 + P3 * HC)
 n = n + 1
 
 If n > 150 Then GoTo 50
Loop Until Abs(XN - XF) < EPS
  
 TCL = 100 * XN - 273:  'surface temperature of the clothing
    
'---HEAT LOSS COMPONENTS---
'heat loss diff. through skin
 Ediff = 3.05 * 0.001 * (5733 - 6.99 * MW - PA)
'heat loss by sweating (comfort)
 If MW > 58.15 Then Esw = 0.42 * (MW - 58.15) Else Esw = 0!
'latent respiration heat loss
 LRES = 1.7 * 0.00001 * m * (5867 - PA)
'dry respiration heat loss
 DRES = 0.0014 * m * (34 - Ta)
'heat loss by radiation
 R = 3.96 * FCL * (XN ^ 4 - (TrA / 100) ^ 4)
'heat loss by convection
 C = FCL * HC * (TCL - Ta)
    
'--- CALCULATE PMV AND PPD ---
'Thermal sensation transfer coefficient
 TS = 0.303 * Exp(-0.036 * m) + 0.028
    
'Predicted Mean Vote
 PMV = TS * (MW - Ediff - Esw - LRES - DRES - R - C)

GoTo 80
50 PMV = 999999!
80 '
End Function

'Predicted PercenTage of Dissatisfied
Function PPD(ByVal PMV As Single)
    PPD = 100 - 95 * Exp(-0.03353 * PMV ^ 4 - 0.2197 * PMV ^ 2)
End Function