Distillation (dééthaniseur, équation d’état Peng-Robinson)

Le module Distillation simule une colonne de distillation multiconstituant (un dééthaniseur d’usine à gaz, 16 constituants C1→C12 + CO2 + N2) par la méthode de Thiele-Geddes avec convergence sur θ (theta). La thermodynamique repose sur une équation d’état de Peng-Robinson réécrite « from scratch » (règles de mélange avec \(k_{ij}\), solveur cubique, coefficients de fugacité et enthalpies de départ codés à la main), thermo n’étant utilisé que pour les propriétés critiques (\(T_c, P_c, \omega\)), les \(C_p^{ig}\) gaz parfait et les chaleurs de vaporisation.

Les deux briques du module sont :

  • PengRobinsonEOS — l’équation d’état de Peng-Robinson (Distillation.PengRobinson) ;

  • Deethanizer — la colonne dééthaniseur complète (Distillation.Deethanizer).

Deux scripts de diagnostic accompagnent le module (diag_enthalpy_models.py, reverse_engineer_enthalpies.py) : ils comparent les modèles d’enthalpie liquide Python aux références Delphi (condenseur partiel, écart \(\Delta H_{vap}\) Watson vs départ PR) et ne font pas partie de l’API.

Référence : chapitres 2, 5 et 6 (application au dééthaniseur).

Équation d’état de Peng-Robinson (PengRobinsonEOS)

Rôle

Fournit à la colonne les coefficients d’équilibre \(K_i\), les enthalpies (vapeur et liquide) et les calculs de flash. La classe PengRobinsonEOS(component_names) récupère depuis thermo (ChemicalConstantsPackage.from_IDs) les \(T_c, P_c, \omega\), les corrélations \(C_p^{ig}\) et \(\Delta H_{vap}\) ; les masses molaires viennent d’une base interne COMPONENT_DB (Mw), et la matrice d’interactions binaires \(k_{ij}\) du module kij_peng_robinson (build_kij_matrix). Constante des gaz : \(R = 8{,}314\)J/(mol·K).

Constantes pures et fonction \(\alpha(T)\)

Pour chaque constituant \(i\), les paramètres purs de Peng-Robinson sont :

\[ \begin{align}\begin{aligned}a_i^0 = 0{,}45724\,\frac{R^2 T_{c,i}^2}{P_{c,i}}, \qquad b_i = 0{,}07780\,\frac{R\,T_{c,i}}{P_{c,i}}\\\kappa_i = 0{,}37464 + 1{,}54226\,\omega_i - 0{,}26992\,\omega_i^2\end{aligned}\end{align} \]

La dépendance en température passe par \(a_i(T) = a_i^0\,\alpha_i(T)\) avec :

\[ \begin{align}\begin{aligned}\alpha_i(T) = \left[\,1 + \kappa_i\left(1 - \sqrt{T/T_{c,i}}\right)\right]^2\\\frac{d\alpha_i}{dT} = -\,\frac{\kappa_i\left[1 + \kappa_i(1 - \sqrt{T/T_{c,i}})\right]} {\sqrt{T\,T_{c,i}}}\end{aligned}\end{align} \]

Règles de mélange (van der Waals, un fluide) avec \(k_{ij}\)

\[ \begin{align}\begin{aligned}a_{ij} = \sqrt{a_i(T)\,a_j(T)}\,\bigl(1 - k_{ij}\bigr)\\a_m = \sum_i\sum_j z_i z_j\,a_{ij}, \qquad b_m = \sum_i z_i\,b_i\end{aligned}\end{align} \]

La dérivée \(da_m/dT\) (nécessaire à l’enthalpie de départ) est obtenue en dérivant \(a_{ij}\) terme à terme (méthode _mix_params, qui renvoie \(a_m, b_m, da_m/dT, a_i(T), a_{ij}\)).

Résolution de la cubique

En posant les grandeurs réduites \(A = a_m P/(RT)^2\) et \(B = b_m P/(RT)\), le facteur de compressibilité \(Z\) est racine de la forme réelle de Peng-Robinson :

\[Z^3 - (1 - B)\,Z^2 + \bigl(A - 3B^2 - 2B\bigr)\,Z - \bigl(AB - B^2 - B^3\bigr) = 0\]

La méthode _solve_cubic résout par numpy.roots et renvoie \((Z_V, Z_L)\) : la plus grande racine réelle positive est la vapeur, la plus petite (bornée à \(B\)) le liquide. Lorsqu’une seule racine réelle existe (mélange supercritique), la partie réelle de la paire complexe conjuguée sert de racine « virtuelle » pour la phase manquante.

Coefficients de fugacité

\[\ln \phi_i = \frac{b_i}{b_m}(Z - 1) - \ln(Z - B) - \frac{A}{2\sqrt{2}\,B} \left[\frac{2\sum_j z_j a_{ij}}{a_m} - \frac{b_i}{b_m}\right] \ln\!\left[\frac{Z + (1 + \sqrt{2})B}{Z + (1 - \sqrt{2})B}\right]\]

fugacity_coefficients(T, P, z, phase) prend \(Z = Z_V\) (phase='vapor') ou \(Z = Z_L\) (phase='liquid'). On en déduit les coefficients d’équilibre :

\[K_i = \frac{\phi_i^L}{\phi_i^V}\]

(méthode K_values). Une estimation initiale par la corrélation de Wilson est disponible (K_values_wilson) : \(K_i = (P_{c,i}/P)\,\exp\!\left[5{,}373\,(1+\omega_i)(1 - T_{c,i}/T)\right]\).

Enthalpies de départ et de mélange

Enthalpie de départ Peng-Robinson (J/mol) :

\[H^{dep} = R\,T\,(Z - 1) + \frac{T\,\dfrac{da_m}{dT} - a_m}{2\sqrt{2}\,b_m} \ln\!\left[\frac{Z + (1 + \sqrt{2})B}{Z + (1 - \sqrt{2})B}\right]\]

L’enthalpie gaz parfait intègre les \(C_p^{ig}\) depuis \(T_{ref} = 298{,}15\)K (ideal_gas_enthalpy). Deux modèles d’enthalpie de mélange sont disponibles :

  • mixture_enthalpy : \(H = H^{ig}(z) + H^{dep}\) (départ PR complet) ;

  • pure_component_enthalpy : approche « constituants purs / mélange idéal » (comme ProSim, Mixprops) — vapeur \(H = \sum_i z_i H^{ig}_i(T)\) ; liquide \(h = \sum_i z_i\bigl[H^{ig}_i(T) - \Delta H_{vap,i}(T)\bigr]\) (le \(\Delta H_{vap}\) DIPPR est retranché pour les constituants sous-critiques, nul si \(T \ge T_{c,i}\)).

Flash et point de bulle

  • flash_TP(T, P, z) — flash TP : boucle de substitution successive sur les \(K_i\) avec résolution de Rachford-Rice (Newton sur le taux de vapeur \(VF\)). Renvoie \((VF, x, y, K)\), avec détection des cas monophasiques (\(\sum z_i K_i \le 1\) → liquide ; \(\sum z_i / K_i \le 1\) → vapeur).

  • bubble_T(P, x, T_guess=None) — température de bulle par Newton-Raphson sur \(\sum_i K_i x_i - 1 = 0\). Renvoie \((T, K, y)\).

Constructeur et méthodes

Élément

Description

Unité

component_names

Liste des constituants (clés courtes : "C1", "C2", "iC4", "CO2", "N2"…)

—

fugacity_coefficients(T, P, z, phase)

\(\ln\phi_i\) (phase = 'vapor' / 'liquid')

—

K_values(T, P, x, y)

Coefficients d’équilibre \(K_i\)

—

K_values_wilson(T, P)

Estimation initiale de Wilson

—

enthalpy_departure(T, P, z, phase)

Enthalpie de départ PR

J/mol

mixture_enthalpy(T, P, z, phase)

\(H^{ig} + H^{dep}\)

J/mol

pure_component_enthalpy(T, P, z, phase)

Enthalpie « constituants purs »

J/mol

flash_TP(T, P, z)

Flash TP (Rachford-Rice)

\((VF, x, y, K)\)

bubble_T(P, x, T_guess)

Température de bulle

K

Exemple

import numpy as np
from Distillation.PengRobinson import PengRobinsonEOS

eos = PengRobinsonEOS(["C1", "C2", "C3", "nC4"])

T = 250.0        # K
P = 27.49e5      # Pa
z = np.array([0.30, 0.25, 0.25, 0.20])

# Flash TP : taux de vapeur, compositions liquide/vapeur, K
VF, x, y, K = eos.flash_TP(T, P, z)
print("VF =", VF)
print("K  =", K)

# Enthalpie de mélange (départ PR complet), J/mol
H_vap = eos.mixture_enthalpy(T, P, y, "vapor")
h_liq = eos.pure_component_enthalpy(T, P, x, "liquid")

# Température de bulle d'un liquide donné
Tb, Kb, yb = eos.bubble_T(P, x)

Colonne dééthaniseur (Deethanizer)

Rôle

Simule un dééthaniseur : colonne qui sépare, en tête, méthane + éthane (distillat \(D\)) et, en pied, les hydrocarbures plus lourds C3+ (résidu \(B\)). La résolution suit la méthode de Thiele-Geddes : bilans matière par plateau résolus par l”algorithme de Thomas modifié (Boston & Sullivan), correction des débits de sortie par convergence sur θ, mise à jour des températures par la méthode :math:`K_b` (constante 5,42), puis rebouclage par un bilan d’énergie bidirectionnel.

Numérotation des étages

Étage \(j\)

Rôle

j = 0

Condenseur (partiel — le distillat sort en vapeur)

j = 1 … N

Plateaux (\(N = 28\) par défaut)

j = N+1 = 29

Rebouilleur (partiel)

La colonne reçoit deux alimentations diphasiques : une froide (FC, ~-18 °C) et une chaude (FH, ~115 °C). Chaque alimentation est flashée puis répartie sur deux étages (éq. 5.24) : la vapeur entre à l’étage \(j_F-1\), le liquide à l’étage \(j_F\).

Bilans matière par plateau (MESH) et algorithme de Thomas

Pour chaque constituant \(i\), le facteur d’entraînement (stripping factor) de chaque étage est :

\[S_{ji} = \frac{K_{ji}\,V_j}{L_j}\]

Le système tridiagonal des débits liquides par constituant \(l_{ji}\) est résolu par l’algorithme de Thomas modifié (éq. 5.26), avec le vecteur d’alimentation \(F_j\) construit par le partage diphasique. Les débits vapeur en découlent : \(v_{ji} = S_{ji}\,l_{ji}\).

Convergence sur θ (correction des débits de sortie)

Les débits calculés \(d_i^{ca} = S_{0i}\,l_{0i}\) (distillat, éq. 5.13) et \(b_i^{ca} = l_{N+1,i}\) (résidu) sont corrigés pour respecter la spécification \(\sum_i d_i = D\) :

\[d_i^{co} = \frac{f_i}{1 + \theta\,(b_i/d_i)^{ca}}, \qquad b_i^{co} = f_i - d_i^{co}\]

où \(f_i = VFC_i + LFC_i + VFH_i + LFH_i\) est l’alimentation totale du constituant \(i\). Le facteur \(\theta\) est trouvé par Newton-Raphson sur \(g(\theta) = \sum_i d_i^{co} - D = 0\) (éq. 5.30–5.33), puis les débits internes \(l_{ji}, v_{ji}\) sont remis à l’échelle (éq. 5.34–5.35).

Mise à jour des températures — méthode \(K_b\) (constante 5,42)

La corrélation d’équilibre du document (éq. 2.13) emploie la même constante 5,42 pour tous les constituants (pas de dépendance en \(\omega\)) :

\[K_{ji} = \frac{P_{c,i}}{P}\,\exp\!\left[5{,}42\left(1 - \frac{T_{c,i}}{T_j}\right)\right]\]

À chaque étage, on calcule les volatilités relatives par rapport au constituant de référence \(b\) = C3 (propane), puis le nouveau \(K_b\) et la température par inversion (éq. 5.36–5.37) :

\[ \begin{align}\begin{aligned}K_{b,j}^{new} = \frac{1}{\sum_i \alpha_{ji}\,x_{ji}}, \qquad \alpha_{ji} = \frac{K_{ji}}{K_{b,j}}\\T_j^{new} = \frac{5{,}42\;T_{c,b}} {5{,}42 + \ln\!\bigl[P_{c,b}/(K_{b,j}^{new}\,P)\bigr]}\end{aligned}\end{align} \]

Bilan d’énergie bidirectionnel

Les débits \(V_j, L_j\) sont recalculés par un bilan enthalpique (éq. 5.53–5.67) descendant depuis le condenseur (section haute : \(L\) calculé, \(V\) par bilan matière) et remontant depuis le rebouilleur (section basse : \(V\) calculé, \(L\) par bilan matière), avec les termes d’alimentation aux étages de feed. Le condenseur est partiel :

\[ \begin{align}\begin{aligned}V_1 = L_0 + D, \qquad Q_C = V_1 H_1 - L_0 h_0 - D\,H_D\\Q_R = -\,h_{F}^{tot} + B\,h_B + D\,H_D + Q_C\end{aligned}\end{align} \]

où \(H_D = H_0\) (enthalpie vapeur au condenseur), \(h_B = h_{N+1}\) (enthalpie liquide au rebouilleur) et \(h_F^{tot}\) la somme des enthalpies d’alimentation. La mise à jour des débits est amortie (damping=0.5).

Boucle de résolution (solve)

À chaque itération externe :

  1. Thomas modifié avec partage d’alimentation (éq. 5.24, 5.26) → \(l_{ji}, v_{ji}\) ;

  2. convergence sur θ (éq. 5.27–5.33) → \(d_i, b_i\) ;

  3. mise à jour des températures \(K_b\) / 5,42 (éq. 5.36–5.37) ;

  4. recalcul des \(K_{ji}\) (éq. 2.13) à la nouvelle température ;

  5. bilan d’énergie bidirectionnel (éq. 5.53–5.67) → \(V_j, L_j, Q_C, Q_R\).

Convergence quand \(\max_j |T_j - T_j^{old}| < \text{tol}\).

Paramètres

Paramètre / attribut

Description

Défaut

n_trays

Nombre de plateaux (hors condenseur et rebouilleur)

28

components

Liste des constituants

16 corps : C1, C2, C3, iC4, nC4, iC5, nC5, C6–C12, CO2, N2

P

Pression de colonne

2,74912e6 Pa (27,49 bar)

D_total

Débit distillat (spécification)

562,1 kmol/h

B_total

Débit résidu (spécification)

724,443 kmol/h

L0

Débit de reflux (liquide au condenseur)

306,2 (cas 1/2) ; 236,1 (cas 3, reflux optimal)

j_FC / j_FH

Étages du liquide d’alimentation froide / chaude

5 / 21 (cas 1) ; 6 / 19 (cas 2 et 3)

Trois jeux de données prédéfinis (chapitre 6) sont fournis via solve(case=…) :

  • Cas 1 (set_default_case1) : FC@5, FH@21, \(L_0 = 306{,}2\) (fonctionnement courant) ;

  • Cas 2 (set_case2) : FC@6, FH@19, \(L_0 = 306{,}2\) ;

  • Cas 3 (set_case3) : FC@6, FH@19, \(L_0 = 236{,}1\) (reflux optimal).

Les alimentations flashées (VFC, LFC, VFH, LFH, en kmol/h) sont chargées à partir des tables de référence 6-3/6-4 (_load_reference_flash) et les profils initiaux \(T, V, L\) des tables 6-6/6-9/6-12.

Sorties

solve renvoie un dictionnaire (converged, iterations, T, V, L, K, d, b, QC, QR, components) et construit le DataFrame col.df (colonnes T (K), T (C), V (kmol/h), L (kmol/h), une ligne par étage). Grandeurs de référence (cas 1, table 6-7) : \(T_{tête} = 260{,}2\)K (-12,95 °C), \(T_{pied} = 416{,}6\)K (143,45 °C), \(Q_C = 4\,133\,587{,}91\)kJ/h, \(Q_R = 22\,854\,385{,}49\)kJ/h.

Exemple

from Distillation.Deethanizer import Deethanizer

col = Deethanizer()                 # 28 plateaux, 16 constituants
res = col.solve(case=1, max_outer=20, tol=1.0, verbose=True)

print("Convergé :", res["converged"], "en", res["iterations"], "itérations")
print("QC =", abs(res["QC"]), "kJ/h   QR =", abs(res["QR"]), "kJ/h")
print(col.df)                       # profils T, V, L par étage

# Débits de distillat (d) et de résidu (b) par constituant, kmol/h
for name, di, bi in zip(res["components"], res["d"], res["b"]):
    print(f"{name:>4}  d={di:12.4f}  b={bi:12.4f}")