import math
import numpy as np
# ============================================================
# 1. CONSTANTES DU PROBLEME
# ============================================================
# --- Air sec ---
R_AIR = 287.04 # r : constante spécifique de l'air [J/(kg.K)]
# --- Loi de Sutherland pour la viscosité dynamique de l'air ---
MU_REF = 18.2e-6 # mu0 : viscosité dynamique à T_REF [Pa.s]
T_REF = 293.15 # T_ref : température de référence = 20 °C [K]
SUTHERLAND_C = 113.0 # C_suth : constante de Sutherland [K]
# --- Matériaux retenus pour l'exemple ---
# Hypothèse pédagogique :
# - diaphragme en acier inoxydable 316L (EN 1.4404)
# - conduite en acier au carbone
# Le document DGA ne prescrit pas ces matériaux : ils doivent être adaptés
# au montage réel.
MATERIAU_DIAPHRAGME = "Acier inoxydable 316L (EN 1.4404)"
MATERIAU_CONDUITE = "Acier au carbone"
# --- Paramètres numériques de convergence ---
TOL_C = 1.0e-8 # tolérance sur C [-]
TOL_T = 1.0e-6 # tolérance sur T statique [K]
MAX_ITER_C = 100
MAX_ITER_T = 100
# ============================================================
# 2. GEOMETRIE FIXE DE L'INSTALLATION
# d0 et D0 sont mesurés à T_REF = 20 °C
# ============================================================
d0 = 0.060 # d0 : diamètre de référence du diaphragme [m]
D0 = 0.100 # D0 : diamètre intérieur de référence de la conduite [m]
# ============================================================
# 3. GRANDEURS MESUREES PENDANT UN POINT D'ESSAI
# ============================================================
p = 5.00e5 # p : pression statique absolue en amont [Pa]
delta_p = 1.00e4 # Delta_p = p_amont - p_aval [Pa]
Ti = 300.0 # Ti : température totale de l'air [K]
Td = 302.0 # Td : température du diaphragme [K]
TD = 301.0 # TD : température de la conduite [K]
# ============================================================
# 4. DILATATION THERMIQUE DES MATERIAUX
# ============================================================
# IMPORTANT : dans les relations du document :
#
# d = d0 * [1 + alpha_d(Td) * (Td - T_REF)]
# D = D0 * [1 + alpha_D(TD) * (TD - T_REF)]
#
# alpha_d(T) et alpha_D(T) sont ici des coefficients MOYENS de dilatation
# entre 20 °C et la température T. C'est cohérent avec cette écriture.
def alpha_diaphragme(T):
"""
alpha_d(T) : coefficient moyen de dilatation du diaphragme 316L [1/K].
Valeurs moyennes entre 20 °C et T :
20-100 °C : 16.0e-6 /K
20-200 °C : 16.5e-6 /K
20-300 °C : 17.0e-6 /K
20-400 °C : 17.5e-6 /K
20-500 °C : 18.0e-6 /K
Entre les points, on effectue une interpolation linéaire.
"""
theta = T - 273.15 # température en °C
# A proximité de 20 °C, on utilise la première valeur moyenne tabulée.
if 20.0 <= theta <= 100.0:
return 16.0e-6
temperatures_C = [100.0, 200.0, 300.0, 400.0, 500.0]
alpha_moyen = [16.0e-6, 16.5e-6, 17.0e-6, 17.5e-6, 18.0e-6]
if theta < temperatures_C[0] or theta > temperatures_C[-1]:
raise ValueError(
f"Température du diaphragme hors du domaine tabulé : {theta:.1f} °C"
)
# Interpolation linéaire avec la fonction préprogrammée NumPy
return float(np.interp(theta, temperatures_C, alpha_moyen))
def dilatation_relative_acier_carbone(T):
"""
Delta L / L0 d'un acier au carbone, L0 étant défini à 20 °C.
Relation utilisée pour 20 °C <= theta <= 750 °C.
"""
theta = T - 273.15 # °C
if not (20.0 <= theta <= 750.0):
raise ValueError(
"Modèle de dilatation de l'acier au carbone utilisé ici : "
"température attendue entre 20 °C et 750 °C."
)
return (
-2.416e-4
+ 1.2e-5 * theta
+ 0.4e-8 * theta**2
)
def alpha_conduite(T):
"""
alpha_D(T) : coefficient MOYEN de dilatation de la conduite en acier
au carbone entre 20 °C et T [1/K].
On le déduit de la dilatation relative DeltaL/L0 :
alpha_moyen = (DeltaL/L0) / (T - T_REF)
"""
theta = T - 273.15
if abs(theta - 20.0) < 1.0e-12:
# Limite locale de la loi au voisinage de 20 °C.
return 1.216e-5
elongation = dilatation_relative_acier_carbone(T)
return elongation / (theta - 20.0)
def dilatation_diaphragme(d_reference, temperature_diaphragme):
"""Calcule d, diamètre du diaphragme à la température Td."""
alpha_d = alpha_diaphragme(temperature_diaphragme)
return d_reference * (
1.0 + alpha_d * (temperature_diaphragme - T_REF)
)
def dilatation_conduite(D_reference, temperature_conduite):
"""Calcule D, diamètre intérieur de conduite à la température TD."""
alpha_D = alpha_conduite(temperature_conduite)
return D_reference * (
1.0 + alpha_D * (temperature_conduite - T_REF)
)
# ============================================================
# 5. FONCTIONS PHYSIQUES DE L'AIR
# ============================================================
def viscosite_air(T):
"""mu(T) : viscosité dynamique de l'air sec [Pa.s]."""
return (
MU_REF
* (T / T_REF) ** 1.5
* (SUTHERLAND_C + T_REF)
/ (SUTHERLAND_C + T)
)
def cp_air(T):
"""Cp(T) : capacité thermique massique à pression constante [J/(kg.K)]."""
exp_term = math.exp(3090.0 / T)
return R_AIR * (
3.5
- 2.8e-5 * T
+ 2.24e-8 * T**2
+ (3090.0 / T) ** 2
* exp_term
/ (exp_term - 1.0) ** 2
)
def facteur_compressibilite_Z(p_pa, T):
"""Z(p,T) : facteur de compressibilité de l'air [-], annexe II."""
p_bar = p_pa / 1.0e5
T_c = T - 273.15
if not (1.0 <= p_bar <= 65.0):
raise ValueError("Z : pression hors du domaine 1 à 65 bar de l'annexe II")
if not (-20.0 <= T_c <= 800.0):
raise ValueError("Z : température hors du domaine -20 à 800 °C de l'annexe II")
if p_bar <= 10.0:
if T_c <= 50.0:
return 1.0 - 1.0e-7 * p_bar * (T_c - 77.0) ** 2
return 1.0
if T_c <= 45.0:
return 0.9974 + 5.0e-6 * p_bar * (T_c - 47.0)
if T_c <= 180.0:
return (
0.9974
+ 1.3e-6 * (p_bar - 35.0) ** 2
+ 3.9375e-6 * p_bar * (T_c - 47.0)
)
return 1.0 + 4.0e-4 * p_bar
def coefficient_detente(beta, p_amont, delta_p, gamma):
"""epsilon : coefficient de détente du gaz [-]."""
p_aval = p_amont - delta_p
return 1.0 - (
0.351 + 0.256 * beta**4 + 0.93 * beta**8
) * (
1.0 - (p_aval / p_amont) ** (1.0 / gamma)
)
def coefficient_decharge(beta, Re_D):
"""C : coefficient de décharge pour prises de pression dans les angles [-]."""
a1 = (19000.0 * beta / Re_D) ** 0.8
return (
0.5961
+ 0.0261 * beta**2
- 0.216 * beta**8
+ 0.000521 * (beta * 1.0e6 / Re_D) ** 0.7
+ (0.0188 + 0.0063 * a1)
* beta**3.5
* (1.0e6 / Re_D) ** 0.3
)
# ============================================================
# 6. ALGORITHME - MEME ORDRE QUE L'ORGANIGRAMME DU DOCUMENT
# ============================================================
if delta_p <= 0.0:
q = 0.0
print("Delta_p <= 0 : débit massique q = 0 kg/s")
raise SystemExit
# ------------------------------------------------------------
# Etape A à F : géométrie à la température de l'essai
# ------------------------------------------------------------
d = dilatation_diaphragme(d0, Td) # (A) diamètre diaphragme [m]
D = dilatation_conduite(D0, TD) # (B) diamètre conduite [m]
beta = d / D # (C) rapport de diamètres [-]
s = math.pi * d**2 / 4.0 # (D) section de l'orifice [m²]
S = math.pi * D**2 / 4.0 # (E) section de la conduite [m²]
E = 1.0 / math.sqrt(1.0 - beta**4) # (F) coefficient d'approche [-]
# ------------------------------------------------------------
# Initialisation de l'organigramme
# ------------------------------------------------------------
C = 0.60 # C0 : première estimation de C
T_est = Ti # première estimation de T statique
# ------------------------------------------------------------
# Boucle extérieure : convergence de la température statique T
# ------------------------------------------------------------
for iteration_T in range(1, MAX_ITER_T + 1):
T_precedente = T_est
# Ordre exact de l'organigramme : mu, Cp, gamma, epsilon, Z, rho
mu = viscosite_air(T_est) # (5)
Cp = cp_air(T_est) # (7)
gamma = Cp / (Cp - R_AIR) # (6)
epsilon = coefficient_detente(beta, p, delta_p, gamma) # (3)
Z = facteur_compressibilite_Z(p, T_est) # (13)
rho = p / (Z * R_AIR * T_est) # (12)
# --------------------------------------------------------
# Boucle intérieure : convergence du coefficient C
# Ordre de l'organigramme : q, Re_D, C
# --------------------------------------------------------
for iteration_C in range(1, MAX_ITER_C + 1):
C_precedent = C
q = C * epsilon * E * s * math.sqrt(2.0 * delta_p * rho) # (1)
Re_D = 4.0 * q / (math.pi * D * mu) # (4)
C = coefficient_decharge(beta, Re_D) # (2)
if abs(C - C_precedent) < TOL_C:
break
else:
raise RuntimeError("La boucle sur C n'a pas convergé")
# Recalcul de q et Re_D avec C convergé
q = C * epsilon * E * s * math.sqrt(2.0 * delta_p * rho)
Re_D = 4.0 * q / (math.pi * D * mu)
# Ordre de l'organigramme : V, a, M, T
V = q / (rho * S) # (11)
a = math.sqrt(gamma * R_AIR * T_est) # (10)
M = V / a # (9)
T = Ti / (1.0 + (gamma - 1.0) * M**2 / 2.0) # (8)
if abs(T - T_precedente) < TOL_T:
T_est = T
break
T_est = T
else:
raise RuntimeError("La boucle sur T n'a pas convergé")
# ============================================================
# 7. RECALCUL FINAL AVEC T CONVERGEE
# Toujours dans l'ordre de l'organigramme
# ============================================================
T = T_est
mu = viscosite_air(T)
Cp = cp_air(T)
gamma = Cp / (Cp - R_AIR)
epsilon = coefficient_detente(beta, p, delta_p, gamma)
Z = facteur_compressibilite_Z(p, T)
rho = p / (Z * R_AIR * T)
# On recalcule C et q jusqu'à convergence finale avec cette T.
for iteration_C_finale in range(1, MAX_ITER_C + 1):
C_precedent = C
q = C * epsilon * E * s * math.sqrt(2.0 * delta_p * rho)
Re_D = 4.0 * q / (math.pi * D * mu)
C = coefficient_decharge(beta, Re_D)
if abs(C - C_precedent) < TOL_C:
break
else:
raise RuntimeError("La boucle finale sur C n'a pas convergé")
q = C * epsilon * E * s * math.sqrt(2.0 * delta_p * rho)
Re_D = 4.0 * q / (math.pi * D * mu)
V = q / (rho * S)
a = math.sqrt(gamma * R_AIR * T)
M = V / a
alpha_d = alpha_diaphragme(Td)
alpha_D = alpha_conduite(TD)
# ============================================================
# 8. AFFICHAGE DES RESULTATS
# Dans l'ordre de l'organigramme
# ============================================================
print("\n" + "=" * 86)
print("MESURE DE DEBIT PAR DIAPHRAGME - RESULTATS")
print("=" * 86)
print("\nDONNEES FIXES / MODELE")
print(f"r : constante spécifique de l'air = {R_AIR:.5f} J/(kg.K)")
print(f"mu0 : viscosité de référence de l'air = {MU_REF:.8e} Pa.s")
print(f"T_ref : température de référence = {T_REF:.2f} K")
print(f"C_suth : constante de Sutherland = {SUTHERLAND_C:.2f} K")
print(f"mat. dia : matériau du diaphragme = {MATERIAU_DIAPHRAGME}")
print(f"mat. tube: matériau de la conduite = {MATERIAU_CONDUITE}")
print(f"d0 : diamètre de référence du diaphragme = {d0:.8f} m")
print(f"D0 : diamètre intérieur de référence de la conduite = {D0:.8f} m")
print("\nGRANDEURS MESUREES")
print(f"p : pression statique absolue en amont = {p:.3f} Pa")
print(f"Delta_p : différence de pression statique amont-aval = {delta_p:.3f} Pa")
print(f"Ti : température totale de l'air = {Ti:.6f} K")
print(f"Td : température du diaphragme = {Td:.6f} K")
print(f"TD : température de la conduite = {TD:.6f} K")
print("\nRESULTATS - ORDRE DE L'ORGANIGRAMME")
print("\n[A à F] Géométrie corrigée de la dilatation")
print(f"alpha_d : coefficient moyen de dilatation du diaphragme = {alpha_d:.8e} 1/K")
print(f"d (A) : diamètre du diaphragme à Td = {d:.9f} m")
print(f"alpha_D : coefficient moyen de dilatation de la conduite = {alpha_D:.8e} 1/K")
print(f"D (B) : diamètre intérieur de la conduite à TD = {D:.9f} m")
print(f"beta(C) : rapport des diamètres d/D = {beta:.9f} -")
print(f"s (D) : section de l'orifice du diaphragme = {s:.9e} m²")
print(f"S (E) : section intérieure de la conduite = {S:.9e} m²")
print(f"E (F) : coefficient d'approche = {E:.9f} -")
print("\n[5, 7, 6, 3, 13, 12] Propriétés de l'air")
print(f"mu (5) : viscosité dynamique de l'air = {mu:.9e} Pa.s")
print(f"Cp (7) : capacité thermique massique à pression constante = {Cp:.6f} J/(kg.K)")
print(f"gamma(6) : rapport Cp/Cv = {gamma:.9f} -")
print(f"eps (3) : coefficient de détente = {epsilon:.9f} -")
print(f"Z (13) : facteur de compressibilité = {Z:.9f} -")
print(f"rho(12) : masse volumique de l'air en amont = {rho:.9f} kg/m³")
print("\n[1, 4, 2] Boucle sur le coefficient de décharge")
print(f"q (1) : débit massique = {q:.9f} kg/s")
print(f"ReD (4) : nombre de Reynolds basé sur D = {Re_D:.3f} -")
print(f"C (2) : coefficient de décharge = {C:.9f} -")
print("\n[11, 10, 9, 8] Fermeture de la boucle sur la température")
print(f"V (11) : vitesse moyenne de l'air dans la conduite = {V:.9f} m/s")
print(f"a (10) : vitesse du son en amont = {a:.9f} m/s")
print(f"M (9) : nombre de Mach en amont = {M:.9f} -")
print(f"T (8) : température statique en amont = {T:.9f} K")
print("\nINFORMATIONS DE CONVERGENCE")
print(f"iter_T : nombre d'itérations sur la température statique = {iteration_T}")
print(f"iter_C : nombre d'itérations finales sur C = {iteration_C_finale}")
To embed this project on your website, copy the following code and paste it into your website's HTML: