École Nationale Supérieure en Génie des Technologies

Industrielles

Filière : Génie des procédés

Option : Conception de Procédés Assistés par Ordinateur

— Projet de modélisation des opérations unitaires II —

Modélisation d’une distillation multi-constituants

multi-étagées en régime permanent

Réalisé par : DI PASQUALE Timothé

Encadré par : M. MARIAS

Année universitaire
2021 / 2022
Table des matières
Glossaire

Introduction 1

1 Modèle physique 2
1.1 Équations MESH . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.1.1 Équations aux condenseur . . . . . . . . . . . . . . . . . . . . . . . . 2
1.1.2 Équations aux plateaux intermédiaires . . . . . . . . . . . . . . . . . . 2
1.1.3 Équations au bouilleur . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.2 Équations constitutives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.2.1 Modèle thermodynamique . . . . . . . . . . . . . . . . . . . . . . . . 3
1.2.2 Modèle enthalpique . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

2 Modèle mathématique 5
2.1 Variables du problème . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.2 Méthode de résolution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.2.1 La méthode MRSL01 . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.2.2 La méthode MRSL21 . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

3 Organisation du programme 8
3.1 Structure globale du programme principale . . . . . . . . . . . . . . . . . . . 8
3.2 Organisation des vecteurs de variables . . . . . . . . . . . . . . . . . . . . . . 9
3.3 Initialisation du vecteur des variables . . . . . . . . . . . . . . . . . . . . . . . 9
3.3.1 Les bilans matières et d’énergie . . . . . . . . . . . . . . . . . . . . . 9
3.3.2 Débit liquide L initial . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
3.3.3 Débit vapeur V initial . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
3.3.4 Température et compositions liquides et vapeurs initiales . . . . . . . . 10

4 Analyse du programme 12
4.1 Résultats avec les conditions de base . . . . . . . . . . . . . . . . . . . . . . . 12
4.2 Influence du nombre de plateaux . . . . . . . . . . . . . . . . . . . . . . . . . 13
4.3 Influence de la pression . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
4.4 Influence de la puissance au bouilleur . . . . . . . . . . . . . . . . . . . . . . 15
4.5 Influence de la puissance au condenseur . . . . . . . . . . . . . . . . . . . . . 16

Conclusion 18

Annexe 1 - Exécution du programme 19

Annexe 2 - Exécution du programme 19

Table des figures

1 Représentation schématique de la colonne de distillation étudiée . . . . . . . . 1
2 Schéma de l’algorithme de Newton-Raphson . . . . . . . . . . . . . . . . . . 5
3 Matrice tridiagonale par blocs . . . . . . . . . . . . . . . . . . . . . . . . . . 6
4 Logigramme non-exhaustif du programme principal . . . . . . . . . . . . . . . 8
5 Diagramme ternaire avec les conditions de base pour deux types de composition 13
6 Diagramme ternaire montrant l’influence du nombre de plateaux sur la distillation 14
7 Diagramme ternaire montrant l’influence de la pression opératoire sur la distil-
lation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
8 Diagramme ternaire montrant l’influence de la puissance au bouilleur . . . . . 16
9 Diagramme ternaire montrant l’influence de la puissance au condenseur . . . . 17
Glossaire

Lettre Unité Description

Indices
i [-] Référence au numéro du constituant
j [-] Référence au numéro du plateau

Variables
𝑉𝑗 [mol/s] Débit molaire vapeur sortant du plateau j
𝑦𝑖, 𝑗 [-] Fraction molaire vapeur en constituant i au plateau j
𝑇𝑗 [K] Température au plateau j
𝑥𝑖, 𝑗 [-] Fraction molaire liquide en constituant i au plateau j
𝐿𝑗 [mol/s] Débit molaire liquide sortant du plateau j

Équations
𝐵𝑀𝑇 𝑗 [-] Bilan matière total au plateau j
𝐵𝑀 𝑃𝑖, 𝑗 [-] Bilan matière partiel en constituant i au plateau j
𝐵𝐸 𝑗 [-] Bilan d’énergie au plateau j
𝐸𝑄 𝑖, 𝑗 [-] Équation d’équilibre en constituant i au plateau j
𝑆𝑂 𝑀 𝑗 [-] Équation de sommation au plateau j

Paramètres
𝑁𝐶 [-] Nombre de constituants
𝑁𝑃 [-] Nombre de plateaux dans la colonne
𝑃 [Pa] Pression opératoire
𝐹𝑗 [mol/s] Débit molaire d’alimentation au plateau j
𝐷𝐿 𝑗 [mol/s] Débit molaire de soutirage liquide au plateau j
𝐷𝑉 𝑗 [mol/s] Débit molaire de soutirage vapeur au plateau j
𝑧𝑖, 𝑗 [-] Fraction molaire en constituant i dans l’alimentation au plateau j
𝑚ℎ 𝐿𝑗 [J/mol] Enthalpie molaire de la phase liquide au plateau j
𝑚𝐻𝑉𝑗 [J/mol] Enthalpie molaire de la phase vapeur au plateau j
𝑣𝑎 𝑝
𝐻𝑖 (𝑇) [J/mol] Enthalpie de vaporisation du constituant i
𝐻𝐹 𝑗 [J/mol] Enthalpie molaire de l’alimentation du plateau j
𝑄𝑗 [J/s] Énergie échangée au plateau j
𝑄𝑐 [J/s] Énergie échangée au condenseur
𝑄𝑏 [J/s] Énergie échangée au bouilleur
𝑇𝑏𝑢𝑙𝑙𝑒 [K] Température de bulle
𝑇𝑟𝑜𝑠𝑒𝑒 ´ [K] Température de rosée
𝑃𝑖𝑠𝑎𝑡 (𝑇) [Pa] Pression de saturation du constituant i
𝛾𝑖 (𝑇) [-] Coefficient d’activité du constituant i
Introduction
Le projet de modélisation des opérations unitaires II qui est exposé dans ce rapport avait pour
but de modéliser une distillation. Cette distillation est particulière puisqu’elle est modulable.
Premièrement, le nombre de plateaux ainsi que le nombre de constituants sont variables. Ensuite,
on peut faire intervenir plusieurs alimentations, des soutirages liquides et vapeurs aux plateaux
que l’on souhaite. Enfin, on peut jouer sur les quantités de chaleur au bouilleur, au condenseur
ainsi qu’à différents plateaux. On retrouve les différents paramètres au niveau de la figure 1.

Outre le fait de modéliser simplement cette distillation, il y a un aspect important puisque

cela est utile dans le monde industriel afin de simuler différents scénarios virtuellement. En
effet, cela permet de gagner du temps et de l’argent car le comportement de la distillation sera
connue sans devoir construire la colonne.

Figure 1 – Représentation schématique de la colonne de distillation étudiée

La résolution du problème se fait par une méthode de Newton-Raphson en langage Fortran

90. L’utilisateur aura le choix de la méthode de résolution qui intervient dans le calcul de la
résolution d’un système linéaire. Ainsi l’utilisateur pourra visualiser l’influence de la méthode
de résolution. Ces deux méthodes sont détaillées dans la suite du projet.

1
1 Modèle physique
Le modèle physique sert à poser le problème en reposant sur les phénomènes physiques
mis en jeu lors de l’opération. Pour cela, il faut poser une certains nombre d’équations afin de
calculer toutes les inconnues.

Notre colonne est modélisée de façon à distinguer trois zones distinctes : le bouilleur (dernier
plateau), le condenseur (premier plateau) et les plateaux intermédiaires. De plus, l’hypothèse
principale est que chaque étage est à l’équilibre, c’est-à-dire que la vapeur et le liquide sortants
de chaque plateau sont en équilibre thermodynamique.

1.1 Équations MESH

1.1.1 Équations aux condenseur

[𝐵𝑀𝑇1 ] : −(𝐿 1 + 𝐷 𝐿 1 ) − 𝑉1 + 𝑉2 + 𝐹1 = 0 (1)

[𝐵𝑀 𝑃𝑖,1 ] : −(𝐿 1 + 𝐷 𝐿 1 ) · 𝑥𝑖,1 − 𝑉1 · 𝑦𝑖,1 + 𝑉2 · 𝑦𝑖,2 + 𝐹1 · 𝑧𝑖,1 = 0 (2)
[𝐸𝑄 𝑖,1 ] : 𝑦𝑖,1 − 𝑚𝐾𝑖,1 · 𝑥𝑖,1 = 0 (3)
𝑁𝐶
∑︁
[𝑆𝑇1 ] : 𝑥 𝑘,1 − 𝑦 𝑘,1 = 0 (4)
𝑘=1
[𝐵𝐸 1 ] : −(𝐿 1 + 𝐷 𝐿 1 ) · 𝑚ℎ1𝐿 − 𝑉1 · 𝑚𝐻1𝑉 + 𝑉2 · 𝑚𝐻2𝑉 + 𝐹1 · 𝐻 𝐹1 − 𝑄 𝐶 = 0 (5)

1.1.2 Équations aux plateaux intermédiaires

[𝐵𝑀𝑇 𝑗 ] : 𝐿 𝑗−1 − (𝐿 𝑗 + 𝐷 𝐿 𝑗 ) − (𝑉 𝑗 + 𝐷𝑉 𝑗 ) + 𝑉 𝑗+1 + 𝐹 𝑗 = 0 (6)

[𝐵𝑀 𝑃𝑖, 𝑗 ] : 𝐿 𝑗−1 · 𝑥𝑖, 𝑗−1 − (𝐿 𝑗 + 𝐷 𝐿 𝑗 ) · 𝑥𝑖, 𝑗 − (𝑉 𝑗 + 𝐷𝑉 𝑗 ) · 𝑦𝑖, 𝑗 + 𝑉 𝑗+1 · 𝑦𝑖, 𝑗+1 + 𝐹 𝑗 · 𝑧𝑖, 𝑗 = 0 (7)

[𝐸𝑄 𝑖, 𝑗 ] : 𝑦𝑖, 𝑗 − 𝑚𝐾𝑖, 𝑗 · 𝑥𝑖, 𝑗 = 0 (8)

𝑁𝐶
∑︁
[𝑆𝑇 𝑗 ] : 𝑥 𝑘, 𝑗 − 𝑦 𝑘, 𝑗 = 0 (9)
𝑘=1

[𝐵𝐸𝑖 ] : 𝐿 𝑗−1 ·𝑚ℎ 𝐿𝑗−1 −(𝐿 𝑗 +𝐷 𝐿 𝑗 )·𝑚ℎ 𝐿𝑗 −(𝑉 𝑗 +𝐷𝑉 𝑗 )·𝑚𝐻𝑉𝑗 +𝑉 𝑗+1 ·𝑚𝐻𝑉𝑗+1 +𝐹 𝑗 ·𝐻 𝐹 𝑗 −𝑄 𝑗 = 0 (10)

1.1.3 Équations au bouilleur

[𝐵𝑀𝑇𝑁 𝑃 ] : 𝐿 𝑁 𝑃−1 − 𝐿 𝑁 𝑃 − (𝑉𝑁 𝑃 + 𝐷𝑉𝑁 𝑃 ) + 𝐹𝑁 𝑃 = 0 (11)

[𝐵𝑀 𝑃𝑖,𝑁 𝑃 ] : 𝐿 𝑁 𝑃−1 · 𝑥𝑖,𝑁 𝑃−1 − 𝐿 𝑁 𝑃 · 𝑥𝑖,𝑁 𝑃 − (𝑉𝑁 𝑃 + 𝐷𝑉𝑁 𝑃 ) · 𝑦𝑖,𝑁 𝑃 + 𝐹𝑁 𝑃 · 𝑧𝑖,𝑁 𝑃 = 0 (12)

[𝐸𝑄 𝑖,𝑁 𝑃 ] : 𝑦𝑖,𝑁 𝑃 − 𝑚𝐾𝑖,𝑁 𝑃 · 𝑥𝑖,𝑁 𝑃 = 0 (13)

𝑁𝐶
∑︁
[𝑆𝑇 𝑗 ] : 𝑥 𝑘,𝑁 𝑃 − 𝑦 𝑘,𝑁 𝑃 = 0 (14)
𝑘=1

[𝐵𝐸𝑖 ] : 𝐿 𝑁 𝑃−1 · 𝑚ℎ 𝑁𝐿 𝑃−1 − 𝐿 𝑁 𝑃 · 𝑚ℎ 𝑁𝐿 𝑃 − (𝑉𝑁 𝑃 + 𝐷𝑉𝑁 𝑃 ) · 𝑚𝐻𝑉𝑁 𝑃 + 𝐹𝑁 𝑃 · 𝐻 𝐹𝑁 𝑃 + 𝑄 𝐵 = 0 (15)

2
1.2 Équations constitutives
1.2.1 Modèle thermodynamique
Le modèle thermodynamique va nous permettre de déterminer les équilibres entre le liquide
et la vapeur pour chaque constituant. Pour cela, nous utilisons une approche dissymétrique 𝛾 − 𝜙.
En partant de cela, nous avons une égalité des fugacités :

𝑓𝑖𝐿 (𝑇, 𝑃, 𝑥𝑖 ) = 𝑓𝑖𝑉 (𝑇, 𝑃, 𝑦𝑖 ) (16)

En détaillant chacun des termes par leur expression :

𝑥𝑖 · 𝛾𝑖 (𝑇, 𝑃, 𝑥𝑖 ) · 𝑓𝑖∗,𝐿 (𝑇, 𝑃) = 𝑦𝑖 · 𝜙𝑉𝑖 (𝑇, 𝑃, 𝑦𝑖 ) · 𝑃 (17)

Nous pouvons ensuite détailler la fugacité du corps pur, 𝑓𝑖∗,𝐿 (𝑇, 𝑃) comme :

𝑣 𝑖∗,𝐿
∫ 𝑃
𝑓𝑖∗,𝐿 (𝑇, 𝑃) = 𝑓𝑖∗,𝐿 (𝑇, 𝑃𝑖𝑠𝑎𝑡 (𝑇)) · exp d𝑃 (18)
| {z } 𝑃𝑖𝑠𝑎𝑡 (𝑇) 𝑅 · 𝑇

𝑖 (𝑇,𝑃𝑖 ) ·𝑃𝑖 (𝑇)

𝜙𝑉 𝑠𝑎𝑡 𝑠𝑎𝑡 | {z }
𝐹𝑎𝑐𝑡𝑒𝑢𝑟 𝑑𝑒 𝑃𝑜𝑦𝑛𝑡𝑖𝑛𝑔

L’égalité devient finalement :

𝑣 𝑖∗,𝐿
∫ 𝑃
𝑥𝑖 · 𝛾𝑖 (𝑇, 𝑃, 𝑥𝑖 ) · 𝜙𝑉𝑖 (𝑇, 𝑃𝑖𝑠𝑎𝑡 ) · 𝑃𝑖𝑠𝑎𝑡 (𝑇) · exp d𝑃 = 𝑦𝑖 · 𝜙𝑉𝑖 (𝑇, 𝑃, 𝑦𝑖 ) · 𝑃 (19)
𝑃𝑖𝑠𝑎𝑡 (𝑇) 𝑅 ·𝑇

Cependant nous allons poser des hypothèses qui permettent de simplifier l’équation 19.
— Le facteur de Poynting est égal à 1 si la pression opératoire est de l’ordre de grandeur
de la pression de saturation.

𝑥𝑖 · 𝛾𝑖 (𝑇, 𝑃, 𝑥𝑖 ) · 𝜙𝑉𝑖 (𝑇, 𝑃𝑖𝑠𝑎𝑡 ) · 𝑃𝑖𝑠𝑎𝑡 (𝑇) = 𝑦𝑖 · 𝜙𝑉𝑖 (𝑇, 𝑃, 𝑦𝑖 ) · 𝑃 (20)

— Le gaz est supposé être un gaz parfait (P<40 bars) :

𝜙𝑉𝑖 (𝑇, 𝑃, 𝑦𝑖 ) = 1 (21)

— La dernière hypothèse est que la solution est idéale :

𝜙𝑖𝑠𝑎𝑡 (𝑇, 𝑃𝑖𝑠𝑎𝑡 ) = 1 (22)

On obtient alors l’expression finale :

𝑥𝑖 · 𝛾𝑖 (𝑇, 𝑃, 𝑥𝑖 ) · 𝑃𝑖𝑠𝑎𝑡 (𝑇) = 𝑦𝑖 · 𝑃 (23)

1.2.2 Modèle enthalpique

Enthalpie liquide
Afin de calculer l’enthalpie liquide, une référence doit être posée. La référence choisie ici est
l’état corps pur vapeur à une température 𝑇𝑟𝑒 𝑓 = 298, 15𝐾. De plus, une hypothèse est ajoutée au
modèle : les capacités thermiques massiques de chaque constituant sont supposées constantes.

Pour que notre modèle puisse représenter au mieux la réalité, nous allons ajouter le terme
qui permet de caractériser l’énergie de mélange. En effet, nous ne possédons pas de corps pur
directement. Ce terme est alors l’enthalpie d’excès. L’enthalpie peut être calculée comme suit :

𝑚ℎ 𝐿 (𝑇, 𝑃, 𝑥𝑖 ) = 𝑚ℎ0,𝐿 (𝑇, 𝑃, 𝑥𝑖 ) + ℎ 𝐸 (𝑇, 𝑃, 𝑥𝑖 ) (24)

3
 Avec l’enthalpie des corps purs :

𝑁𝐶
" ∫ ∫ !#
∑︁ 𝑇𝑒𝑏,𝑖 𝑇
0,𝐿 𝑟𝑒 𝑓 𝑣𝑎 𝑝 𝑣𝑎 𝑝 𝑙𝑖𝑞
𝑚ℎ (𝑇, 𝑃, 𝑥𝑖 ) = ℎ𝑖 + 𝑥𝑖 · 𝐶 𝑝𝑖 d𝑇 − 𝐻𝑖 (𝑇𝑒𝑏,𝑖 ) + 𝐶 𝑝𝑖 d𝑇 (25)
𝑖=1 𝑇𝑟 𝑒 𝑓 𝑇𝑒𝑏,𝑖

Et l’enthalpie d’excès :
𝑁𝐶  
𝐸 2
∑︁ 𝜕 (𝑙𝑛𝛾𝑖 )
ℎ (𝑇, 𝑃, 𝑥𝑖 ) = −𝑅 · 𝑇 · 𝑥𝑖 · (26)
𝑖=1
𝜕𝑇 𝑃,𝑥𝑖

Enthalpie vapeur
𝑁𝐶
" ∫ !#
∑︁ 𝑇𝑒𝑏,𝑖
− (𝐻𝑉 − 𝐻 ∗ )𝑇,𝑃,𝑦𝑖
𝑟𝑒 𝑓 𝑣𝑎 𝑝
𝑚𝐻𝑉 (𝑇, 𝑃, 𝑦𝑖 ) = ℎ𝑖 + 𝑥𝑖 · 𝐶 𝑝𝑖 d𝑇 (27)
𝑖=1 𝑇𝑟 𝑒 𝑓

Cependant, l’hypothèse posée pour l’étude est que l’écart à l’idéalité, (𝐻𝑉 − 𝐻 ∗ )𝑇,𝑃,𝑦𝑖 , est
négligé.

4
2 Modèle mathématique
2.1 Variables du problème
Précedemment, nous avons vu le système d’équations permettant de résoudre notre pro-
blème. Voyons les dimensions de notre système. Nous avons 2NC+3 équations pour 2NC+5
variables. Le système possède alors deux degrés de liberté. Pour saturer le système, nous fixons
deux variables : la puissance au condenseur 𝑄 𝐶 et la puissance au bouilleur 𝑄 𝐵 . De cette ma-
nière, nous retombons sur un système à 2CN+3 équations avec 2NC+3 inconnues. Ce système
va ensuite se répéter sur les 𝑁 𝑃 plateaux composant la colonne de distillation. Les inconnues
sont alors :
— 𝑉 𝑗 [mol/s] : Débit molaire vapeur sortant du plateau j ;
— 𝑦𝑖, 𝑗 [-] : Fraction molaire vapeur en constituant i au plateau j ;
— 𝑇 𝑗 [K] : Température au plateau j ;
— 𝑥𝑖, 𝑗 [i] : Fraction molaire liquide en constituant i au plateau j ;
— 𝐿 𝑗 [mol/s] : Débit molaire liquide sortant du plateau j.

2.2 Méthode de résolution

Les équations du problème, MESH et constitutives ne sont pas toutes linéaires. Pour cela,
il est nécessaire d’utiliser une méthode de résolution adapter aux équations non linéaires. La
méthode de Newton-Raphson a été choisie. Cette dernière permet de linéariser des équations
non linéaires. Son principe repose sur la recherche de la solution en prenant une tangente nulle.
Sa mise en oeuvre est assez simple. Son algorithme peut être représenté comme un cycle comme
sur la figure 2.

Figure 2 – Schéma de l’algorithme de Newton-Raphson

On peut remarquer que la méthode fait appel au vecteur des variables ainsi qu’au vecteur
des équations à chaque itération. De plus, ces derniers permettent le calcul du jacobien qui se
présente comme :

5
 𝜕𝐹1 𝜕𝐹1
© 𝜕 𝑋1 ... ... ... 𝜕 𝑋 (2∗𝑁 𝐶+3)∗𝑁 𝑃 ª
­ .. .. .. ®
­ . . . ®
­ ®
.. 𝜕𝐹𝑖 ..
𝑀=­
­ ®
. 𝜕𝑋𝑗 . ® (28)
­ ®
­ .. .. .
.. ®
­ . . ®
­ 𝜕𝐹(2∗𝑁 𝐶+3)∗𝑁 𝑃 𝜕𝐹(2∗𝑁 𝐶+3)∗𝑁 𝑃 ®
𝜕 𝑋1 ... ... . . . 𝜕 𝑋 (2∗𝑁 𝐶+3)∗𝑁 𝑃
« ¬

∀𝑖 𝜖 [1; (2 ∗ 𝑁𝐶 + 3) ∗ 𝑁 𝑃] et ∀ 𝑗 𝜖 [1; (2 ∗ 𝑁𝐶 + 3) ∗ 𝑁 𝑃]

La matrice de la jacobienne peut très vite avoir une taille conséquente. Par exemple, pour un
nombre de plateaux NP=30 et pour trois constituants NC=3, nous avons (2 ∗ 𝑁𝐶 + 3) ∗ 𝑁 𝑃 =
(2 ∗ 3 + 3) ∗ 30 = 270 lignes et colonnes soit 270 ∗ 270 = 72900 éléments.

Cependant, une grande partie de la matrice contient des éléments vides. En effet, elle se
présente comme une matrice tridiagonale par blocs (voir figure 3). Le comportement tridiagonale
par blocs de la matrice est du à l’organisation particulière des variables et des équations dans
leur vecteur. Pour l’exemple précédent, nous obtenons un système creux avec environ 97%
d’éléments nuls. Le système de remplissage de la matrice de la jacobienne et la méthode de
résolution du système linéaire vont alors avoir une influence sur la résolution globale. Deux
méthodes sont proposées : MRSL01 et MRSL21, les deux étant fournies.

Figure 3 – Matrice tridiagonale par blocs

6
2.2.1 La méthode MRSL01
Cette méthode nécessite la matrice de la jacobienne en paramètre d’entrée. De cette manière,
la matrice est remplie de façon classique, c’est-à-dire que toutes les lignes et toutes les colonnes
sont remplies, y compris les éléments nuls. Par conséquent, l’opération coûte du temps plus que
ce qu’il ne faudrait.

La méthode renvoie ensuite le pas 𝛿𝑋 𝑘 pour chaque variable à chaque plateau permettant
les variations.

2.2.2 La méthode MRSL21

Cette méthode doit être plus rapide que MRSL01. En effet, nous ne retrouvons pas la matrice
jacobienne dans son entièreté mais seulement les blocs où les éléments sont non nuls. Ces blocs
sont représentés en violet, vert et orange sur la figure 3. Il faut alors avoir une certaine rigueur sur
le remplissage mais cette étape est alors plus rapide puisque nous effectuons moins d’opérations.

Voilà les dimensions des différentes matrices :

— Matrice A, en violet : A(NLA,N1A,M) avec NLA = N1A =NC+2 et M = NP ;
— Matrice B, en vert : B(N,N,M) avec N = (2*NC+3) ;
— Matrice C, en orange : C(N,NCC,M) avec NCC = NC+2.
La méthode renvoie ensuite le pas 𝛿𝑋 𝑘 pour chaque variable à chaque plateau permettant
les variations des variables pour la recherche de la solution.

7
3 Organisation du programme
3.1 Structure globale du programme principale

Figure 4 – Logigramme non-exhaustif du programme principal

8
3.2 Organisation des vecteurs de variables
La méthode de Newton-Raphson est utilisée dans cette résolution. Les variables du problème
ainsi que toutes les équations sont alors stockées sous forme de vecteur. Le vecteur X contient
les 𝑁 𝑃 · (2𝑁𝐶 + 3) variables. Puis le vecteur F contient les 𝑁 𝑃 · (2𝑁𝐶 + 3) équations.

L’ordre dans lequel sont rangées les variables et les équations dans leur vecteur est important
pour obtenir une matrice jacobienne tridiagonale par blocs.

𝑉1 𝐵𝑀𝑇1
© 𝑦 ª © 𝐵𝑀 𝑃 ª
­ 1,1 ® ­ 1,1 ®
­ . ® ­ .. ®
­ .. ® ­ . ®
­ ® ­ ®
­ 𝑦 𝑁𝐶,1 ® ­ 𝐵𝑀 𝑃 𝑁𝐶,1 ®
­ ® ­ ®
­ 𝑇1 ® ­ 𝐵𝐸 1 ®
­ ® ­ ®
­ 𝑥 1,1 ® ­ 𝐵𝐸𝑄 1,1 ®
­ ® ­ ®
­ .. ® ­ .. ®
­ . ® ­ . ®
­ ® ­ ®
­ 𝑥 𝑁𝐶,1 ® ­ 𝐵𝐸𝑄 𝑁𝐶,1 ®
­ ® ­ ®
­ 𝐿1 ® ­ 𝑆𝑂 𝑀1 ®
­ . ® ­ .. ®
­ . ® ­ ®
­ . ® ­ . ®
­ . ® ­ . ®
­ .. ® ­ .. ®
­ ® ­ ®
­ 𝑉 ® ­ 𝐵𝑀𝑇 ®
­ 𝑗 ® ­ 𝑗 ®
­ 𝑦 ® ­ 𝐵𝑀 𝑃 ®
­ 1, 𝑗 ® ­ 1, 𝑗 ®
­ .. ® ­ .. ®
­ . ® ­ . ®
­ ® ­ ®
­ 𝑦 𝑁𝐶, 𝑗 ® ­ 𝐵𝑀 𝑃 𝑁𝐶, 𝑗 ®
­ ® ­ ®
𝑋 = ­­ 𝑇 𝑗 ®® (29) 𝐹 = ­­ 𝐵𝐸 𝑗 ®
® (30)
­ 𝑥 1, 𝑗 ® ­ 𝐵𝐸𝑄 1, 𝑗 ®
­ ® ­ ®
­ .. ® ­ .. ®
­ . ® ­ . ®
­ ® ­ ®
­ 𝑥 𝑁𝐶, 𝑗 ® ­ 𝐵𝐸𝑄 𝑁𝐶, 𝑗 ®
­ ® ­ ®
­ 𝐿𝑗 ® ­ 𝑆𝑂 𝑀 𝑗 ®
­ . ® ­ .. ®
­ . ® ­ ®
­ . ® ­ . ®
­ . ® ­ .. ®
­ .. ® ­ . ®
­ ® ­ ®
­ 𝑉 ® ­ 𝐵𝑀𝑇 ®
­ 𝑁 𝑃 ® ­ 𝑁 𝑃 ®
­ 𝑦 ® ­ 𝐵𝑀 𝑃 ®
­ 1,𝑁 𝑃 ® ­ 1,𝑁 𝑃 ®
­ .. ® ­ .. ®
­ . ® ­ . ®
­ ® ­ ®
­ 𝑦 𝑁𝐶,𝑁 𝑃 ® ­ 𝐵𝑀 𝑃 𝑁𝐶,𝑁 𝑃 ®
­ ® ­ ®
­ 𝑇𝑁 𝑃 ® ­ 𝐵𝐸 𝑁 𝑃 ®
­ ® ­ ®
­ 𝑥 1,𝑁 𝑃 ® ­ 𝐵𝐸𝑄 1,𝑁 𝑃 ®
­ ® ­ ®
­ .. ® ­ .. ®
­ . ® ­ . ®
­ ® ­ ®
­𝑥 𝑁𝐶,𝑁 𝑃 ® ­ 𝐵𝐸𝑄 𝑁𝐶,𝑁 𝑃 ®
« 𝐿𝑁𝑃 ¬ « 𝑆𝑂 𝑀𝑁 𝑃 ¬

3.3 Initialisation du vecteur des variables

3.3.1 Les bilans matières et d’énergie
Dans toute résolution, il est nécessaire d’initialiser les inconnues. C’est-à-dire indiquer au
programme les valeurs initiales sur lesquelles il doit s’appuyer pour lancer la résolution.

Pour réaliser notre initialisation, nous devons nous appuyer sur trois bilans effectués sur la

9
colonne.

Le bilan matière total

𝐹 = 𝐷 +𝑊 (31)
𝐷 = Distillat ; 𝑊 = Résidus

Le bilan matière partiel sur un constituant

𝐹 · 𝑧 𝑎𝑐 𝑒𝑡𝑜𝑛𝑒,1
´ = 𝐷 · 𝑦 𝑎𝑐 𝑒𝑡𝑜𝑛𝑒
´ + 𝑊 · 𝑥 𝑎𝑐 𝑒𝑡𝑜𝑛𝑒
´ (32)

Le bilan d’énergie au condenseur

𝑄 𝐶 = (𝐷 + 𝐿 1 ) · Δ𝐻𝑣𝑎 𝑝 (33)
À partir de ces trois équations, les débits liquides et vapeur peuvent être déterminés pour
chaque plateau.

3.3.2 Débit liquide L initial

Les contraintes sur la distillation sont variables et diverses. En effet, le système permet à
l’utilisateur d’avoir plusieurs alimentations à différents plateaux, divers soutirages liquides et/ou
vapeurs, des échanges de chaleur à chaque plateau. De cette façon, il faut adapter notre modèle
pour recevoir l’ensemble des contraintes.

Cependant, en ce qui concerne l’alimentation, nous résolvons notre problème avec une seule
alimentation. Cette alimentation peut toutefois être adaptée en terme de composition, de plateau
d’injection et de débit. Seul son état sera fixé à l’état de liquide bouillant. Cet état indique alors
que la température d’alimentation est la température de bulle à la composition de l’alimentation.
— Le débit liquide pour les plateaux inférieurs au plateau d’alimentation ( 𝑗 < 𝑁𝑎𝑙𝑖𝑚 )
  
𝑄𝐶 𝑧 𝑎𝑐 𝑒𝑡𝑜𝑛𝑒,1
´ − 𝑦 1,𝑐𝑜𝑛𝑑𝑒𝑛𝑠𝑒𝑢𝑟
𝐿𝑗 = −𝐹 · 1− (34)
Δ𝐻𝑣𝑎 𝑝 𝑥 1,𝑏𝑜𝑢𝑖𝑙𝑙𝑒𝑢𝑟 − 𝑦 1,𝑐𝑜𝑛𝑑𝑒𝑛𝑠𝑒𝑢𝑟

— Le débit liquide pour les plateaux supérieurs au plateau d’alimentation ( 𝑗 >= 𝑁𝑎𝑙𝑖𝑚 )
 
𝑄𝐶 𝑧 𝑎𝑐 𝑒𝑡𝑜𝑛𝑒,1
´ − 𝑦 1,𝑐𝑜𝑛𝑑𝑒𝑛𝑠𝑒𝑢𝑟
𝐿𝑗 = +𝐹· (35)
Δ𝐻𝑣𝑎 𝑝 𝑥1,𝑏𝑜𝑢𝑖𝑙𝑙𝑒𝑢𝑟 − 𝑦 1,𝑐𝑜𝑛𝑑𝑒𝑛𝑠𝑒𝑢𝑟

3.3.3 Débit vapeur V initial

Le débit vapeur va être considéré constant sur l’ensemble des plateaux.
𝑄𝐶
𝑉𝑗 = ∀ 𝑗𝜖 [1, 𝑁 𝑃] (36)
Δ𝐻𝑣𝑎 𝑝

3.3.4 Température et compositions liquides et vapeurs initiales

Pour les compositions liquides et vapeurs ainsi que pour la température, les conditions aux
limites sont connues. Grâce à cela, nous pouvons calculer ces variables pour les plateaux restants
par régression linéaire. Voilà comment sont calculées les températures et les compositions selon
le plateau :

10
Au bouilleur
— 𝑦𝑖,𝑏𝑜𝑢𝑖𝑙𝑙𝑒𝑢𝑟 = 𝑧𝑖,𝑎𝑙𝑖𝑚 ;
— 𝑥𝑖,𝑏𝑜𝑢𝑖𝑙𝑙𝑒𝑢𝑟 ;
— 𝑇 = 𝑇𝑟𝑜𝑠𝑒𝑒
´ .
Le calcul de la température de rosée ne peut pas s’effectuer par une méthode classique de
dichotomie puisque celle-ci est reliée directement aux compositions liquides. Cela forme alors
un système de NC+1 équations à NC+1 inconnues. Dans notre cas général, nous avons trois
constituants ce qui donne un système de la forme :

𝛾1 · 𝑃1𝑠𝑎𝑡
0= · 𝑥 1,1 − 𝑦 1,1 (37)
𝑃
𝛾2 · 𝑃2𝑠𝑎𝑡
0= · 𝑥 2,1 − 𝑦 2,1 (38)
𝑃
𝛾3 · 𝑃3𝑠𝑎𝑡
0= · 𝑥 3,1 − 𝑦 3,1 (39)
𝑃
𝑁𝐶
∑︁
0=1− 𝑥𝑖,1 (40)
𝑖=1
Pour résoudre ce système, il a été choisi d’utiliser une méthode de Newton-Raphon avec la
subroutine MRSL01 puisque le système est de petite taille. La solution de la résolution donne
alors directement les compositions liquides au bouilleur ainsi que la température de rosée.

Au condenseur
— 𝑥𝑖,𝑐𝑜𝑛𝑑𝑒𝑛𝑠𝑒𝑢𝑟 = 𝑧𝑖,𝑎𝑙𝑖𝑚 ;
— 𝑦𝑖,𝑏𝑜𝑢𝑖𝑙𝑙𝑒𝑢𝑟 : déterminés grâce à l’équilibre liq/vap ;
— 𝑇 = 𝑇𝑏𝑢𝑙𝑙𝑒 .
La température d’ébullition peut être calculée par une méthode de dichotomie. En utilisant
l’équilibre liquide/vapeur, les compositions vapeurs au condenseur peuvent être calculées avec
l’équation 41.

𝑥𝑖 · 𝛾𝑖 (𝑇, 𝑃, 𝑥𝑖 ) · 𝑃𝑖𝑠𝑎𝑡 (𝑇)

𝑦𝑖 = (41)
𝑃

11
4 Analyse du programme
L’analyse du programme a pour but de tester le programme en modifiant des paramètres
d’entrée via le fichier texte contenant les conditions opératoires. De plus, les résultats nous
permettent aussi d’étudier le comportement de la distillation.

Plusieurs paramètres seront mis en évidence :

— Le nombre de plateaux ;
— La pression opératoire ;
— La puissance au bouilleur ;
— La puissance au condenseur ;
— L’ajout de soutirages liquides et vapeurs.

Pour chaque changement de paramètre, outre la pression opératoire et la puissance au

condenseur, deux compositions en entrée seront étudiées (Acétone / Benzène / Chloroforme) :
— Composition de base : 0,6 / 0,3 / 0,1
— Autre composition : 0,1 / 0,2 / 0,7

4.1 Résultats avec les conditions de base

Établir le diagramme ternaire dans les conditions de base ne permet pas seulement d’observer
le comportement général de la distillation, mais permet aussi d’avoir une base de comparaison
lors des variations des paramètres qui sont effectuées par la suite.

Voici les conditions de simulation de base :

— Nombre de plateaux = 30 ;
— Puissance au bouilleur = 79420 J/s ;
— Puissance au condenseur = 62700 J/s ;
— Pression opératoire = 1 atm ;
— Soutirages = Aucun.

12
 Figure 5 – Diagramme ternaire avec les conditions de base pour deux types de composition

En rouge sont représentées les compositions vapeurs et en bleu les compositions liquides.
Nous pouvons remarquer que pour la composition de base, de l’acétone quasi-pur sort en tête
de colonne. En revanche, nous retrouvons une majorité de benzène au bouilleur. Pour l’autre
composition, nous avons une majorité de chloroforme en tête de colonne est un mélange des
trois constituants au bouilleur.

4.2 Influence du nombre de plateaux

L’hypothèse de base est que le liquide et la vapeur quittant un plateau sont en équilibre
thermodynamique. Cependant, ce n’est pas réellement le cas, de ce fait le nombre de plateaux
n’est que théorique et pas réel. Dans la réalité, si l’ont voulait une séparation parfaite, nous
devrions disposer d’une colonne à plateaux infinies. L’attente ici est alors de voir une précision
dans les résultats lorsque le nombre de plateaux augmente.

Les conditions de simulation sont :

— Nombre de plateaux = 30 / 50 ;
— Puissance au bouilleur = 79420 J/s ;
— Puissance au condenseur = 62700 J/s ;
— Pression opératoire = 1 atm ;
— Soutirages = Aucun.

13
 Figure 6 – Diagramme ternaire montrant l’influence du nombre de plateaux sur la distillation

De la même manière que précédemment, la vapeur est représentée en rouge et le liquide est
représenté en bleu. De plus, les triangles représentent la simulation pour 30 plateaux.

L’observation qui peut être faite est que le nombre de plateaux n’influe pas sur la finalité de
la distillation. En effet, les points se superposent et possèdent le même point initial et final. La
conclusion est qu’elle permet d’affiner les résultats tout au long de la colonne.

4.3 Influence de la pression

Le mélange ternaire que nous avons fourni un azéotrope binaire entre l’acétone et le chlo-
roforme. Théoriquement, il est possible de déplacer le point azéotropique afin d’augmenter une
pureté en jouant sur la pression. Cette simulation va permettre d’observer si c’est effectivement
le cas ici. Par ailleurs, nous pourrons observer le comportement général d’un changement de
pression sur la distillation.

Les conditions de simulation sont :

— Nombre de plateaux = 30 ;
— Puissance au bouilleur = 79420 J/s ;
— Puissance au condenseur = 62700 J/s ;
— Pression opératoire = 1 / 2 / 3 atm ;
— Soutirages = Aucun.

14
 Figure 7 – Diagramme ternaire montrant l’influence de la pression opératoire sur la distillation

La vapeur est représentée en rouge et le liquide en rouge.

— Croix : P = 3 atm ;
— Triangle : P = 2 atm ;
— Carré : P = 1 atm.
La principale observation est que la pression a une forte influence sur le domaine de
distillation. En effet, les courbes balayent un plus grand domaine lorsque la pression augmente.
Cependant, les compositions finales sont peu influencées. Cela se remarque au condenseur avec
de l’acétone quasi-pur. Au bouilleur, il y une légère augmentation de la composition en benzène
avec l’augmentation de la pression.

4.4 Influence de la puissance au bouilleur

La puissance au bouilleur va jouer sur la quantité de liquide vaporisé. Plus la quantité de
chaleur apportée est grande, plus la quantité de liquide vaporisée doit être grande. Cela va ainsi
jouer sur les compositions selon les volatilités des composants.

Les conditions de simulation sont :

— Nombre de plateaux = 30 ;
— Puissance au bouilleur = 75240 / 79420 / 83600 J/s ;
— Puissance au condenseur = 62700 J/s ;
— Pression opératoire = 1 atm ;

15
 — Soutirages = Aucun.

Figure 8 – Diagramme ternaire montrant l’influence de la puissance au bouilleur

La puissance au bouilleur a une forte influence sur les résultats de la distillation. La première
remarque est que plus la puissance apportée est faible (courbe du bas = puissance la plus faible),
plus la différence de compositions entre le bouilleur et le condenseur est faible. Cela s’explique
par la quantité de liquide vaporisée. Plus la quantité de chaleur est importante, plus la quantité
de liquide vaporisée augmente et crée donc des écarts significatifs entre le bouilleur et le
condenseur.

4.5 Influence de la puissance au condenseur

La puissance soutirée au niveau du condenseur va jouer sur la quantité de liquide en reflux.
Cela va ainsi jouer sur les compositions selon les volatilités des composants.
— Nombre de plateaux = 30 ;
— Puissance au bouilleur = 79420 J/s ;
— Puissance au condenseur = 54340 / 62700 J/s ;
— Pression opératoire = 1 atm ;
— Soutirages = Aucun.

16
 Figure 9 – Diagramme ternaire montrant l’influence de la puissance au condenseur

— Trait bleu = Liquide ; Qc = 62700 J/s ;

— Trait rouge = Vapeur ; Qc = 62700 J/s ;
— Trait rose = Vapeur ; Qc = 54340 J/s ;
— Trait vert = Liquide ; Qc = 54340 J/s.
Ce diagramme met en évidence l’importance du reflux sur les compositions finales. En
effet, plus la puissance retirée au condenseur est élevée, plus la quantité de vapeur condensée est
grande. De cette manière, le reflux se rapproche d’un condenseur total et donc plus le produit
récupéré en tête de colonne est pur en acétone (courbe bleue et rouge). En revanche, si l’on
condense moins de vapeur, la majorité de l’acétone est évacuée et le mélange s’enrichit en
benzène.

17
Conclusion
La mise en oeuvre d’une simulation de distillation montre à quel point cet outil est puissant.
Il permet d’étudier en détails de nombreux paramètres de la distillation en quelques instants.
Cela inclut alors un fort potentiel économique et un gain de temps non négligeable.

Malgré la puissance de cet outil, il est nécessaire d’avoir un code performant et complet
pouvant accueillir toutes les modifications que l’utilisateur veut apporter. Les axes d’amélio-
ration seraient alors de prendre en compte plusieurs alimentations à des états différents. Il est
possible aussi de prendre en compte les phénomènes de mécanique des fluides en compte, par
exemple les pertes de charges entre la tête et le bouilleur. En effet, l’hypothèse d’une pression
constante et uniforme dans toute la colonne reste assez forte.

D’un point de vue personnel, ce projet m’a réconforté dans mon choix de parcours pro-
fessionnel. La technologie et les outils numériques pour la simulation occupent une place
importante de l’industrie. Bien qu’il faille de l’expérience pour réaliser des modélisations com-
plètes et sûres, ce projet m’a permis de me familiariser encore plus avec les réflexes nécessaires
lors de la mise en place d’un code informatique. Enfin, ce projet m’a donné l’envie de réaliser
d’autres modélisations permettant de découvrir des opérations unitaires, de monter un code et
observer le comportement avec les résultats de la simulation.

18
Annexe 1 - Les différents fichiers utilisateur
Dimension de la distillation
L’utilisateur choisi le nombre de plateau et le nombre de constituants. Dans le fichier «di-
mensions.txt».

Les données thermodynamiques

Le fichier «DonneesThermo.txt» est là pour que l’utilisateur puisse renseigner l’ensemble
des données thermodynamiques concernant ses constituants.

Les conditions opératoires

Le fichier «Conditions Opératoires.txt» a pour objet de contenir les informations des condi-
tions opératoires. On retrouve notamment les pression ainsi que les différentes puissances.

Annexe 2 - Exécution du programme

Le programme comporte une certaines quantité de fichiers texte et de sous-programme.
Dans ce cas, il est compliqué de compiler correctement l’ensemble des fichiers. Afin de faciliter
la compilation de l’ensemble du programme, l’utilisation d’un MakeFile a été choisi. Il permet
de compiler l’ensemble des fichiers à l’aide d’une seule commande.

La commande à utiliser pour compiler le programme est «make» dans le terminal.

19