# Thermodynamique Moléculaire Appliquée au Génie des Procédés (HDR)

Patrice Paricaud

► To cite this version:

Patrice Paricaud. Thermodynamique Moléculaire Appliquée au Génie des Procédés (HDR). Génie des procédés. Université Paris 13 - Sorbonne Paris Cité, 2014. tel-04540270

HAL Id: tel-04540270

https://hal.science/tel-04540270

Submitted on 10 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
<figure><img src="image_1.png" /></figure>

ENSTA
ParisTech

Laboratoire de Chimie et Procédés
ENSTA-ParisTech

UNIVERSITÉ

PARIS 13
NORD

Curriculum Vitae

Mémoire de Travaux de Recherche

Présenté en vue de l'obtention de
l'Habilitation à Diriger des Recherches

Thermodynamique Moléculaire Appliquée
au Génie des Procédés

Patrice Paricaud

Mars 2014
A ma famille
# Table des Matières

## Curriculum Vitae

<table><tr><td>1.</td><td>Expérience Professionnelle</td><td>4</td></tr><tr><td>2.</td><td>Formation</td><td>4</td></tr><tr><td>3.</td><td>Enseignements</td><td>5</td></tr><tr><td>4.</td><td>Encadrements</td><td>6</td></tr><tr><td>5.</td><td>Publications</td><td>7</td></tr><tr><td>5.</td><td>Communications</td><td>9</td></tr><tr><td>6.</td><td>Activités Administratives</td><td>11</td></tr><tr><td>7.</td><td>Parcours et Synthèse des Travaux de Recherche</td><td>12</td></tr></table>

## Mémoire de Travaux de Recherche

<table><tr><td>1.</td><td>Introduction</td><td>15</td></tr><tr><td>2.</td><td>Equations d'Etat et Théories de Perturbation</td><td>17</td></tr><tr><td>2.1.</td><td>L'équation de van der Waals</td><td>18</td></tr><tr><td>2.1.1.</td><td>Terme répulsif de van der Waals</td><td>19</td></tr><tr><td>2.1.2.</td><td>Terme attractif de van der Waals</td><td>20</td></tr><tr><td>2.2.</td><td>L'équation d'état des sphères dures (Carnahan-Starling et BMCSL)</td><td>22</td></tr><tr><td>2.3.</td><td>Théories de Perturbation</td><td>23</td></tr><tr><td>2.3.1.</td><td>Théorie de Zwanzig</td><td>24</td></tr><tr><td>2.3.2.</td><td>Théorie de Barker-Henderson</td><td>24</td></tr><tr><td>2.3.3.</td><td>Théorie de perturbation pour les fluides multipolaires</td><td>25</td></tr><tr><td>2.3.4.</td><td>Théorie de Wertheim pour les fluides associatifs</td><td>26</td></tr><tr><td>2.3.5.</td><td>Théorie de Wertheim TPT1 pour les chaînes polymériques</td><td>28</td></tr><tr><td>2.4.</td><td>Travaux Réalisés dans le Développement d'Equations d'Etat</td><td>29</td></tr><tr><td>2.4.1.</td><td>Modification de la théorie de Barker - Henderson</td><td>29</td></tr><tr><td>2.4.2.</td><td>Développement d'une équation d'état pour les sphères dures non additives</td><td>30</td></tr><tr><td>2.5.</td><td>Modification du terme associatif de SAFT. Application aux acides carboxyliques</td><td>31</td></tr><tr><td>3.</td><td>Thermodynamique des Solutions de Polymère</td><td>35</td></tr><tr><td>3.1.</td><td>Types de comportement de phase</td><td>35</td></tr><tr><td>3.1.1.</td><td>Comportement de type I</td><td>35</td></tr><tr><td>3.1.2.</td><td>Comportement de type II</td><td>36</td></tr><tr><td>3.1.3.</td><td>Comportement de type V</td><td>37</td></tr><tr><td>3.1.4.</td><td>Comportement de type IV</td><td>38</td></tr><tr><td>3.1.5.</td><td>Comportement de type III</td><td>40</td></tr><tr><td>3.1.6.</td><td>Comportement de type VI</td><td>41</td></tr><tr><td>3.1.7.</td><td>Comportement de phase des solutions de polymères</td><td>43</td></tr><tr><td>3.2.</td><td>Modélisation des équilibres de phases de solutions de polymère</td><td>44</td></tr><tr><td>3.2.1.</td><td>Solutions de polyéthylène</td><td>44</td></tr><tr><td>3.2.2.</td><td>Solubilités de gaz dans le polyéthylène</td><td>45</td></tr><tr><td>4.</td><td>Systèmes Polydispersés et Thermodynamique Continue</td><td>47</td></tr><tr><td>4.1.</td><td>Distributions discrètes et pseudoconstituants</td><td>48</td></tr><tr><td>4.2.</td><td>Distributions continues</td><td>49</td></tr><tr><td>4.3.</td><td>Description de la méthode des moments</td><td>49</td></tr><tr><td>5.</td><td>Propriétés Thermodynamiques des Solutions d'Electrolytes</td><td>54</td></tr><tr><td>5.1.</td><td>Modélisation avec le modèle SAFT-VRE</td><td>55</td></tr><tr><td>5.1.1.</td><td>Systèmes eau-sel-hydrocarbures</td><td>56</td></tr><tr><td>5.1.2.</td><td>Modélisation des équilibres de phases liés au procédé SI</td><td>57</td></tr><tr><td>5.1.3.</td><td>Modélisation des semi-clathrates de gaz</td><td>57</td></tr><tr><td>5.2.</td><td>Développement d'une version électrolyte de PC-SAFT</td><td>62</td></tr></table>
# Curriculum Vitae (Mars 2014)

## Patrice Paricaud

Enseignant-chercheur depuis le 1er Décembre 2004, à l'Unité d'enseignement et de recherche de Chimie et Procédés (UCP), Ecole Nationale Supérieure de Techniques Avancées (ENSTA-ParisTech).
Adresse professionnelle: 828 Boulevard des Maréchaux, 91762 Palaiseau Cedex, France.

Téléphone professionnel: +33 (0)1 81 87 20 26.

Page web: http://ucp.ensta-paristech.fr/equipe/paricaud.php

Email: patrice.paricaud@ensta-paristech.fr

Né le 1er Janvier 1976 à Guéret.

Nationalité Française, Marié.

## 1. Expérience Professionnelle

Depuis Déc. 2004. Enseignant-chercheur à l'Unité d'enseignement et de recherche Chimie et Procédés (UCP), Ecole Nationale Supérieure de Techniques Avancées (ENSTA-ParisTech).

Fév. 2003-Oct. 2004. Postdoc à l'université de Vanderbilt, département de génie chimique, Nashville, USA, sous la direction de Peter T. Cummings. Sujet: "Développement d'un modèle moléculaire polarisable pour l'eau. Simulation de Monte-Carlo (GEMC) et dynamique moléculaire".

Juillet-Nov. 1998. Stage ingénieur à l'IFP Energies Nouvelles, Rueil-Malmaison, sous la direction d'Alain Barreau. "Modélisation des propriétés thermodynamiques des gaz à condensât".

Juillet 1997. Stage au laboratoire SEPA, AREVA, Bessines sur Gartempe. Gestion des produits du laboratoire sous une base de données Microsoft Access.

## Prix & Distinctions

2006. Prix récompensant les cinq publications les plus significatives de l'année de l'Université de Vanderbilt, pour l'article : Paricaud et al., J. Chem. Phys., 122, 244511, (2005).

2003. Dudley Newitt Prize : prix annuel récompensant la meilleure thèse en modélisation du département de génie chimique de l'Imperial College de Londres.

## 2. Formation

2000-2003. Doctorat à l'Impérial College de Londres (R.U.), département de génie chimique. Titre: "Understanding the fluid phase behaviour of polymer systems with the SAFT theory". Thèse sous la direction de George Jackson, et soutenue le 23 Janvier 2003.

1996-1999. Ingénieur et DEA en génie chimique (mention bien), École Nationale Supérieure des Industries Chimiques (ENSIC), Nancy, France. Stage de DEA en coopération avec IFP Energies Nouvelles et encadré par Jean-Noël Jaubert et Roland Solimando. Sujet : "Représentation de la fraction lourde des gaz à condensât à l'aide de fonctions de distribution".

1994-1996. Classes préparatoires mathématiques supérieures et spéciales P', Lycée Gay-Lussac, Limoges.
## Langues

Anglais courant lu, écrit, parlé : trois ans au Royaume-Uni, et un an et demi aux USA.
Allemand: notions (cinq années d'études).

## Informatique

Langages de programmation: Fortran, Matlab, notions de C, Turbo Pascal, et Basic, et de programmation en parallèle avec MPI.

Logiciels de simulation de procédés: PROSIM, ASPEN, CHEMCAD, G-PROMS.

Logiciels de Mathématiques: Maple, Mathcad, GAMMS.

Logiciel de calculs *ab initio*: DMOL3 (Materials Studio), Turbomole, Gaussian03.

## Connaissances en simulation moléculaire et en calculs *ab initio*

Utilisation et développement de code fortran de simulations de Monte-Carlo (NVT, NPT, et ensemble de Gibbs), de la dynamique moléculaire à l'équilibre ; utilisation et développement de champs de forces classiques polarisables pour l'eau. Calculs *ab initio* sous Materials Studio (DMOL3), turbomole et Gaussian03.

* Formations suivies en simulation moléculaire :

Sept. 2013. Vienne, Autriche. Formation aux logiciels de Materials Design (Medea, Gibbs, VASP...).
Nov. 2006. Workshop au CECAM, Lyon, France. Sujet: "simulation de la matière à l'échelle nano en utilisant la théorie de la fonctionnelle de la densité, les pseudopotentiels et les ondes planes".

Août 2003. Workshop au laboratoire national d'Oak Ridge. Sujets: "techniques avancées en simulation moléculaire: *ab-initio* MD (Car-Parrinello Molecular Dynamics), DFT (Density Functional Theory,), and DPD (Dispersive Particle Dynamics)". "Introduction aux logiciels NWChem (*ab-initio* MD) et LAMMS (MD classique)".

Juillet 2002. Ecole d'été CCP5, Kings College, Londres, R.U. "Introduction à la simulation moléculaire (Monte-Carlo et MD)".

# 3. Enseignements

Les heures d'enseignement indiquées entre parenthèses correspondent aux heures cumulées depuis 2000 (thèse) jusqu'à l'année 2013-2014.

* Cours en première année de cursus ingénieur ENSTA

**IN103 : Introduction à Matlab** (201heures). Le but de ce cours est d'apprendre à programmer sous Matlab. Le cours est divisé en deux parties : une partie cours + TD permettant d'apprendre les principales fonctions de Matlab, et une partie dédiée à la réalisation d'un projet. Ce projet consiste à la programmation du modèle thermodynamique NRTL, à l'optimisation des paramètres binaires NRTL, puis au dimensionnement d'une colonne de distillation par minimisation des coûts.

* Cours en deuxième année de cursus ingénieur ENSTA

**CB203 : Introduction au génie des procédés** (32heures). Responsable du cours. Cours qui permet d'introduire les bases du génie chimique aux étudiants. Il est constitué de cours magistraux et de travaux dirigés. Participation aux TDs et au projet (exercices + projet sous aspen plus).
**CB203 : Introduction à la simulation moléculaire** (6heures). Ce cours présente les principales techniques de : simulation moléculaire classique (Monte-Carlo, dynamique moléculaire), les algorithmes associés et les applications.

* Cours en troisième année de cursus ingénieur ENSTA et masters.

**A5-1 : Equations d'état et thermodynamique moléculaire** (15heures). Ce cours retrace l'historique du développement des équations d'état: de l'équation de van der Waals à la théorie SAFT. Les différents types d'interactions intermoléculaires et la classification de van Konynenburg et Scott des six principaux types de comportement de phase en mélanges binaires sont également décrits.

**C4-2. La filière hydrogène:** procédés industriels de production de l'hydrogène (34h). Vaporeformage du méthane, procédés par oxydation partielle. Gazéification du charbon. Pyrolyse et gazéification de la biomasse. Thermolyse de l'eau. Electrolyse alcaline. Electrolyse à haute température. Cycles thermochimiques.

**B3-1. Simulation de Procédés Batch** (34heures). Modélisation des opérations unitaires batch en génie chimique. Mise en équations de systèmes dynamiques. Modélisation des réacteurs et les colonnes de distillation fonctionnant en mode batch. Travaux dirigés et projet de simulation de procédés batch sous matlab.

**Semaine de découverte des milieux industriels:** responsable depuis 2005 de l'organisation d'une semaine de visites de sites industriels chimiques pour les étudiants ENSTA en 3ème année. Plusieurs grands groupes visités depuis 2005, un peu partout en France: Arkema, Ineos, Exxon, Total, LyondellBasell, BASF, Arcelor-Mittal, Areva, EDF (centrale nucléaire), Air Liquide, Sanofi-Aventis, Veolia, Dalkia, Rouillier (engrais) ... L'organisation de ces visites m'a permis d'avoir un excellent aperçu des usines chimiques présentes en France.

**Master des énergies renouvelables (REST)** de l'école Polytechnique (Palaiseau). Cours sur "la filière hydrogène" : procédés de production de H₂ (32h).

* Autres enseignements

2011. Semaine BEST (Board of European Students of Technology), 3heures de cours sur les "procédés industriels de fabrication de l'hydrogène".
2002-2003 Cours de Français scientifique donnés à l'Imperial College de Londres (40h).

# 4. Encadrements

* Thèses

Oct. 2013 à Oct. 2016. Participation à l'encadrement de la thèse de Alfonso Gonzalez (30%), sous la direction de C. Coquelet, CTP, Mines ParisTech et de A. Chapoy, université Heriot-Watt d'Edinbourg.

Oct. 2008 à Jan. 2012. Co-encadrement de la thèse de Thi Thu Suong Ho (50%). Sujet : "Développement d'une équation d'état pour les solutions d'électrolytes". Thèse dirigée par Walter Fürst (ENSTA). Un article a été soumis en 2014. Un autre est en cours de rédaction. Thèse soutenue le 10 Janvier 2012.

2011. Participation à l'encadrement de la thèse de Wassila Bouchaffaa (10%) pour la partie modélisation. Sujet : "Etude des conditions de formation des semi-clathrates hydrates de gaz". Thèse dirigée par Didier Dalmazzone et soutenue en 2011. Un article a été récemment soumis.
➤ Postdocs

Sept. 2013 à Aout 2014. Co-encadrement (50%) de Ayoko Fukumoto, postdoc financé par la chaire Areva. Coopération avec Didier Dalmazzone et Walter Fürst (UCP, ENSTA-ParisTech). Sujet : études expérimentales et de modélisation des hydrates et semi-clathrates d'hydrogène. Application au stockage de l'hydrogène et à la séparation de gaz. Un article a été soumis en 2014.

Mars 2013 à Février 2014. Co-encadrement de Jiri Janecek (50%), Postdoc financé par Mines ParisTech. Coopération avec C. Coquelet de Mines-ParisTech et M. Dicko de l'université Paris XIII. Sujet : développement d'équations d'état de type crossover, vérifiant les lois d'échelles autour des points critiques. Un article a été soumis en 2014.

Sept. 2011 à Février 2013. Encadrement à 100% de Jiri Janecek. Postdoc financé par l'ANR (ANR-09-CP2D-10-03). Sujet : "Modélisation des équilibres de phases des mélanges impliquant les acides carboxyliques, avec la prise en compte de la formation de dimères en phase gazeuse". Trois publications ont été faites.

Jan. 2010-Jan. 2011 Encadrement à 100% de Jinlong Li. Postdoc financé par l'ANR (ANR-09-CP2D-10-03). Sujet : "Modélisation des équilibres de phases de molécules extraites des bio-huiles, en utilisant les approches COSMO-RS et COSMO-SAC". Deux articles sur ce sujet ont été publiés.

➤ Stages de Master et de fin d'études

Fév. à Juillet. 2008. Encadrement à 100% de Laetitia Tazi. Stage de fin d'étude financé par le CEA. Sujet : "modélisation des équilibres de phases du mélange Hix lié au procédé iode-soufre pour la production industrielle d'hydrogène". Ce stage a débouché à une publication et à plusieurs communications.

Fév. à Juillet. 2006. Co-encadrement de Mohammed Ellaite (50%). Sujet : "modélisation des coefficients de partages de molécules d'arome dans le mélange eau-éthanol". Ce stage de master a débouché à une publication et un proceeding (congrès SFGP 2009).

**5. Publications**

26 articles ; facteur h de 13; 567 citations dont 537 sans autocitation (web of knowledge, Mars 2014).

➤ Articles dans des revues internationales à comité de lecture

28) A. Fukumoto, P. Paricaud, D. Dalmazzone, W. Bouchafaa, T. T. -S. Ho, W. Fürst, Modeling the dissociation conditions of carbon dioxide + TBAB, TBAC, TBAF, and TBPB semiclathrate hydrates, J. Chem. Eng. Data, soumis (2014).

27) J. Janecek, P. Paricaud, M. Dicko, C. Coquelet, Modified Kiselev crossover approach applied to the CPA equation of state, Chem. Eng. Res. Des., soumis (2014.)

26) J. Janecek, P. Paricaud, Size distribution of associated clusters in liquid alcohols: Interpretation of simulation results in the frame of SAFT approach, J. Chem Phys., (2013).

25) J. Janecek, P. Paricaud, Influence of Cyclic Dimer Formation on the Phase Behavior of Carboxylic Acids. II Cross associating Systems, J. Phys. Chem. B., 117, 9430-9438 (2013).
24) M. Nala, E. Auger, I. Gedik, N. Ferrando, M. Dicko, P. Paricaud, F. Volle, J. P. Passarello, J.-C. de Hemptinne, P. Tobaly, P. Stringari, C. Coquelet, D. Ramjugernath, P. Naidoo, R. Lugo, Vapour Liquid Equilibrium (VLE) for the Systems Furan + n-Hexane and Furan + Toluene. Measurements, data treatment and modelling using molecular models, Fluid. Phase. Equilib., 337, 234-245 (2013).

23) E. El Ahmar, A. Valtz, P. Paricaud, C. Coquelet, L. Abbas, W. Rached, Vapour-Liquid Equilibrium of Binary Systems Containing Pentafluorochemicals from 363 to 413 K: Measurement and Modelling with Peng-Robinson and three SAFT-like Equations of States, Int. J. Refrig., 35, 2297, 2310 (2012).

22) J. Janecek, P. Paricaud, Influence of Cyclic Dimer Formation on the Phase Behavior of Carboxylic Acids, J. Phys. Chem. B., 116, 7874-7882 (2012).

21) J. Li, and P. Paricaud, Application of the COSMO-RS and COSMO-SAC models for the predictions of partition coefficients, vapor-liquid and liquid-liquid equilibria of bio-oil related mixtures, Energy & Fuels, 26, 3756-3768 (2012).

20) J. Li, C. He, C. Peng, H. Liu, Y. Hu, P. Paricaud, "Modeling of the thermodynamic properties of aqueous ionic liquid solutions with an equation of state for square well chain fluid with variable range", Ind. Eng. Chem. Res., 50, 7027-7040 (2011).

19) D. Fan, J. Li, J. Shi, C. Peng, H. Liu, Y. Hu, and P. Paricaud, "Vapor-Liquid Equilibria in the Propyl Acetate + Ethanoic Acid Binary System from (323.15 to 353.15) K: Measurement with a Static Method and Modeling with the NRTL, Wilson, UNIQUAC, and COSMO-SAC Approaches", J. Chem. Eng. Data, 56, 1323-1329 (2011).

18) P. Paricaud, "Modeling of the dissociation conditions of salt hydrates and gas semi-clathrate hydrates. Application to lithium bromide, hydrogen iodide, and tetra-n-butylammonium bromide + carbon dioxide systems", J. Phys. Chem. B, 115, 288-299 (2011).

17) P. Paricaud, L. Tazi, J.M. Borgard, "Modeling the phase equilibria of the H1x mixture using the SAFT-VRE equation of state: binary systems" Int. J. Hydrogen Energy, 35, 978-991 (2010).

16) X. Courtial, C.-B Soo, C. Coquelet, P. Paricaud, D. Ramjugernath and D. Richon, "Vapor-liquid equilibrium in the n-butane + methanol system, measurement and modeling from 323.2 to 443.2 K ", Fluid Phase Equilibr., 277, 152-161 (2009).

15) P. Paricaud, "Phase equilibria in polydisperse nonadditive hard-sphere systems", Phys. Rev. E., 78, 021202 (2008).

14) V. Athès, P. Paricaud, M. Ellaite, I. Souchon, W. Fürst, "Vapour-liquid equilibria of aroma compounds in hydroalcoholic solutions: Measurements with a recirculation method and modelling with the NRTL and COSMO-SAC approaches", Fluid Phase Equilib., 265, 139-154 (2008).

13) P. Paricaud, A. Galindo, G. Jackson, "Examining the effect of chain length polydispersity on the phase behavior of polymer solutions with the statistical associating fluid theory (Wertheim TPT1) using discrete and continuous distributions", J. Chem. Phys., 127, 154906 (2007).

12) J. L. Rivera, Francis W. Starr, P. Paricaud, and Peter T. Cummings, "Polarizable Contributions to the Surface Tension of Liquid Water", J. Chem. Phys., 125, in press (2006).

11) A.J. Haslam, N. von Solms, C.S. Adjiman, A. Galindo, G. Jackson, P. Paricaud, M.L. Michelsen and G.M. Kontogeorgios, "Predicting enhanced absorption of light gases in polyethylene using simplified PC-SAFT and SAFT-VR", Fluid Phase Equilib. 243, 74-91 (2006).

10) P. Paricaud, "A General Perturbation Approach for Equation of State Development: Applications to Simple Fluids, Ab initio Potentials and Fullerenes. ", J. Chem. Phys, 124, 154505 (2006).

9) P. Paricaud, M. Predota, A. A. Chialvo, P. T. Cummings, "From dimer to condensed phases at extreme conditions: Accurate predictions of the properties of water by a Gaussian charge polarizable model", J. Chem. Phys., 122, 244511 (2005).

8) A. Valtz, A. Chapoya, C. Coquelet, P. Paricaud, D. Richon, "Vapour-liquid equilibria in the carbon dioxide-water system, measurement and modelling from 278.2 to 318.2K", Fluid Phase Equilib., Fluid Phase Equilib., 226, 333-344 (2004).

7) P. Paricaud, S. Varga, P. T. Cummings, G. Jackson, "Effect of polymer chain-length polydispersity on the phase behavior of model athermal mixtures of colloids and flexible self-excluding polymers", Chem. Phys. Lett., 398, 489-494 (2004).

6) P. Paricaud, A. Galindo, G. Jackson, "Modeling the Cloud Curves and the Solubility of Gases in Amorphous and Semicrystalline Polyethylene with the SAFT-VR Approach and Flory Theory of Crystallization", Ind. Eng. Chem. Res., 43, 6871-6889 (2004).
5) P. Paricaud, A. Galindo, and G. Jackson, "Understanding liquid-liquid immiscibility and LCST behaviour in polymer solutions with a Wertheim TPT1 description", Mol. Phys., 101, 2575-2600 (2003).

4) B. H. Patel, P. Paricaud, A. Galindo, and G. C. Maitland, "Prediction of the salting out effect of strong electrolytes on water + oil systems", Ind. Eng. Chem. Res., 42, 3809-3823 (2003).

3) P. Paricaud, S. Varga, and G. Jackson, "Study of the demixing transition in model athermal mixtures of colloids and flexible self-excluding polymers using the thermodynamic perturbation theory of Wertheim", J. Chem. Phys., 118, 8525-8536 (2003).

2) P. Paricaud, A. Galindo, and G. Jackson, "Recent advances in the use of the SAFT approach in describing electrolytes, interfaces, liquid crystals and polymers", Fluid Phase Equilib., 194, 87-96 (2002).

1) J.N. Jaubert, R. Solimando, P. Paricaud, and A. Barreau, "Use of distribution functions: A useful tool to calculate the properties of condensate gases", Ind. Eng. Chem. Res., 39, 5029-5036 (2000).

Chapitre de livre

M. Dicko, C. B. Soo, C. Coquelet, P. Paricaud, Etude et détermination des points critiques : mesures expérimentales et modélisation, Récents Progrès en Génie des Procédés - Numéro 98, 2007 ISBN 2-910239-72-1, Ed. SFGP, Paris, France.

**5. Communications**

Conférences internationales avec comité scientifique

26) J. Janecek, P. Paricaud, M. Dicko, C. Coquelet, Simplified crossover CPA equation of state for binary mixtures, 27th European Symposium on Applied Thermodynamics (ESAT), 6-9 Juillet 2014, presentation orale acceptée.

25) J. Janecek, P. Paricaud, O. Baudouin, E. Auger, C. Coquelet, J.-P. Passarello, R. Lugo, J.-C. de Hemptinne, Influence of Cyclic Dimer Formation on the Equilibrium Behavior of Carboxylic Acids. From Theory to Process Simulation. 27th European Symposium on Applied Thermodynamics (ESAT), 6-9 Juillet 2014, presentation orale acceptée.

24) A. Fukumoto, P. Paricaud, D. Dalmazzone, W. Bouchafaa, T. T. -S. Ho, W. Fürst, Experimental measurements and modeling of the dissociation conditions of gas semiclathrate hydrates formed with TBAB, TBAC, TBAF, and TBPB salts, 27th European Symposium on Applied Thermodynamics (ESAT), 6-9 Juillet 2014, presentation orale acceptée.

23) R. Lugo, P. Paricaud, C. Coquelet, J. P. Passarello, J. C. de Hemptinne, Molecular-based predictive models applied to high-added value oxygen-bearing molecules and their mixtures; XIX International Conference on Chemical Thermodynamics in Russia, presentation orale, Moscow, Russie, 24-28 Juin 2013.

22) J. Jacenek, P. Paricaud, "Influence of dimer formation on the phase behaviour of carboxylic acids", Inmother 2012, Lyon, 19-20 Mars 2012, poster.

21) M. Dicko, P. Paricaud, C. Coquelet, Thermodynamic 2011 "Choosing a crossover approach for industrial applications", Thermodynamics 2011, Athène, Grèce, 1-3 Sept. 2011, présentation orale.

20) P. Paricaud, W. Bouchafaa, D. Dalmazzone, "Modeling of the dissociation conditions of salt hydrates and gas semi-clathrate hydrates. Application to lithium bromide, hydrogen iodide hydrates, and mixed hydrates of gas + Tetra-n-butylammonium bromide", 25th European Symposium on Applied Thermodynamics (ESAT) 2011. Saint-Pétersbourg, Russie, 24-27 Juin, présentation orale.

19) E. Auger, J-P. Passarello, P. Paricaud, P. Tobaly, F. Volle, "Including polarizability effects in GC-PPC-SAFT: application to alkyl-ether containing systems", 25th European Symposium on Applied Thermodynamics (ESAT) 2011. Saint-Pétersbourg, Russie. 24-27 Juin, poster.

18) P. Paricaud, L. Tazi, J.-M. Borgard, "Modeling the phase equilibria of the HIx mixture using the SAFT-VRE equation of state", SAFT 2010, Barcelone, 19-21 Sept. 2010, Poster.
17) X. Courtial, C.-B. Soo, C. Coquelet, P. Paricaud, D. Ramjugernath, D. Richon, "Phase diagram of methanol and alkane binary systems in the field of natural gas fractionators: experimental work and modeling", AIChE annual meeting, Philadelphia, USA, 20th Nov. 2008. Présentation orale.

16) P. Paricaud, M. Predota, A. A. Chialvo, P. T. Cummings. "From dimer to condensed phases at extreme conditions: accurate prediction properties of water by a Gaussian charge polarisable model", Joint European Thermodynamics Conference (JETC) IX, Saint-Étienne, France, 12-15 June 2007. Présentation orale.

15) P. Paricaud, B. Patel, A. Galindo, G. Jackson: "Recent advances in the use of the SAFT approach to describe the phase behaviour of associating molecules, electrolytes and polymers". (Présentation orale + poster).

14) P. Paricaud. "A General Perturbation Approach for Equation of State Development: Applications to Simple Fluids, Ab initio Potentials and Fullerenes", Thermodynamics 2007, Institut Français du Pétrole (IFP), Rueil-Malmaison, France. Présentation orale.

13) J. L. Rivera, Francis W. Starr, P. Paricaud, and Peter T. Cummings, "Polarizable Contributions to the Surface Tension of Liquid Water". AIChE annual meeting, San Francisco, Nov. 2006. Présentation orale.

12) P. Paricaud, Y. Peng, C. McCabe, "Modelling of the solubility of gases in semi-crystalline polyethylene", 4th Workshop of molecular simulation of amorphous polymers and materials, Clermont-Ferrand, France, 16 Juin 2005. Présentation orale.

11) N. von Solms, A. J. Haslam, P. Paricaud, C. S. Adjiman, A. Galindo, G. Jackson, M. L. Michelsen, and G. M. Kontogeorgis, "Phase equilibrium in the polyethylene manufacturing process - a Comparison of PC-SAFT and SAFT-VR approaches". Thermodynamics 2005, Sesimbra, Portugal, Avril 2005. Présentation orale.

10) Y. Peng, C. McCabe, P. Paricaud. "Modeling the Solubility of Gases in Semi-crystalline Polyethylene with the SAFT-VR Theory". Polymer workshop, CECAM, ENS Lyon. Sept. 2005. Présentation orale.

9) P. Paricaud, A. Galindo, G. Jackson. "Modeling the solubility of gas in amorphous and semi-crystalline polyethylene with SAFT-VR", AIChE annual Meeting, San Francisco, Nov. 2003, Poster.

8) P. Paricaud, A. Galindo, G. Jackson. "Modeling the solubility of hydrocarbons in amorphous and semi-crystalline polyethylene with SAFT-VR", FOMMS 2003, Keystone, Colorado, USA, July 2003. Poster.

7) P. Paricaud, S. Varga, A. Galindo, and G. Jackson., "Understanding the thermodynamics and fluid phase behaviour of polymer-colloid and polymer-solvent systems". Thermodynamics 2003, Cambridge, R. U., Avril 2003. Présentation orale.

6) P. Paricaud, A. Galindo, and G. Jackson. "Understanding the thermodynamics and fluid phase behaviour of polymer solutions". 16th European Conference on Thermophysical Properties, Imperial College, Londres, R.U., Sept 2002. Présentation orale.

5) G. Jackson, P. Paricaud, A. Galindo, "Modeling the thermodynamics and phase equilibria of hydrocarbons in polyethylene", 6th Liblice conference on the Statistical Mechanics of Liquids, Liblice, République Tchèque, Juin 2002. Présentation orale.

4) B. H. Patel, P. Paricaud, A. Galindo, G. Jackson. "Predicting the salting out of alkanes in aqueous solutions with strong electrolytes", 6th Liblice conference on the Statistical Mechanics of Liquids, Liblice, République Tchèque, Juin 2002. Poster.

3) G. Jackson, P. Paricaud, A. Galindo, "Recent advances in the use of the SAFT approach in describing electrolytes interfaces, liquid crystals and polymers", PPEPPD (Properties and Phase Equilibria for Product and Process Design), Kurashiki, Japon, May 2001. Présentation orale.

2) P. Paricaud, A. Galindo, G. Jackson, "Understanding the fluid phase equilibria in n-alkanes - polyethylene systems". Thermodynamics 2001, Université de Bristol, R.U., avril 2001. Présentation orale.

1) B. Patel, P. Paricaud, A. Galindo, G. Jackson, "Predicting the anomalous temperature dependence of n-alkanes in aqueous solution using SAFT-VR approach", Thermodynamics 2001, Université de Bristol, R.U., avril 2001. Poster.
➤ Conférences internationales invitées

P. Paricaud, "Challenges and future directions for the SAFT approach", SAFT 2010, Barcelone, 19-21 Sept. 2010, Présentation orale.

P. Paricaud, L. Tazi, J.-M. Borgard, "Modeling the phase equilibria of the H1x mixture using the SAFT-VRE equation of state", talk. Colloquium to promote experimental work in thermophysical properties for science research and industry. 60th birthday of Prof. D. Richon. Mines-ParisTech, Paris, Sept. 2009. Présentation orale.

➤ Conférences nationales avec comité scientifique

V. Athès, W. Khalil, P. Paricaud, W. Furst, I. Souchon, "Effect of thermodynamic properties estimation on simulation of batch distillation: a case study using aroma compounds in hydro-alcoholic solutions", 12ème congrès de la Société Française de Génie des Procédés (SFGP2009), Marseille, Oct. 2009.

M. Dicko, C.-B. Soo, C. Coquelet, P. Paricaud, "Etude et détermination des points critiques: mesures expérimentales et modélisation", 12ème congrès de la Société Française de Génie des Procédés (SFGP2009), Marseille, Oct. 2009. Présentation orale.

P. Paricaud, G. Jackson, A. Galindo, "Modeling the fluid phase equilibria of polyethylene solutions". 8ème congrès de la Société Française de Génie des Procédés, Nancy, Oct. 2001. Présentation orale.

➤ Colloques, séminaires

P. Paricaud "Développements récents et perspectives en équations d'état. Journée du GDR thermodynamique, Paris, 11 Déc. 2012.

P. Paricaud, "Récents développements en modélisation des solutions d'électrolytes ", Colloque SFGP, Paris, 13 juin 2012. Présentation orale.

P. Paricaud, W. Bouchafaa, D. Dalmazzone, Modélisation Thermodynamique des Conditions de Formation des Hydrates de Sels et des Semi-clathrates de Gaz, Journée SFG sur l'utilisation simultanée de la thermodynamique des équilibres de phases et la thermodynamique énergétique, Nancy, 6 Jan 2011.

T. Thu-Suong, P. Paricaud, W. Fürst, Développement d'une Equation d'Etat pour les Solutions d'Electrolytes, Journée SFG sur l'utilisation simultanée de la thermodynamique des équilibres de phases et la thermodynamique énergétique, Nancy, 6 Jan 2011.

J. Li, P. Paricaud, Modélisation des équilibres de phases de mélanges de constituants issus de biohuiles en utilisant les approches COSMO-RS et COSMO-SAC, Journée SFG sur l'utilisation simultanée de la thermodynamique des équilibres de phases et la thermodynamique énergétique, Nancy, 6 Jan 2011.

P. Paricaud, J.-P. Passarello, "Apport de la simulation moléculaire au développement et la validation des équations d'état". Présentation. Journée SFGP sur l'utilisation de la thermodynamique en génie des procédés, IFP Energies Nouvelles, Rueil-Malmaison, Jan 2009. Présentation orale.

## 6. Activités Administratives

Membre du comité de pilotage du groupe "thermodynamique" de la Société Française de Génie des Procédés (SFGP). Responsable de l'organisation d'une journée thématique SFGP dédiée à "l'enseignement de la thermodynamique en génie des procédés" (8 Juillet 2009, ENSTA-ParisTech). Participation à l'organisation d'une journée SFGP dédiée à la thermodynamique des électrolytes (2012).

Membre du comité d'organisation de la conférence internationale "Thermodynamics 2007".

2004-2010. Responsable du site internet du laboratoire de Chimie et Procédés (conception et mises à jour du site web).
Membre du jury / rapporteur des thèses suivantes:

5) Rapporteur de la thèse de Hoda Tahani, soutenue le 12 Mars 2012, Institute of Petroleum Engineering, Heriot-Watt.

4) Membre du jury de la thèse de Peyman Memari, soutenue le 16 Février 2011. Laboratoire de Chimie Physique d'Orsay, Université de Paris-Sud. Thèse dirigée par Bernard Rousseau.

3) Membre du jury de la thèse de Moussa Dicko, soutenue le 26 Novembre 2010. Laboratoire de Thermodynamique des Equilibres de Phases (TEP), Mines-ParisTech. Thèse dirigée par Christophe Coquelet.

2) Membre du jury de la thèse de Romain Privat, soutenue le 27 Novembre 2008. LTMP - Laboratoire de Thermodynamique des Milieux Polyphasés, ENSIC, Nancy. Thèse dirigée par Jean-Noël Jaubert.

1) Membre du jury de la thèse de Xavier Courtial, soutenue le 12 Décembre 2008. Laboratoire de Thermodynamique des Equilibres de Phases (TEP), Mines-ParisTech. Thèse dirigée par Christophe Coquelet.

2009 Membre de la commission pour le recrutement d'un maître de conférence à l'ENSIC de Nancy.

Révisions régulières d'articles scientifiques: Fluid Phase Equilib., J. Phys. Chem. B and C, Mol. Phys., J. Chem. Phys., Ind. Eng. Chem. Res., J. Chem. Eng. Data, ... Evaluations de projets de recherche.

# 7. Parcours et Synthèse des Travaux de Recherche

Mes activités de recherche sont consacrées à la modélisation des équilibres de phases de mélanges rencontrés dans les procédés industriels chimiques et énergétiques. Les applications sont la pétrochimie et les biocarburants, la production d'hydrogène, les polymères, les hydrates, les fluides frigorigènes, et les solutions d'électrolytes.

Mes recherches en modélisation thermodynamique ont commencé lors de mon stage de DEA à l'ENSIC (Nancy) sous la responsabilité de J.N. Jaubert (ENSIC), de R. Solimando (ENSIC) et d'A. Barreau (IFPEN). Le sujet était la modélisation des équilibres de phases et des propriétés thermodynamiques des gaz à condensat avec l'équation d'état Peng-Robinson, en représentant la coupe lourde par des fonctions de distribution. Ce stage fut particulièrement enrichissant car il m'a permis d'approfondir mes connaissances en équilibre de phases précédemment acquises à l'ENSIC, et d'apprendre la plupart des algorithmes de flash appliqués aux fluides pétroliers. Ce travail a permis l'amélioration de la méthode développée précédemment par M. Sportisse et a conduit à un premier article.

Ma thèse au département de génie chimique de l'Imperial College de Londres a porté sur la modélisation des équilibres de phases des solutions de polymère, en utilisant l'équation d'état SAFT. Elle fut encadrée par G. Jackson, l'un des développeurs de la théorie SAFT, et financée par BP Chemicals. Nous avons modélisé les équilibres liquide-vapeur et liquide-liquide de mélanges d'hydrocarbures + polyéthylène, et étudié l'effet de polydispersité des chaînes de polymère sur les équilibres de phase. Ces travaux ont été présentés dans plusieurs congrès internationaux, ont conduit à la publication de huit articles et ont été récompensés par le prix annuel de la meilleure thèse en calculs du département. Cette thèse m'a permis d'acquérir de solides connaissances sur les théories des fluides et la thermodynamique statistique ainsi de solides compétences dans l'équation d'état SAFT.

Au cours d'un postdoc de 19 mois au département de génie chimique de l'université de Vanderbilt (USA), j'ai développé, sous la responsabilité de P. T. Cummings, un modèle polarisable de l'eau capable de décrire la plupart des propriétés physiques de l'eau. Nous avons développé des codes de simulation moléculaire dynamique et de simulation de Monte-Carlo dans divers espaces (NVT, NPT, Gibbs) pour déterminer les propriétés de l'eau. Ce travail a conduit à deux articles dont un a été récompensé par un prix.
Depuis mon arrivé à l'unité de chimie et procédés (UCP) de l'ENSTA-ParisTech en tant qu'enseignant-chercheur, je poursuis mes activités de recherche en modélisation et en simulation moléculaire des équilibres de phases de mélanges, pour des applications au génie chimique. L'objectif principal de mes activités est le développement de modèles thermodynamiques, principalement d'équations d'état, basés sur une description à l'échelle microscopique des constituants et sur les théories issues de la thermodynamique statistique. Ces développements s'accompagnent de simulations moléculaires afin de valider les modèles thermodynamiques. Parallèlement, nous développons des modèles thermodynamiques pour des systèmes spécifiques comme les solutions d'électrolytes, en coopération avec W. Fürst, directeur de l'UCP. Nous travaillons sur la modélisation des semi-clathrates de gaz, et nous avons été les premiers à proposer une approche thermodynamique pour prédire leurs conditions de formation. Ces hydrates de gaz sont actuellement très étudiés pour le stockage et la séparation de gaz. Nous avons étudié une grande variété de systèmes, allant des systèmes aqueux en présence d'électrolytes aux mélanges multiconstituants : hydrocarbures, biocarburants, fluides frigorigènes ... Un autre objectif est la valorisation des modèles que nous développons, par l'utilisation de l'interface cape-open ou par une implémentation directe de nos codes dans des logiciels commerciaux de simulation de procédés. D'autre part, mes activités de recherche sont en cohérence avec mes enseignements (ENSTA, Master), lesquels sont orientés vers la thermodynamique, la simulation moléculaire, la programmation, et la simulation de procédés.

Une grande partie de mes activités de recherche a été financée par des financements publics (projets ANR) et des contrats. J'ai participé au cours de ces projets à l'encadrement de trois thèses et de plusieurs stagiaires masters, et encadré trois postdocs. Ces travaux ont conduit à des publications dans des journaux internationaux et ont été présentés dans différentes conférences internationales. Nous travaillons étroitement avec le laboratoire CTP de Mines-ParisTech dans différents projets ainsi qu'avec d'autres laboratoires (Institut du pétrole de l'université Heriot-Watt (Ecosse), Paris 13). Je compte dans les prochaines années étendre mes activités dans le domaine de la simulation moléculaire. Nous avons d'ailleurs obtenu récemment un projet ANR de 4 ans sur cette thématique avec une application aux fluides frigorigènes de nouvelle génération. En parallèle, nous comptons maintenir une activité forte dans l'étude des solutions d'électrolytes qui est l'une des principales thématiques du laboratoire. Voici ci-dessous le détail des contrats de recherche et des coopérations.

## Contrats de Recherche de 2007 à 2013.

Nov. 2013-Nov. 2017 Project ANR PrediRef. Le projet PrediRef porte sur la modélisation thermodynamique et la simulation moléculaire des équilibres de phases de mélanges des réfrigérants (fluides frigorigènes) de nouvelle génération. L'aide financière par ANR est de 208 k€ pour l'ENSTA-ParisTech. Ce projet est dirigé par l'entreprise Arkema, et implique les partenaires suivants: Arkema, le laboratoire UCP de l'ENSTA-ParisTech, le laboratoire CTP de Mines-ParisTech, et l'entreprise PROSIM.

Oct. 2013 – Sept. 2016 Contrat de 15 k€. Coopération avec l'institut du pétrole de l'université Heriot-Watt d'Edimbourg. Co-encadrement d'une thèse sur la modélisation des mélanges hydrocarbure-eaux avec SAFT-γ-Mie.

Sep. 2013- Août. 2014 Project chaire Areva. Contrat de 43k€. Co-encadrement d'un postdoc travaillant sur l'étude des conditions de formation des hydrates et semi-clathrates d'hydrogène.

Mars 2013 – Mars 2014 Contrat de 15 k€. Coopération le CTP de Mines-ParisTech et l'université de Paris 13. Co-encadrement d'un postdoc travaillant sur le développement d'une équation d'état crossover vérifiant les lois d'échelle autour des points critiques.

Oct. 2009- Mars 2013. Project ANR MeMoBiol.

Le projet memobiol porte sur la modélisation thermodynamique et la simulation moléculaire des équilibres de phases de mélanges de molécules valorisables, issues d'huiles de pyrolyse. Le projet a débuté le 1er Octobre 2009 et dure trois ans. L'aide financière a été de 113.5k€ pour l'ENSTA-ParisTech. Le projet a impliqué les partenaires suivants: IFP Energies Nouvelles, le laboratoire CTP
de Mines-ParisTech, le laboratoire UCP de l'ENSTA-ParisTech, le laboratoire LSPM (LIMHP) de l'université Paris 13, et les entreprises Materials design et PROSIM.

2008. Contrat de 10k€ avec le CEA de Saclay. Sujet: étude des équilibres de phases du mélange H1x (H₂O+HI+I₂+H₂) rencontré dans le procédé de production de l'hydrogène par cycle thermochimique.

## Coopérations avec d'autres laboratoires

### Coopérations nationales

Christophe Coquelet, directeur du laboratoire CTP (Centre Thermodynamique des Procédés) de Mines-ParisTech, Fontainebleau.

Moussa Dicko, Fabien Volle, Jean-Philippe Passarello, Laboratoire LSPM de l'université de Paris 13.

Rafael Lugo, Jean-Charles de Hemptinne, IFP Energies Nouvelles, Rueil-Malmaison.

Violaine Athès, unité de génie et microbiologie des procédés alimentaires, AgroParisTech, Grignon.

### Coopérations internationales

Antonin Chapoy, Chemical Engineering Department, Heriot-Watt University, Edinburgh, Ecosse.

Jinlong Li, Honglai Liu, department of chemical engineering, East China University of Science and Technology, Shanghai, China.

George Jackson, Amparo Galindo, Department of Chemical Engineering, Imperial College de Londres.

J. L. Rivera, Chemical Engineering Dept., Université de de San Nicolas de Hidalgo, Mexique.

Peter T. Cummings, Département de génie chimique, Vanderbilt University, Nashville, USA.

Szabolcs Varga, Département de Physique, Université de Veszprém, Hongrie.

Deresh Ramjugernath, School of Chemical Engineering, University KwaZulu Natal, South Africa.
# Mémoire de Travaux de Recherche

## 1. Introduction

La connaissance des propriétés thermodynamiques et des équilibres de phases des mélanges est essentielle pour le développement de procédés industriels chimiques. Les mesures expérimentales de ces propriétés et leur prédiction par des modèles mathématiques sont nécessaires pour le design et l'optimisation d'opérations unitaires telles que les colonnes de distillation, les séparateurs par extraction liquide-liquide et les réacteurs polyphasiques. Par exemple, le dimensionnement précis d'une colonne de distillation ne peut se faire que par une connaissance très précise des équilibres liquide-vapeur et de l'enthalpie du mélange à distiller, et par l'utilisation d'un simulateur de procédé dans lequel est incorporé un modèle thermodynamique.

Les propriétés des phases liquides, gazeuses ou solides dépendent de la nature des espèces chimiques impliquées, en particulier de la nature des interactions entre les molécules. Par exemple, le fait que la température critique de l'eau soit beaucoup plus élevée que celle du méthane s'explique par l'existence de fortes interactions dipôle - dipôle et de liaisons hydrogène dans l'eau. Afin de pouvoir prédire les propriétés de mélange, on peut représenter les forces d'interactions entre les molécules par un modèle mathématique que l'on appelle champ de force. Un champ de force est une somme d'interactions de différentes natures (répulsion, dispersion, dipôle - dipôle, charge - charge,...), chacune étant représentée par un potentiel intermoléculaire. Le potentiel intermoléculaire le plus connu est probablement le potentiel de Lennard-Jones $\phi^{LJ}$ (Figure 1-1) dont l'expression est donnée par

$$ \phi^{LJ}(r) = 4\epsilon \left( \left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^{6} \right). \qquad (1.1) $$

Ce potentiel permet de représenter les forces répulsives à courte distance dues au recouvrement des nuages électroniques des molécules, et les forces attractives de dispersion à moyenne portée. Les forces attractives de dispersion sont dues aux fluctuations des nuages électroniques autour des noyaux, créant des dipôles instantanés qui interagissent. Le potentiel de dispersion peut être utilisé pour les gaz rares ou les molécules quasi-sphériques (méthane). Pour les molécules plus compliquées, on utilise des interactions à sites multiples, des potentiels électrostatiques (charge-charge, dipôle-dipôle, effet de polarisation, ...), des potentiels intramoléculaires (vibration, torsion, ...).

La thermodynamique statistique permet de faire le lien direct entre les interactions modélisées par un champ de force et les propriétés thermodynamiques macroscopiques à l'équilibre d'un fluide ou d'un solide. Ce lien est fait via la relation fondamentale¹

$$ A = -kT \ln Q, \qquad (1.2) $$

où A est l'énergie libre de Helmholtz du système étudié, Q la fonction de partition dans l'ensemble canonique, et T la température. Si l'expression analytique de l'énergie libre A du système étudié est connue sur tous les domaines de température, volume et composition, toutes les autres propriétés thermodynamiques (pression, enthalpie, entropie, potentiels chimiques, capacité calorifique, ...) peuvent être prédites en utilisant les relations de la thermodynamique classique. Cependant, il est très difficile voire impossible d'établir une expression mathématique exacte pour Q pour la plupart des modèles de champs de force. L'expression analytique exacte de Q en fonction de la température et du volume est impossible à déterminer même pour les fluides simples tel que le fluide de Lennard-Jones, mis à part quelques cas idéaux, comme par exemple le gaz parfait.

La thermodynamique moléculaire tend à rendre la thermodynamique statistique applicable au génie des procédés, et consiste au développement de modèles thermodynamiques pour l'ingénieur, basés sur une description microscopique des mélanges. En général, plus la représentation des molécules est détaillée, plus l'expression mathématique du modèle thermodynamique est compliquée et donc coûteuse en temps de calcul. Lors du développement d'un modèle thermodynamique, un compromis doit donc être trouvé entre la complexité du modèle thermodynamique, et ses capacités de représentation et de prédiction. Des théories issues de la thermodynamique statistique permettent d'élaborer des expressions approximées de l'énergie libre A. Les équations d'état théoriques se
réfèrent donc à un modèle moléculaire simple et bien déterminé, et peuvent être développées à partir de deux approches: les théories de perturbation et les théories des équations intégrales. Dans ce mémoire, nous parlerons essentiellement des théories de perturbation.

Aujourd'hui, la vitesse des ordinateurs et l'utilisation de clusters sous linux permettent la simulation de systèmes larges, sur des temps plus longs, avec une description détaillée des interactions intermoléculaires. Les simulations moléculaires dynamiques et de Monte-Carlo sont de plus en plus utilisées pour prédire les propriétés et les équilibres de phases de mélanges, notamment grâce au développement de logiciels. Ces techniques sont utiles pour comprendre des phénomènes à l'échelle des molécules et pour prédire les propriétés de mélanges dans des conditions extrêmes, pour lesquelles les mesures expérimentales seraient dangereuses et coûteuses (haute température, haute pression, composés toxiques ou explosifs).

<figure><img src="image_2.png" /><figcaption>Figure 1-1. Potentiel intermoléculaire de Lennard-Jones. $\sigma$ est tel que $\phi^{LJ}(\sigma) = 0$ et $\epsilon$ est la profondeur du potentiel.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 1-2. Liens entre les mesures expérimentales, la simulation moléculaire et les modèles thermodynamiques. Application à la simulation de procédé.</figcaption></figure>
Cependant, ces techniques ont toujours besoin de données expérimentales pour la validation des résultats ou le calage de paramètres de champ de force, et la détermination des propriétés est coûteuse en temps de calculs. L'utilisation de modèles thermodynamiques nécessitant peu de temps de calcul est toujours nécessaire pour simuler un procédé industriel.

D'un point de vue industriel, la qualité d'un modèle thermodynamique est souvent jugée selon ses capacités à décrire les données expérimentales et à prédire en dehors des conditions expérimentales, avec un nombre réduit de paramètres. La clef pour le développement un tel modèle est que celui-ci repose sur une description microscopique des interactions et sur des fondements solides, afin que les paramètres du modèle aient une signification physique précise. La simulation moléculaire peut alors être très utile pour valider, tester un tel modèle.

La Figure 1-2 montre les liens étroits entre les mesures expérimentales, la simulation moléculaire, et les modèles thermodynamiques, et rappelle le schéma de Gray et Gubbins ². Les mesures expérimentales seront toujours indispensables pour le développement et l'optimisation de nouveaux procédés. Cependant, la tendance actuelle dans l'industrie est de minimiser le nombre de mesures expérimentales, afin de réaliser un premier dimensionnement du procédé, puis de réaliser d'autres mesures accompagnées systématiquement d'une analyse détaillée des erreurs expérimentales, pour valider et affiner les modèles et le dimensionnement ³. La simulation moléculaire peut alors jouer un rôle clef, en validant et en complétant les mesures.

Ce mémoire est une synthèse de mes travaux de recherche depuis mon stage de recherche en DEA en 1999, jusqu'à 2013 en tant qu'enseignant-chercheur à l'ENSTA-ParisTech. Nos recherches ont été consacrées à la modélisation thermodynamique, au développement d'équations d'état, ainsi qu'à la simulation moléculaire, avec comme objectif l'amélioration des modèles existants, la compréhension des équilibres de phases, et l'optimisation et le design de procédés industriels chimiques. Ce mémoire est organisé par thème et ne suit pas l'ordre chronologique de nos recherches. Un rappel sur les théories de perturbation et une description des principaux types de comportement de phase sont données afin d'introduire nos recherches. Les différents thèmes sont décrits, puis nos principaux résultats dans chacun de ces thèmes sont exposés. Dans un dernier chapitre, les principaux projets de recherche pour les cinq prochaines années sont présentés.

## 2. Equations d'Etat et Théories de Perturbation

Les équations d'état dérivées de l'équation de van der Waals⁴ ont suscité un très fort intérêt, à la fois scientifique et industriel. Le développement de nouvelles équations d'état a été motivé par la forte demande industrielle notamment par celle des industries pétrolière et pétrochimique. Les équations d'états cubiques dérivées de l'équation de van der Waals ont été très utilisées pour traiter des mélanges d'hydrocarbures multiconstituants (fluides pétroliers) et les mélanges à hautes pressions (fluides supercritiques). Les équations d'états cubiques utilisant des règles de mélanges combinant l'équation d'état avec un modèle à coefficient d'activité (Wong-Sandler ⁵, MHV2 ⁶, ...) ont permis de représenter les équilibres de phases de mélanges fortement non idéaux. Les équations d'état plus sophistiquées issues de la thermodynamique statistique telles que les équations de type SAFT ⁷ (Statistical Associating Fluid Theory) sont de plus en plus utilisées dans l'industrie car elles peuvent être facilement combinées à des méthodes de contribution de groupe, et permettent de représenter une grande variété de mélanges avec un nombre limité de paramètres binaires: mélanges d'hydrocarbures, de molécules polaires et apolaires, solutions de polymères ou d'électrolytes, etc. D'excellentes revues ont été faites sur les modèles thermodynamiques utilisés dans l'industrie ⁸,⁹. Afin de présenter nos travaux sur le développement de nouvelles équations d'état et de montrer quelles ont été nos contributions et améliorations, nous ferons quelques rappels en thermodynamique statistique appliquée aux théories des liquides, en particulier sur les principales avancées dans le développement des équations d'état, et montrer l'évolution de l'équation de van de Waals à SAFT.
## 2.1. L'équation de van der Waals

L'équation de van der Waals⁴ de 1873 repose sur une description des fluides à l'échelle moléculaire. En 1873, l'approche moléculaire était acceptée par la plupart des scientifiques et servait d'hypothèse pour le développement de modèles. Van der Waals fut le premier à développer une équation d'état permettant de décrire à la fois les phases gazeuses et liquides, et les équilibres liquide-vapeur se terminant au point critique du fluide. Il a ainsi montré que la phase liquide d'un corps pur n'était pas fondamentalement différente de la phase vapeur et que le passage d'une phase à l'autre était continu au point critique. Van der Waals fut également le premier à écrire une équation d'état comme étant la somme d'un terme répulsif et d'un terme attractif. Le terme répulsif de van der Waals est dérivé de la théorie cinétique des gaz et tient compte du fait que les molécules ont une certaine taille, tandis que le terme attractif tient compte des interactions attractives et a été obtenu à partir de la théorie de capillarité de Laplace. Nous allons ici déterminer les termes attractifs et répulsifs de l'équation en utilisant une approche différente de celle utilisée par van der Waals, et basée sur la thermodynamique statistique, afin de montrer les améliorations théoriques qui ont été faites par la suite sur chaque terme de l'équation d'état.

Considérons un système de N particules sphériques identiques contenues dans un volume V. On suppose que les particules interagissent via un potentiel d'interaction $\phi^S$ de type Sutherland (Figure 2-1). Ce potentiel comporte une partie répulsive dure et une partie attractive évoluant en $1/r^6$. Il est défini par

$$ \phi^S(r) = \begin{cases} +\infty & \text{if } r < \sigma \\ -\epsilon(\sigma/r)^6 & \text{if } r \ge \sigma, \end{cases} \qquad (2.1) $$

où $\sigma$ est le diamètre des particules et $\epsilon > 0$ la profondeur du potentiel.

<figure><img src="image_4.png" /><figcaption>Figure 2-1. Potentiel intermoléculaire de Sutherland en $1/r^6$.</figcaption></figure>

On peut écrire $\phi^S$ selon la somme de deux potentiels: $\phi^S(r) = \phi_\sigma^{hs}(r) + \phi^{att}(r)$, où le potentiel d'interaction $\phi_\sigma^{hs}$ des sphères dures de diamètre $\sigma$ est donné par

$$ \phi_{\sigma}^{hs}(r) = \begin{cases} +\infty & \text{if } r < \sigma \\ 0 & \text{if } r \ge \sigma, \end{cases} \qquad (2.2) $$

Le symbole "hs" signifie "hard spheres". La fonction de partition $Q$ dans l'ensemble canonique (NVT) du système de N particules sphériques (molécules atomiques) est donnée par

$$ Q = \frac{Z_N}{N!\Lambda^{3N}}, \qquad (2.3) $$

où $\Lambda = h\sqrt{2\pi mkT}$ est la longueur de de Broglie; $Z_N$ est l'intégrale de configuration et est donnée par
$$
Z_N = \int \dots \int \exp \left( - \frac{U_p(\mathbf{r}_1, \dots, \mathbf{r}_N)}{kT} \right) d\mathbf{r}_1 \dots d\mathbf{r}_N. \qquad (2.4)
$$

On détermine la moyenne d'une propriété X dans l'ensemble canonique selon,

$$
\langle X \rangle = \frac{\int \dots \int X \exp \left( -\frac{U_p}{kT} \right) d\mathbf{r}_1 \dots d\mathbf{r}_N}{Z_N}. \qquad (2.5)
$$

L'énergie potentielle d'interaction $U_p = \sum_{i,j>i} \phi(r_{ij})$ est la somme d'une partie répulsive et
attractive :

$$
U_p = U^{hs} + U^{att} = \sum_{i,j>i} \phi_{\sigma}^{hs}(r_{ij}) + \sum_{i,j>i} \phi^{att}(r_{ij}). \quad (2.6)
$$

En combinant les Eqs. (2.4) et (2.6), on obtient

$$
Z_N = \int \dots \int \exp \left( -\frac{U^{att}}{kT} \right) \exp \left( -\frac{U^{hs}}{kT} \right) d\mathbf{r}_1 \dots d\mathbf{r}_N. \quad (2.7)
$$

En remplaçant $X$ par le terme $\exp(-\frac{U^{att}}{kT})$ dans l'Eq. (2.5), on montre que

$$
Z_N = Z_N^{hs} \left\langle \exp \left( -\frac{U^{att}}{kT} \right) \right\rangle_{hs}, \qquad (2.8)
$$

où $Z_N^{hs}$ est l'intégrale de configuration du système de sphères dures. Le symbole $\langle \dots \rangle$ dans l'Eq. (2.8) signifie une moyenne statistique du terme $\exp(-\frac{U^{att}}{kT})$ sur l'ensemble $NVT$ du système de sphères dures (moyenne sur l'ensemble des positions possibles des sphères dures). L'équation d'état du système est obtenue en fonction de l'énergie libre de Helmholtz $A$, en utilisant la relation

$$
A = -kT \ln Q = A^{hs} - kT \ln \left\langle \exp \left( -\frac{U^{att}}{kT} \right) \right\rangle_{hs}, \quad (2.9)
$$

où $A^{hs}$ est l'énergie libre de Helmholtz du système de sphères dures sans attraction. L'équation (2.9) est exacte et permet de faire apparaître l'équation d'état selon la somme d'une partie répulsive et d'une partie attractive. Elle est également la base des théories de perturbation. L'équation de van der Waals correspond à certaines approximations pour les termes répulsif et attractif.

2.1.1. Terme répulsif de van der Waals

Van der Waals avait développé le terme répulsif de manière intuitive en considérant le volume
d'exclusion des molécules. Afin d'introduire les améliorations de ce terme puis nos travaux de
recherche, nous considérons ici un autre développement du terme répulsif basé sur l'équation du Viriel.
Considérons un système de N sphères dures contenues dans un volume V. Le développement du
Viriel en fonction du facteur de compressibilité des sphères dures Z^hs = PV/(NkT) est de la forme

$$
Z^{hs} = 1 + B\rho + C\rho^2 + D\rho^3 + \dots, \qquad (2.10)
$$

où $\rho = N/V$ est la densité du système, et $B$, $C$, et $D$ sont les 2ème, 3ème et 4ème coefficients du Viriel.
La terme répulsif de van der Waals consiste à faire l'hypothèse $C = B^2$, $D = B^3$, etc., c'est-à-dire que
le développement du Viriel peut s'écrire selon la suite géométrique
$Z^{hs} = 1 + B\rho + B^2\rho^2 + B^3\rho^3 + \dots$ En conséquence, pour $B\rho < 1$, $Z^{hs}$ devient

$$
Z^{hs} = \frac{1}{1 - B\rho}. \tag{2.11}
$$

Notez que tout le long de ce mémoire, nous employons le terme « densité » pour $\rho$ (en référence au
terme anglo-saxon « density ») plutôt que masse volumique, car il s'agit du nombre de molécules par
unité de volume, et non pas de la masse par unité de volume. Les coefficients du Viriel peuvent être
déterminés de manière exacte en fonction du potentiel d'interaction $\phi(r)$, d'après la thermodynamique
statistique. Pour une particule sphérique, $B$ est donné de manière exacte par

$$
B = 2\pi \int_{0}^{\infty} \left( 1 - \exp \left( -\frac{\phi(r)}{kT} \right) \right) r^2 dr. \qquad (2.12)
$$
Le potentiel d'interaction des sphères dures est donné par l'Eq. (2.2). On montre à partir des équations (2.12) et (2.2) que $B = \frac{2}{3}\pi\sigma^3$. Notez que $B$ est égal à la moitié du volume d'exclusion d'une particule par rapport à une autre.

On définit la compacité selon $\eta = \frac{\pi}{6}\rho\sigma^3$ ; $\eta$ est égale à la fraction de volume occupé par les particules. On peut exprimer $Z^{hs}$ en fonction de $\eta$ selon

$$Z^{hs} = \frac{1}{1 - 4\eta} \qquad (2.13)$$

La partie répulsive de la pression, $P^{hs} = \rho k T Z^{hs}$, est donnée par

$$P^{hs} = \frac{\rho k T}{1 - \frac{2\pi}{3}\rho\sigma^3} \qquad (2.14)$$

Si l'on définit le covolume $b$ par $b = N_{Av}B$, où $N_{Av}$ est le nombre d'Avogadro, on obtient le terme répulsif de van der Waals $P^{hs} = \frac{RT}{v-b}$ en fonction du volume molaire $v = N_{Av}/\rho$ et de la constante des gaz parfaits $R = kN_{Av}$. En intégrant par rapport à $V$ la relation $P = -\left(\frac{\partial A}{\partial V}\right)_T$, il est possible d'exprimer l'énergie libre selon $A^{hs} = A^{ideal} + A^{res,hs}$, où $A^{ideal}$ est l'énergie libre d'un gaz idéal, et $A^{res,hs}$ l'énergie libre résiduelle des sphères dures, qui est donnée par

$$\frac{A^{res,hs}}{NKT} = \int_{0}^{\rho} \frac{Z^{hs} - 1}{\rho} d\rho = -\ln\left(1 - \frac{2\pi}{3}\rho\sigma^3\right) \qquad (2.15)$$

### 2.1.2. Terme attractif de van de Waals

Considérons maintenant un système de particules interagissant via le potentiel de Sutherland (Eq. (2.1)). Déterminons l'énergie potentielle d'interaction attractive $dU_i^{att}$ entre une particule $i$, et les particules $j$ comprises dans le volume $dv = 4\pi r^2 dr$ délimité par deux sphères concentriques de rayons $r$ et $r + dr$, avec $r \ge \sigma$ (Figure 2-2).

<figure><img src="image_5.png" /><figcaption>Figure 2-2. Energie potentiel d'interaction entre une particule donnée $i$ et les particules $j$ telles que la distance soit comprise entre $r$ et $r + dr$.</figcaption></figure>

Si l'on suppose que la densité $\rho$ est parfaitement homogène, le nombre de particules $j$ comprises dans le volume $dv$ est $\rho 4\pi r^2 dr$, donc l'énergie d'interactions attractive entre la particule $i$ et les particules $j$ vaut $dU_i^{att} = \phi^S(r)\rho 4\pi r^2 dr$. Si l'on intègre sur tout le volume $V$, on obtient l'énergie potentielle d'interaction attractive $U_i^{att}$ entre la particule $i$ et toutes les autres $N-1$ particules,

$$U_i^{att} = \int_{\sigma}^{+\infty} \phi^S(r) \rho 4\pi r^2 dr \qquad (2.16)$$
Notez que l'on a pris "+\infty" pour la borne supérieure de l'intégrale, ce qui correspond à la limite thermodynamique ($N \to \infty$, $V \to \infty$). D'après l'Eq. (2.16), l'énergie potentielle d'interaction attractive $U_i^{att}$ entre la particule $i$ et toutes les autres $N-1$ particules ne dépend pas de la position des particules, et est proportionnelle à la densité $\rho$. On appelle "théorie du champs moyen" une telle approximation.

L'énergie potentiel totale d'interaction attractive ($U^{att}$) dans le système de $N$ particules est égale à $\frac{N}{2}U_i^{att}$, c'est-à-dire

$$U^{att} = 2\pi N \rho \int_{\sigma}^{+\infty} \phi^S(r) r^2 dr = -\frac{2}{3}\pi N \rho \sigma^3 \epsilon. \qquad (2.17)$$

Le facteur $\frac{1}{2}$ permet de ne pas compter deux fois les paires d'interactions. L'expression de l'énergie libre de l'équation de van der Waals est obtenue en remplaçant l'expression de $U^{att}$ (Eq. (2.17)) dans l'Eq. (2.9). Etant donné que $U^{att}$ ne dépend pas des positions des particules, le symbole $\langle...\rangle$ dans l'Eq. (2.9) peut être enlevé, et l'on obtient l'équation de van der Waals sous forme de l'énergie libre :

$$\frac{A}{NkT} = \frac{A^{ideal}}{NkT} - \ln \left( 1 - \frac{2\pi}{3}\rho\sigma^3 \right) - \frac{2}{3}\pi\rho\sigma^3 \frac{\epsilon}{kT}. \qquad (2.18)$$

On rappelle que $A^{ideal}$ pour les molécules atomiques est donnée par

$$\frac{A^{ideal}}{NkT} = \ln(\rho\Lambda^3) - 1. \qquad (2.19)$$

Notez que $U^{att}$ donnée par l'Eq. (2.17) est égale à l'énergie interne résiduelle de l'équation de van der Waals car la partie répulsive de l'énergie interne résiduelle est nulle. On peut comparer l'Eq. (2.17) avec l'expression exacte de l'énergie interne résiduelle $U^{res}$ d'un système de particules sphériques¹,

$$U^{res} = 2\pi N \rho \int_{0}^{+\infty} \phi(r) g(r) r^2 dr, \qquad (2.20)$$

où $g(r)$ est la fonction de distribution radiale du système; $g(r)$ caractérise la structure d'un fluide et les variations locales de densité. On voit $U^{att}$ peut donc être obtenue en faisant l'hypothèse que $g(r) = 0$ pour $r < \sigma$ et $g(r) = 1$, pour $r \ge \sigma$. L'équation de van der Waals est une équation d'état dite cubique car la résolution de la densité à température et pression fixées consiste à résoudre analytiquement un polynôme de degré 3. D'autre part, il est possible de relier les paramètres de l'équation d'état $(\epsilon, \sigma)$ aux propriétés critiques $(T_c, P_c, V_c)$ calculées: en appliquant les conditions aux points critiques $(\partial P/\partial V = 0$ et $\partial^2 P/\partial V^2 = 0)$, on montre que

$$T_c = \frac{8}{27} \frac{\epsilon}{k}, P_c = \frac{\epsilon}{18\pi\sigma^3}, Z_c = \frac{3}{8}. \qquad (2.21)$$

Il existe alors deux méthodes de détermination des paramètres $\epsilon$ et $\sigma$ de l'équation d'état pour un corps pur donné. L'une d'elles consiste à optimiser les paramètres $\epsilon$ et $\sigma$ sur les données expérimentales d'équilibre L-V ou de type PVT. L'autre méthode consiste à faire en sorte que les données critiques calculées correspondent aux données critiques expérimentales. Les paramètres $\epsilon$ et $\sigma$ sont alors déterminés à partir de l'Eq. (2.21). Par exemple pour le méthane ($T_c = 190.6$ K, $P_c = 4.61$MPa), on trouve que $\epsilon/k = 643.275$K, $\sigma = 3.242$Å via l'Eq. (2.21). L'équation de vdw permet de décrire qualitativement un équilibre liquide-vapeur et l'allure des isothermes d'un corps pur. Cependant, quelque soit la méthode de détermination de paramètre utilisée, la représentation des données expérimentales L-V est plutôt mauvaise. Le terme répulsif de vdw suretime très largement la pression répulsive aux densités élevées. En pratique, la compacité $\eta$ d'un liquide se situe entre 0.3 et 0.5. Cependant, le terme répulsif de vdw n'est valable que pour des compacités inférieures à 0.25 (voir l'Eq. (2.13)).

Afin d'améliorer la représentation des données, de nombreux auteurs ont essayé de modifier de manière semi-empirique le terme attractif de vdw tout en s'imposant les contraintes de pouvoir résoudre analytiquement le volume molaire à T et P fixées, et de pouvoir relier analytiquement les paramètres de l'équation d'état aux propriétés critiques. Ils ont gardé le terme répulsif de vdw bien que celui-ci est mauvais d'un point de vue quantitatif, car ce terme permet la résolution analytique des racines de l'équation d'état. Ainsi, Redlich et Kwong¹¹ (RK) ont proposé une équation d'état cubique qui a été très utilisée dans l'industrie, notamment l'industrie pétrolière. L'équation RK donne de bonnes prédictions pour les seconds coefficients du Viriel. De nombreux auteurs ont ensuite essayé
d'améliorer l'équation RK. On peut citer les équations cubiques SRK<sup>12</sup>, Peng-Robinson<sup>13</sup>, et Patel-Teja<sup>14</sup> parmi celles les plus utilisées dans l'industrie. Les équations cubiques sont particulièrement précises pour le calcul des compositions des phases pour les mélanges d'hydrocarbures tels que les fluides pétroliers. Les équations cubiques peuvent être utilisées avec des règles de mélanges spécifiques faisant appel à un modèle à énergie de Gibbs d'excès, pour représenter les mélanges fortement non idéaux (règles de mélanges d'Huron Vidal<sup>15</sup> et de Wong-Sandler<sup>5</sup>). Des méthodes de contributions de groupe telles que PPR78<sup>16-18</sup> ont été développées pour les équations cubiques et donnent d'excellentes prédictions des diagrammes de phases et des lignes de points critiques.

## 2.2. L'équation d'état des sphères dures (Carnahan-Starling et BMCSL)

Van der Waals savait que son terme répulsif (Eq. (2.11)) n'était pas bon quantitativement et qu'il était uniquement valable aux faibles densités, mais il ne disposait pas à l'époque de toutes les lois de la thermodynamique statistique pour élaborer une expression très précise. Le volume exclu effectif pour chaque particule a tendance à diminuer lorsque la densité augmente. A des densités plus élevées, les molécules sont en moyenne plus proches les unes des autres. Les volumes d'exclusion ne sont pas additifs car ils se recouvrent. En conséquence, le terme de van der Waals (Eq. (2.13)) est trop répulsif car il surestime le volume d'exclusion effectif par molécule aux densités élevées. L'effet de recouvrement des volumes d'exclusion est illustré dans la Figure 2-3. En 1896, Van der Waals proposa une expression tenant compte de la diminution du volume d'exclusion effectif lorsque la densité augmente. Cette expression permet d'estimer le 3ème coefficient du Viriel C, mais celle-ci s'est avérée fausse ($C = \frac{15}{32}B^2$), bien qu'assez proche de l'expression exacte $C = \frac{5}{8}B^2$ déterminée plus tard par Jäger et Boltzmann<sup>4</sup>.

<figure><img src="image_6.png" /><figcaption>Figure 2-3 Sphères dures de diamètres $\sigma$ et volumes d'exclusion. Le volume d'exclusion de la particule $i$ par rapport aux deux particules $j$ et $k$ est inférieur à la somme des deux volumes d'exclusion (sphères de diamètres $2\sigma$ délimitées par les pointillés) car ceux-ci se recouvrent en partie si les particules $j$ et $k$ sont proches. Le volume de recouvrement est illustré par la zone hachurée et s'annule si la distance entre les particules $j$ et $k$ est supérieure à $2\sigma$.</figcaption></figure>

Suite aux progrès en thermodynamique statistique et au développement de la simulation moléculaire, plusieurs expressions du facteur de compressibilité des sphères dures en fonction de la compacité $\eta$ ont été proposées. L'expression qui a connu le plus de succès à cause de sa grande précision et de sa simplicité est celle de Carnahan et Starling (CS)<sup>19</sup>. Elle est donnée par

$$Z^{hs} = \frac{1 + \eta + \eta^2 - \eta^3}{(1 - \eta)^3}. \qquad (2.22)$$

L'équation CS a été développée à partir de l'équation du Virial exacte des sphères dures, qui est donnée par<sup>20</sup>

$$Z^{hs} = 1 + 4\eta + 10\eta^2 + 18.365\eta^3 + 28.225\eta^4 + 39.74\eta^5 + 53.5\eta^6 + 70.8\eta^7 + \dots \quad (2.23)$$

L'équation CS donne un développement très proche,

$$Z^{hs} = 1 + 4\eta + 10\eta^2 + 18\eta^3 + 28\eta^4 + 40\eta^5 + 54\eta^6 + 70\eta^7 + \dots \qquad (2.24)$$
Comme le montre la Figure 2-4, le modèle CS permet de prédire avec une excellente précision le facteur de compressibilité des sphères dures déterminé par simulation moléculaire. L'énergie libre résiduelle correspondante au terme répulsif CS est donnée par

$$ \frac{A^{res,hs}}{NkT} = \frac{4\eta - 3\eta^2}{(1 - \eta)^2}. \qquad (2.25) $$

On peut combiner le terme répulsif CS avec le terme attractif de van der Waals pour donner une nouvelle équation d'état que l'on nommera l'équation d'état de "van der Waals améliorée".

<figure><img src="image_7.png" /><figcaption>Figure 2-4 Facteur de compressibilité des sphères dures en fonction de la compacité. Les cercles représentent les données de simulation moléculaire d'Alder et Wainwright 21. La ligne continue correspond au modèle CS, tandis que la courbe en pointillé correspond au terme vdW.</figcaption></figure>

Pour cette équation, l'énergie libre de Helmholtz est donnée par

$$ \frac{A}{NkT} = \frac{A^{ideal}}{NkT} + \frac{4\eta - 3\eta^2}{(1 - \eta)^2} - 4\frac{\epsilon}{kT}\eta. \qquad (2.26) $$

Une telle équation d'état est applicable aux molécules quasi-sphériques interagissant uniquement via des forces de dispersion (méthane, argon, ...). L'optimisation des paramètres $\epsilon$ et $\sigma$ sur les données liquide-vapeur (densités liquides saturantes et pressions de vapeur) permet d'obtenir une équation d'état plus précise (voir le cas du méthane, thèse de Paricaud 22, page 213).

Plus tard, Boublík 23 et indépendamment Mansoori et al. 24 ont développé une équation d'état (modèle BMCSL) pour les mélanges de sphères dures dites additives, correspondant à des diamètres de collision croisés additifs ($\sigma_{ij} = \frac{\sigma_i+\sigma_j}{2}$). Le modèle BMCSL est très précis pour les sphères de tailles assez proches, et a été très utilisé car son expression est simple et compatible avec l'équation d'état de Carnahan et Starling (Eq. (2.25)). L'énergie libre résiduelle du modèle BMCSL est donnée par

$$ \frac{A^{res,hs}}{NkT} = \frac{6}{\pi\rho} \left[ \left( \frac{\zeta_2^3}{\zeta_3^2} - \zeta_0 \right) \ln(1 - \zeta_3) + \frac{3\zeta_1\zeta_2}{1 - \zeta_3} + \frac{\zeta_2^3}{\zeta_2^3(1 - \zeta_3)^2} \right], \qquad (2.27) $$

où $\zeta_l = \frac{\pi}{6}\rho \sum_{i}^{n} x_i \sigma_i^l$. L'extension de l'Eq. (2.27) que nous avons proposée pour les mélanges de sphères non additives 25 (c'est-à-dire pour les mélanges tels que $\sigma_{ij} \neq \frac{\sigma_i+\sigma_j}{2}$) est présentée dans la Section 2.4.2.

## 2.3. Théories de Perturbation
### 2.3.1. Théorie de Zwanzig

Une théorie de perturbation permet de développer une équation d'état correspondant à un potentiel d'interaction donné $\phi(r)$ en utilisant un potentiel d'interaction de référence $\phi^{ref}(r)$ dont l'équation d'état est connue. Le lecteur peut trouver plusieurs revues sur ce sujet dans la littérature¹,²,²⁶⁻²⁸. L'expression de l'équation d'état est obtenue à partir d'un développement de Taylor de la fonction de partition par rapport à un ou plusieurs paramètres. Ces paramètres permettent de passer du potentiel de référence $\phi^{ref}(r)$ au potentiel $\phi(r)$. Par exemple, il est commun d'utiliser un paramètre $\lambda$ tel que $\phi(r, \lambda) = \phi^{ref}(r) + \lambda\phi^{att}(r)$ et $\phi(r, \lambda = 1) = \phi(r)$. On peut alors développer l'équation d'état correspondant à $\phi(r)$ à partir d'un développement de Taylor de la fonction de partition autour de $\lambda = 0$ et calculée en $\lambda = 1$²⁹. Il existe de nombreuses théories de perturbation, correspondant à des choix différents de potentiels de référence ou à des techniques de perturbation spécifique. Les premiers ayant proposé de telles théories sont Pople³⁰ et Zwanzig²⁹. Voici un rappel rapide de la théorie de perturbation de Zwanzig²⁹ qui permet de déterminer de manière exacte le terme d'ordre 1 de la plupart des théories de perturbation. On utilise généralement comme système de référence un système de sphères dures ayant un certain diamètre $d$ (c'est-à-dire, on prend $\phi^{ref}(r) = \phi_d^{hs}(r)$), car le facteur de structure du fluide de sphères dures est très similaire à celui des fluides simples (potentiel d'interaction de type Lennard-Jones), ce qui permet de faire converger rapidement la série de Taylor. Le diamètre $d$ est choisi de manière à faire converger au plus vite la théorie de perturbation ou à annuler certains termes. En appliquant un développement de Taylor par rapport à $U^{att}$ au terme attractif de l'Eq. (2.9), on montre que²⁹

$$ 
\begin{aligned}
\frac{A}{NkT} &= \frac{A^{hs}}{NkT} + \frac{A_1}{NkT} + \frac{A_2}{NkT} + \dots \\
&= \frac{A^{hs}}{NkT} + \left\langle \frac{U^{att}}{NkT} \right\rangle_{hs} + \frac{1}{2N} \left( \left\langle \left( \frac{U^{att}}{kT} \right) \right\rangle_{hs}^2 - \left\langle \left( \left( \frac{U^{att}}{kT} \right) \right) \right\rangle_{hs}^2 \right) + \dots, \quad (2.28)
\end{aligned}
$$

Lorsque $\phi(r)$ est un potentiel d'interaction par paire, on montre que le premier terme de l'expansion est donné par¹

$$ \left\langle \frac{U^{att}}{NkT} \right\rangle_{hs} = \left\langle \sum_{i<j} \phi^{att}(r_{ij}) \right\rangle_{hs} = \frac{2\pi\rho}{kT} \int \phi^{att}(r) g^{hs}(r) r^2 dr, \quad (2.29) $$

où $g^{hs}(r)$ est la fonction de distribution du fluide de référence (sphères dures); $g^{hs}(r)$ caractérise la microstructure du fluide de référence, c'est-à-dire le fait que la densité de ce fluide ne soit pas homogène au niveau microscopique. Pour les sphères dures, il existe dans la littérature plusieurs expressions analytiques de $g^{hs}(r)$ et l'intégrale dans l'Eq. (2.29) se calcule analytiquement ou numériquement selon la complexité de $g^{hs}(r)$ et $\phi^{att}(r)$. On remarque une similitude entre les Eqs. et (2.29) : le terme attractif de van der Waals est basé sur l'approximation $g^{hs}(r) = 1$, c'est-à-dire que l'on suppose que la densité du fluide est homogène au niveau microscopique.

### 2.3.2. Théorie de Barker-Henderson

La théorie de Barker-Henderson (BH) est utilisée dans la plupart des versions de l'équation d'état SAFT. Il existe deux versions de la théorie BH: une version pour les potentiels d'interaction à répulsion dure³¹ (square-well, Sutherland, Yukawa, ..), et une autre pour les potentiels à répulsion souple³² (Lennard-Jones). La théorie de Barker-Henderson à l'ordre 2 est utilisée dans la plupart des versions de l'équation d'état SAFT. Le terme d'ordre 1 de la théorie BH est le même que celui de la théorie de Zwanzig et est donné par l'Eq. (2.29). Le terme exact d'ordre 2 est compliqué car il implique des intégrales multiples et les fonctions de distributions radiales à 3 et 4 particules. Barker et Henderson ont proposé plusieurs approximations de ce terme. Le lecteur peut se reporter au livre de McQuarrie¹ pour un développement complet du terme d'ordre 2. L'approximation couramment utilisée dans les versions de l'équation d'état SAFT est l'approximation dite de "compressibilité locale". Le terme d'ordre 2 ($A_2$) est donné par
$$ \frac{A_2}{NkT} = -\frac{\pi \rho}{(kT)^2} \int_{\sigma}^{\infty} [\phi^{att}(r)]^2 kT \left[ \frac{\partial}{\partial P} (\rho g^{hs}(r)) \right] r^2 dr. \quad (2.30) $$

La théorie BH pour le potentiel de Lennard-Jones ³² $\phi^{LJ}(r)$ fait appel à un potentiel d'interaction intermédiaire $\phi_{\alpha,\gamma,d}^{pert}(r)$ tel que $\phi_{\alpha=1,\gamma=1,d}^{pert}(r) = \phi_d^{LJ}(r)$ et $\phi_{\alpha=0,\gamma=0,d}^{pert}(r) = \phi_d^{hs}(r)$ ($\phi_d^{hs}$ est le potentiel pour sphères dures de diamètre $d$). L'expression de l'équation d'état est obtenue à partir d'un développement de Taylor par rapport aux deux paramètres $\alpha$ et $\gamma$. L'expression de $d$ a été choisée de manière à annuler certains termes de l'expansion ³². Elle est donnée par

$$ d = \int_{0}^{\sigma} \left(1 - \exp\left[-\frac{\phi^{LJ}(r)}{(kT)}\right]\right) dr. \quad (2.31) $$

Le diamètre effectif de sphère dure $d$ diminue lorsque $T$ augmente (Eq. (2.31)): les molécules à haute température ont une énergie cinétique plus importante et ont tendance à se rentrer d'avantage les unes dans les autres.

### 2.3.3. Théorie de perturbation pour les fluides multipolaires.

Dans les années 50 à 70, les chercheurs ont essayé de développer des équations d'état permettant de différencier les forces de dispersion des forces électrostatiques entre les molécules, l'idée étant d'améliorer les prédictions des règles de mélange et de combinaison pour les mélanges de composés polaires et apolaires, donc d'éviter l'utilisation de paramètres d'interactions binaires. L'approche consiste à modéliser les forces électrostatiques entre molécules polaires et quadripolaires par des interactions entre dipôles et quadripôles ponctuels. Cette approche est basée sur l'expansion multipolaire d'une distribution de charges dont la somme est nulle : la somme des énergies d'interactions de type charge-charge entre deux molécules électriquement neutres peut être exprimée comme la somme d'interactions de type dipôle-dipôle, quadripôle-quadripôle, dipôle-quadripôle, octopôle-octopôle, etc ... Le lecteur est dirigé vers le livre de Gray et Gubbins pour d'avantage de détails sur ce sujet ², ainsi qu'à une récente revue de Paricaud et Passarello ²⁷. Pople ³⁰ fut le premier à déterminer les termes de perturbation dans l'expression de l'énergie libre, pour tenir compte des interactions multipolaires, en séparant le potentiel intermoléculaire en une partie ne dépendant que de la distance entre deux molécules, et en une partie dépendant de la distance et de l'orientation relative entre les deux molécules. Plusieurs équipes de recherche comme les groupes de Stell, de Gubbins ou de Henderson ont par la suite travaillé sur ce sujet. Le fluide de Stockmayer est probablement le système qui a été le plus étudié, ou par simulation moléculaire, ou par les théories de perturbation. Il s'agit d'un corps pur dont les interactions sont représentées par un potentiel de Lennard-Jones, auquel on a ajouté un dipôle ponctuel centré autour du potentiel de Lennard-Jones. Un tel modèle moléculaire est très bien adapté pour des petites molécules polaires comme HCl.

Le potentiel d'interaction dipôle-dipôle entre deux molécules i et j est donné par ²:

$$ u_{\mu\mu}(\mathbf{r}_i\mathbf{r}_j) = \mu_i\mu_j \left( \frac{-2 \cos \theta_i \cos \theta_j + \sin \theta_i \sin \theta_j \cos(\phi_i - \phi_j)}{r_{ij}^3} \right) \quad (2.32) $$

où $\mu_i$ et $\mu_j$ sont les dipôles des molécules i et j, et $r_{ij}$ la distance entre les centres des deux particules. Les angles $\theta_i$, $\theta_j$, $\phi_i$ et $\phi_j$ sont définis dans la Figure 2-5.

<figure><img src="image_8.png" /><figcaption>Figure 2-5. Deux molécules sphériques LJ dipolaires (Fluide de Stockmayer).</figcaption></figure>

On peut écrire le potentiel de Stockmayer $u_{ij}$ comme étant la somme d'un terme isotrope et d'un terme anisotrope proportionnel à un paramètre $\lambda$, selon $u_{ij,\lambda} = u_{LJ}(r_{ij}) + \lambda u_{\mu\mu}(\mathbf{r}_i\mathbf{r}_j)$. Pour $\lambda = 0$, le potentiel correspond au fluide de Lennard-Jones, et pour $\lambda = 1$ on obtient le potentiel de Stockmayer.
On peut alors développer l'équation d'état du fluide de Stockmayer comme étant une perturbation par rapport au fluide Lennard-Jones (fluide de référence). L'idée consiste à faire une expansion de Taylor de l'énergie libre par rapport à λ (idée similaire aux théories de perturbation de Zwanzig et de Barker-Henderson). Sachant que

$$ \left\langle \frac{\partial u_{ij}}{\partial \lambda} \right\rangle_{\omega_i \omega_j} = 0, \qquad (2.33) $$

Le symbole $\langle \rangle$ signifie une moyenne selon l'orientation des deux molécules. Le terme d'ordre 1 de la perturbation est nul. La théorie de perturbation commence donc à partir du terme d'ordre 2, lequel est donné par

$$ \frac{A_2}{NkT} = -\frac{\rho}{4(kT)^2} \int \langle u_{\mu\mu}^2(r_{ij}, \omega_i, \omega_j) \rangle_{\omega_i \omega_j} g^{LJ}(r_{ij}) d\mathbf{r}_{ij}. \qquad (2.34) $$

$g^{LJ}$ est la fonction de distribution radiale du fluide de référence (Lennard-Jones ici). Le terme d'ordre 3 fait intervenir la fonction de distribution à trois particules $^{26}$. Contrairement au cas des fluides simples (Lennard-Jones), les termes d'ordre supérieur pour les fluides polaires ne sont pas négligeables. Stell et col. $^{33}$ ont proposé de les estimer via une approximation de Pade supposant une suite géométrique. Ainsi l'énergie libre du fluide de Stockmayer s'exprime selon

$$ \frac{A}{NkT} = \frac{A^{LJ}}{NkT} + \frac{A_2}{NkT} \left( \frac{1}{1 - A_2/A_3} \right). \qquad (2.35) $$

Des expressions analytiques de type polynomial ou approximation de Pade, et fonctions de la température et de la densité $\rho$, ont été proposées pour calculer les intégrales intervenant dans les termes de la perturbation. Stell et col. $^{33}$ et Gubbins et col. $^{34}$ ont proposé de telles expressions pour les sphères dures multipolaires et les fluides de Lennard-Jones multipolaires, respectivement. Différents termes multipolaires de l'énergie libre ont été repris par la suite pour développer des versions multipolaires de l'équation d'état SAFT $^{35-42}$. Tandis que certaines versions ne considèrent qu'un seul dipôle par molécule, Jog et Chapman $^{42,43}$ ont proposé une version multipolaire de SAFT pouvant considérer plusieurs dipôles locaux sur une même molécule (approche segment). Cette version est bien adaptée au développement de méthodes de contributions de groupes. En utilisant l'approche segment de Jog et Chapman $^{42,43}$ et les expressions des intégrales proposées par Gubbins et Twu, $^{34}$, Passarello et col $^{36,44-46}$ ont développé la version polaire PPC-SAFT de l'équation d'état PC-SAFT, ainsi qu'une méthode de contribution de groupe combinée avec cette équation d'état, appelée GC-PPC-SAFT. Nous avons utilisé l'équation d'état PPC-SAFT dans le cadre du projet ANR memobiol et dans le cadre d'une coopération avec Mines Paristech et Arkema sur des mélanges de fluides frigorigènes $^{47}$. Nous avons également utilisé une autre version polaire de PC-SAFT (terme de Jog et Chapman $^{43}$) dans le cadre d'une thèse (voir la section 5.2).

### 2.3.4. Théorie de Wertheim pour les fluides associatifs

Au début des années 80, Wertheim $^{48-50}$ a développé une théorie pour décrire les propriétés thermodynamiques des fluides associatifs (fluides comportant des liaisons hydrogène). Dans l'approche de Wertheim, le potentiel intermoléculaire par paire considéré est la somme d'un potentiel de sphère dure $\phi_R$ (potentiel de référence) et d'un certain nombre de potentiels à puits carrés $\phi_A$ (square-well) à courte portée et centrés près de la surface sphérique décrite par le potentiel $\phi_R$. Le potentiel intermoléculaire s'écrit selon

$$ \phi(\mathbf{r}_i, \mathbf{r}_j) = \phi_R(r_{ij}) + \sum_{a_i, b_j} \phi_A(\mathbf{r}_{a_i}, \mathbf{r}_{b_j}), \qquad (2.36) $$

où $\mathbf{r}_{a_i}$ et $\mathbf{r}_{b_j}$ sont les vecteurs positions des sites $a_i$ et $b_j$ situés sur les molécules $i$ et $j$, et $r_{ij} = |\mathbf{r}_i - \mathbf{r}_j|$. L'obtention de la théorie de Wertheim est compliquée car celle-ci est basée sur la thermodynamique statistique, l'expansion du Viriel et la théorie des graphes. Nous allons ici réobtenir l'expression finale de Wertheim avec un site associatif, à partir d'une approche basée sur la loi d'action de masse (d'après Paricaud $^{22}$). Considérons un corps pur $A$ associatif, comportant $N$ sphères dures, chacune comportant un site associatif (Figure 2-6).
<figure><img src="image_9.png" /><figcaption>Figure 2-6. Sphères dures associatives ayant un site associatif (modèle 1A) représenté par un potentiel à puits carré (cercles roules).</figcaption></figure>

On peut alors considérer l'équilibre chimique $A + A \rightleftharpoons AA$. Soit $\rho_0$ la densité en nombre des monomères et $\rho_D$ la densité des dimères à l'équilibre. La densité en nombre total du composé ($\rho = N/V$) est donc égale à $\rho = 2\rho_D + \rho_0$ (conservation de la masse). L'équilibre chimique induit que $\rho_D/\rho_0^2 = K(\rho, T)$, où K est la constante d'équilibre apparente définie par rapport aux densités en nombre. Dans la limite du mélange de gaz idéaux ($\rho \to 0$), la constante d'équilibre devient uniquement fonction de la température (coefficients de fugacité = 1), et peut s'exprimer de manière exacte d'après la thermodynamique statistique du gaz idéal dimérisant, selon⁵¹

$$ \lim_{\rho \to 0} K = \frac{1}{2} \int e_R(12) f(12) d(12). \qquad (2.37) $$

où la notation (12) représente la position relative de deux particules 1 et 2 ; $e_R(12) = \exp(-\beta\phi_R(12))$ avec $\beta = 1/kT$ et $f(12) = \exp\phi_A(12) - 1$ est la fonction de Mayer. On peut donc montrer pour un gaz dimérisant à faible pression que

$$ \rho = \rho_0 + \rho_0^2 \int e_R(12) f(12) d(12). \qquad (2.38) $$

Wertheim a montré que pour une phase dense, on obtient

$$ \rho = \rho_0 + \rho_0^2 \int g_R(12) f(12) d(12). \qquad (2.39) $$

où $g_R$ est la fonction de distribution radiale du système de référence non associé (sphères dures non associatives). L'Eq. (2.39) correspond à la théorie de perturbation à l'ordre 1 (théorie TPT1) de Wertheim. Déterminons maintenant la contribution de l'association à l'énergie libre résiduelle. Considérons un gaz pur idéal non associatif contenant N molécules dans un volume V, à la température T. L'énergie libre est donnée par l'Eq. (2.19) et son potentiel chimique est donné par

$$ \mu^{\text{ideal}} = kT \ln (\rho\Lambda^3). \qquad (2.40) $$

Considérons maintenant un mélange de gaz idéaux constitués de monomères et de dimères en équilibre. L'équilibre chimique donne $2\mu_0 = \mu_D$. La pression d'un tel système est $P = (\rho_0 + \rho_D)kT$. L'énergie libre de ce mélange $A_D^{\text{ideal}}$ est donnée par $A_D^{\text{ideal}} = N_D\mu_D + N_0\mu_0 - PV$. Soit N le nombre total de monomère associés ou pas. On a $N = N_0 + 2N_D$. Le terme $A_D^{\text{ideal}}$ réduit par $NkT$ s'écrit donc

$$ \frac{A_D^{\text{ideal}}}{NkT} = \frac{(2N_D + N_0)\mu_0}{2N_D + N_0} - \frac{\rho_0 + \rho_D}{\rho} = \ln(\rho\Lambda^3) - \frac{\rho_0 + \rho_D}{\rho}. \qquad (2.41) $$

La contribution nette due à l'association ($A^{\text{assoc}}$) est donc égale à

$$ \frac{A^{\text{assoc}}}{NkT} = \frac{A_D^{\text{ideal}} - A^{\text{ideal}}}{NkT} = \ln \frac{\rho_0}{\rho} - \frac{\rho_0 + \rho_D}{\rho} + 1. \qquad (2.42) $$

En posant $X_A = \rho_0/\rho$, on montre que

$$ \frac{A^{\text{assoc}}}{NkT} = \ln X_A - \frac{X_A}{2} + \frac{1}{2}, \qquad (2.43) $$

avec

$$ X_A = \frac{1}{1 + \rho X_A \int e_R(12) f(12) d(12)}. \qquad (2.44) $$

$X_A$ est la fraction de molécules non associées au site A dans le système, et est ici égale à la fraction de monomère. Chapman et col⁷,⁵²,⁵³ ont généralisé la contribution associative pour les mélanges de composés ayant plusieurs sites d'association. L'expression générale est donnée par
$$ \frac{A^{\text{assoc}}}{NkT} = \sum_{i=1}^{n} x_i \left[ \sum_{a=1}^{s_i} \left( \ln X_{a,i} - \frac{X_{a,i}}{2} \right) + \frac{s_i}{2} \right], \quad (2.45) $$

où $X_{a,i}$ est la fraction de molécules $i$ non liée au site $a$, et est telle que

$$ X_{a,i} = \frac{1}{1 + \sum_{j=1}^{n} \sum_{b=1}^{s_j} \rho_{x_j} X_{b,j} \Delta_{a,i,b,j}}. \quad (2.46) $$

### 2.3.5. Théorie de Wertheim TPT1 pour les chaînes polymériques

Il existe plusieurs manières d'obtenir la contribution de chaîne utilisée dans l'équation d'état SAFT. Nous présentons ici l'approche qui consiste à prendre comme limite une énergie d'association infinie (association complète), et à considérer des molécules comportant un site donneur et un site accepteur (modèle associatif 2B). Une équation d'état pour les chaînes polymériques est alors obtenue comme étant la limite de l'association complète de la théorie de Wertheim pour les fluides associatifs. Le système étudié est illustré par la Figure 2-7. Il contient des sphères associatives de type 2B ($S_{2B}$), mais aussi des sphères associatives de type 1A ($S_{1A}$) qui serviront de fin de chaînes.

<figure><img src="image_10.png" /><figcaption>Figure 2-7. Mélanges de sphères dures associatives avec un site ou deux sites associatifs (modèles associatifs 1A et 2B), formant une chaîne de $m$ sphères en prenant la limite de l'association complète ($X \to 0$).</figcaption></figure>

On considère le facteur de compressibilité suivant

$$ Z_s^{\text{assoc}} = \frac{PV}{N_s kT} = \sum_{i=1}^{m} \rho_s x_i \sum_{a=1}^{s_i} \left[ \frac{\partial X_{a,i}}{\partial \rho_s} \left( \frac{1}{X_{a,i}} - \frac{1}{2} \right) \right], \quad (2.47) $$

où $N_s = Nm$ est le nombre total de segments sphériques dans le mélange et $\rho_s = N_s/V$. $X_{a,i}$ est la fraction de molécules $i$ non liées au site $a$. En tenant compte de la stœchiométrie lors de la formation des chaînes $C_m$ de $m$ segments ($2S_{1A} + (m-2)S_{2B} \to C_m$), on montre que

$$ Z_s^{\text{assoc}} = \rho_s \frac{2(m-1)}{m} \left[ \frac{\partial X}{\partial \rho_s} \left( \frac{1}{X} - \frac{1}{2} \right) \right], \quad (2.48) $$

où $X$ est donné par

$$ X = \frac{1}{1 + \frac{\rho_s}{m} X \Delta_{AB}}. \quad (2.49) $$

Lorsque l'énergie d'association tend vers l'infini, la fraction de molécules non associées $X$ tend vers 0. En substituant l'Eq. (2.49) dans l'Eq. (2.48), et en prenant la limite $X \to 0$, on obtient l'expression suivante

$$ \lim_{X \to 0} Z_s^{\text{assoc}} = Z_s^{\text{chain}} = - \frac{(m-1)}{m} \left[ 1 + \frac{\rho_s}{\Delta_{AB}} \frac{\partial \Delta_{AB}}{\partial \rho_s} \right], \quad (2.50) $$

$\Delta_{AB}$ est une intégrale pouvant être simplifiée en supposant que la fonction de distribution radiale $g_R$ du système de référence est constante le long de l'intervalle d'intégration, et est égale à la valeur de contact $g_R(\sigma)$. $\Delta_{AB}$ est alors estimée par

$$ \Delta_{AB} = \kappa_{AB} \left( \exp \left( \frac{\epsilon_{AB}}{kT} \right) - 1 \right) g_R(\sigma). \quad (2.51) $$
$\kappa_{AB}$ et $\epsilon_{AB}$ sont respectivement le volume et l'énergie d'association. Comme seul $g_R$ dépend de la densité du système dans l'expression (2.51), l'Eq. (2.50) devient

$$Z_s^{\text{chain}} = - \frac{(m-1)}{m} \left[ 1 + \rho_s \frac{\partial \ln g_R(\sigma)}{\partial \rho_s} \right]. \qquad (2.52)$$

La contribution résiduelle nette due à la formation de chaîne est donnée par

$$Z^{\text{chain}} = m \left( Z_s^{\text{chain}} - \lim_{\rho \to 0} Z_s^{\text{chain}} \right) = -(m-1) \rho \frac{\partial \ln g_R(\sigma)}{\partial \rho}. \qquad (2.53)$$

En intégrant par rapport à $\rho$ l'Eq. (2.53), on obtient la contribution de chaîne $A^{\text{chain}}$ à l'énergie libre résiduelle,

$$\frac{A^{\text{chain}}}{NkT} = -(m-1) \ln g_R(\sigma). \qquad (2.54)$$

L'Eq. (2.54) est l'expression que l'on retrouve dans la plupart des versions de l'équation d'état SAFT. La théorie TPT1 de Wertheim décrit bien les propriétés des chaînes flexibles pour des segments interagissant via des potentiels de sphères dures ou via des potentiels souples tels que le potentiel de Lennard-Jones. Cependant, cette théorie a certaines limitations : d'une part elle ne tient pas compte du positionnement des sites associatifs les uns par rapport aux autres sur une même molécule. Cette limitation peut être résolue en utilisant la théorie de Wertheim à l'ordre 2 (théorie TPT2). La théorie TPT1 ne peut pas décrire non plus la formation de clusters de molécules associées en boucle ou de dimères cycliques, qui jouent pourtant un grand rôle pour des molécules comme l'acide fluorhydrique ou les acides carboxyliques. Nous avons, au cours de nos travaux, utilisé une extension de la théorie TPT1 réalisée par Sear et Jackson$^{54,55}$ permettant de tenir compte de la formation de groupes de molécules associées en cycle, et nous l'avons étendue aux mélanges et appliquée aux cas des acides carboxyliques$^{56}$ (voir la Section 2.5).

## 2.4. Travaux Réalisés dans le Développement d'Equations d'Etat

### 2.4.1. Modification de la théorie de Barker - Henderson

Nous avons apporté plusieurs modifications à la théorie de Barker-Henderson afin d'améliorer la prédiction des propriétés thermodynamiques pour les potentiels d'interaction à répulsion dures et souples$^{57}$. Nous avons ainsi proposé un nouveau diamètre effectif $d$ dépendant de la température, et pris en compte les termes d'ordre supérieur de la théorie de perturbation afin que la théorie de Barker-Henderson prédise de manière exacte le second coefficient du Viriel correspondant au potentiel intermoléculaire choisi. Dans la nouvelle approche de perturbation, l'énergie libre est donnée par$^{57}$

$$\frac{A}{NkT} = \frac{A^{hs}}{NkT} + \sum_{i=1}^{\infty} \frac{A_i}{NkT}, \qquad (2.55)$$

où la somme des termes de perturbation est approximée par

$$\sum_{i=1}^{\infty} \frac{A_i}{NkT} = - \frac{2\pi\rho}{K^{hs}(1+\chi)} \int_{\sigma}^{\infty} \left( \exp\left(-\frac{K^{hs}(1+\chi)\phi(r)}{kT}\right) - 1 \right) g^{hs}(r) r^2 dr. \qquad (2.56)$$

$K^{hs}$ est la compressibilité isotherme du mélanges référence de sphères dures et $\chi$ est un facteur de correction dépendant de la densité réduite. L'équation d'état théorique (2.56) a été testée sur plusieurs potentiels d'interaction (puits carré, Yukawa, Sutherland, Mie, Exp-6, ...) et donne d'excellents prédictions des données de simulation moléculaire, comme le montre la Figure 2-8 pour le potentiel Exp-6.
<figure><img src="image_11.png" /></figure>

Figure 2-8 Equilibres liquide-vapeur d'un corps pur dont les interactions sont représentées par un potentiel de Buckingham Exp-6. Effet de la dureté du potentiel (partie répulsive du potentiel) (a) densités réduites des phases liquides et vapeurs, (b) pressions de vapeur réduites. Comparaison entre la théorie de perturbation et les données de simulation moléculaire (GEMC). Ref. ⁵⁷.

L'idée de modifier le terme d'ordre 2 a été reprise par la suite par le groupe de G. Jackson dans le cadre du développement de la nouvelle version de SAFT-VR Mie (SAFT gamma Mie) ⁵⁸⁻⁶¹, en particulier l'idée d'introduction le terme correctif $\chi$ qui permet une nette amélioration dans la prédiction des densités liquides et des pressions de vapeur. G. Jackson et col. ⁵⁸⁻⁶¹ ont rajouté le terme d'ordre 3 qui lui permet d'améliorer la prédiction du point critique.

### 2.4.2. Développement d'une équation d'état pour les sphères dures non additives

Nous avons proposé une extension du modèle BMCSLS ²³,²⁴ pour mélanges de sphères dures additives aux mélanges de sphères dures non additives ²⁵. Le mot "additif" signifie que les diamètres de collision croisée $\sigma_{ij}$ vérifient la règle de combinaison de Lorentz, c'est à dire $\sigma_{ij} = \frac{\sigma_i+\sigma_j}{2}$. Le cas des sphères non - additives correspond au cas plus général $\sigma_{ij} \neq \frac{\sigma_i+\sigma_j}{2}$. En effet, la règle de Lorentz n'a pas de réel fondement théorique et repose uniquement sur l'intuition que les molécules ou segments de molécule supposés sphériques se comporteraient comme des billes de billard et auraient des chocs élastiques. La loi de Lorentz (moyenne arithmétique) est appliquée la plupart du temps que ce soit avec l'équation d'état SAFT ou en simulation moléculaire, bien qu'en pratique, absolument rien de prouve que les paramètres croisées devraient suivre cette loi. Il a été montré par plusieurs auteurs, notamment par Y. Nezbeda⁶²,⁶³, que le fait d'autoriser que $\sigma_{ij} \neq \frac{\sigma_i+\sigma_j}{2}$ permet d'ajuster beaucoup plus facilement certaines propriétés comme les grandeurs de mélanges en phase liquide. Il est certain que les propriétés des phases liquides sont très sensibles au paramètre $\sigma_{ij}$ car les molécules en phase liquide sont très proches les unes des autres, donc la partie répulsive du potentiel joue un rôle important en phase liquide. L'extension que nous avons proposée est compatible avec le modèle BMCSLS ²³,²⁴ et est basée sur une modification du modèle MIX1. L'expression de l'énergie libre est donnée par ²⁵:
$$ \frac{A}{NkT} = \frac{A^{ahs}}{NkT} + \frac{2\pi}{3} \rho \sum_{i,j} x_i x_j (\sigma_{ij}^3 - \delta_{ij}^3) g_{ij}^{ahs}, \quad (2.57) $$

où $\delta_{ij} = \frac{\sigma_i + \sigma_j}{2}$; $A^{ahs}$ est l'énergie libre résiduelle du mélange de sphères dures additives, donnée par l'Eq. (2.27), et $g_{ij}^{ahs}$ est la valeur de contact de la fonction de distribution radiale entre les sphères i et j dans ce mélange de sphères additives. L'équation (2.57) est basée sur une perturbation par rapport au mélange de sphères additives. Ce modèle peut-être appliqué pour prédire les équilibres de phases dans les mélanges polymères + colloïdes. En effet, ce type de mélanges peut être représenté par un mélange binaire de sphères non additives, lequel permet de décrire les forces dites de "déplétion" entre les particules colloïdales crée par la présence chaînes de polymère.

Le mélange de sphères dures non additives peut aussi être utilisé comme système de référence dans une théorie de perturbation (donc une équation d'état) afin de décrire la propriétés de cas réels où $\sigma_{ij} \neq \frac{\sigma_i + \sigma_j}{2}$. C'est par exemple le cas des mélanges de gaz à très haute température et très haute pression$^{64}$. On observe pour ce type de mélange des équilibres dits gaz - gaz, qui sont en fait des équilibres de phases de type fluide-fluide, dus essentiellement aux répulsions croisées (entre molécules différentes), lesquelles sont défavorables par rapport aux répulsions entre molécules de même type (cas où $\sigma_{ij} > \frac{\sigma_i + \sigma_j}{2}$). D'autre part, le paramètre $\sigma_{ij} \neq \frac{\sigma_i + \sigma_j}{2}$ peut être très utile pour mieux représenter les équilibres liquide-liquide, car les propriétés des phases denses sont très sensibles à ce paramètre. Nous envisageons de développer une telle équation d'état dans un futur proche.

<figure><img src="image_12.png" /><figcaption>Figure 2-9 Equilibres fluide-fluide de mélanges binaires de sphères dures non-additives, avec $\sigma_2 = 0.1\sigma_1$. Les symboles correspondent aux données de simulation moléculaire obtenus pour différentes valeurs de non-additivité $\Delta_{12}$, tandis que les courbes correspondent aux prédictions de deux modèles : lignes en pointillés : théorie MIX1. Lignes continues : nouveau modèle$^{25}$.</figcaption></figure>

Le nouveau modèle (Eq. (2.57)) permet de prédire très précisément les équilibres fluide-fluide des mélanges de sphères dures non-additives, comme le montre la Figure 2-9.

## 2.5. Modification du terme associatif de SAFT. Application aux acides carboxyliques.

Dans le cadre d'un projet ANR (projet memobiol 2009-20012) portant sur l'étude des biocarburants, en particulier sur l'étude des propriétés physiques de molécules organiques oxygénées issues de biohuiles, nous avons travaillé sur la modélisation de mélanges comportant des acides carboxyliques. Ce travail a été réalisé en grande partie par Jiri Janecek, qui a été postdoc sous ma responsabilité durant le projet memobiol. La modélisation des équilibres de phases en présence d'acides carboxyliques n'est
pas simple car ces acides, notamment les plus courts, forment des dimères cycliques en phase vapeur ou à l'état dilué en phase organique non polaire. La Figure 2-10 montre le schéma d'un tel dimère cyclique pour l'acide acétique. Le terme cyclique signifie que les sites associatifs donneurs et accepteurs des deux molécules sont tous associés. Il existe différents types de modèles associatifs SAFT permettant de représenter les acides carboxyliques. Le plus simple est le modèle associatif appelé modèle "1A" suppose que l'acide possède un seul site associatif de type A. La dimérisation de l'acide est alors prise en compte par des associations A-A. Un tel modèle est réaliste en phase vapeur où les dimères prédominent. Cependant, le modèle associatif 1A ne permet pas de tenir compte de la formation de groupes de molécules associées au-delà du dimère (trimère, ...). Breil et al.⁶⁵ ont utilisé le modèle 1A dans l'approche CPA (équation d'état basée sur la combinaison d'une équation d'état cubique et du terme associatif de SAFT), et ont obtenu de très bons résultats pour les corps purs et les mélanges avec composés polaires.

Le deuxième modèle associatif est le modèle appelé 2B, comportant deux sites : un site accepteur A et un site donneur B (voir la Figure 2-7). Ce modèle est communément utilisé pour décrire les alcools simples. Concernant les acides carboxyliques, le modèle 2B décrit bien les phases liquides, mais ne permet pas de tenir compte de la formation de dimère cyclique en phase vapeur. Bien que ce modèle permet de décrire assez bien les équilibres de phases des mélanges binaires avec composés polaires, il est moins précis pour les mélanges d'acide carboxyliques + alcanes, et certaines propriétés thermodynamiques des acides pures, comme l'enthalpie de vaporisation ou le volume molaire en phase vapeur, est très mal représenté par le modèle 2B : l'enthalpie de vaporisation est largement surestimée⁵⁶.

Sear et Jackson⁵⁴,⁵⁵ ont étendu la théorie de Wertheim aux molécules formant des structures en boucle (n-mers cycliques). Cette théorie appliquée aux dimères cycliques est une extension du modèle 2B, et est appelé le modèle DBD (doubly bonded dimers). Ce modèle associatif possède également un site accepteur A et un site donneur B, cependant les deux sites sur la molécule sont suffisamment rapprochés (et non pas opposés comme sur la Figure 2-7) pour que des dimères cycliques puissent se former par deux liaisons A-B. Une fois le dimère formé, les deux molécules associées ne peuvent plus s'associer d'avantage car leurs sites sont déjà tous liés, et le dimère se comporte comme une nouvelle espèce moléculaire. C'est pourquoi la formation de dimères en phase vapeur entraîne une forme diminution du volume molaire en phase vapeur car l'acide occupe moins d'espace que prévu pour une température et une pression donnée si des dimères se forment: selon la loi des gaz parfaits, $V = NkT/P$, si N est divisé par deux, V est divisé par deux à T et P données. Supposons que les paramètres des modèles associatif 2B et DBD soient ajustés sur les données liquide-vapeur d'un acide carboxylique, les deux approches peuvent très bien décrire les volumes molaires liquides et les pressions de vapeur. Cependant, seul le modèle DBD peut décrire les volumes molaires en phase vapeur et les enthalpies de vaporisation. Le fait que le modèle classique 2B a tendance à surévaluer l'enthalpie de vaporisation peut s'expliquer par la relation de Clausius-Clapeyron qui s'applique le long de la courbe d'équilibre liquide-vapeur et qui est donnée par

$$ \Delta H_{vap} = T (v_{gaz} - v_{liq}) \frac{dP}{dT}. \qquad (2.58) $$

En effet, si le modèle 2B décrit très bien les données $v_{liq}$ et $dP/dT$, alors il surestime $\Delta H_{vap}$ si $v_{gaz}$ est surestimé (Figure 2-10).

Nous avons incorporé l'approche de Sear et Jackson dans l'équation d'état PC-SAFT⁶⁶,⁶⁷ et nous l'avons étendue aux mélanges (mélanges d'acides carboxyliques, mélanges d'un acide + un composé polaire associatif, mélanges d'un acide + un composé non polaire). Dans le cas d'un acide pure, on peut montrer que le terme associatif DBD s'écrit selon⁵⁴,⁵⁶

$$ \frac{A^{\text{assoc}}}{NkT} = \left( \ln X_0 - X_G + \frac{1}{2} + \frac{X_G^2}{2X_0} \right). \qquad (2.59) $$

$X_0$ est la fraction de molécules non associées (monomères), tandis que $X_G$ est la fraction de molécules non liée au site A. $X_0$ et $X_G$ sont déterminées par la résolution du système d'équations suivant

$$ -1 + \frac{X_G}{X_0} = \rho X_G \Delta_{AB}, \qquad (2.60) $$

$$ \frac{1}{X_0} - \frac{X_G^2}{X_0} = \rho X_0 \Phi_{AB/BA}. \qquad (2.61) $$
Les fonctions $\Delta_{AB}$ et $\Phi_{AB/BA}$ sont données par

$$ \Phi_{AB} = \Delta_{AB} \left[ \exp\left(\frac{\epsilon_{AB}}{kT}\right) - 1 \right] f_v, \qquad (2.62) $$

avec

$$ \Delta_{AB} = \kappa_{AB} \left[ \exp\left(\frac{\epsilon_{AB}}{kT}\right) - 1 \right] g^{hs}(d), \qquad (2.63) $$

<figure><img src="image_13.png" /><figcaption>Figure 2-10 Acide acétique en équilibre liquide-vapeur. La figure de droite correspond au modèle associatif 2B, tandis que la figure de gauche correspond au modèle associatif DBD impliquant la formation de dimère cyclique. Dans la figure du milieu, on compare les données expérimentales d'enthalpie de vaporisation ⁶⁸ aux prédictions des modèles PC-SAFT-2B (traits) et PC-SAFT-DBD (lignes).</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure 2-11 Facteur de compressibilité de la phase vapeur en équilibre avec la phase liquide, en fonction de la température, pour l'acide acétique. Les losanges correspondent aux données expérimentales de Miyamoto et al. ⁶⁹. Les cercles correspondent aux données expérimentales issues du livre de Timmermans ⁷⁰. Les lignes correspondent aux prédictions de trois équations d'état : PC-SAFT-1A, PC-SAFT-2B et PC-SAFT-DBD.</figcaption></figure>

$\Delta_{AB}$ est similaire à l'expression (2.51). Le paramètre $f_v$ permet de caractériser l'aptitude de l'acide à former des dimères : plus $f_v$ est grand, plus l'acide aura tendance à se dimériser. En prenant $f_v = 0$, on montre que le modèle DBD devient équivalent au modèle 2B. Le nouveau terme associatif étendu aux cas de mélanges d'acides et de composés avec associations croisés, a été proposé pour la première fois par Janecek et Paricaud ⁷¹. Le lecteur est invité à consulter l'article pour plus de détails. Les paramètres des acides purs ont été déterminés pour le nouveau modèle PC-SAFT-DBD en ajustant les volumes molaires des phases liquides et vapeurs en équilibre, ainsi que sur les pressions de vapeur. Les enthalpies de vaporisation et facteurs de compressibilité $Z = PV/NkT$ ont alors été prédits par le modèle. D'excellents résultats sont obtenus pour l'acide acétique comme le montre les Figure 2-10 et Figure 2-11. Le modèle associatif DBD permet non seulement d'améliorer la représentation des acides
purs, mais aussi les équilibres de phases des mélanges. La prise en compte de la dimérisation en phase vapeur permet d'améliorer sensiblement la représentation des courbes de rosée, pour les mélanges avec des alcanes (Figure 2-12) ou avec des composés polaires (Figure 2-13).

<figure><img src="image_15.png" /><figcaption>Figure 2-12 Equilibres liquide-vapeur du mélange n-hexane + acide acétique, à T=313.15K. Les symbols correspondent aux données expérimentales (Ref. 72) et les lignes sont les courbes calculées avec différents modèles : (a) modèle 1A, $k_{ij} = 0$ (lignes en pointillés) et $k_{ij} = 0.05$ (ligne noire). (b) modèles 2B ($k_{ij} = 0 / k_{ij} = 0.02$) et (c) modèle DBD ($k_{ij} = 0 / k_{ij} = 0.035$).</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure 2-13 Equilibre liquide-vapeur du mélange éthanol(1) + acide propénoïque (2) à T=323.2K. La ligne rouge en pointillés correspond au modèle 2B tandis que la ligne noire correspond au modèle DBD. Pour les deux modèles, $k_{ij} = 0$.</figcaption></figure>
# 3. Thermodynamique des Solutions de Polymère

Les équations d'état de type SAFT sont particulièrement adaptées pour décrire le comportement de phase des mélanges impliquant les polymères. Il est important de connaître les équilibres de phases liquide-vapeur et liquide-liquide de ces systèmes pour le design et l'optimisation de réacteurs de polymérisation<sup>22,73</sup>. Lors de la polymérisation en phase gazeuse, les particules de polymères sont suspendues par un débit gazeux (supercritique) de monomères. Les monomères sont absorbés dans les zones amorphes des grains de polymère, et la vitesse de réaction de polymérisation est directement liée à la solubilité des monomères. Dans le cas d'une copolymérisation, la composition finale des chaînes en copolymère est alors reliée aux solubilités des deux monomères. Des effets de synergie d'absorption sont observés et il est important que le modèle puisse en tenir compte<sup>74</sup>.

Lors d'une polymérisation en masse en phase liquide, le solvant est le monomère. Pour certaines conditions de température et de pression, une immiscibilité liquide-liquide peut être observée entre une phase riche en monomère et une phase riche en chaînes polymérique. Cette immiscibilité est due à des effets entropiques s'expliquant par la différence de taille entre la molécule de solvant et les chaînes, et peut engendrer un polymère de mauvaise qualité<sup>75</sup>. Il convient donc de connaître à l'avance les conditions pour lesquelles cette immiscibilité est observée afin de l'éviter.

Avant de présenter nos travaux sur ce sujet, nous rappelons brièvement le comportement de phase des mélanges binaires polymère-solvant en les présentant comme étant un cas extrême d'un mélange d'une molécule courte et d'une chaîne longue. L'immiscibilité liquide-liquide due à la différence de taille est observée pour les mélanges d'alcanes linéaires de type V et IV, tel que le mélange méthane + n-hexane. Nous avons expliqué l'immiscibilité liquide-liquide dans les mélanges d'alcanes, et montré qu'elle était due à des effets entropiques. Nous verrons ensuite quelques résultats de modélisation d'équilibres de phases pour des systèmes polymère-solvant. Ces travaux furent en grande partie réalisés dans le cadre de ma thèse et financés par BP Chemicals. La dernière section de ce chapitre sera consacrée à l'application de la thermodynamique continue aux systèmes polydispersés, en particuliers aux polymères polydispersés en terme de masse moléculaire (ou longueur de chaîne).

## 3.1. Types de comportement de phase

A cause de la grande variété de comportements de phase observés expérimentalement, il est très utile de classer ces comportements en différents types. Scott et van Konynenburg<sup>76</sup> ont proposé une classification simple des comportements de phase fluides des mélanges binaires, en six types principaux. Ils ont montré que l'équation de van der Waals est capable de décrire cinq de ces six types (Figure 3-1). Une discussion générale sur les types de comportement de phase a été proposée par plusieurs auteurs<sup>9,17,22,77</sup>. Nous allons voir ici les principaux comportements de phase de ces six types afin d'expliquer celui des solutions de polymère. Le comportement de phase en présence d'azéotrope n'est pas abordé. Le lecteur peut se reporter à l'excellente étude de Privat<sup>17</sup> sur les six types de diagrammes de phase avec prise en compte des lignes azéotropiques.

### 3.1.1. Comportement de type I

Dans le cas d'un mélange binaire de type I, la ligne critique liquide-vapeur est continue et joint les deux points critiques des corps purs. Il n'y a pas de zone d'équilibre liquide-liquide. Comme exemple de mélanges de type I, on peut citer les mélanges argon + krypton, méthane + éthane, méthane + diazote, dioxyde de carbone + dioxygène. Les mélanges de type I sont tels que les molécules présentent des tailles et des interactions similaires. Pour les mélanges méthane + alcanes linéaires, un comportement de type I est observé pour les mélanges méthane + éthane, méthane + propane, méthane + n-butane et méthane + n-pentane. Cependant, un comportement de type V est observé pour le mélange méthane + n-hexane<sup>77</sup>. Pour les mélanges d'éthane + n-alcane, le passage du type I au type V se produit à partir du n-nonadecane (C19). Pour les mélanges avec le propane, la transition se
produit autour du n-C40, cependant, la zone d'équilibre liquide-liquide à basse température devient difficile à repérer car elle devient métastable par rapport aux phases solides.

<figure><img src="image_17.png" /><figcaption>Figure 3-1 Les six principaux types de comportements de phase des mélanges binaires ⁷⁶, représentés dans les diagrammes pression-température. Les lignes continues correspondent aux pressions de vapeur des corps pur 1 et 2 et aux lignes triphasiques liquide-liquide-vapeur. Les courbes en pointillés sont les lignes de points critiques. Les triangles représentent points critiques terminaux supérieurs et inférieurs. Source: Paricaud ²²</figcaption></figure>

### 3.1.2. Comportement de type II

Les mélanges de type II se différentient des mélanges de type I par l'existence d'une région d'équilibre liquide-liquide à basse température (Figure 3-2). Une ligne triphasique est alors observée proche de la courbe de pression de vapeur du composé le plus volatile, et se termine en un point triple appelé point critique terminal supérieur (Upper Critical End Point). La région d'équilibre liquide-liquide apparaît dans les diagrammes de phase pression composition et température-composition. Dans le diagramme T-x, le point critique L-L est appelé la température critique supérieure de solution (Upper Critical Solution Temperature). Le comportement de type II est observé lorsque les deux composés ont des interactions de natures différentes : l'un est polaire ou quadripolaire, tandis que l'autre est apolaire. Quelques exemples de comportements de type II sont les mélanges CO₂ + n-octane, CO₂ + n-décane. Certains pensent que la plupart des mélanges de type I sont en fait des mélanges de type II dont la région liquide-liquide à basse température est masquée par les phases solides. Les mélanges liquides de type II ont souvent des enthalpies de mélange positives, caractérisant le fait que les interactions croisées (interactions 1-2) sont défavorables par rapport aux interactions entre molécules similaires (interactions 1-1 et 2-2). L'immiscibilité est due à des effets enthalpiques. Considérons l'enthalpie libre de mélange, $\Delta G_m = \Delta H_m - T\Delta S_m$ à température et pression fixées. A haute température, l'entropie de mélange est positive car dominée l'entropie idéale de mélange, $\Delta S_m^{ideal} = -R\sum_i n_i \ln x_i$. A haute température, le terme $-T\Delta S_m$ est négatif et domine le terme $\Delta H_m$. $\Delta G_m$ est donc négative et concave pour tout le domaine de composition $x_1$, et le mélange est miscible. A basse température, $\Delta H_m$ domine le terme $-T\Delta S_m$. $\Delta G_m$ change de courbure selon certaines compositions $x_1$, voire devient positive. Le mélange devient alors immiscible. Cela permet d'expliquer l'existence d'une température critique supérieure de solution pour ce type de mélanges. Dans le cas des mélanges CO₂-alcanes, le CO₂ possède un moment quadripolaire important. Donc les interactions CO₂-CO₂ sont beaucoup plus attractives que les interactions CO₂-alcane. De plus les forces de dispersion attractives alcane-alcane sont globalement plus importantes que les interactions CO₂-
alcanes, si l'alcane est suffisamment long (interactions site-site entre groupements CH₂). Ceci explique l'existence d'une immiscibilité.

<figure><img src="image_18.png" /></figure>

Figure 3-2 Diagramme de phase pression-température d'un mélange binaire de type II. Les lignes continues correspondent aux pressions de vapeur des corps purs 1 et 2 et à la ligne triphasique liquide-liquide-vapeur. Les courbes en pointillés sont les lignes de points critiques. Le triangle représente le point critique terminal supérieur. Les cercles représentent les points critiques des corps purs. Source: Paricaud ²².

<figure><img src="image_19.png" /></figure>

Figure 3-3 Diagrammes de phase pression-composition (a) et température-composition (b) d'un mélange binaire de type II, correspondant respectivement aux coupes $T = T_1$ et $P = P_1$ de la Figure 3-2. Les lignes continues correspondent aux courbes de coexistence liquide-vapeur et liquide-liquide. Le cercle dans la Figure (b) représente un point critique liquide-liquide de type UCST. Source: Paricaud ²².

Il faut noter que les mélanges de type II ont souvent un comportement hétéroazéotropique, qui n'est pas ici représenté dans les diagrammes de phase. L'existence d'azéotropes est favorisée lorsque les courbes de pression de vapeur des corps purs sont proches. Lorsque ces courbes se croisent à un point dit point de Bancroft, un azéotrope existe nécessairement pour certaines conditions. Notez que les azéotropes peuvent aussi être observés pour les mélanges de type I (eau-ethanol).

### 3.1.3. Comportement de type V

Le diagramme pression- température du comportement de type V est représenté dans la Figure 3-4. Le type V est caractérisé par l'apparition d'une zone d'immiscibilité liquide-liquide autour du point critique du composé le plus volatile. La ligne continue de points critiques liquide-vapeur observée dans les types I et II devient discontinue pour le type V, et une ligne triphasique liquide-liquide-vapeur apparaît, délimitée par un point critique terminal supérieur et un point critique terminale inférieur (UCEP, LCEP).
<figure><img src="image_20.png" /></figure>

Figure 3-4 Diagramme de phase pression-température d'un mélange binaire de type V. Les lignes continues correspondent aux pressions de vapeur des corps purs 1 et 2 et à la ligne triphasique liquide-liquide-vapeur. Les courbes en pointillés sont les lignes de points critiques. Les triangles représentent les points critiques terminaux inférieur et supérieur. Les cercles représentent les points critiques des corps purs. Source: Paricaud 22.

Pour montrer le comportement de phase du type V, nous avons considéré deux isothermes : une isotherme en T₁ située au point critique terminal inférieur (LCEP) et une isotherme en T₂ située entre les points LCEP et UCEP. On observe clairement la zone d'immiscibilité L-L dans le diagramme P-X (b)) qui se termine à haute pression en un point critique (UCSP). On peut citer comme exemple de mélanges de type V le mélange méthane + n-hexane.

### 3.1.4. Comportement de type IV

La différence entre les comportements de types V et IV est l'existence d'une seconde région d'équilibre liquide-liquide à basse température pour le type IV, similaire à celle du type II (Figure 3-5).

<figure><img src="image_21.png" /></figure>

Figure 3-5 Diagramme de phase pression-température d'un mélange binaire de type IV. Les lignes continues correspondent aux pressions de vapeur des corps purs 1 et 2 et aux deux lignes triphasiques liquide-liquide-vapeur. Les courbes en pointillés sont les lignes de points critiques. Les triangles représentent les points critiques terminaux inférieurs et supérieurs. Les cercles représentent les points critiques des corps purs. Source: Paricaud 22.
<figure><img src="image_22.png" /></figure>

Figure 3-6 Diagrammes pression-composition d'un mélange binaire de type IV aux températures $T_1$ (a) et $T_2$ (b) correspondant à la Figure 3-6. Les lignes continues correspondent aux courbes de coexistence et aux lignes tri-phasiques liquide-liquide-vapeur. Les cercles représentent les points critiques V-L et L-L de mélange. Le triangle représente le point critique terminal.

<figure><img src="image_23.png" /></figure>

Figure 3-7 Diagramme température-composition d'un mélange binaire de type IV à la pression $P_1$ correspondant à la Figure 3-6. Les lignes continues correspondent aux courbes de coexistence et à la ligne tri-phasique liquide-liquide-vapeur. Les cercles représentent les points critiques V-L et L-L de mélange.

Des exemples de diagrammes de phase pression-composition (Figure 3-6) et température compositions sont donnés ci-dessous. Deux régions d'équilibres liquide-vapeur apparaissent dans le diagramme de phase isotherme à $T = T_1$. Ces deux régions fusionnent à la température $T_2$ correspondant à l'un des points terminaux UCEP. D'autre part, deux régions d'immiscibilité liquide-liquide peuvent être observées dans le diagramme température-composition à la pression $P_1$, l'une se terminant par un maximum en température (point critique de type UCST), et l'autre se terminant par un minimum (point critique de type LCST). On peut citer le mélange $CO_2 + n$-tridécane comme exemple de mélange de type IV. Certains auteurs, comme pour le type II, prétendent qu'en pratique il n'existe que les mélanges de type IV et non pas de type V car pour les mélanges supposés de type V, la zone liquide-liquide à basse température existerait mais serait métastable par rapport aux phases solides.

L'immiscibilité liquide-liquide avec point critique de type LCST observée dans les mélanges de types V et IV autour du point critique du composé le plus volatile, peut s'expliquer par le fait que les enthalpies et entropies de mélange soient négatives ⁷⁵. Selon $\Delta G_m = \Delta H_m - T\Delta S_m$, à basse température, le terme enthalpique $\Delta H_m$ domine et est concave et négatif. Donc il n'y a pas d'immiscibilité. Lorsque $T$ augmente, le terme entropique $-T\Delta S_m$ est positif et convexe, et commence à dominer. $\Delta G_m$ devient convexe pour certaines compositions, voire positive, et le mélange devient immiscible. Cette immiscibilité est donc due à des effets entropiques. Paricaud et al. ⁷⁵ ont expliqué le comportement de type V et IV par le fait que le terme $-T\Delta S_m$ soit négatif dans les mélanges de type IV ou V, lorsqu'il y a une forte différence de taille entre les deux composés (cas typique du mélange n-hexane + méthane ou des systèmes polymère-solvant). Considérons un mélange
d'un solvant représenté par une sphère et d'une chaîne de sphères représentant un alcane linéaire long. La négativité de l'entropie de mélange est due au fait qu'une chaîne entourée de petites molécules de solvant, a tendance à fortement attirer des molécules de solvant par les forces attractives de dispersion, ce qui entraîne une contraction du système (le volume molaire diminue). Or on sait qu'une augmentation de la densité (N/V) est suivie d'une diminution d'entropie (le système est plus ordonné). Si cependant on augmente la concentration en chaîne dans le mélange, les forces répulsives entre les chaînes vont dominer les forces attractives, ce qui va entrainer une ré-augmentation du volume molaire. On observe alors un maximum de la densité en fonction de la concentration en chaîne. L'entropie du mélange à ce maximum en densité (minimum en volume molaire) est alors minimale et l'entropie de mélange est négative. La Figure 3-8 illustre ce phénomène à l'échelle microscopique.

<figure><img src="image_24.png" /><figcaption>Figure 3-8 Schéma représentant le mélange binaire solvants (sphères grises) + chaîne (chaîne de dix sphères noires), pour expliquer l'évolution de la densité du mélange en fonction de la composition en chaîne à pression et température fixées. D'après Paricaud et al. ²²,⁷⁵.</figcaption></figure>

### 3.1.5. Comportement de type III

Lorsque les deux composés sont de nature très différente (mélange d'un composé fortement polaire comme l'eau et d'un alcane tel que le méthane) et sont immiscibles sur un grand domaine de température et de pression, le comportement de phase est la plupart du temps de type III. On peut observer une progression du type IV au type III à l'aide d'une équation d'état et en augmentant l'immiscibilité entre les deux composés via un paramètre binaire d'interaction croisée ($k_{12}$). Le type III ne présente qu'une seule ligne de points critiques, correspondant à la fusion des deux lignes critiques observées dans le type IV. Les points terminaux UCEP et LCEP ont fusionné et finalement une seule ligne triphasique se terminant en un point UCEP est observée dans le diagramme de phase pression-température.

La courbe de points critiques pour les mélanges de type III peuvent présenter un maximum ou un minimum comme sur la Figure 3-9, mais ce n'est pas forcément le cas: la ligne critique peut monter directement vers les hautes pressions dans extremum, lorsque les deux composés sont très immiscibles. Comme exemples de mélanges de type III on peut citer les mélanges méthane + n-alcane très long, les mélanges $CO_2 + \text{tetradécane}$, le mélange eau+$CO_2$⁷⁸.
<figure><img src="image_25.png" /></figure>

Figure 3-9 Diagramme de phase pression-température d'un mélange binaire de type III. Les lignes continues correspondent aux pressions de vapeur des corps purs 1 et 2 et à la ligne triphasique liquide-liquide-vapeur. La courbe en pointillés est la ligne de points critiques. Le triangle représente le point critique terminal supérieur. Les cercles représentent les points critiques des corps purs. Source: Paricaud ²².

<figure><img src="image_26.png" /></figure>

Figure 3-10 Diagrammes de phase température-composition d'un mélange binaire de type III aux pressions P₁ (a) et P₂ (b) correspondant à la Figure 3-10. Les lignes continues correspondent aux courbes de coexistence fluide-fluide. Les cercles représentent les points critiques V-L et L-L de mélange.

### 3.1.6. Comportement de type VI

Le comportement de type VI est assez proche du comportement de type II. La différence avec les types II et VI est l'existence pour le type VI d'une seconde ligne de points critiques liquide-liquide délimitant la zone d'immiscibilité à basses températures. Ces deux lignes peuvent se rejoindre à haute pression comme sur l'exemple donné par la Figure 3-11. Il existe une seule ligne triphasique délimitée par des points critiques terminaux supérieur et inférieur (LCEP et UCEP).

Les mélanges de type VI ont la caractéristique d'avoir des régions d'immiscibilité liquide-liquide en boucle fermée dans les diagrammes température - composition (Figure 3-12). L'eau est en général d'un des composés dans ce type de mélanges et l'autre composé en général présente des liaisons hydrogènes croisées avec l'eau. Comme exemples de type VI, on peut citer les mélanges eau + 2-butanone, water + 2-butanol, water + nicotine. Ces mélanges présentent des points critiques de type UCST et LCST. La nature des points critiques LCST est différente de celle des points LCST dans les mélanges de type IV et V, tandis que la nature du point UCST est similaire à celle d'un point UCST d'un type II. Nous avons comparé les deux types de comportement VI et V en étudiant les propriétés de mélange ⁷⁵ de mélanges de sphères dures associatives. Dans le type VI, l'existence d'un LCST est
due à la présence de liaisons hydrogènes croisées: à basses températures en dessous de la température LCST, les interactions croisées par liaison hydrogène dominent les autres types d'interaction, et le mélange est miscible.

<figure><img src="image_27.png" /></figure>

Figure 3-11 Diagramme de phase pression-température d'un mélange binaire de type VI. Les lignes continues correspondent aux pressions de vapeur des corps purs 1 et 2 et à la ligne triphasique liquide-liquide-vapeur. La courbe en pointillés est la ligne de points critiques. Le triangle représente le point critique terminal supérieur. Les cercles représentent les points critiques des corps purs. Source: Paricaud ²².

<figure><img src="image_28.png" /></figure>

Figure 3-12 Diagrammes de phase température-composition d'un mélange binaire de type VI aux pressions $P_1$ et $P_2$ correspondant à la Figure 3-11. Les lignes continues correspondent aux courbes de coexistence liquide-vapeur et liquide-liquide. Les cercles représentent les points critiques V-L et L-L de mélange.

Lorsque la température augmente, ces liaisons hydrogène croisées commencent à se casser à cause de l'agitation thermique. Une immiscibilité L-L apparaît due au fait que les interactions en molécules de même espèce chimique (interactions 1-1 et 2-2) sont plus favorables que les interactions croisées. A hautes températures, l'entropie idéale de mélange domine à nouveau le terme d'enthalpie de mélange, et le mélange redevient miscible.
<figure><img src="image_29.png" /></figure>

Figure 3-13 Propriétés thermodynamiques de mélange d'un mélange symétrique binaire de type VI, de deux composés 1 et 2 similaires sphériques avec un site associatif. (a) Enthalpie libre de mélange, (b) enthalpie de mélange, (c) entropie de mélange, (d) volume molaire. Seules les associations croisées sont considérées, et l'énergie de dispersion $ϵ_{12}$ est telle que $ϵ_{12} < \sqrt{ϵ_{11}ϵ_{22}}$ (forces de dispersion croisées défavorables). Trois températures sont considérées, telles que $T_3^* < T_{LCST}, T_{LCST} < T_4^* < T_{UCST}$ et $T_5^* > T_{UCST}$. A basses températures en dessous du point LCST, $ΔH_m$ est négative (à cause des associations croisées) et domine le terme $TΔS_m$ lui-aussi négatif (due à la diminution de l'ordre à cause des associations croisées). Le mélange est donc miscible. Quand $T$ augmente, $TΔS_m$ devient positif mais est toujours dominé par $ΔH_m$ qui est cette fois-cis positive : le mélange devient immiscible. A hautes températures au dessus du point UCST, $TΔS_m$ est largement positive car dominée par l'entropie idéale de mélange, et domine $ΔH_m$ : le mélange redevient miscible. D'après Paricaud et al.⁷⁵.

### 3.1.7. Comportement de phase des solutions de polymères

En supposant que les chaînes polymériques sont monodispersées, les solutions de polymères se comportent comme des mélanges binaires de deux composés de tailles très différentes (petite molécule + chaîne longue). Le comportement de phase de tels mélanges est généralement de type IV ou III et la présence d'immiscibilité est due à la grande différence de taille entre le polymère et les molécules de solvant. Le type VI peut également être observé pour certains systèmes aqueux : pour les chaînes courtes de polyéthylène glycol en solution aqueuse, on peut observer des courbes de coexistence liquide-liquide en boucle fermée, typique du comportement de type VI.

<figure><img src="image_30.png" /></figure>

Figure 3-14 Diagrammes de phase température-composition d'un mélange binaire de type IV (a) et d'un mélange polymère-solvant (b). Dans le cas (b), la première région d'équilibre L-V est confondue avec l'axe des ordonnées ($W_{polymer}$). D'après Paricaud ²².

Dans la Figure 3-14, le diagramme de phase $T-x$ d'un mélange IV est comparé à celui d'une solution de polymère (typiquement, un hydrocarbure + polyéthylène). On observe, pour les deux
mélanges, deux régions d'immiscibilité liquide-liquide. Cependant les courbes L-L sont très asymétriques pour les solutions de polymère, et l'une des régions d'équilibre liquide-vapeur n'est pas visible car confondue avec l'axe des ordonnées. Pour la deuxième région liquide-vapeur, la courbe de rosée est elle-aussi confondue avec l'axe des ordonnées.

## 3.2. Modélisation des équilibres de phases de solutions de polymère

Nous avons utilisé l'équation d'état SAFT-VR combinée avec le potentiel à puits carré (square-well), développée par Gil-Villegas et al.⁷⁹, pour modéliser les équilibres de phases liquide-vapeur et liquide-liquide de systèmes polymères + solvant. Une grande partie de ce travail a été effectuée durant ma thèse et j'ai continué de coopérer quelques temps avec l'Imperial College sur ce sujet.

### 3.2.1. Solutions de polyéthylène

L'obtention des paramètres SAFT-VR des corps purs est généralement faite par régression sur les données expérimentales liquide-vapeur (masses volumiques des phases liquides et vapeurs et pressions de vapeur). Cependant, ces données expérimentales n'existent pas pour le polyéthylène car ce dernier n'est évidemment pas vaporisable. Les paramètres SAFT-VR du polyéthylène peuvent être alors obtenus par extrapolation de corrélations faites sur les paramètres des n-alcanes, à des masses moléculaires élevées ⁷³.

<figure><img src="image_31.png" /><figcaption>Figure 3-15 Paramètres SAFT-VR de corps purs pour les n-alcanes ⁷³ (nombre de segments *m*, et paramètres *σ*, *ε*, *λ* du potentiel à puits carré). Corrélation linéaire des paramètres en fonction de la masse molaire.</figcaption></figure>
<figure><img src="image_32.png" /><figcaption>Figure 3-16 Modélisation des équilibres liquides-vapeur et liquide-liquide des solutions de n-pentane + polyéthylène. Les symboles correspondent aux données expérimentales<sup>80</sup> tandis que les lignes correspondent aux prédictions du modèle SAFT-VR avec $k_{ij}=0$.</figcaption></figure>

Nous avons ensuite prédit les équilibres de phases de plusieurs mélanges (solvant + polyéthylène). La représente le diagramme de phase $T-x$ du mélange polyéthylène + n-pentane. On observe clairement les régions d'équilibres liquide-vapeur et liquide-liquide. Le modèle SAFT-VR prédit très bien les données expérimentales avec $k_{ij} = 0$.

### 3.2.2. Solubilités de gaz dans le polyéthylène

L'entreprise BP Chemicals, qui a financé en partie mon travail de thèse, était particulièrement intéressée par des modèles thermodynamiques capables de prédire les solubilités de gaz (monomères et co-solvants) dans les polymères polyoléfines semi-cristallins, pour deux applications différentes. D'une part, il est important de connaître la solubilité des monomères dans les grains de polymère lors de la polymérisation en phase gazeuse / supercritique, car celle-ci influence les vitesses de réaction de polymérisation. D'autre part, il est important de connaître la quantité de monomères restant dans les grains de polymères stockés, pour des raisons de sécurité (atmosphères explosives). Nous avons travaillé essentiellement sur le polyéthylène haute et basse densité (HDPE et LDPE). Pour prédire la quantité absorbée dans les grains, on a fait trois hypothèses principales :

- on suppose que les molécules de monomères ou de co-solvants (hydrocarbures) s'absorbent au niveau des zones amorphes du polymère semi-cristallin.

- on suppose que l'absorption dans les zones amorphes s'apparente à un équilibre liquide-vapeur, car la structure des zones amorphes ressemble à celle d'un polymère fondue.

- En réalité, la solubilité est limitée pour deux raisons : d'une part, les zones amorphes ne sont pas parfaitement homogènes car elles présentent certaines zones plus structurées que d'autre. D'autre part, la solubilité du monomère est limitée car le gonflement du polymère suite à l'absorption de monomères est limité par un effet élastique : les chaînes dites « liantes » qui sont piégées par deux zones cristallisées consécutives sont étirées lors du gonflement des zones amorphes, et se comportent comme des élastiques qui s'opposent à l'expansion de ces zones amorphes. En pratique, on tient compte du fait que la solubilité dans les zones amorphes d'un polymère semi-cristallin est limitée de manière effective, par l'utilisation d'un paramètre binaire $k_{ij}$ positif (interactions monomère-polymère plus faible que prévues, entrainant une solubilité plus faible).
<figure><img src="image_33.png" /><figcaption>Figure 3-17 Cristallinité de plusieurs types de polyéthylène en fonction de la température. Les symboles correspondent aux données expérimentales des références. Les lignes correspondent aux prédictions du modèle (a) polyéthylène produits à partir de catalyseurs de type métallocène ; (b) polyéthylènes produits à partir de catalyseurs de type Ziegler-Natta.</figcaption></figure>

Nous avions dans un premier temps développé un modèle basé sur une théorie de Flory, pour prédire le taux de cristallinité des polyéthylènes à haute et basse densité ⁷³. Ce modèle décrit de manière très précise les données expérimentales de cristallinité, comme le montre la Figure 3-17. La solubilité d'un gaz dans un polyéthylène semi-cristallin, définie comme étant la masse en kg de gaz absorbée par kg de polyéthylène semi-cristallin, est estimée à l'aide de l'équation ci-dessous :

$$ Sol. = (1 - w_{cris}) \frac{w_{gaz}^{(am)}}{w_{PE}^{(am)}}, \qquad (3.1) $$

où $w_{cris}$ est la cristallinité à la température T, et $w_{gaz}^{(am)}$, $w_{PE}^{(am)}$ les fractions massiques en gaz et en polyéthylène dans les zones amorphes du polyéthylène. On suppose en effet que les molécules de gaz ne sont absorbées que dans les zones amorphes. Celles-ci sont supposées se comporter comme un polymère fondu (liquide), et l'on estime $w_{gaz}^{(am)}$ et $w_{PE}^{(am)}$ en résolvant un équilibre liquide-vapeur avec l'équation d'état SAFT-VR. Il faut noter que la densité et la micro-structure des zones amorphes ne sont pas homogènes et n'ont pas localement la même densité ni la même structure que celle d'un polymère fondue : la transition entre zone amorphe et zone cristallisé n'est pas nette mais plutôt continue. De plus, les effets élastiques dus à la présence de chaînes liantes, s'oppose au gonflement du polymère donc à l'absorption de molécule de gaz. Aussi, le modèle tend à surévaluer la solubilité de gaz dans ces zones ⁸¹,⁸²,⁸³. Nous avons modélisé cet effet de manière implicite via un paramètre d'interaction binaire positif $k_{ij}$ entre le polyéthylène et le gaz ⁷³,⁷⁴ lequel permet de réduire les attractions entre les molécules de gaz et les chaînes, donc de prédire des solubilités de gaz plus faibles. Une bonne description des données expérimentales de solubilité est obtenue avec cette approche, comme le montre la Figure 3-18.
<figure><img src="image_34.png" /></figure>

Figure 3-18 Solubilité du 1-butène et du 1-hexène dans le polyéthylène semi-cristallin, en fonction de la pression et de la température. Comparaison entre les données expérimentales et les prédictions du modèle SAFT-VR + cristallinité.

Afin de tenir compte de manière plus réaliste de l'effet d'élasticité et du fait que la masse volumique des zones amorphes situées entre deux zones cristallisées soient supérieures à celle du polymère supposé fondu, il peut être envisagé de prendre en compte une surpression⁸¹,⁸², ou un terme d'énergie libre propre aux effets élastiques et qui se rajoute à l'équation d'état (terme de Langevin⁵¹).

## 4. Systèmes Polydispersés et Thermodynamique
Continue

Les mélanges dits "polydispersés" sont rencontrés dans de nombreuses applications industrielles. Quelques exemples de mélanges polydispersés sont les fluides pétroliers, les solutions de polymères (polydispersité des masses moléculaires) et les systèmes colloïdaux (polydispersité des tailles, des charges, ...). Dans le cas des polymères, la polydispersité a une grande influence sur les équilibres liquide-liquide, et il est donc important d'en tenir compte. Deux principales difficultés sont rencontrées pour traiter ce genre de systèmes. D'une part, ces systèmes impliquent de nombreux constituants, ce qui nécessite de résoudre un grand nombre d'équations non linéaires lors de la résolution des équilibres de phases. D'autre part, la composition de ces systèmes est souvent caractérisée par une distribution continue (et non pas discrète) des fractions molaires. On peut par exemple mentionner la distribution des points d'ébullition pour une fluide pétrolier (courbe TBP), mesurée avec une colonne de distillation, ou la distribution des masses molaires (c'est à dire des longueurs de chaîne) d'un polymère, mesurée par chromatographie d'exclusion stérique. Ainsi, la thermodynamique continue constitue un ensemble de méthodes permettant de traiter ces mélanges dont la composition est donnée par une fonction de distribution continue. Le lecteur est dirigé vers l'excellente revue de Browarzik et Kehlen⁸⁴ sur la thermodynamique continue. Il existe deux
méthodes distinctes pour traiter les mélanges multi-constituants polydisersés. Nous avons utilisé ces
deux méthodes au cours de nos travaux de recherche.

**4.1. Distributions discrètes et pseudoconstituants.**

Une première méthode dite "lumping" consiste à représenter le mélange par un certain nombre de
constituants discrets ou pseudoconstituants. Cette méthode est très utilisée dans l'industrie pétrolière.
La thermodynamique classique et les algorithmes de flash peuvent alors être utilisés pour les mélanges
de pseudoconstituants. La distribution continue TBP des points d'ébullitions est traditionnellement
discrétisée en coupes pétrolières de telle sorte que la coupe $C_k$ contienne les constituants dont les
points d'ébullition sont compris entre les points d'ébullition des n-alcanes $C_{k-1}$ et $C_k$. Pedersen et col.
⁸⁵⁻⁸⁸ ont utilisé l'équation d'état SRK et l'approche par pseudoconstituants pour décrire les équilibres
de phases des fluides pétroliers (huiles, gaz à condensat).

<figure><img src="image_35.png" /><figcaption>Figure 4-1 Courbe de dépôt liquide d'un gaz à condensat à composition globale constante calculée avec différents modèles : + données expérimentales. (a) T=353.15K, T=378.15K. Lignes continues: courbe calculées avec notre model ; - - - courbe calculée prédite avec le modèle d'Avaullée et al.⁸⁹ ; ... courbe calculée avec le modèle de Pedersen et al.⁹⁰.</figcaption></figure>

<figure><img src="image_36.png" /><figcaption>Figure 4-2. Composition de la coupe lourde du gaz à condensat, correspondant aux mesures de dépôt liquide de la Figure 4-1.</figcaption></figure>

Lors de mon stage de DEA effectué au laboratoire de thermodynamique de l'ENSIC en coopération
l'IFP, j'ai amélioré une approche développée précédemment par Sportisse et col. ⁹¹ et basée sur des
fonctions de distribution discrètes, pour décrire la composition des gaz à condensat. Chaque coupe Ck
au-delà du $C_{20}$ était représentée par un alcane, un alkylbenzène et un polyaromatique. Ce modèle qui
permet de décrire la composition de la coupe lourde a été combiné avec l'équation d'état de Peng- Robinson ⁹². Nous avons obtenu une très bonne description des enveloppes de phases et des courbes de dépôts liquides de gaz à condensat (Jaubert et al. ⁹³). Un exemple d'application de ce modèle est donné dans la Figure 4-1.

## 4.2. Distributions continues

La seconde approche consiste à considérer une distribution continue des fractions molaires des constituants du mélange. On parle alors de thermodynamique continue. L'équation d'état (la pression, l'énergie libre) devient une fonctionnelle de la fonction de distribution (une fonctionnelle est une fonction de fonction). Le principal avantage d'une telle approche est qu'elle permet de diminuer considérablement le nombre d'équations à résoudre lors de la résolution d'équilibres de phases (calcul de flash). La thermodynamique continue a été utilisée pour la première fois par Katz and Brown ⁹⁴ pour calculer un équilibre liquide-vapeur d'une mélange polydispersé en utilisant la loi de Raoult. Selon le choix du modèle thermodynamique et la forme analytique de la fonction de distribution (fonction gaussienne, gamma, ..), il peut être possible d'obtenir des solutions analytiques à l'équilibre de phase en particulier pour le calcul des points de bulles et de rosée : dans ce cas, les compositions des deux phases en équilibre sont également décrites par des fonctions de distribution dont les expressions analytiques peuvent être déterminées. Cependant, les calculs dits de Flash pour une composition globale du mélange donnée, ne permettent pas en général d'obtenir des solutions analytiques. Il existe alors des méthodes numériques dites des moments pour déterminer des distributions des fractions molaires dans les phases à l'équilibre. Cotterman et Prausnitz ⁹⁵ ont proposé une telle méthode pour les calculs de flash de mélanges polydispersés. Hu et Prauznitz ⁹⁶⁻⁹⁸ ont étendu la méthode aux systèmes polydispersés dits semi-continus (distribution continue de composés + solvants discrets).

Le cas spécifique des systèmes polymères-solvants a été traité par de nombreux auteurs. Schulz ⁹⁹ fut l'un des premiers à utiliser la thermodynamique continue pour ses systèmes, en représentant la composition des chaînes polymériques par une fonction de distribution de type gamma. La fonction de distribution gamma des masses moléculaires est d'ailleurs également appelée la distribution de Schulz-Flory. On peut mentionner également les travaux de Koningsveld ¹⁰⁰, Huggins ¹⁰¹, Solk ¹⁰², Heidemann et col. ¹⁰³⁻¹⁰⁷. Jog et Chapman ¹⁰⁸ ont proposé un algorithme combiné avec l'équation d'état SAFT pour déterminer les courbes de polyspersité (courbes cloud et shadow) qui correspondent aux compositions d'une phase en son point de bulle et de l'autre phase en équilibre, pour un système polydispersé. Beheme et al. ¹⁰⁹ ont utilisé une méthode similaire combinée avec l'équation d'état PC-SAFT, pour prédire les équilibres liquide-liquide de systèmes polymère-solvant.

## 4.3. Description de la méthode des moments

Nous avons utilisé la thermodynamique continue à divers systèmes de polymères ¹¹⁰ et de colloïdes ²⁵,¹¹¹. Voici une description rapide de la méthode pour calculer les courbes cloud et shadow d'un mélange polydispersé. Prenons l'exemple d'un mélange contenant un solvant (composé 1) et un polymère polydispersé (mélanges de chaînes de longueur différente) dont la composition globale est déterminée par une fonction de distribution connue. Nous avons étudié un cas particulier pour lequel une chaîne de polymère de longueur *m* est représentée par une chaîne flexible de *m* segments sphériques et adjacents, tandis que le solvant est ici simplement représenté par une sphère. On suppose que la molécule de solvant et les segments de chaînes ont tous le même diamètre *σ*. Prenons par exemple une distribution de Schulz-Flory (fonction gamma) pour décrire la distribution globale des longueurs de chaîne. Cette fonction est définie selon

$$X^{(0)}(m) = \frac{b^a}{\Gamma(a)} m^{a-1} e^{-bm}. \qquad (4.1)$$

où $X^{(0)}$ est la fraction globale en nombre de chaînes ayant comme longueur *m* ; *b* et *a* sont deux paramètres caractérisant la fonction de distribution et sont reliés à la longueur moyenne des chaînes
$\langle m \rangle_n$ et l'indice de polydispersité $I_p$, tels que $\langle m \rangle_n = a/b$ et $I_p = (a+1)/a$. Un index $I_p =1$ correspond à un polymère monodispersé (toutes les chaînes ont la même longueur). Plus $I_p$ est grand, plus le polymère a une distribution large des masses moléculaires. L'idée principale de la thermodynamique continue est d'exprimer l'équation d'état, c'est-à-dire l'énergie libre résiduelle de Helmholtz et la pression, en fonction de la fonction de distribution, puis de faire apparaître dans les expressions mathématiques les moments de la fonction de distribution. Lors d'un équilibre entre deux phase, on doit considérer alors trois fonctions de distribution différentes: la fonction de distribution globale $X^{(0)}$; et deux autres caractérisant la composition des deux phases en équilibre (comme lors d'un calcul de Flash biphasique, trois vecteurs de composition doivent être considérés : la composition globale et celles des deux phases en équilibre).

Considérons la distribution des densités (nombre de molécules par unité de volume) des
chaînes de longueur m dans une phase donnée, selon

$$
\rho_m(m) = \langle \rho_m \rangle_0 X(m). \tag{4.2}
$$

où $\langle \rho_m \rangle_0$ est le moment de densité d'ordre zéro de cette distribution. $\langle \rho_m \rangle_0$ correspond à la densité totale des chaînes de polymère. Les moments de densité d'ordre $k$ sont définis par

$$
\langle \rho_m \rangle_k = \int dm \rho_m(m) m^k . \qquad (4.3)
$$

La méthode des moments consiste à exprimer la densité d'énergie libre résiduelle, en fonction des moments de densité. La partie idéale de la densité d'énergie libre est donnée par

$$
\frac{A^{\text{ideal}}}{VkT} = \frac{f^{\text{ideal}}}{kT} = \rho_1 [\ln \rho_1 - 1] + \int dm \rho_m(m) [\ln \rho_m(m) - 1]. \quad (4.4)
$$

où $\rho_1$ est le nombre de densité du solvant. Le terme idéal donné par l'Eq. (4.4) est une fonctionnelle de la fonction de distribution. La contribution idéale au potentiel chimique $\mu^{\text{ideal}}$ de la chaîne de longueur $m$ est donnée par la dérivée de cette fonctionnelle par rapport à la fonction $\rho_m(m)$, c'est-à-dire¹¹⁰

$$
\mu^{\text{ideal}}(m) = kT \left( \frac{\partial \left( \int dm \rho_m \left[ \ln \rho_m(m) - 1 \right] \right)}{\partial \rho_m(m)} \right)_{\rho_1} = kT \ln \rho_m(m). \quad (4.5)
$$

La longueur d'onde de de Broglie a été enlevée des expressions car elle n'a pas d'influence sur les équilibres de phase. Lorsqu'une équation d'état de type SAFT est utilisée pour des systèmes polydispersés, la densité d'énergie libre résiduelle $f^{res}$ peut en général s'écrire uniquement en fonction de la température, du nombre de densité du solvant ($\rho_1$), et de certains moments de densité $\langle \rho_m \rangle_k$ définis par l'Eq. (4.3). Pour une polydispersité des longueurs de chaîne $m$, seuls les moments d'ordres 0 et 1 apparaissent dans les expressions analytiques. Ainsi on peut monter que $f^{res}$ est uniquement une fonction de $T$, $\rho_1$, $\langle \rho_m \rangle_0$, et $\langle \rho_m \rangle_1$. Par exemple, si l'on suppose que les segments sphériques ont tous le même diamètre $\sigma$, le terme de sphère dure $f^{res,hs}$ peut s'écrire selon¹¹⁰

$$
\frac{f^{\text{res,hs}}}{kT} = \left( \sum_{i=1}^{n} \rho_i m_i \right) \frac{4\eta - 3\eta^2}{(1-\eta)^2} = (\rho_1 + \langle \rho_m \rangle_1) \frac{4\eta - 3\eta^2}{(1-\eta)^2}, \quad (4.6)
$$

avec $\eta = \frac{\pi}{6}\sigma^3 (\rho_1 + \sum_{i=2}^{n} \rho_m(m_i)m_i) = \frac{\pi}{6}\sigma^3 (\rho_1 + \langle \rho_m \rangle_1)$. La partie résiduelle $\mu^{res}$ du potentiel chimique de la chaîne de longueur $m$ est donnée par

$$
\mu^{\text{res}}(m) = \left( \frac{\partial f^{\text{res}}}{\partial \langle \rho_m \rangle_0} \right)_{\rho_1, \langle \rho_m \rangle_1} \cdot \left( \frac{\partial \langle \rho_m \rangle_0}{\partial \rho_m(m)} \right)_{\rho_1} \\
+ \left( \frac{\partial f^{\text{res}}}{\partial \langle \rho_m \rangle_1} \right)_{\rho_1, \langle \rho_m \rangle_0} \cdot \left( \frac{\partial \langle \rho_m \rangle_1}{\partial \rho_m(m)} \right)_{\rho_1}. \qquad (4.7)
$$

Sachant que $(\frac{\partial \langle \rho_m \rangle_k}{\partial \rho_m(m)})_{\rho_1} = m^k$, on peut montrer que la partie résiduelle $\mu^{res}$ est linéaire par rapport à $m$. On peut donc écrire $\mu^{res}$ selon

$$
\mu^{res}(m) = \langle \mu \rangle_0^{res} + m \langle \mu \rangle_1^{res}. \quad (4.8)
$$

avec $\langle \mu \rangle_k^{\text{res}} = \left( \frac{\partial f^{\text{res}}}{\partial \langle \rho_m \rangle_k} \right)_{\rho_1, \langle \rho_m \rangle_{j \neq k}}$. Pour le calcul d'un point de bulle, la distribution des longueurs de chaîne dans une des deux phases à l'équilibre (phase cl "cloud") est égale à la distribution globale, donc est connue. La distribution dans l'autre phase (phase sh "shadow") est obtenue en égalant le
potentiel chimique de la chaîne de longueur m dans chacune des phases, c'est-à-dire
$\mu^{(Cl)}(m) = \mu^{(Sh)}(m)$. On en déduit que

$$
\rho_m^{(Sh)}(m) = \rho_m^{(Cl)}(m) \exp \left[ \sum_{k=0,1} \left( \langle \mu \rangle_k^{res(Cl)} - \langle \mu \rangle_k^{res(Sh)} \right) \frac{m^k}{kT} \right] . \quad (4.9)
$$

Les Eqs. (4.9) et (4.2) permettent de déterminer la fonction de distribution des chaînes dans la phase Sh. Pour déterminer les moments de densité $\langle\rho_m\rangle_0$, et $\langle\rho_m\rangle_1$ et les masses volumiques des deux phases en équilibre, on doit résoudre plusieurs conditions : l'égalité des pressions dans les deux phases, et deux autres équations correspondant à la définition de $\langle\rho_m\rangle_0$, et $\langle\rho_m\rangle_1$ : il suffit pour cela de considérer l'Eq. (4.9) et de l'intégrer. On obtient ainsi les deux conditions suivantes ¹¹⁰ :

$$
\langle \rho_m \rangle_0^{(Sh)} = \int dm \rho_m^{(Cl)}(m) \exp \left[ \sum_{k=0,1} \left( \langle \mu \rangle_k^{\text{res(Cl)}} - \langle \mu \rangle_k^{\text{res(Sh)}} \right) \frac{m^k}{kT} \right], \quad (4.10)
$$

$$
\langle \rho_m \rangle_1^{(Sh)} = \int dm \rho_m^{(Cl)}(m) m \exp \left[ \sum_{k=0,1} \left( \langle \mu \rangle_k^{\text{res(Cl)}} - \langle \mu \rangle_k^{\text{res(Sh)}} \right) \frac{m^k}{kT} \right] \quad (4.11)
$$

Si la distribution globale est une distribution de Schulz-Flory, alors la distribution dans la phase Sh est également une distribution de Schulz-Flory, d'après l'Eq. (4.9). Cependant, les moments de cette distribution sont changés (par exemple la longueur de chaîne moyenne sera modifiée).

<figure><img src="image_37.png" /><figcaption>Figure 4-3. Equilibre liquide-liquide dans un système solvant + polymère polydispersé, prédit avec la théorie SAFT-HS et la méthode des moments (thermodynamique continue), pour des distributions de masse molaire de</figcaption></figure>
Schulz-Flory ayant une longueur moyenne de chaîne $\langle m \rangle_n^{(0)} = 100$ et des indices de polydispersité différents. (a) Distribution des chaînes moléculaires étudiées. (b) Courbes cloud et shadow prédites pour différents indices de polydispersité. (c) Evolution de la longueur de chaîne moyenne dans la phase Sh, en fonction de la composition totale en polymère. Les cercles représentent les points critiques de mélange. Le cas $I_p = 1$ correspond à un mélange binaire solvant + chaîne de longueur 100. D'après Paricaud et al.¹¹⁰.

On peut ainsi changer la masse molaire moyenne d'un polymère en provoquant un équilibre de phase. La Figure 4-3 montre les résultats obtenus pour une mélange solvant + polymère polydispersé : les courbes cloud (Cl) et shadow (Sh) (équilibres liquide-liquide) pour des distributions de Schulz-Flory ayant la même longueur de chaîne moyenne, mais des indices de polydispersité différents. On observe clairement que la polydispersité des masses moléculaires élargit la zone d'immiscibilité liquide-liquide (Figure 4-3 (b)). D'autre part, la longueur de chaîne moyenne dans la phase shadow peut être ou inférieure ou supérieure à la longueur de chaîne moyenne globale, selon la composition en solvant.

<figure><img src="image_38.png" /></figure>

Figure 4-4. Equilibre liquide-liquide dans un système particules colloïdales monodispersées + polymère polydispersé, prédit avec la théorie SAFT-HS et la méthode des moments (thermodynamique continue), pour des distributions de masse molaire de Schulz-Flory ayant une longueur moyenne de chaîne (nombre de segments) $\langle m \rangle_n^{(0)} = 50$ et des indices de polydispersité différents. (a) Courbes cloud et shadow prédites pour différents indices de polydispersité. (b) Evolution de la longueur de chaîne moyenne dans la phase Sh, en fonction de la
composition totale en polymère. Les symboles représentent les points critiques des mélanges. Le cas $I_p = 1$ correspond à un mélange binaire solvant + chaîne de longueur. D'après Paricaud et al. ¹¹¹.

<figure><img src="image_39.png" /></figure>

Figure 4-5. Equilibre liquide-liquide dans un système de particules colloïdales monodispersées (composé 1) mélangées avec un autre type de particules colloïdales polydispersées en taille (diamètre). (a) Courbes cloud et shadow prédites pour différents indices de polydispersité avec le modèle développé par Paricaud ²⁵ pour les sphères dures non-additives. (b) Evolution du diamètre moyen des particules polydispersées dans la phase Sh, en fonction de la composition totale en composé 1. (c). Fonction de distribution décrivant la composition globale des particules polydispersées.

On peut ainsi créer un polymère de longueur de chaîne supérieure ou inférieure à celle issue d'une polymérisation, simplement en réalisant une équilibre liquide-liquide et en jouant sur la composition en solvant et les conditions de température et de pression. Dans ce cas particulier où l'indice de polydispersité de la distribution globale de Schulz-Flory est inchangé par l'équilibre de phase et seule la longueur moyenne change. Nous avons obtenu des résultats similaires pour les mélanges de polymères polydisersés + colloïdes monodispersés ¹¹¹ (Figure 4-4). Ce n'est cependant pas le cas pour un système de particules colloïdales polydispersées selon leur taille (diamètre) ²⁵, car dans ce cas, l'indice de polydispersité varie beaucoup d'une phase à l'autre (Figure 4-5).
# 5. Propriétés Thermodynamiques des Solutions d'Electrolytes

Les solutions d'électrolytes sont rencontrées dans de nombreux procédés industriels (capture des gaz acides, fabrication des engrais, traitements des minerais, extractions liquide-liquide dans l'industrie nucléaire, ...), et la connaissance de leurs propriétés thermodynamiques est très importante pour le dimensionnement des unités opératoires (réacteurs, cristalliseurs, extracteurs, ...). Les solutions d'électrolytes ne peuvent pas être traitées directement par les équations d'état classiques pour les mélanges non électrolytiques, car elles présentent des interactions très spécifiques entre espèces chimiques (ions). Les ions en solutions interagissent par des forces coulombiennes de grande portée dont l'amplitude dépend de la constante diélectrique du solvant. Les interactions ion-solvant, dites interactions de solvatation, sont également très spécifiques à ces systèmes par leur amplitude et leur courte portée : les ions fortement solvatés sont en permanence entourés d'une ou de plusieurs couches de molécules de solvant (couches de solvatation). Une description adéquate de ces forces nécessite l'utilisation de modèles ou de termes spécifiques aux solutions d'électrolytes. Lors d'une journée SFGP (société française de génie des procédés) dédiée aux solutions d'électrolytes, Fürst ¹¹² et Paricaud ¹¹³ ont fait une revue détaillée des approches thermodynamiques utilisées pour décrire les propriétés de ces systèmes. D'excellentes revues des modèles utilisés pour les solutions d'électrolytes peuvent être trouvées dans la littérature ⁸,¹¹⁴,¹¹⁵. Nous allons présenter ici uniquement le cas des équations d'état de type SAFT, et exposer nos travaux dans ce domaine.

La plupart des modèles et termes d'énergie libre pour les solutions d'électrolytes supposent que le solvant est un milieu continu ayant une certaine constante diélectrique $D_m$, et que les ions sont des sphères chargées. Le modèle de Debye-Hückel ¹¹⁶⁻¹¹⁸ a été obtenue à la suite de la résolution de l'équation de Poisson-Boltzman permettant de déterminer le potentiel électrostatique moyen. Le modèle DH n'est valable qu'à très faibles concentrations des ions (faible force ionique). Le modèle MSA (Mean Spherical Approximation), issu de la résolution analytique des équations intégrales d'Ornstein-Zernike, est applicable sur un large domaine de concentration. Il existe plusieurs niveaux de représentation d'une solution d'électrolytes (voir la Figure 5-1) via un terme MSA. Par exemple, le modèle primitif restreint (RPM) ¹¹⁹ suppose que le solvant est un milieu diélectrique continu, et que les ions sont des sphères dures charges ayant toutes le même diamètre. Le modèle primitif non restreint (PM) ¹²⁰ est similaire, mais autorise le fait que les ions puissent avoir des diamètres différents. Les modèles non primitifs ¹²¹ considèrent les molécules de solvant de manières explicites, celles-ci étant modélisées par des sphères dipolaires. Enfin, le modèle BiMSA ¹²²,¹²³ considère des associations site-site (terme similaire à la théorie de Wertheim) entre ions, et permet de tenir compte de la formation de paires d'ions, lesquelles sont observées pour les ions larges, ou à hautes températures et à fortes concentrations. Des termes MSA ont été développés pour chacun de ces modèles (Figure 5-1).

<figure><img src="image_40.png" /><figcaption>Figure 5-1. Classification des modèles utilisés pour représenter les solutions d'électrolytes. RPM : ions de même taille, diélectrique continu; PM : ions de tailles différentes, milieu continu; NPM : ions de tailles différentes, solvant explicite (sphère dipolaire); BiMSA : ions de tailles différentes, diélectrique continu, formation de paires d'ions. Source : Paricaud ¹¹³.</figcaption></figure>
Harvey et Prausnitz ¹²⁴, Jin et Donohue ¹²⁵, Aasberg-Petersen et al. ¹²⁶ puis Fürst et Renon ¹²⁷, ont été parmi les premiers à proposer des équations d'état pour solutions d'électrolytes. Le modèle de Fürst et Renon ¹²⁷ a été utilisé sur des systèmes industriels complexes ¹²⁸⁻¹³⁴. La particularité de ce modèle est l'utilisation de corrélations dépendant des diamètres ioniques expérimentaux, permettant de réduire le nombre de paramètres et de rendre le modèle prédictif. L'avantage des équations d'état pour électrolytes est que ces modèles sont aussi applicables pour les mélanges non-électrolytiques sur l'ensemble du domaine de température et de pression, tandis que l'utilisation d'un modèle pour électrolytes à coefficients d'activité nécessite la combinaison de plusieurs modèles thermodynamiques et corrélations (pressions de vapeur des corps pures, modèles à énergies de Gibbs d'excès, loi de Henri, etc.), ce qui entraine souvent l'utilisation de nombreux paramètres empiriques ajustés et non extrapolables.

L'énergie libre de Helmholtz d'une équation d'état pour solutions d'électrolytes s'écrit de la forme

$$A = A^{\text{eos-ne}} + A^{\text{ion}} + A^{\text{solv}}, \quad (5.1)$$

où $A^{\text{eos-ne}}$ est l'énergie libre de l'équation d'état pour les mélanges non-électrolytiques (équations cubiques, CPA, SAFT, ...), $A^{\text{ion}}$ le terme électrostatique décrivant les interactions colombiennes entre les ions (terme DH ou MSA), et $A^{\text{solv}}$ le terme de solvatation décrivant des interactions spécifiques entre les ions et les molécules de solvant. Certains modèles comme SAFT-VRE ¹³⁵ n'incluent pas de terme de solvatation, celle-ci étant traité dans le terme $A^{\text{eos-ne}}$. Les différents modèles pour électrolytes se différencient non seulement par rapport aux modèles choisis pour le terme MSA, le terme de solvatation et la partie non électrolytique de l'équation d'état, mais aussi par les différentes hypothèses à faire sur les diamètres des ions (diamètres solvatés ou diamètres cristallins) ainsi que sur la constante diélectrique utilisées dans le terme MSA (constante diélectrique du solvant pure, ou prise en compte de l'effet des ions, ...).

## 5.1. Modélisation avec le modèle SAFT-VRE

Une partie de mes travaux de recherche sur les solutions d'électrolytes ont été réalisés avec l'équation d'état SAFT-VRE développée par Galindo et al. ¹³⁵, qui est l'extension de SAFT-VR aux solutions d'électrolytes. Dans l'approche SAFT-VR, les molécules sont représentées par des chaînes flexibles de segments sphériques, et les interactions entre espèces chimiques sont décrites par des potentiels à puits carrés (square-well). L'eau est représentée par une sphère interagissant via un potentiel à puits carré et des sites associatifs. Les paramètres de l'eau ont été optimisés sur les données liquide-vapeur de l'eau pure ⁷⁸,¹³⁶,¹³⁷. Un potentiel à puits carré est également utilisé pour décrire les interactions de solvatation.

L'équation d'état SAFT-VRE s'exprime selon la somme de différentes contributions de l'énergie libre :

$$\frac{A}{NkT} = \frac{A^{\text{ideal}}}{NkT} + \frac{A^{\text{seg}}}{NkT} + \frac{A^{\text{chain}}}{NkT} + \frac{A^{\text{assoc}}}{NkT} + \frac{A^{\text{MSA}}}{NkT} . \quad (5.2)$$

où $A^{\text{ideal}}$ l'énergie libre d'un mélange de gaz idéaux, $A^{\text{seg}}$ la contribution due aux interactions entre segments, $A^{\text{chain}}$ la contribution due à la formation des chaînes, $A^{\text{chain}}$ la contribution associative due aux liaisons hydrogène. Le terme $A^{\text{MSA}}$ (Mean Spherical Approximation) le terme développé par Waisman et Lebowitz¹¹⁹ pour le modèle restreint primitif (RPM), et la constante diélectrique utilisée est celle de l'eau pure liquide saturante. Le potentiel chimique $\mu_i$ du constituant $i$ est donné par $\mu_i = (\partial A / \partial n_i)_{T,V,n_j, j \neq i}$, et son coefficient d'activité est donné par

$$\gamma_i = \frac{\rho}{\rho_{\text{ref}}} \exp \left( \frac{\mu_i^{\text{res}} - \mu_i^{\text{res(ref)}}}{RT} \right) \quad (5.3)$$

où $\rho$ et $\rho_{\text{ref}}$ sont la densité de la solution à T, P et celle à l'état de référence (eau pure liquide à T, P), respectivement, et $\mu_i^{\text{res}}$, $\mu_i^{\text{res(ref)}}$ les potentiels chimiques résiduels molaires correspondant. Le coefficient d'activité dans l'échelle des molalités est donné par¹¹⁵ $\gamma_i^{(m)} = \gamma_i x_w$, où $x_w$ est la fraction
molaire de l'eau. Le coefficient d'activité moyen pour un électrolyte 1-1 est donné
$\gamma^{+−} = \sqrt{(\gamma_C^{(m)} \gamma_A^{(m)})}$, le coefficient osmotique est donné par$^{115}$ $\phi = -1000 \ln (x_w \gamma_w) / (2m_s M_w)$.

**5.1.1. Systèmes eau-sel-hydrocarbures**

Nous avons appliqué ce modèle aux mélanges ternaire eau + alcane + sel, et prédit les effets de sel (salting out) sur les équilibres liquide-vapeur et liquide-liquide. Ces mélanges sont d'une grande importance dans l'industrie pétrolière. Les modèles moléculaires SAFT sont décrits dans la Figure 5-2. L'alcane est représenté par une chaîne de sphères, tandis que les ions sont des sphères chargées. L'eau possède quatre sites associatifs.

<figure><img src="image_41.png" /><figcaption>Figure 5-2. Modèles moléculaires SAFT utilisés pour l'eau (a), les alcanes (b) et les ions (c) $^{138}$.</figcaption></figure>

<figure><img src="image_42.png" /><figcaption>Figure 5-3. Effet de la présence de sel sur la solubilité du méthane dans l'eau à T =443.15K. Comparaisons entre les points expérimentaux et la prédiction du modèle SAFT-VRE <sup>138</sup>. Les cercles correspondent à une composition nulle en sel (mélange binaire eau-méthane) et les carrés correspondent à une molalité en NaCl de 5.2mol/kg.</figcaption></figure>
Une application est par exemple la prédiction des conditions de formations des hydrates de méthane en présence de sels. Nous avons prédit la solubilité du méthane dans l'eau à différentes températures et en présence de sel (NaCl): on voit que la présence de sel diminue sensiblement la solubilité du méthane, pour une pression et une température données (Figure 5-3).

### 5.1.2. Modélisation des équilibres de phases liés au procédé SI.

Le procédé par cycle thermochimique soufre-iode (SI) (Figure 5-4) est considéré comme l'un des possibles procédés de production de l'hydrogène à l'échelle industrielle à partir d'eau et de chaleur. Ce cycle est une série de réactions chimiques et d'étapes de séparation des produits dont la somme donne la décomposition de l'eau en hydrogène et oxygène. Il comprend trois réactions chimiques dont la décomposition de l'acide HI en H₂ et I₂ (partie HIx du procédé) et la réaction de Bunsen (I₂ + SO₂ + 2H₂O → 2HI + H₂SO₄) ayant lieu autour de 100°C et conduisant à une séparation liquide-liquide entre deux phases acides. Il est prévu de combiner le procédé iode-soufre avec un réacteur nucléaire à très haute température de type (réacteur VHTR).

<figure><img src="image_43.png" /><figcaption>Figure 5-4. Schéma du procédé à cycle soufre-iode. Zoom sur la partie HIx du procédé. Modélisation des équilibres liquide-vapeur du mélange HI+H₂O avec l'équation d'état SAFT-VRE ¹³⁷.</figcaption></figure>

Afin de dimensionner et d'optimiser le procédé SI, il est nécessaire de développer un modèle thermodynamique capable de représenter à la fois les équilibres de phases et les enthalpies des mélanges. Dans le cadre d'un contrat avec le CEA, nous avons modélisé avec l'équation d'état SAFT-VRE les équilibres de phases liés à la partie HIx du procédé de production de l'hydrogène par cycle thermochimique iode-soufre, en particulier les propriétés des mélanges HI + H₂O, I₂ + H₂O et HI + I₂. Ce travail a été fait en partie par Laetitia Tazi, stagiaire sous ma responsabilité, et a conduit à la publication d'un article scientifique ¹³⁷. La modélisation du mélange HI + H₂O à haute température est particulière car elle nécessite la prise en compte de l'équilibre chimique 2HI ⇌ H⁺ + I⁻. Cet équilibre est résolu en phase liquide, en plus de l'équilibre liquide-vapeur.

### 5.1.3. Modélisation des semi-clathrates de gaz

Les semi-clathrates sont des hydrates constitués d'un réseau cristallin de molécules d'eau et d'ions. Ils se comportent comme des complexes solides car leur formation est stœchiométrique (composition
fixe). L'anion est un halogénure (Cl⁻, Br⁻, F⁻) ou un nitrate (NO₃⁻), tandis que le cation est un alkyl ammonium quaternaire. Ces hydrates sont stables à pression ambiante et à des températures proches de l'ambiant. Leur température de dissociation dépend de la nature des ions. Ces hydrates présentent des cavités formées par les molécules d'eau, qui peuvent éventuellement encapsuler des molécules de gaz (Figure 5-5). Contrairement aux hydrates de gaz classiques, les semi-clathrates sont stables même sans la présence de gaz. Les semi-clathrates de bromure de tétrabutylammonium (TBAB) ont été très étudiés car le TBAB est très utilisé dans l'industrie et donc assez bon marché par rapport aux autres sels d'alkymmonium.

<figure><img src="image_44.png" /><figcaption>Figure 5-5. Structure d'un semi-clathrate hydrate de TBAB de type B, contenant deux molécules de gaz (sphères jaunes). D'après Shimada *et al.*¹³⁹.</figcaption></figure>

Le diagramme de phase du mélange TBAB+H₂O system (Figure 5-6) présente différentes phases hydrates ayant différentes compositions: l'hydrate de type A a un nombre d'hydratation égal à 26 (26 molécules d'eau pour une molécule de TBAB) et l'hydrate de type B a un nombre d'hydratation égal à 38.

<figure><img src="image_45.png" /><figcaption>Figure 5-6 Diagramme de phase solide-liquide du mélange binaires eau + TBAB. Les symboles représentent les points expérimentaux (voir les références¹⁴⁰,¹⁴¹ pour une revue) : astérisques: Ref. ¹⁴²; losanges vides: type A, Ref. ¹⁴³; cercles: type B, Ref. ¹⁴³; + Ref. ¹⁴²; losanges pleins : Ref. ¹⁴¹; triangles vides: Ref. ¹⁴⁴.Les lignes épaisses sont les courbes de coexistence solide-liquide calculées avec SAFT-VRE et l'équation (7).</figcaption></figure>

En présence de molécules de gaz de petite taille et sous pression (H₂, CO₂, N₂, CH₄,...), les solutions de TBAB en cristallisant peuvent former des semi-clathrates de gaz. Les hydrates de TBAB peuvent ainsi servir de matériaux absorbant les gaz et permettant la séparation de gaz par différence d'absorption. Les applications potentielles de ces systèmes sont l'absorption du CO₂ des fumées de
combustion, et la séparation des composés du gaz de synthèse (H₂, CO₂, et CO). De nombreuses mesures expérimentales des conditions de formation des hydrates d'alkylammonium¹⁴⁰,¹⁴¹,¹⁴⁵⁻¹⁵⁶ ont été récemment réalisées. Nous avons été les premiers à proposer un tel modèle pour déterminer les conditions de dissociation des semi-clathrates de gaz¹⁵⁷. Ce modèle permet de prédire la température et l'enthalpie de dissociation des hydrates pour une pression et une composition en sel de TBAB données, et sera particulièrement utile pour le dimensionnement et l'optimisation d'unité de séparation de gaz. Herri et col.¹⁵⁸,¹⁵⁹ a développé une approche similaire basée sur le modèle e-NRTL et a obtenu de bons résultats pour les semi-clathrates de méthane.

Nous rappelons ici les principales équations du modèle thermodynamique développé, puis nous présentons les résultats de modélisation appliqués aux hydrates de TBAB + CO₂.

## Equilibre solide-liquide entre la phase aqueuse et un hydrate de TBAB sans présence de gaz

On considère un système contenant une solution de TBAB de composition massique $w_{TBAB}$ (fraction massique en TBAB), à la température $T$ et à la pression $P$, en équilibre avec un hydrate semi-clathrate de TBAB, ayant un nombre d'hydratation $\nu_w$ défini comme étant le nombre de molécule d'eau par molécule de TBAB dans la phase hydrate. On peut montrer que la condition d'équilibre est donnée par¹⁵⁷

$$ \Delta g_{dis} = \nu_w \mu_w^{liq} + \mu_C^{liq} + \mu_A^{liq} - \nu_w \mu_w^{hyd} - \mu_C^{hyd} - \mu_A^{hyd} = 0. \quad (5.4) $$

où $\mu_w^{ref(x)}$, $\mu_C^{ref(x)}$, $\mu_A^{ref(x)}$, sont les potentiels chimiques de l'eau, du cation (TBA⁺) et de l'anion (Br⁻). Les indices « liq » et « hyd » correspondent aux phases liquide et hydrate. L'équation (5.4) est également valide pour les systèmes TBAB + eau + gaz. Les potentiels chimiques $\mu_i^{liq}$ en phase liquide peuvent s'exprimer selon

$$ \mu_i^{liq} = \mu_i^{ref(x)} + RT \ln(x_i \gamma_i). \quad (5.5) $$

où $\mu_i^{ref(x)}$ est le potentiel de référence dans l'échelle des fractions molaires, $\gamma_i$ le coefficient d'activité de l'espèce i dans l'échelle des fractions molaires, et $x_i$ la fraction molaire du constituant i. L'état de référence est la dilution infinie dans l'eau liquide pure à T et P. La condition d'équilibre (Eq. (5.4)) s'écrit alors

$$ \frac{\Delta g^0}{RT} + \ln(x_C \gamma_C) + \ln(x_A \gamma_A) + \nu_w \ln(x_w \gamma_w). \quad (5.6) $$

avec $\Delta g^0 = \nu_w \mu_w^* + \mu_C^{ref(x)} + \mu_A^{ref(x)} - \mu_{hyd}$ dépendant uniquement de T et P. En faisant l'hypothèse que les capacités calorifiques des phases liquides et hydrates ne dépendent pas de la température, on peut exprimer $\Delta g^0$ selon

$$ \frac{\Delta g^0}{RT} = \frac{\Delta h^0}{RT} \left(1 - \frac{T}{T_0}\right) - \frac{\Delta c_p^0}{R} \left(\frac{T_0 - T}{T}\right) + \frac{\Delta c_p^0}{R} \ln\left(\frac{T_0}{T}\right) + \frac{\Delta v^0}{RT} (P - P_0) + \frac{\Delta g^0 (T_0, P_0)}{RT} \quad (5.7) $$

où $\Delta h^0$, $\Delta c_p^0$ et $\Delta v^0$ sont des paramètres caractérisant la phase hydrate. $\Delta h^0$ peut être estimé par l'enthalpie de fusion de la phase hydrate en son point de congruence ou directement ajusté sur la courbe d'équilibre solide-liquide, tandis que l'on peut négliger les paramètres $\Delta c_p^0$ et $\Delta v^0$. Soit $T_0$ le point de fusion de l'hydrate à sa composition stœchiométrique à la pression atmosphérique $P_0=0.101325MPa$; $\Delta g^0(T_0, P_0)$ est la valeur de $\Delta g^0$ à $T_0$ et $P_0$. En appliquant l'Eq. (5.6) à la composition stœchiométrique $(x_w^{(st)}, x_C^{(st)}, x_A^{(st)})$ et à $T_0$ et $P_0$, on peut déterminer $\Delta g^0(T_0, P_0)$ selon

$$ \frac{\Delta g^0 (T_0, P_0)}{RT} = -\ln(x_C^{(st)} \gamma_C^{(st)}) - \ln(x_A^{(st)} \gamma_A^{(st)}) - \nu_w \ln(x_w^{(st)} \gamma_w^{(st)}). \quad (5.8) $$

où $x_w^{(st)}/\nu_w = x_C^{(st)} = x_A^{(st)} = 1/(\nu_w + 2)$. Le nombre d'hydratation $\nu_w$ est fixé à sa valeur expérimentale tandis que $T_0$, $\Delta h^0$, $\Delta c_p^0$ et $\Delta v^0$ sont ajustées sur les courbes de coexistence solide-liquide et les enthalpies de dissociation. Notez que $\Delta h^0$ a des valeurs très proches de l'enthalpie de
dissociation $\Delta H_{dis}$, donnée par $\Delta H_{dis} = -RT^2\partial(\Delta g_{dis}/(RT))/\partial T$. Les paramètres du TBAB ont été optimisés pour représenter les données expérimentales de coefficients osmotiques. L'approche thermodynamique développée ci-dessus est utilisée pour décrire les équilibres solide-liquide du système eau + TBAB. Les paramètres caractérisant les phases hydrates de types A (ou H26) et B (ou H38) sont présentés dans le Tableau 5-1. Comme le montre la Figure 5-6, le nouveau modèle décrit très bien le diagramme de phase.

<table><caption>Tableau 5-1 Paramètres des phases semi-clathrates de TBAB (Ref. ¹⁵⁷): $v_w$ est le nombre d'hydratation, $T_0$ la température de fusion de la phase hydrate à la composition stœchiométrique, $\Delta h^0$ le paramètre enthalpique exprimé en kJ par mole de TBAB; $\Delta v^0$ est exprimé en cm³ par mole de TBAB; $n_i$ est le nombre de cavité par molécule de TBAB dans l'hydrate. Les paramètres de potentiel de cavité pour le CO2, $\epsilon_{ij}^{cell}$ et $V_{ij}^{cell}$, sont également reportés.</caption><thead><tr><th>Type</th><th>v<sub>w</sub></th><th>T<sub>0</sub> (K)</th><th>Δh<sup>0</sup> (kJ.mol<sup>-1</sup>)</th><th>Δv<sup>0</sup> (cm<sup>3</sup>.mol<sup>-1</sup>)</th><th>n<sub>i</sub></th><th>V<sub>ij</sub><sup>cell</sup> (Å<sup>3</sup>)</th><th>ε<sub>ij</sub><sup>cell</sup>/k (K)</th></tr></thead><tbody><tr><td>B</td><td>38</td><td>283.50</td><td>200.7</td><td>-30</td><td>3</td><td>0.012415</td><td>2750</td></tr><tr><td>A</td><td>26</td><td>285.15</td><td>150.2</td><td>-30</td><td>2</td><td>0.006208</td><td>2750</td></tr></tbody></table>

## Equilibre solide-liquide entre la phase aqueuse et un semi-clathrate hydrate de TBAB en présence de gaz

On considère maintenant un système contenant une solution de TBAB de composition initiale massique $w_{TBAB}$ (fraction massique en TBAB avant la mise sous pression), à la température $T$ et à la pression $P$, en équilibre avec un hydrate semi-clathrate de TBAB ayant un nombre d'hydratation $v_w$, en présence de $N_g$ différents gaz (H₂, CO₂, N₂, ...). Ce système est en équilibre triphasique (solide-liquide-gaz). L'approche que nous avons proposée est une combinaison du modèle décrit dans la section précédente avec le modèle de van der Waals-Plattteeuw (vdW-P) model¹⁶⁰⁻¹⁶² utilisés pour modéliser les hydrates de gaz. L'hypothèse principale du nouveau modèle est de supposer que la phase métastable vide $\beta$ du modèle vdW-P correspond aux phases hydrates de TBAB qui se forment sans présence de gaz. Il est alors possible de montrer que la condition d'équilibre s'écrit¹⁵⁷

$$ 
\begin{aligned}
& - \frac{\Delta h^0}{RT} \left(1 - \frac{T}{T_0}\right) + \frac{\Delta c_p^0}{R} \left(\frac{T_0 - T}{T}\right) - \frac{\Delta c_p^0}{R} \ln \left(\frac{T_0}{T}\right) - \frac{\Delta v^0}{RT} (P - P_0) - \frac{\Delta g^0 (T_0, P_0)}{RT} \\
& = \ln (x_C \gamma_C) + \ln (x_A \gamma_A) + v_w \ln (x_w \gamma_w) + \sum_{i=1}^{N_{cav}} n_i \ln \left( 1 + \sum_{j=1}^{N_g} C_{ij} f_j \right),
\end{aligned}
\quad (5.9) $$

Dans l'Eq. (5.9), $N_{cav}$ est le nombre de types de cavité différents, $n_i$ est le nombre de cavités de type $i$, $f_i$ la fugacité du gaz $j$, et $C_{ij}$ la constante de Langmuir donné par

$$ C_{ij} = \frac{4\pi}{kT} \int_{0}^{R_{cav,i}} \exp\left(\frac{-w_{ij}(r)}{kT}\right) r^2 dr. \quad (5.10) $$

Dans l'Eq. (5.10), $w_{ij}(r)$ est le potentiel d'interaction entre la molécule de gaz $j$ et les molécules constituant la cavité (molécules d'eau + influence des ions). En supposant que le potentiel $w_{ij}(r)$ a la forme d'un puits carrés, on peut exprimer $C_{ij}$ selon

$$ C_{ij} = \frac{4\pi}{kT} V_{ij}^{cell} \exp \left( \frac{\epsilon_{ij}^{cell}}{kT} \right), \quad (5.11) $$

où $V_{ij}^{cell}$ est le volume libre de la molécule de gaz $j$ à l'intérieur de la cavité $i$, et $\epsilon_{ij}^{cell}$ est l'énergie d'attraction (profondeur du puits de potentiel de cavité) entre la cavité $i$ et la molécule de gaz $j$. Les paramètres $\epsilon_{ij}^{cell}$ et $V_{ij}^{cell}$ sont ajustés sur les données expérimentales (température de dissociation de la phase hydrate à pression et composition en TBAB données).
<figure><img src="image_46.png" /></figure>

Figure 5-7 Equilibres de phase liquide-vapeur-hydrate (L-V-H) pour le système CO2 + TBAB + H2O. Les symboles représentent les données expérimentales et les lignes continuent sont les lignes triphasiques prédites par le nouveau modèle, pour $w_{TBAB}$ = 0.03, 0.045, 0.05, 0.1 et 0.427. Triangles vides : $w_{TBAB}$ = 0.03, Ref. ¹⁵⁰; cercles pleins: $w_{TBAB}$ = 0.045, Ref. ¹⁵⁰; astérisques : $w_{TBAB}$ = 0.05, Ref. ¹⁶³; triangles pleins, $w_{TBAB}$ = 0.1, Ref. ¹⁶³; cercles vides, $w_{TBAB}$ = 0.0442, Ref. ¹⁴⁹; losanges vides : $w_{TBAB}$ = 0.07, Ref. ¹⁴⁹; carrés vides: $w_{TBAB}$ = 0.0901, Ref. ¹⁴⁹; losanges pleins : $w_{TBAB}$ = 0.1, Ref. ¹⁴⁰; carrés pleins: $w_{TBAB}$ = 0.427, Ref. ¹⁴⁰; croix : $w_{TBAB}$ = 0.4, Ref. ¹⁴¹. La ligne continue fine correspond à $w_{TBAB}$ = 0.32 (composition stœchiométrique de l'hydrate de type B) et la ligne en pointillées correspond à $w_{TBAB}$ = 0.55. La ligne en traits fins correspond à $w_{TBAB}$ = 0.6.

<figure><img src="image_47.png" /></figure>

Figure 5-8 Diagramme température-composition des systèmes TBAB + H2O et TBAB + H2O + CO2. Les lignes sont les prédictions du nouveau modèle: les lignes continues correspondent au diagramme de phase T-x du mélange TBAB + H2O à pression atmosphérique, et les lignes en traits mixtes fins correspond à la ligne triphasique L-V-H38 à P=2MPa. Le cercle gris est le point quadruple prédit par le modèle, où coexistent 4 phases : liquide, vapeur et les hydrates de types A et B. Les symboles noirs sont les points expérimentaux pour $w_{TBAB}$ = 0.045 (Ref. ¹⁵⁰) et $w_{TBAB}$ = 0.1 (Ref. ¹⁶³). Les croix sont les données expérimentales de Deschamps et Dalmazzone¹⁴¹.

Selon Shimada et al.¹³⁹, la structure de l'hydrate de TBAB de type B est composed d'un cation TBA⁺, d'un anion Br⁻ et de anions, et de 38 molécules d'eau. Elle a trois cavités de forme dodécaédrique
par molécule de TBAB, pouvant contenir des molécules de gaz. En conséquence, le parameter $n_i$ vaut 3 pour l'hydrate de type B. On fait l'hypothèse que l'hydrate de type A n'a que deux cavités par molécule de TBAB. Comme le montre la Figure 5-7, notre modèle peut décrire de façon très précise les pressions de dissociation à température et composition en TBAB données. En particulier, l'effet de la composition en TBAB est parfaitement reproduit, sans aucun ajustement de paramètre. Comme le montre la Figure 5-8, le diagramme température-composition à la pression P = 2MPa est très bien prédit par le nouveau modèle, sans ajustement de paramètre. Le modèle prédit un changement de structure du type A au type B à haute pression et pour de fortes concentrations en TBAB. D'autre part, le modèle prédit une déstabilisation de l'hydrate lorsque que la fraction massique en TBAB dépasse 0.32, c'est-à-dire la composition stœchiométrique de l'hydrate de type B.

## 5.2. Développement d'une version électrolyte de PC-SAFT.

Au cours de sa thèse dirigée par Walter Fürst, à laquelle j'ai participé en tant que co-encadrant, Thi Thu Suong Ho ¹⁶⁴ a développée une nouvelle équation d'état pour les solutions d'électrolytes. La nouvelle équation d'état (iPC-SAFT) combine le terme MSA, un terme de solvation, et une version polaire du modèle PC-SAFT EOS (terme polaire de Jog et Chapman ⁴²,⁴³). Le terme MSA utilisé correspond au modèle non primitif pour les interactions coulombiennes entre ions. L'expression de Pottel ¹⁶⁵ a été utilisée pour calculer la constante diélectrique $D_m$ du milieu : $D_m = (D_s - 1) \frac{1-\xi_{3,ions}}{1+\xi_{3,ions}/2}$, où $D_s$ est la constante diélectrique du solvant pure, et $\xi_{3,ions}$ est la compacité les ions dans le milieu. Cette expression tient compte du fait que la constante diélectrique diminue lorsque la concentration en ion augmente. La prise en compte de l'effet des sels sur la constante diélectrique du milieu est originale car rarement utilisée dans la littérature.

<figure><img src="image_48.png" /></figure>

Figure 5-9 Coefficients d'activité moyens ($\gamma_{\pm}$) de cinq solutions aqueuses contenant un sel de bromure, à T=298.15K et à pression atmosphérique, en fonction de la molalité du sel. Les symboles correspondent aux données expérimentales ¹¹⁵,¹⁶⁶, tandis que les lignes ont été calculées avec le modèle iPC-SAFT-SR2.
<figure><img src="image_49.png" /></figure>

Figure 5-10 Masses volumiques de cinq solutions aqueuses contenant un sel de bromure, à T=298.15K et pression atmosphérique, en fonction de la molalité du sel. Les symboles correspondent aux données expérimentales<sup>115</sup>, tandis que les lignes ont été calculées avec le modèle iPC-SAFT-SR2.

<figure><img src="image_50.png" /></figure>

Figure 5-11 Coefficients d'activités moyens de solutions aqueuses de KCl, à différentes températures, en fonction de la molalité en KCl. Comparaison entre les données expérimentales<sup>166</sup> et le modèle iPC-SAFT-SR2.

Deux versions du modèle iPC-SAFT ont été développées, correspondant à deux différents termes de solvatation. Le premier modèle, iPC-SAFT-SR2, utilise le terme SR2 comme terme de solvatation. Ce terme a été développé par Planche et Renon<sup>167</sup>, puis simplifié par Fürst et col.<sup>127,168</sup>. L'autre modèle, iPC-SAFT-Born, fait intervenir le terme de Born. Le modèle est caractérisé par un certain nombre de paramètres intervenant dans les différents termes. Trois paramètres caractérisent chaque ion : le diamètre ionique, le paramètre d'énergie dispersive de l'ion et un paramètre d'interaction binaire ion-solvant. Pour certains sels, un paramètre binaire cation-anion est utilisé. Ces paramètres ont été optimisés en ajustant les données expérimentales de coefficients d'activités moyens, de coefficients osmotiques, de densités liquides et de volumes molaires relatifs apparents à 25°C et 1 bar. Nous avons considéré dans cet ajustement 71 systèmes d'électrolytes forts de type 1:1, 1:2, 2:1 et 2:2, dans un domaine de molalités allant jusqu'à 6 kg.mol<sup>-1</sup>. Les coefficients d'activité et les volumes apparents sont des propriétés très sensibles par rapport aux paramètres du modèle. Des résultats comparables ont été obtenus pour les deux modèles avec un léger avantage pour iPC-SAFT-SR2. Comme le montre les Figure 5-9 et, le modèle iPC-SAFT-SR2 permet de bien décrire les données de coefficients d'activités moyens et les masses volumiques des solutions d'électrolytes forts
En utilisant des idées proches de celles de Fürst et Renon ¹²⁷, nous avons corrélé certains paramètres du modèle iPC-SAFT-SR2 (comme l'énergie de solvation) aux diamètres des ions, dans le but de réduire le nombre total de paramètres. Des essais ont également été faits afin de vérifier les capacités du modèle à représenter les données à forte concentration en sel. De bonnes prédictions ont été obtenues pour les solutions de LiBr. Enfin, ce modèle a été étendu en température pour plusieurs électrolytes forts, le domaine étudié allant de 298.15K à 498.15 K. Les tests ont porté sur les chlorures et sulfates, deux séries d'électrolytes pour lesquelles nous disposions de grandeurs expérimentales sur un large intervalle de température. La représentation des données par le modèle iPC-SAFT-SR2 est très satisfaisante (Figure 5-11).

# 6. Modélisation avec les Approches COSMO-RS et COSMO-SAC: Applications aux Biocarburants.

Les modèles de type COSMO sont des modèles prédictifs à coefficients d'activité (ou à énergie de Gibbs d'excès), comme l'est également le modèle UNIFAC. Ces modèles peuvent être très utiles pour prédire les propriétés de mélange de molécules oxygénées et polyfonctionnelles (par exemple celles issues de la biomasse végétale) dont les propriétés thermodynamiques et leur comportement dans différents solvants ne sont pas connus. Les modèles prédictifs peuvent alors être utilisés pour prédire par exemple la solubilité (solide-liquide, liquide-vapeur ou liquide-liquide) d'une molécule dans un solvant donné, et les coefficients de partage (liquide-vapeur, liquide-liquide) d'une molécule en présence de deux solvants.

Klamt et col. ¹⁶⁹-¹⁷¹,¹⁷²,¹⁷³ ont tout d'abord proposé le modèle COSMO-RS (Conductor-Like Screening MOdel for Real Solvent). Un calcul COSMO consiste à faire un calcul ab initio (on utilise en général la méthode DFT (Density functional Theory) pour un calcul cosmo) sur une seule molécule entourée d'un conducteur idéal censé représenter le solvant. La cavité qui entoure la molécule et définit la surface du conducteur idéal, est construite comme la somme de surfaces sphériques centrées autour des atomes de la molécule. Le calcul COSMO détermine les charges de polarisation sur la surface du conducteur. Le signe de ces charges surfaciques est localement opposé à celui des charges partielles de la molécule. La Figure 6-1 illustre un calcul COSMO pour la molécule de formaldéhyde.

<figure><img src="image_51.png" /><figcaption>Figure 6-1 Description d'un calcul ab initio de type COSMO, et découpage de la surface de la cavité en segments de surface.</figcaption></figure>

La surface de la cavité est découpée en segments surfaciques de surface $s_i$, et une charge surfacique $q_i$ est attribuée à chaque segment $i$, lesquels sont repérés par le vecteur position $\mathbf{r}_i$ (Figure 6-1). On détermine alors le "profile-sigma", c'est à dire la distribution de charge surfacique sur la surface de la cavité. Dans les approches COSMO-RS et COSMO SAC, chaque molécule est caractérisée par sa surface $A_i$, son volume $v_i$, et son profile sigma $p_i(\sigma)$. Nous allons maintenant donner les expressions des modèles COSMO-RS et COSMO-SAC sans rentrer dans le détail, puis nous donnerons des exemples d'applications. Enfin, nous montrerons un exemple de simulation d'un décanteur sous Aspen plus, faisant intervenir nos propres codes de calculs pour le modèle thermodynamique.
## 6.1. COSMO-RS

Dans l'approche COSMO-RS<sup>169,170,172</sup>, le potentiel chimique $\mu_{i/S}$ du composé $i$ dans la solution $S$ est donné par

$$ \mu_{i/S} = \int A_i p_i(\sigma) \mu_S(\sigma) d\sigma + \mu_{i/S}^C + kT \ln x_i. \quad (6.1) $$

$\mu_{i/S}^C$ est la partie combinatoire du potentiel chimique, $x_i$ la fraction molaire du constituant $i$ ; $\mu_S(\sigma)$ est appelé le potentiel-sigma de la solution<sup>173</sup> et est donné par

$$ \mu_S(\sigma) = -kT \ln \left( \int p_S(\sigma') \exp \left( \frac{a_{eff}(\mu_S(\sigma') - e(\sigma, \sigma'))}{kT} \right) d\sigma' \right). \quad (6.2) $$

où $a_{eff}$ est une constante (paramètre universel), et $e(\sigma, \sigma')$ l'énergie d'interaction entre les segments de surface $\sigma$ et $\sigma'$, qui comprend un terme électrostatique, un terme de dispersion, et un terme de liaison hydrogène. $p_S(\sigma)$ est le profile-sigma de la solution, obtenu par une règle de mélange par rapport aux profiles-sigma des constituants du mélange<sup>173</sup>. Le coefficient d'activité est ensuite obtenu avec l'expression

$$ \gamma_i = \exp \left( \frac{\mu_{i/S} - \mu_{i/i}}{kT} \right), \quad (6.3) $$

où $\mu_{i/i}$ est le potentiel chimique du constituant $i$ pure. Une description plus détaillée du modèle est donnée dans le livre de Klamt<sup>173</sup>.

## 6.2. COSMO-SAC

Lin et Sandler<sup>174</sup> ont proposé une version modifiée de COSMO-RS, appelée COSMO-SAC. Les expressions de COSMO-RS ont été simplifiées et réarrangées (intégrales remplacées par des sommes), et les termes de combinaison ainsi que les constantes universelles des deux modèles sont différents. D'autre part, COSMO-SAC ne contient pas de terme de dispersion dans l'énergie d'interaction entre segment, contrairement à COSMO-RS. Le coefficient d'activité du constituant $i$ est donné par

$$ \ln \gamma_i = \sum_{\sigma} \frac{A_i p_i(\sigma)}{a_{eff}} [\ln \Gamma_S(\sigma) - \ln \Gamma_i(\sigma)] + \ln \gamma_{i/S}^C, \quad (6.4) $$

où $\ln \gamma_{i/S}^C$ est le terme de combinaison de Staverman-Guggenheim. Les termes $\Gamma_S(\sigma)$ et $\Gamma_i(\sigma)$ sont les coefficients d'activités des segments de charges $\sigma$ et sont déterminés en résolvant

$$ \ln \Gamma_l(\sigma) = -\ln \left( \sum_{\sigma'} p_l(\sigma') \Gamma_l(\sigma') \exp \left( \frac{-\Delta W(\sigma, \sigma')}{kT} \right) \right), \quad (6.5) $$

où l'indexe $l$ vaut pour tous les constituants ou pour la solution $S$. L'énergie d'interaction $W(\sigma, \sigma')$ entre deux segments de charges $\sigma$ et $\sigma'$ est la somme d'un terme électrostatique (appelé énergie « misfit » par Klamt) et d'un terme décrivant les liaisons hydrogène.

## 6.3. Résultats

### 6.3.1. Coefficients de partage d'arômes

Dans le cadre d'une coopération avec AgroParisTech en 2007, nous avons étudié les coefficients de partage liquide-vapeur de molécules organiques oxygénées (arômes) dans des solutions hydro-alcooliques, l'application industrielle étant la distillation d'alcools forts. Ce travail a été réalisé en grande partie par Mohammed Ellaite lors de son stage de Master que j'ai co-encadré, et a conduit à une publication<sup>175</sup>. Nous avons mesuré les coefficients de partage (mesures faites à AgroParisTech), et modélisé les équilibres de phases avec NRTL<sup>176</sup> et COSMO-SAC<sup>174</sup>. Comme le montre la Figure 6-2,
les prédictions de COSMO-SAC sont très satisfaisantes pour ces systèmes, sur l’ensemble du domaine
de composition en éthanol.

<figure><img src="image_52.png" /><figcaption>Figure 6-2 Coefficient de partage liquide-vapeur ($K_i = y_i/x_i$) d'arômes à dilution infinie dans des solutions hydro-alcooliques, à pression atmosphérique. Les symboles correspondent aux données expérimentales¹⁷⁵ : (a). (○) isobutanol, (*) cis-3-hexenol et (●) 2-phenylethanol; (×) isopentanol and (♦) linalool; (◊) isobutanal et (□) furfural; (d) (Δ) ethyl acetate, (■) isopentyl acetate et (+) 2-phenylethyl acetate. Les lignes en pointillés ont été calculées avec le modèle NRTL ajusté sur les données expérimentales des binaires et des ternaires. Les lignes continues sont les prédictions du modèle COSMO-SAC.</figcaption></figure>

6.3.2. Etude des propriétés des constituants de biohuile

Le projet memobiol (2008-2013), financé par l'ANR et dirigé par l'IFPEN, concernait l'étude des
propriétés thermodynamiques de mélanges impliquant des molécules oxygénées extraites des bio-
huiles. Les bio-huiles sont des mélanges liquides contenant des centaines de constituants (eau, alcools,
cétones, furanes, esters, acides carboxyliques,...), et sont issues de la pyrolyse flash de biomasses
végétales. Certaines de ces molécules sont à haute valeur ajoutée, et peuvent être utilisées comme
matières premières pour la chimie verte. Il est donc intéressant de pouvoir extraire ces molécules des
bio-huiles.
<figure><img src="image_53.png" /><figcaption>Figure 6-3 Coefficients de partage octanol-eau à 298.15K, prédits par COSMO-RS <sup>177</sup> (carrés), COSMO-SAC <sup>174</sup> (triangles) et mCOSMO-SAC <sup>178</sup> (cercles), et comparés aux données expérimentales <sup>179</sup> (la ligne correspond à la droite $Log_{10}P_i^{exp} = Log_{10}P_i^{cal}$).</figcaption></figure>

<figure><img src="image_54.png" /><figcaption>Figure 6-4 Equilibre solide-liquide du mélange binaire eau-sorbitol à pression atmosphérique. Les cercles correspondent aux données expérimentales et les lignes continues correspondent aux prédictions du modèle COSMO-SAC. Les lignes en pointillés sont les prédictions en utilisant l'hypothèse de la solution idéale.</figcaption></figure>

L'extraction de ces molécules est cependant compliquée, et nécessite l'utilisation de différents solvants, et de plusieurs techniques de séparations (extraction liquide-liquide, distillation, solide-liquide, ...). Le design et l'optimisation des procédés d'extraction nécessitent une bonne connaissance des équilibres de phase. Les mesures expérimentales d'équilibre de phase sont parfois inexistantes pour les mélanges de ces molécules avec des solvants. Les modèles prédictifs tels que les contributions de groupes <sup>18,44,180</sup> ou les approches type COSMO peuvent alors être très utiles.

Au cours du projet, Jinlong Li, postdoc sous ma responsabilité, a utilisé plusieurs modèles de type COSMO (COSMO-RS, COSMO-SAC, m-COSMO-SAC, ..) pour prédire les équilibres de phases (liquide-vapeur, liquide-liquide et solide-liquide) ainsi que les coefficients de partage octanole-eau de molécules issues des bio-huiles. Deux articles ont été publiés sur ce sujet <sup>181,182</sup>. Les approches de type COSMO donnent des prédictions assez satisfaisantes comme le montre les Figures 6.3 et 6.4.