import math
# ============================================================
# 1. CONSTANTES DU PROBLEME
# ============================================================
# Constante spécifique de l'air
R_AIR = 287.04
# [J/(kg.K)]
# ------------------------------------------------------------
# Loi de Sutherland
# ------------------------------------------------------------
MU_REF = 18.2e-6
# viscosité dynamique de référence [Pa.s]
T_REF = 293.15
# température de référence = 20 °C [K]
SUTHERLAND_C = 113.0
# constante de Sutherland [K]
# ------------------------------------------------------------
# Matériaux
# ------------------------------------------------------------
MATERIAU_DIAPHRAGME = "Acier inoxydable"
MATERIAU_CONDUITE = "Acier au carbone"
# ------------------------------------------------------------
# Paramètres numériques de convergence
# ------------------------------------------------------------
TOL_C = 1.0e-8
# tolérance sur le coefficient de décharge C
TOL_T = 1.0e-6
# tolérance sur la température statique T [K]
MAX_ITER_C = 100
MAX_ITER_T = 100
# ============================================================
# 2. GEOMETRIE FIXE DE L'INSTALLATION
# ============================================================
d0 = 0.060
# diamètre de référence du diaphragme [m]
D0 = 0.100
# diamètre intérieur de référence de la conduite [m]
# ============================================================
# 3. GRANDEURS MESUREES PENDANT L'ESSAI
# ============================================================
p = 5.00e5
# pression statique absolue en amont du diaphragme [Pa]
delta_p = 1.00e4
# différence de pression statique amont - aval [Pa]
Ti = 300.0
# température totale de l'air [K]
Td = 302.0
# température du diaphragme [K]
TD = 301.0
# température de la conduite [K]
# ============================================================
# 4. COEFFICIENTS DE DILATATION THERMIQUE
# ============================================================
def alpha_d(T):
theta = T - 273.15
# température en °C
if theta < 20.0:
return 16.0e-6
elif theta <= 100.0:
return (
16.0e-6
+ (16.5e-6 - 16.0e-6)
* (theta - 20.0)
/ (100.0 - 20.0)
)
elif theta <= 200.0:
return (
16.5e-6
+ (17.0e-6 - 16.5e-6)
* (theta - 100.0)
/ (200.0 - 100.0)
)
elif theta <= 300.0:
return (
17.0e-6
+ (17.5e-6 - 17.0e-6)
* (theta - 200.0)
/ (300.0 - 200.0)
)
else:
return 17.5e-6
def alpha_D(T):
theta = T - 273.15
# température en °C
if theta < 20.0:
return 12.0e-6
elif theta <= 100.0:
return (
12.0e-6
+ (12.4e-6 - 12.0e-6)
* (theta - 20.0)
/ (100.0 - 20.0)
)
elif theta <= 200.0:
return (
12.4e-6
+ (12.8e-6 - 12.4e-6)
* (theta - 100.0)
/ (200.0 - 100.0)
)
elif theta <= 300.0:
return (
12.8e-6
+ (13.2e-6 - 12.8e-6)
* (theta - 200.0)
/ (300.0 - 200.0)
)
else:
return 13.2e-6
# ============================================================
# 5. FONCTIONS PHYSIQUES DE L'AIR
# ============================================================
def viscosite_air(T):
mu = (
MU_REF
* (T / T_REF) ** 1.5
* (SUTHERLAND_C + T_REF)
/ (SUTHERLAND_C + T)
)
return mu
def cp_air(T):
exp_term = math.exp(3090.0 / T)
Cp = R_AIR * (
3.5
- 2.8e-5 * T
+ 2.24e-8 * T**2
+ (3090.0 / T) ** 2
* exp_term
/ (exp_term - 1.0) ** 2
)
return Cp
def facteur_compressibilite_Z(p, T):
p_bar = p / 1.0e5
T_C = T - 273.15
if not (1.0 <= p_bar <= 65.0):
raise ValueError(
"Pression hors du domaine de calcul de Z : 1 à 65 bar."
)
if not (-20.0 <= T_C <= 800.0):
raise ValueError(
"Température hors du domaine de calcul de Z : "
"-20 à 800 °C."
)
# --------------------------------------------------------
# Pression entre 1 et 10 bar
# --------------------------------------------------------
if p_bar <= 10.0:
if T_C <= 50.0:
Z = (
1.0
- 1.0e-7
* p_bar
* (T_C - 77.0) ** 2
)
else:
Z = 1.0
# --------------------------------------------------------
# Pression entre 10 et 65 bar
# --------------------------------------------------------
else:
if T_C <= 45.0:
Z = (
0.9974
+ 5.0e-6
* p_bar
* (T_C - 47.0)
)
elif T_C <= 180.0:
Z = (
0.9974
+ 1.3e-6
* (p_bar - 35.0) ** 2
+ 3.9375e-6
* p_bar
* (T_C - 47.0)
)
else:
Z = (
1.0
+ 4.0e-4
* p_bar
)
return Z
def coefficient_detente(beta, p, delta_p, gamma):
p_aval = p - delta_p
epsilon = (
1.0
- (
0.351
+ 0.256 * beta**4
+ 0.93 * beta**8
)
* (
1.0
- (p_aval / p) ** (1.0 / gamma)
)
)
return epsilon
def coefficient_decharge(beta, Re_D):
A = (
19000.0
* beta
/ Re_D
) ** 0.8
C = (
0.5961
+ 0.0261 * beta**2
- 0.216 * beta**8
+ 0.000521
* (beta * 1.0e6 / Re_D) ** 0.7
+ (
0.0188
+ 0.0063 * A
)
* beta**3.5
* (1.0e6 / Re_D) ** 0.3
)
return C
# ============================================================
# 6. VERIFICATION DE LA MESURE DE PRESSION DIFFERENTIELLE
# ============================================================
if delta_p <= 0.0:
print(
"delta_p <= 0 : débit massique q = 0 kg/s"
)
raise SystemExit
# ============================================================
# 7. CORRECTION DES DIAMETRES PAR DILATATION
# ============================================================
# ------------------------------------------------------------
# Coefficients de dilatation
# ------------------------------------------------------------
alpha_d_val = alpha_d(Td)
alpha_D_val = alpha_D(TD)
# ------------------------------------------------------------
# (A) Diamètre du diaphragme à Td
# ------------------------------------------------------------
d = d0 * (
1.0
+ alpha_d_val
* (Td - T_REF)
)
# ------------------------------------------------------------
# (B) Diamètre de la conduite à TD
# ------------------------------------------------------------
D = D0 * (
1.0
+ alpha_D_val
* (TD - T_REF)
)
# ============================================================
# 8. GRANDEURS GEOMETRIQUES
# ============================================================
# ------------------------------------------------------------
# (C) beta : rapport des diamètres
# ------------------------------------------------------------
beta = d / D
# ------------------------------------------------------------
# (D) s : section de l'orifice du diaphragme
# ------------------------------------------------------------
s = (
math.pi
* d**2
/ 4.0
)
# ------------------------------------------------------------
# (E) S : section intérieure de la conduite
# ------------------------------------------------------------
S = (
math.pi
* D**2
/ 4.0
)
# ------------------------------------------------------------
# (F) E : coefficient d'approche
# ------------------------------------------------------------
E = 1.0 / math.sqrt(
1.0 - beta**4
)
# ============================================================
# 9. INITIALISATION DES ITERATIONS
# ============================================================
# Première estimation du coefficient de décharge
C = 0.60
# T est la température statique inconnue.
#
# Ti est la température totale mesurée.
#
# Pour commencer :
# T = Ti constitue la première estimation de T.
T = Ti
# Erreurs initiales
erreur_T = 1.0
# Compteur de la boucle extérieure
iteration_T = 0
# Cette variable permettra de conserver
# le nombre d'itérations de la dernière boucle sur C.
iteration_C = 0
# ============================================================
# 10. BOUCLE EXTERIEURE : CONVERGENCE DE T
# ============================================================
while erreur_T > TOL_T:
iteration_T += 1
if iteration_T > MAX_ITER_T:
raise RuntimeError(
"La boucle sur T n'a pas convergé."
)
# ========================================================
# PROPRIETES DE L'AIR POUR LA VALEUR ACTUELLE DE T
# ========================================================
# --------------------------------------------------------
# (5) mu : viscosité dynamique
# --------------------------------------------------------
mu = viscosite_air(T)
# --------------------------------------------------------
# (7) Cp : capacité thermique massique
# --------------------------------------------------------
Cp = cp_air(T)
# --------------------------------------------------------
# (6) gamma : rapport Cp/Cv
# --------------------------------------------------------
gamma = (
Cp
/ (Cp - R_AIR)
)
# --------------------------------------------------------
# (3) epsilon : coefficient de détente
# --------------------------------------------------------
epsilon = coefficient_detente(
beta,
p,
delta_p,
gamma
)
# --------------------------------------------------------
# (13) Z : facteur de compressibilité
# --------------------------------------------------------
Z = facteur_compressibilite_Z(
p,
T
)
# --------------------------------------------------------
# (12) rho : masse volumique
# --------------------------------------------------------
rho = (
p
/ (
Z
* R_AIR
* T
)
)
# ========================================================
# 11. BOUCLE INTERIEURE : CONVERGENCE DE C
# ========================================================
erreur_C = 1.0
iteration_C = 0
while erreur_C > TOL_C:
iteration_C += 1
if iteration_C > MAX_ITER_C:
raise RuntimeError(
"La boucle sur C n'a pas convergé."
)
# ----------------------------------------------------
# Ancienne valeur de C
# ----------------------------------------------------
C_ancien = C
# ----------------------------------------------------
# (1) q : débit massique
# ----------------------------------------------------
q = (
C
* epsilon
* E
* s
* math.sqrt(
2.0
* delta_p
* rho
)
)
# ----------------------------------------------------
# (4) Re_D : nombre de Reynolds
# ----------------------------------------------------
Re_D = (
4.0
* q
/ (
math.pi
* D
* mu
)
)
# ----------------------------------------------------
# (2) Nouveau coefficient de décharge
# ----------------------------------------------------
C = coefficient_decharge(
beta,
Re_D
)
# ----------------------------------------------------
# Erreur sur C
# ----------------------------------------------------
erreur_C = abs(
C - C_ancien
)
# ========================================================
# FIN DE LA BOUCLE SUR C
#
# C est maintenant convergé pour la valeur actuelle de T.
# ========================================================
# --------------------------------------------------------
# Recalcul du débit avec le C convergé
# --------------------------------------------------------
q = (
C
* epsilon
* E
* s
* math.sqrt(
2.0
* delta_p
* rho
)
)
# --------------------------------------------------------
# Nombre de Reynolds correspondant
# --------------------------------------------------------
Re_D = (
4.0
* q
/ (
math.pi
* D
* mu
)
)
# --------------------------------------------------------
# (11) V : vitesse moyenne de l'air
# --------------------------------------------------------
V = (
q
/ (
rho
* S
)
)
# --------------------------------------------------------
# (10) a : vitesse du son
# --------------------------------------------------------
a = math.sqrt(
gamma
* R_AIR
* T
)
# --------------------------------------------------------
# (9) M : nombre de Mach
# --------------------------------------------------------
M = V / a
# --------------------------------------------------------
# (8) Nouvelle température statique
# --------------------------------------------------------
T_nouveau = (
Ti
/ (
1.0
+ (
gamma - 1.0
)
* M**2
/ 2.0
)
)
# --------------------------------------------------------
# Erreur sur T
# --------------------------------------------------------
erreur_T = abs(
T_nouveau - T
)
# --------------------------------------------------------
# Mise à jour de T pour l'itération suivante
# --------------------------------------------------------
T = T_nouveau
# ============================================================
# FIN DES DEUX BOUCLES
# ============================================================
# ============================================================
# 12. MISE A JOUR FINALE DES GRANDEURS
# ============================================================
#
# T est maintenant convergée.
#
# On recalcule simplement les propriétés finales avec cette T.
#
# Il n'y a PAS de troisième boucle.
# ============================================================
# ------------------------------------------------------------
# (5) Viscosité dynamique finale
# ------------------------------------------------------------
mu = viscosite_air(T)
# ------------------------------------------------------------
# (7) Capacité thermique finale
# ------------------------------------------------------------
Cp = cp_air(T)
# ------------------------------------------------------------
# (6) Rapport Cp/Cv final
# ------------------------------------------------------------
gamma = (
Cp
/ (Cp - R_AIR)
)
# ------------------------------------------------------------
# (3) Coefficient de détente final
# ------------------------------------------------------------
epsilon = coefficient_detente(
beta,
p,
delta_p,
gamma
)
# ------------------------------------------------------------
# (13) Facteur de compressibilité final
# ------------------------------------------------------------
Z = facteur_compressibilite_Z(
p,
T
)
# ------------------------------------------------------------
# (12) Masse volumique finale
# ------------------------------------------------------------
rho = (
p
/ (
Z
* R_AIR
* T
)
)
# ------------------------------------------------------------
# (1) Débit massique final
# ------------------------------------------------------------
q = (
C
* epsilon
* E
* s
* math.sqrt(
2.0
* delta_p
* rho
)
)
# ------------------------------------------------------------
# (4) Nombre de Reynolds final
# ------------------------------------------------------------
Re_D = (
4.0
* q
/ (
math.pi
* D
* mu
)
)
# ------------------------------------------------------------
# (11) Vitesse moyenne finale
# ------------------------------------------------------------
V = (
q
/ (
rho
* S
)
)
# ------------------------------------------------------------
# (10) Vitesse du son finale
# ------------------------------------------------------------
a = math.sqrt(
gamma
* R_AIR
* T
)
# ------------------------------------------------------------
# (9) Nombre de Mach final
# ------------------------------------------------------------
M = V / a
# ============================================================
# 13. AFFICHAGE DES RESULTATS
# ============================================================
print("\n" + "=" * 80)
print(
"MESURE DE DEBIT D'AIR PAR DIAPHRAGME"
)
print("=" * 80)
# ============================================================
# DONNEES FIXES
# ============================================================
print("\nDONNEES FIXES / MODELE")
print(
f"r : constante spécifique de l'air "
f"= {R_AIR:.5f} J/(kg.K)"
)
print(
f"mu_ref : viscosité dynamique de référence "
f"= {MU_REF:.8e} Pa.s"
)
print(
f"T_ref : température de référence "
f"= {T_REF:.2f} K"
)
print(
f"C_suth : constante de Sutherland "
f"= {SUTHERLAND_C:.2f} K"
)
print(
f"mat. dia : matériau du diaphragme "
f"= {MATERIAU_DIAPHRAGME}"
)
print(
f"mat. tube : matériau de la conduite "
f"= {MATERIAU_CONDUITE}"
)
print(
f"d0 : diamètre de référence du diaphragme "
f"= {d0:.8f} m"
)
print(
f"D0 : diamètre de référence de la conduite "
f"= {D0:.8f} m"
)
# ============================================================
# GRANDEURS MESUREES
# ============================================================
print("\nGRANDEURS MESUREES")
print(
f"p : pression statique absolue en amont "
f"= {p:.3f} Pa"
)
print(
f"Delta_p : différence de pression amont-aval "
f"= {delta_p:.3f} Pa"
)
print(
f"Ti : température totale de l'air "
f"= {Ti:.6f} K"
)
print(
f"Td : température du diaphragme "
f"= {Td:.6f} K"
)
print(
f"TD : température de la conduite "
f"= {TD:.6f} K"
)
# ============================================================
# GEOMETRIE
# ============================================================
print("\n[A à F] GEOMETRIE")
print(
f"alpha_d : coefficient de dilatation du diaphragme "
f"= {alpha_d_val:.8e} 1/K"
)
print(
f"d (A) : diamètre du diaphragme à Td "
f"= {d:.9f} m"
)
print(
f"alpha_D : coefficient de dilatation de la conduite "
f"= {alpha_D_val:.8e} 1/K"
)
print(
f"D (B) : diamètre intérieur de la conduite à TD "
f"= {D:.9f} m"
)
print(
f"beta (C) : rapport des diamètres d/D "
f"= {beta:.9f}"
)
print(
f"s (D) : section de l'orifice du diaphragme "
f"= {s:.9e} m²"
)
print(
f"S (E) : section intérieure de la conduite "
f"= {S:.9e} m²"
)
print(
f"E (F) : coefficient d'approche "
f"= {E:.9f}"
)
# ============================================================
# PROPRIETES DE L'AIR
# ============================================================
print(
"\n[5, 7, 6, 3, 13, 12] "
"PROPRIETES DE L'AIR"
)
print(
f"mu (5) : viscosité dynamique de l'air "
f"= {mu:.9e} Pa.s"
)
print(
f"Cp (7) : capacité thermique massique à p constante "
f"= {Cp:.6f} J/(kg.K)"
)
print(
f"gamma (6) : rapport des capacités thermiques Cp/Cv "
f"= {gamma:.9f}"
)
print(
f"eps (3) : coefficient de détente "
f"= {epsilon:.9f}"
)
print(
f"Z (13) : facteur de compressibilité "
f"= {Z:.9f}"
)
print(
f"rho (12) : masse volumique de l'air en amont "
f"= {rho:.9f} kg/m³"
)
# ============================================================
# BOUCLE SUR C
# ============================================================
print(
"\n[1, 4, 2] "
"BOUCLE SUR LE COEFFICIENT DE DECHARGE"
)
print(
f"q (1) : débit massique "
f"= {q:.9f} kg/s"
)
print(
f"Re_D (4) : nombre de Reynolds basé sur D "
f"= {Re_D:.3f}"
)
print(
f"C (2) : coefficient de décharge "
f"= {C:.9f}"
)
# ============================================================
# BOUCLE SUR T
# ============================================================
print(
"\n[11, 10, 9, 8] "
"FERMETURE DE LA BOUCLE SUR LA TEMPERATURE"
)
print(
f"V (11) : vitesse moyenne dans la conduite "
f"= {V:.9f} m/s"
)
print(
f"a (10) : vitesse du son en amont "
f"= {a:.9f} m/s"
)
print(
f"M (9) : nombre de Mach en amont "
f"= {M:.9f}"
)
print(
f"T (8) : température statique en amont "
f"= {T:.9f} K"
)
# ============================================================
# CONVERGENCE
# ============================================================
print("\nCONVERGENCE")
print(
f"Nombre d'itérations sur T "
f"= {iteration_T}"
)
print(
f"Nombre d'itérations sur C lors du dernier passage "
f"= {iteration_C}"
)
print(
f"Erreur finale sur T "
f"= {erreur_T:.3e} K"
)
print(
f"Erreur finale sur C "
f"= {erreur_C:.3e}"
)
To embed this project on your website, copy the following code and paste it into your website's HTML: