Chapitre 8

Introduction aux méthodes

rigoureuses de distillation
multicomposant

8.1 Principe de résolution

Les méthodes présentées dans les chapitres précédents ne s’appliquent

qu’à des cas particuliers:

- pour des systèmes binaires, avec la méthode de Mc Cabe et Thiele

lorsque les hypothèses simplificatrices de Lewis sont applicables ou avec
la méthode de Ponchon-Savarit dans le cas contraire.

- pour la détermination d’une solution approximative à une distillation

multicomposant en utilisant l’approche FUG (Fenske-Underwood-Gilliland).

Les méthodes rigoureuses nécessitent la résolution sur tous les étages de

la colonne des équations « MESH »
•de bilans matières (M),
•d’équilibre liquide-vapeur (E), C C
•des contraintes sur les fractions molaires (S) ∑ xi = 1 et ∑ yi = 1
•et de bilan d’enthalpie (H). i =1 i =1

ULaval - GCH2004 - [Link]
Chap 8 - 1

 • Les premières méthodes de résolution du problème (dans les années 30)
se faisaient par itération, à la main, étage par étage, équation par équation.

• Avec les ordinateurs, ces méthodes ont cependant rencontrées des

problèmes de convergence. On a développé des méthodes (itératives elles
aussi) de résolution par composant (et non plus par étage) qui conduisent à
la résolution de système d’équations linéaires dont la solution s’obtient
rapidement par inversion de matrices tridiagonales.

• De nombreux logiciels de CAO (UniSim, Aspen, ProSim) permettent ainsi

la résolution exacte des opérations de distillation.

• Cependant, avant toute approche de résolution exacte, il est préférable de

disposer d’une solution approximative du problème comme celle obtenue
par la méthode FUG.

8.2 Expression matricielle des bilans matières

Numérotation
des étages

Entrées-sorties sur un étage j

ULaval - GCH2004 - [Link]
Chap 8 - 2

 Le bilan matière partiel sur l’étage j par rapport au composé i ( i varie de 1 à C):

L j −1 xi , j −1 + V j +1 yi , j +1 + F j zi , j − L j xi , j − V j yi , j = 0
On ignore les indices i:
L j −1 x j −1 + V j +1 y j +1 + F j z j − L j x j − V j y j = 0 [8.1]
yj = K jxj et y j +1 = K j +1 x j +1 [8.2]
lj , vj les débits molaires de chaque composé i sortant des étages
lj l j +1
xj = et x j +1 = [8.3]
Lj L j +1

En remplaçant 8.3 et 8.2 dans 8.1:

 VK   V j +1 K j +1 
−l j −1 + 1 + j j
  l j +  −  l j +1 = Fj z j [8.4]
 Lj   L j +1 

 Vj K j   V j +1 K j +1 
−l j −1 + 1 +
  l j +  −  l j +1 = Fj z j [8.4]
 Lj   L j +1 

A j l j −1 + B j l j + C j l j +1 = D j [8.5]
Vj K j
A j = −1 Bj = 1+ [8.6]
Lj
V j +1 K j +1
Cj = − D j = Fj z j
L j +1
Équ. 8.4 et 8.5 sont valables pour 2 ≤ j ≤ N et pour les C composés.

Si il n’y a pas d’alimentation alors Fj = D j = 0

ULaval - GCH2004 - [Link]
Chap 8 - 3

 Pour le condenseur total, j=1, le bilan matière s’écrit:

L1 x1 + DxD − V2 y2 = F1 z1 [8.7]
l1
x1 = xD , y2 = K 2 x2 , x1 =
L1

C1l2 + B1l1 = D1 [8.8]

D V2 K 2
B1 = 1 + C1 = − D1 = F1 z1
L1 L2

Pour le rebouilleur partiel, j=N, le bilan s’écrit :

− LN −1 xN −1 + BxB + VN y N = FN z N [8.10]
y N = K N xN , lN = LN xN = BxB

AN lN −1 + BN lN = DN [8.11]
VN K N V K
AN = −1 BN = 1 + = 1+ N N DN = FN z N
LN B

ULaval - GCH2004 - [Link]
Chap 8 - 4

 Pour chaque composé, l’ensemble des équations 8.5, 8.8 et 8.11 s’exprime:

 B1 C1 0 0 0 0 0 0 0   l1   D1 
A B2 C2 0 0 0 0 0 0   l   D 
 2   2   2 
0 A3 B3 C3 0 0 0 0 0   l3   D3 
     
 • = 
0 0 0 0 0 AN − 2 BN − 2 CN − 2 0  lN − 2   DN − 2 
     
0 0 0 0 0 0 AN −1 BN −1 CN −1   lN −1   DN −1 
0 0 0 0 0 0 0 AN BN   lN   DN 


[8.13]
Si on connait les coefficients Aj Bj Cj Dj
On peut calculer des débits molaires de chaque composé lj
La résolution de ce système d’équation linéaire se fait par inversion de
matrice tridiagonale.

Pour résoudre ce problème certaines

variables doivent être spécifiées.
Il est pratique de se fixer
les variables suivantes

F , zi , TF , N , N F
L
p, Treflux , ,D
D

ULaval - GCH2004 - [Link]
Chap 8 - 5

 8.3 Estimés initiaux des flux molaires L j et V j et des températures

Pour initier les itérations, les coefficients A j B j C j doivent être calculés.

Pour estimer Lj et Vj , on pourra utiliser les hypothèses de Lewis de flux
molaires constants.
L
Les spécifications sur et D étant données, on peut alors estimer
D
les flux molaires liquide et vapeur des zones de rectification et d’épuisement
et tenant compte des conditions de l’alimentation

L
L = D V = L+D L = L + qF
D
L − L ( HV − H F )
V = V − (1 − q ) F q= =
F ( HV − H L )

8.3 (suite) Estimés initiaux des températures

Pour les températures, on peut supposer une variation linéaire entre le

pied et la tête entre les températures de bulle du distillat et du résidu.

La connaissance de ces températures va permettre le calcul des

volatilités Ki de chaque composé sur chaque plateau.

Dans le cas où celles-ci seraient dépendantes des compositions, il

faudrait les estimer en utilisant des composés clés (voir méthode FUG).

Commentaire : ces informations devraient connues à l’issue de la résolution

approximative qui doit être faite préalablement.

ULaval - GCH2004 - [Link]
Chap 8 - 6

 8.4 Calculs des points de bulle

Après la résolution du système 8.13, le débit molaire de chaque composé l j

est connu sur chaque étage j et on peut recalculer de nouvelles

lj
fractions molaires xj =
Lj

Le calcul du point de bulle avec ces nouvelles concentrations va permettre

de recalculer la température du plateau

Avec ces nouvelles valeurs, on pourrait ainsi reprendre les calculs des
sections 8.3 et 8.4. Cependant, la convergence du profil de température
dans la colonne n’est pas garantie et la méthode de convergence dite
méthode θ a été développée.

8.5 Méthode de convergence θ

Cette méthode définit le paramètre θ comme solution de l’équation

 
 
C
 Fzi 
Dspécifié = ∑ [8.14]
  
i =1
1 + θ  B xi , B  
 Dx
  i,D cal 
: Cette équation a une racine positive que l’on peut obtenir à l’aide de la
méthode de Newton en cherchant le zéro de la fonction

 
 
C
 Fzi 
f (θ ) = ∑  − Dspécifié = 0 [8.15]
  
i =1
1 + θ  B xi , B  
 Dx
  i,D cal 

ULaval - GCH2004 - [Link]
Chap 8 - 7

 méthode de Newton: intersection de la tangente à la courbe avec l’axe des θ

3
f (θ ) droite de pente f ' (θ k )
2

4 1
θ
θ sol θ k +1 θk

0 − f (θ k ) f (θ k )
f ' (θ k ) = ⇒ θ k +1 = θ k −
θ k +1 − θ k f ' (θ k )

 
f (θ k )  
θ k +1 = θ k − C
 Fzi 
f ' (θ k ) f (θ ) = ∑  − Dspécifié
  
i =1
1 + θ  B xi , B  
 Dx
cal 
'
 1   −u '    i, D
rappel :   =  2 
u  u 
f (θ k )
Et d’une itération à l’autre θ k +1 = θ k +
 
Équation 8.16:  Fz  B xi , B  
i 
C  D x 
 i , D  cal
∑ 
i =1  
2

  B x   

 1 + θ k  i,B
  
   D xi , D cal  

 B xi , B  se calcule à l’aide des solutions de l’équation 8.13

 
L
 D xi , D cal (B x ) = l
i,B i,N (D x ) = l
i,D i ,1  
D

ULaval - GCH2004 - [Link]
Chap 8 - 8

  
  C
D = ∑ ( Dxi , D )
C
 Fzi 
Dspécifié = ∑ [8.14]
 B xi , B  
i =1
1 + θ   
 Dx i =1
  i,D cal 

La nouvelle valeur de θ
(équ. 8.16) permet de corriger les valeurs des
débits molaires de chaque composé au distillat (8.19) et au résidu (8.20)

 
 
 Fzi 
( Dx ) i , D cor =
  
[8.19]
1 + θ  B xi , B  
 Dx
  i,D cal 

 B xi , B 
( Bx ) i , B cor = ( Dxi , D )cor θ 
Dx  [8.20]
 i,D cal

Les équations 8.14, 8.19 et 8.20 donnent :

C C

∑(D x )
i =1
i,D = Dspécifié ∑( B x )
i =1
i,B
cor
= B = F − Dspécifié
cor

La méthode de convergence θ impose le respect des débits

de distillat et de résidu.

Les débits molaire liquide de chaque composé, sur chaque étage sont
alors corrigés, dans la zone de rectification, et dans la zone d’épuisement:

.
(D x ) l (B x ) l
( li, j ) = ( )
i , D cor
(l ) = ( )
i , B cor

cor
(D x ) i , D cal
i , j cal i , j cor
(B x ) i , B cal
i , j cal

xi , j =
(l ) i , j cor
Ce sont ces fractions molaires qui sont utilisées
pour le calcul des températures de bulles Tj sur
C chaque étage. La convergence des températures
∑ (l )
i =1
i , j cor
Tj s’observe lorsque θ → 1 ± 10−5

ULaval - GCH2004 - [Link]
Chap 8 - 9

 8.6 Expression matricielle des bilans d’enthalpie

Après convergence des températures, il faut tenir compte des bilans

d’enthalpie sur chaque étage qui s’écrit :

L j h j + V j H j = L j −1h j −1 + V j +1 H j +1 + Fj hFj + Q j [8.26]

pour 2 ≤ j ≤ N −1
Le bilan matière entre le dessous d’un étage j et le haut de la colonne:

J J −1
L j = V j +1 − D + ∑ Fk L j −1 = V j − D + ∑ Fk [8.27]
k =1 k =1

En remplaçant 8.27 dans 8.26

V j +1 ( h j − H j +1 ) + V j ( H j − h j −1 ) = [8.28]
 j   j −1 
Fj hF j + Q j − D ( h j −1 − h j ) −  ∑ Fk  h j +  ∑ Fk  h j −1
 k =1   k =1 

Pour le condenseur total en tête (j=1) V2 = L1 + D et h1 = hD

le bilan d’enthalpie: Q1 + V2 H 2 + F1hF 1 = ( L1 + D ) h1 = V2 h1

V2 ( h1 − H 2 ) = Q1 + F1hF 1 [8.30]
Pour le rebouilleur partiel en pied (j=N) LN −1 = B + VN − FN
BhB + VN H N = LN −1hN −1 + FN hFN + QN
VN ( H N − hN −1 ) = QN + FN hFN + B ( hN −1 − hN ) − FN −1hFN −1 [8.33]
pour 2 ≤ j ≤ N −1
V j +1 ( h j − H j +1 ) + V j ( H j − h j −1 ) = [8.28]
 j   j −1 
Fj hFj + Q j − D ( h j −1 − h j ) −  ∑ Fk  h j +  ∑ Fk  h j −1
 k =1   k =1 
Les équations 8.33, 8.30 et 8.28 peuvent s’écrire sous la forme matricielle :

ULaval - GCH2004 - [Link]
Chap 8 - 10

  BE1 C E1 0 0 0 0 0 0 0   V1   DE1 
 0 BE 2 CE 2 0 0 0 0 0 0  V   D 
   2   E2 
 0 0 BE 3 CE 3 0 0 0 0 0   V3   DE 3 
     
 • = 
 0 0 0 0 0 0 BEN − 2 CEN − 2 0  VN − 2   DEN − 2 
     
 0 0 0 0 0 0 0 BN −1 CEN −1  VN −1   DEN −1 
 0 0 0 0 0 0 0 0 BEN   VN   DEN 

[8.34]
BE1 = 0 CE1 = h1 − H 2 DE1 = F1 hF 1 + Q1

pour 2 ≤ j ≤ N − 1 BEj = H j − H j −1 CEj = h j − H j +1

 j   j −1 
DE1 = F j hFj + Q j − D ( h j −1 − h j ) −  ∑ Fk  h j +  ∑ Fk  h j −1
 k =1   k =1 

BEN = H N − H N −1 DEN = FN hFN + QN + B ( hN −1 − hN ) − FN hN −1

Le calcul des enthalpies des courants liquide ou vapeur pour des

mélanges idéaux s’obtient à partir des enthalpies des composés purs

h%i H% i calculées à la température de l’étage Tj:

C C
h j (T ) = ∑ xi , j h%i(T ) H j (T ) = ∑ yi , j H% i (T )
j j j j
i =1 i =1
ou

h% T j = h% Tref + CPL (TJ − Tref )

H% T j = H% Tref + C PV (TJ − Tref ) + λi ,Tref

λi,T ref
chaleur latente de vaporisation du
composé pur i

ULaval - GCH2004 - [Link]
Chap 8 - 11

 Dans le cas d’une colonne adiabatique on aurait :

Qj = 0 pour 2 ≤ j ≤ N − 1; QN = QRebouilleur ; Q1 = QCondenseur

Pour le condenseur total, l’énergie à retirer s’exprime par :

 L
Q1 = V2 ( hD − H 2 ) = D 1 +  ( hD − H 2 )
 D

Pour le rebouilleur partiel, l’énergie à fournir s’exprime à partir du bilan

d’enthalpie sur toute la colonne :

 N 
QN = DhD + BhB −  ∑ Fk hFk  − Q1
 k =1 

L’inversion de la matrice 8.34

donne les nouvelles valeurs de Vj .
Des nouvelles valeurs de Lj sont
estimées avec l’équation 8.27.

J
L j = V j +1 − D + ∑ Fk
k =1
On pourra reprendre alors les calculs à
la section 8.3 et ce jusqu'à convergence.
Entre 2 itérations k et k+1

L j , k +1 − L j ,k
<ε
L j ,k
V j ,k +1 − V j , k
<ε
V j ,k

Une valeur de ε entre 10-5 et 10-4 peut être utilisée.

ULaval - GCH2004 - [Link]
Chap 8 - 12