Hadj Kali
Hadj Kali
Hadj Kali
THESE
prsente
pour obtenir
MIJOULE Claude
Prsident
M.
GERBAUD Vincent
RICHON Dominique
UNGERER Philippe
JOULIA Xavier
LETOURNEAU Jean Jacques
Directeur de thse
Rapporteur
Rapporteur
Membre
Membre
Rsum en Franais
La simulation molculaire est un outil prcieux pour la prdiction de proprits thermodynamiques des
molcules pour lesquelles les donnes sont rares. Base sur la thermodynamique statistique, la simulation
molculaire ncessite des techniques dchantillonnage des tats du systme efficaces et une description
prcise des interactions nergtiques au sein du systme au travers dun champ de force. Ce dernier doit
tre prcis pour raliser des prdictions de qualit comparable celle des mesures exprimentales. La
gnricit dun champ de force est souhaitable et est acquise en utilisant des paramtres pour des groupes
chimiques qui permettent de recomposer les molcules chimiques.
Nous contribuons llaboration dun champ de force gnrique en optimisant des paramtres de Lennard
Jones de la contribution de Van der Waals pour le groupe nitrile. Pour cela, des simulations dquilibre
liquide vapeur dans lensemble de Gibbs Monte Carlo et des simulations de la phase liquide condense
dans lensemble isobare isotherme NPT sont ralises. Le caractre polaire des nitriles ncessite de
prendre en compte la contribution lectrostatique au travers de charges ponctuelles issues dune analyse de
population de type MEP ou Mulliken. Les paramtres de Lennard Jones rgresss avec lanalyse MEP
reproduisent trs bien les proprits liquide et vapeur de lactonitrile, mieux quavec lanalyse Mulliken et
les paramtres de Lennard Jones associs. Mais la transfrabilit (prdiction des proprits du propionitrile
et n-butyronitrile sans ajustement des paramtres) est dmontre pour la seule analyse Mulliken. Des
critres de transfrabilit des paramtres pour les molcules polaires sont noncs : cohrence des charges
ponctuelles sur les sries chimiques homologues, pertinence physique des paramtres de Lennard Jones.
Une autre contribution consiste laborer une mthodologie de recherche des ventuels azotropes pour
des fluides de Lennard Jones. Elle combine des critres macroscopiques de Gnie des Procds et un
enchanement de simulations molculaires dans lensemble de Gibbs Monte Carlo et un mouvement
dchange didentit original.
Rsum en Anglais
Application of molecular simulation for the calculation of liquid vapour equilibria of nitriles
and for the prediction of azeotropes
Molecular simulation is invaluable for the prediction of thermodynamic properties for molecules for which
the data are rare. Based on statistical thermodynamics, molecular simulation requires efficient sampling
techniques of the system states and a precise force field that describe the energy interactions within the
system. The precision of a force field conditions that of the predictions which can be closed to that of
experimental measurements. Generic force fields are acknowledged by using parameters for chemical
groups which compose the chemical molecules.
We contribute to the development of a generic force field by optimizing Lennard Jones parameters in the
Van der Waals contribution for the nitrile group. Simulations of liquid - vapor equilibria using Gibbs
ensemble Monte Carlo and simulations of the condensed liquid phase using NPT simulations are carried
out. The polar character of nitriles requires the electrostatic contribution through atomic charges to be
considered, two population analysis MEP and Mulliken are compared. The set of Lennard Jones parameters
regressed with MEP charges reproduces accurately the liquid and vapor properties of acetonitrile, better
than with Mulliken charges and their associated Lennard Jones parameters. But the parameter
transferability (prediction of properties of propionitrile and n-butyronitrile without adjustment of the
parameters) is effective only for Mulliken charges. Criteria of parameter transferability for polar molecules
are stated: coherence of the atomic charges in homologous chemical series, physical relevance of the
Lennard Jones parameters.
The second contribution consists in a methodology for the search of possible azeotropes for Lennard Jones
fluids. It combines macroscopic criteria of chemical engineering and a sequence of Gibbs ensemble and a
novel identity exchange movement developed within the Grand Canonical ensemble formalism.
Discipline
Mots Cls
Laboratoire
Je ddie ce travail :
A ma petite Tasnim, Mouna et Khawla,
Mohamed Amine, Sabir, Isshaq, et Abderrahmane ;
Ainsi qu toute ma famille.
Chapitre I.
Chapitre II.
A.
Introduction ..........................................................................................................II-1
2.
Moyenne ..............................................................................................................II-1
3.
a.
b.
b.
c.
B.
i.
ii.
iii.
iv.
Introduction ........................................................................................................II-13
2.
3.
a.
Introduction..............................................................................................II-13
b.
c.
d.
Echantillonnage uniforme........................................................................II-16
e.
f.
g.
Algorithme de Metropolis.........................................................................II-20
Introduction..............................................................................................II-22
b.
quations du mouvement........................................................................II-23
c.
d.
Thermostats T et P ..................................................................................II-26
4.
Ergodicit ...........................................................................................................II-27
5.
6.
b.
C.
c.
d.
e.
Le biais configurationnel..........................................................................II-33
f.
2.
Diffrentes approximations......................................................................II-35
b.
c.
d.
e.
ii.
iii.
La mcanique molculaire.................................................................................II-42
a.
b.
c.
d.
e.
f.
D.
i.
i.
ii.
ii.
ii.
iii.
ii.
b.
La rgle de phase....................................................................................II-55
c.
d.
ii
2.
b.
3.
i.
ii.
iii.
ii.
iii.
iv.
c.
d.
b.
E.
ii.
iii.
iv.
ii.
iii.
Pseudo-ensembles .......................................................................II-82
iv.
Chapitre III.
A.
B.
C.
1.
2.
3.
4.
Actonitrile................................................................................................III-8
b.
Propionitrile...............................................................................................III-9
iii
c.
D.
2.
3.
2.
Mthodes dchantillonnage...................................................................III-14
b.
ii.
iii.
iv.
La sommation Ewald....................................................................III-16
v.
vi.
c.
d.
e.
b.
E.
n-butyronitrile..........................................................................................III-10
ii.
iii.
iv.
v.
vi.
ii.
iii.
iv.
v.
vi.
iv
Chapitre IV.
A.
B.
2.
3.
C.
D.
Rsultats......................................................................................................... IV-9
E.
1.
2.
3.
Chapitre V.
Rfrences
Annexe I.
Annexe II.
Annexe III.
Annexe IV.
Chapitre I.
Introduction gnrale
La connaissance des quilibres entre phases est importante pour la conception et la simulation de
beaucoup de procds de sparation comme la distillation et l'extraction liquide liquide. Lorsque les
donnes sont absentes, des modles thermodynamiques peuvent les gnrer, comme des modles de
coefficient d'activit et des quations d'tat prsents dans le chapitre II. Cependant, ds linstant o les
interactions sont non ngligeables, les constituants en prsence se comportent de faon non idale et les
modles thermodynamiques ncessitent des paramtres empiriques d'interaction binaire rgresss partir
de donnes exprimentales existantes et souvent incompltes. La consquence de cette rgression
empirique est que la valeur prdictive de tels modles est faible. C'est un inconvnient srieux sachant
lnorme diffrence entre le nombre de constituants chimiques rfrencs dans les Chemical Abstract
Series (plus de 15 millions) et le nombre de donnes d'quilibre de phase ( peine 30 000) disponibles dans
les bases de donnes les plus compltes.
Des modles thermodynamiques macroscopiques qualits prdictives ont aussi t proposs : ASOG,
UNIFAC seuls ou associs une rgle de mlange complexe type PSRK, MHV2 ou Wong Sandler. Ils
emploient le concept intressant d'additionner des contributions chimiques de groupe pour prdire la
valeur des paramtres dinteraction binaire. Mais la prcision des donnes obtenues est ingale, moins
que des sous-groupes secondaires voire ternaires soient dfinis, gchant l'intrt de la mthode du fait de
lintroduction de plus de complexit : Modified UNIFAC, etc, .
Cet tat de fait est une incitation pour l'arrive de nouvelles mthodes permettant de prvoir des
proprits physico-chimiques avec une grande prcision. Parmi ces nouveaux outils, la simulation
molculaire a merg comme un outil complmentaire pour construire les passerelles entre les dtails
microscopiques d'un systme (atomes, interactions nergtiques, distribution des molcules, etc.) et les
proprits macroscopiques d'intrt pour le Gnie des Procds (tat physique, coefficients de transfert,
proprits d'quilibre, etc.). Ces passerelles sont dsormais possibles grce dune part aux concepts de
thermodynamique statistique qui est, par ailleurs, lorigine de la majorit des modles thermodynamiques
employs en Gnie des Procds, et dautre part laugmentation continuelle de la puissance des
ordinateurs qui a permis, en quelques annes, de rduire considrablement les temps de calcul de plusieurs
jours quelques heures.
Comportant beaucoup d'analogies avec la vraie exprimentation de laboratoire, la simulation molculaire
ralise fondamentalement des expriences numriques directement sur les systmes modles qui visent
gnrer un grand nombre de configurations sur lesquelles des proprits sont moyennes et galises
leur valeur macroscopique grce au premier principe de la thermodynamique statistique. Outre cette
fonction de production de donnes physico-chimiques prcises qui est lobjectif principal de la thse, la
simulation molculaire est et sera de plus en plus dans le futur une source de connaissance des
I-1
phnomnes se droulant lchelle molculaire [ Chen et Mathias, 2002 ; De Pablo et Escobedo, 2002 ;
Sandler 2003 ].
Le chapitre II rappelle les concepts essentiels de la thermodynamique statistique, les techniques
dchantillonnage et les modles de mcanique molculaire employs en simulation molculaire. La
relation avec les modles thermodynamiques macroscopiques est aussi prsente ainsi que les principales
classes de modles reprsentant les quilibres liquide vapeur en Gnie des Procds.
On recense gnralement trois verrous scientifiques pour lapplication de la simulation molculaire en
Gnie des Procds [ Gerbaud, 2003 ] :
(i)
(ii)
llaboration de mthodologies thoriques permettant le calcul direct des proprits qui intressent
le Gnie des Procds. En effet, les mthodes de simulation molculaire actuelles permettent
dobtenir des proprits densemble (courbe dquilibre liquide vapeur) mais cest presque
toujours un point particulier (azotrope, temprature de rose, coordonnes critiques, ) qui
manque en Gnie des Procds.
(iii)
Lintgration des mcanismes molculaires dans des modles multichelles applicables pour ltude
de systmes macroscopiques et de leurs proprits.
Le chapitre III sattaque au premier verrou en dfinissant une dmarche pour dterminer des paramtres
nergtiques de Lennard Jones, epsilon () et sigma (), gnriques pour le groupe CN des nitriles. On
appelle nitriles les composs comportant le groupe cyano (CN) li une chane alkyle ou phnyle
(RCN). Ce sont des molcules polaires dintrt industriel certain :
Solvants industriels slectifs dans des procds divers comme lextraction des cires et des graisses,
Intermdiaires de synthse pour des produits pharmaceutiques et des pesticides,
Intermdiaires pour la fabrication de fibres textiles (lacrylonitrile).
Pour cela, nous avons appliqu la mthodologie propose par Ungerer et al. [ 2000 ] et appliqu sur les nalcanes pour lobtention des paramtres epsilon () et sigma () des n-alcanes1. Les rsultats montrent quil
est possible de rgresser ces paramtres pour le groupe CN sur une molcule nitrile de rfrence,
lactonitrile, et de les employer avec succs pour prdire, sans aucun nouvel ajustement, les proprits
dautres nitriles. La qualit des prdictions en terme de prcision est remarquable et confirme le potentiel
de la simulation molculaire pour la prvision de proprits thermodynamiques prcises. Cependant, une
Notons, cependant, que lutilisation du modle AUA (Anisotropic United Atom) pour la dtermination des paramtres
nergtiques des groupements CH2 et CH3 des n-alcanes ncessite loptimisation dun troisime paramtre delta () qui dfinit la
distance entre le carbone et le centre de force du groupement chimique tudi. Pour les nitriles, =0.
I-2
attention particulire doit tre porte lobtention des charges ponctuelles utilises dans lvaluation de la
contribution lectrostatique de lnergie dinteraction. En effet, les deux jeux de charges lectrostatiques
utiliss dans ce travail sont issus dune analyse de population diffrente : lanalyse Mlliken classique et
lanalyse MEP base sur le calcul de la surface de potentiel. Les rsultats obtenus permettent dnoncer les
critres dobtention dun jeu de paramtres de Lennard Jones transfrable et montrent lattention quil faut
porter ltude prliminaire de mcanique quantique.
Le chapitre IV rpond au deuxime verrou en proposant une mthodologie pour la prdiction des
azotropes de mlanges binaires. Au-del de leur signification de non idalit du mlange, limportance des
azotropes est reconnue puisquelle conditionne directement la faisabilit et la conception des procds de
distillation. En particulier, essentielle est la capacit de prvoir si un mlange donn forme un ou plusieurs
azotropes et de calculer les conditions appropries de composition, temprature et pression qui y sont
associes. La mthodologie est base sur la combinaison, dune part, de proprits macroscopiques des
azotropes frquemment utilises en Gnie des Procds par lintroduction de la constante dquilibre Ki
dfinie comme le rapport des fractions molaires du constituant i en phase vapeur et en phase liquide, et
dautre part, de proprits de thermodynamique statistique permettant dchantillonner correctement un
systme molculaire possdant un azotrope en progressant le long de la courbe de coexistence liquide
vapeur. Cette progression seffectue par une succession de simulations dans lensemble de Gibbs Monte
Carlo volume impos et de mouvements de changement didentit. Se basant sur la fonction de partition
de lensemble grand canonique, on dmontre que la probabilit dacceptation du mouvement de
changement didentit est gale au produit des deux mouvements classiques dinsertion et de destruction.
Appliqus dans ce travail des fluides simples de Lennard Jones, les rsultats sont concluants et des pistes
pour ltendre aux fluides rels sont proposes.
Le chapitre V est consacr la Conclusion gnrale.
I-3
Notions thoriques
Chapitre II.
Ce chapitre prsente des notions thoriques sur la thermodynamique statistique et sur la simulation
molculaire (mthodes dchantillonnage et champs de forces). Un accent particulier est port sur tout ce
qui concerne la simulation des quilibres liquide vapeur dont la modlisation lchelle macroscopique
fait lobjet dune partie complmentaire en fin de chapitre. On trouvera des complments dans diverses
rfrences bibliographiques [ Mc Quarry, 1976 ; Allen et Tildesley, 1987 ; Frenkel et Smit, 1996 ; Dugas,
2000 ; Leach, 1999 ] mais aussi dans un site web ddi du Computational Molecular Science and
Engineering
Forum
prsentant
des
matriaux
pdagogiques
en
ligne
(http://www.ecs.umass.edu/che/am3/AIChEcour.html).
A.
Thermodynamique statistique
1.
Introduction
Moyenne
a.
Lorsquon mesure une proprit macroscopique M, par exemple la temprature, la valeur obtenue
est le rsultat des mouvements chaotiques et des collisions dun trs grand nombre de molcules
(lexemple du mouvement brownien). Ainsi, si cette proprit tait observe sur une trs petite chelle de
temps (par exemple 10-8 s), on remarquerait que cest une quantit fluctuante. En pratique, cependant, le
temps de mesure macroscopique est si grand (de lordre de 1 ms 10 s) par rapport lchelle de temps
des fluctuations que celles-ci ne sont pas observes et que la valeur de la proprit parat constante. En
dautres termes, les valeurs des proprits macroscopiques que lon mesure sont des moyennes
temporelles dun trs grand nombre dtats dgnrs du systme [ Prausnitz, 1999 ].
Formalisons ce concept en dsignant par une configuration particulire du systme tudi un instant
donn, et en supposant quon peut crire la valeur instantane de la proprit M comme une fonction
M(). Etant donn que le systme volue pendant le moment de la mesure, la valeur de M change et la
valeur macroscopique de la proprit M observe exprimentalement est gale la moyenne temporelle de
M() prise sur un long intervalle de temps :
II-1
M ( ( t )) dt
(II-1)
O les crochets < > indiquent quil sagit dune moyenne et le temps dobservation.
b.
Dans les annes 1880, Boltzmann et Gibbs ont propos une approche conceptuelle nouvelle du
problme du calcul de la valeur moyenne dune grandeur physique. A cause de la complexit de lvolution
en fonction du temps de la proprit thermodynamique M((t)), Gibbs a suggr de remplacer la
moyenne temporelle par une moyenne sur un ensemble statistique de configurations. Les configurations
qui sont des ensembles de positions et de moments des constituants du systme sont distribus selon une
densit de probabilit () qui reprsente le systme un instant donn.
De son cot, Boltzmann a quantifi cette densit en stipulant que la probabilit p de trouver un systme,
microscopique ou macroscopique, en quilibre une temprature T dans une configuration s dnergie Es,
est donne par la relation :
ps = Cte exp( Es )
(II-2)
O la grandeur = 1/(kBT), appele facteur de Boltzmann, est probablement la faon la plus naturelle
dexprimer le concept de temprature puisquelle souligne linaccessibilit dune temprature absolue
nulle : lorsque T tend vers 0, devient infini, ce qui implique une probabilit nulle, kB tant la constante
universelle de Boltzmann gale 1,38066 10-23 J/K. Nous allons revenir sur lorigine de ces paramtres
et la manire de les obtenir dans le paragraphe concernant les ensembles statistiques usuels. Notons juste
ici que cest cette distribution de probabilit qui est lorigine de la thorie cintique des gaz.
3.
Exemple prliminaire
II-2
Cet exemple illustre que le rsultat dune exprience rpte sur une seule boite est identique celui
obtenu avec seulement une exprience sur un ensemble de plusieurs boites.
Un tel ensemble, nomm ensemble statistique, est une construction intellectuelle cense reprsenter un
instant donn les diffrentes configurations qui, dans le systme rel, napparatraient quau cours du
temps. Cet ensemble est compos dun grand nombre de systmes, tous construits de la mme manire,
chacun deux tant diffrent et tant une rplique du systme rel fig dans lun des tats dgnrs
accessibles. Sil y a n tats accessibles, il y aura n systmes dans lensemble de telle sorte que chaque tat
quantique accessible au cours du temps par le systme rel est reprsent dans lensemble par un systme
dans cet tat quantique stationnaire.
b.
Un ensemble statistique est une collection virtuelle ou mentale dun nombre A trs grand de
systmes construits de faon tre la rplique au niveau thermodynamique du systme tudi.
Considrant un systme quelconque, son tat thermodynamique peut tre dcrit lchelle macroscopique
par quelques proprits macroscopiques comme la temprature T, la composition x et le volume V. A
lchelle microscopique, une description aussi simple est impossible cause du grand nombre dtats
quantiques possibles ayant les mmes proprits macroscopiques fixes : une mole contient 6,023.1023
molcules et donc le nombre de positions et de moments qui dfinissent une configuration est
incommensurable. Pourtant, notre but est de calculer des proprits thermodynamiques comme lnergie
et lentropie de ce systme partir de la connaissance des configurations lchelle molculaire du systme
et partir de la forme fonctionnelle dun potentiel dinteraction molculaire.
Grce aux travaux de Maxwell, Boltzmann et en particulier ceux de Gibbs qui est le premier avoir
introduit le concept dun ensemble statistique, on nonce le postulat suivant :
Afin de calculer la valeur de nimporte quelle proprit thermodynamique mcanique (par exemple la pression), on
calcule sa valeur dans chacun des tats quantiques en cohrence avec les quelques paramtres ncessaires pour spcifier
le systme au sens macroscopique et on fait la moyenne en attribuant le mme poids chaque tat quantique
possible .
Autrement dit, la moyenne sur lensemble statistique de cette proprit mcanique est gale la moyenne
temporelle macroscopique et donc correspond la proprit thermodynamique macroscopique que lon
cherche [ Diu et al., 1989 ; Allen et Tieldesley, 1987 ].
X ((t ))temps = X
ensemble
=X
postulat n1
(II-3)
II-3
1,
2,
3,
m,
Lnergie
E1 ,
E2 ,
E3,
Em,
Nombre doccupation
a1,
a2,
a3,
al,
II-4
=A
(II-4)
a E
j
(II-5)
Puisque lensemble canonique a lui aussi t isol du milieu extrieur, on peut lui appliquer le principe
dquiprobabilit qui stipule alors que toutes les distributions des nombres doccupation a1, a2, a3,
satisfaisant les deux dernires quations ((II-4) et (II-5)) sont quiprobables et doivent tre affectes du
mme poids lors de lvaluation dune grandeur moyenne.
Le nombre de faons quune distribution particulire des aj, W(a) W(a1 , a2 , a3 , ), peut tre ralise et
arrange en groupes, tel que a1 systmes soient dans le premier, a2 dans le second, etc., est gal :
W ( a ) = A! a ! a ! a !... = A!
1 2 3
ak !
(II-6)
Pour une distribution particulire, a j A est la fraction des configurations de lensemble canonique se
trouvant dans ltat nergtique j (c'est--dire ayant lnergie Ej). La probabilit totale pj quun systme
occupe ltat quantique j est obtenue en faisant une moyenne de a j A sur toutes les distributions
permises auxquelles nous avons appliqu le principe de probabilit gale. Ceci signifie quon exige que chacun
des tats quantiques soit reprsent un nombre gal de fois dans lensemble, puisquon ne dispose
daucune information pour prfrer un systme un autre. Une interprtation de ce principe permet
dnoncer que :
Tous les tats quantiques accessibles et distincts dun systme ferm dnergie fixe (NVE, microcanonique) sont
quiprobables
La probabilit pj est donne alors par :
pj =
aj
W (a).a (a)
=1
A
W (a )
j
(II-7)
Connaissant la probabilit de trouver un systme aux variables N, V et T fixes dans le jime tat quantique,
on peut calculer la moyenne de lensemble canonique de nimporte quelle proprit mcanique par :
M = M j pj
(II-8)
II-5
Cependant, les sommations impliques dans lquation (II-7) sont difficiles calculer mathmatiquement
parce quil existe un nombre infini dtats quantiques accessibles, et donc, en pratique, les quations (II-7)
et (II-8) sont trop compliques pour tre utilises.
Mais le fait de pouvoir faire tendre A vers linfini permet dutiliser les rsultats des thormes
mathmatiques qui sappliquent pour les grands nombres. En effet, il est facile de montrer que les
coefficients multinomiaux, comme W(a), sont extrmement aigus (prsentent un pic) autour de leur valeur
maximale si toutes les variables aj sont trs grandes. Ce rsultat dmontr avec une seule contrainte
(lquation (II-4)) dans Mc Quarry [ 1976 ] est gnral. Disposant dune autre contrainte (lquation (II-5)),
le pic maximal de W(a) ne sera plus un unique point o tous les aj sont gaux mais un ensemble de points
aj dont ltendue sera trs petit. Soit a* = {aj*} cette distribution.
Par consquence, la valeur de W(a) dans lquation (II-7) sera compltement ngligeable pour nimporte
quel ensemble de valeur aj autre que a*. Ceci permet de remplacer la sommation de lquation (II-7) pour
toutes les distributions par un seul terme valu pour a* :
pj =
aj
W ( a ).a ( a )
1
A= A
W ( a )
j
( a j )
W ( a * ).a *j
A W (a * )
a *j
(II-9)
O aj* est la valeur de aj dans cette distribution qui maximise W(a) et qui reprsente la distribution la plus
probable.
En comparant les quations (II-7) et (II-9), on obtient :
pj =
aj
A=
a *j
(II-10)
Ainsi, pour calculer les probabilits ncessaires pour lvaluation des moyennes densemble, on a besoin de
dterminer la distribution a* qui maximise W(a) sous les deux contraintes (II-4) et (II-5). On dfinit ainsi
un problme doptimisation dune fonction plusieurs variables soumises aux contraintes (II-4) et (II-5).
En appliquant la mthode des multiplicateurs de Lagrange, lensemble des valeurs aj qui maximise W(a) est
obtenu par la rsolution de lquation :
ln W (a ) a k a k Ek = 0
a j
k
k
(II-11)
ln a *j 1 E j = 0
j = 1, 2,...
(II-12)
II-6
Do :
j = 1, 2,...
(II-13)
Lquation (II-13) nous donne la distribution la plus probable en fonction de et . Cependant, on peut
exprimer en fonction de en faisant la sommation de lquation (II-13) sur tous les j et en appliquant la
contrainte (II-4) pour obtenir :
exp( + 1) = 1 A exp E j
j
(II-14)
Ainsi, la probabilit de trouver le systme dans ltat quantique j donne par lquation (II-9) devient :
pj =
a *j
A
exp E j (N , V )
exp( E (N ,V ))
(II-15)
Par exemple, en substituant ce rsultat dans lquation (II-8), avec lnergie Ej prise comme proprit
mcanique, on obtient lexpression de lnergie moyenne E dans lensemble canonique qui, daprs le
postulat de Gibbs, correspond lnergie thermodynamique E :
E (N ,V ) exp( E
j
E = E( N , V , ) =
exp( E
( N ,V )
( N ,V )
(II-16)
La pression P1 est une autre proprit mcanique importante. Lorsque le systme se trouve dans ltat j,
dE j = Pj dV est le travail fourni par le systme quand son volume augmente de dV, tout en
maintenant constant le nombre de particules. La pression est alors donne par :
E j
Pj =
V
(II-17)
E j
P = p j Pj =
j
V
j
exp E j
N
exp E j
(II-18)
II-7
Le dnominateur des quations (II-16) et (II-18) apparat travers toutes les moyennes dans lensemble
canonique. Cette quantit est appele la fonction de partition de lensemble canonique quon distinguera
par la notation :
Q NVT = Q (N , V , T ) = Q (N , V , ) = exp E j (N , V )
(II-19)
Ainsi, la connaissance de la fonction de partition d'un systme permet d'accder l'ensemble de ses
grandeurs thermodynamiques. Par exemple, en drivant lquation (II-19) par rapport ses variables et en
appliquant la dfinition dune grandeur moyenne donne par lquation (II-8), on obtient :
d (ln Q NVT ) = Ed
(II-20)
On transforme cette relation pour faire apparatre les multiplicateurs de Lagrange plutt que leurs
drives, ce qui donne :
d ln Q NVT + E = d E
(II-21)
Cette expression peut tre compare lexpression gnrale de lentropie S , qui peut tre obtenue par la
relation thermodynamique fondamentale :
d E = Td S Pd V + i d N i
(II-22)
dS =
dE
T
(II-23)
Selon le postulat de Boltzmann, les quations (II-21) et (II-23) sont relies par une constante kB, ce qui
permet dcrire :
S = k B ln Q NVT + E
(II-24)
Lanalogie entre les quations (II-21) et (II-23) permet didentifier le multiplicateur de Lagrange en
fonction de la temprature du systme, en effet :
dS =
dE
et d S = k B d ln Q NVT + E = k B d E
T
1
kBT
(II-25)
II-8
F = kBT ln Q NVT = E T S
(II-26)
La dmonstration que nous venons de prsenter est celle donne par Mc Quarry pour un systme
monophasique monoconstituant [ McQuarry, 1976 ]. Cest ce dveloppement qui est lorigine de toutes
les quations utilises en simulation molculaire. Ungerer a prsent dans sa thse [ Ungerer, 1999 ] une
adaptation de cette approche pour le cas gnral diphasique multiconstituant en fixant toutes les variables
thermodynamiques intensives (temprature, pression et potentiels chimiques) constantes, ce qui a permis
dcrire lexpression gnrale de la fonction de partition et de dduire par la suite les rsultats connus non
seulement pour les ensembles classiques (NVT, NPT) mais aussi pour lensemble de Gibbs introduit par
Panagiotopoulos en 1987 que nous dcrirons plus en dtail plus loin.
Nous passons ensuite rapidement en revue les principaux ensembles utiliss en mcanique statistique.
Lobtention de la fonction de partition correspondante chacun de ces ensembles ncessitant un
dveloppement similaire, nous nous limiterons ainsi citer uniquement lexpression mathmatique et
quelques caractristiques pour chacun des ensembles statistiques classiques. Nous supposons qu la limite
thermodynamique (pour un systme de taille infinie) et loin du voisinage o se produit les transitions entre
phases, les diffrents ensembles conduisent des proprits thermodynamiques moyennes identiques.
Autrement dit, lnergie calcule est gale quelque soit lensemble utilis.
ii.
Soit un systme de volume V, contenant N molcules dun mme constituant et ayant une nergie
E. Cest un systme isol avec N, V et E fixs. Lensemble correspondant de A configurations aura un
volume total AV, contiendra AN molcules et aura pour nergie totale = AE. Bien que toutes les
configurations de lensemble soient identiques du point de vue de la thermodynamique macroscopique
NVE, elles ne le sont pas lchelle molculaire.
Lensemble NVE est appel lensemble microcanonique et il est utile pour des discussions thoriques.
Pour des applications plus pratiques cependant, on ne considre pas les systmes isols, mais ceux pour
lesquelles la temprature est fixe au lieu de lnergie, comme pour lensemble canonique NVT ou
lensemble isotherme-isobare NPT.
Dans lensemble NVE, le potentiel thermodynamique associ la fonction de partition NVE est
lentropie S :
S = kB ln NVE
(II-27)
II-9
iii.
Lensemble NPT est largement utilis dans les simulations de type Monte Carlo notamment parce
quil permet une comparaison directe avec des rsultats issus dexprimentations relles qui sont
frquemment ralises temprature et pression contrles. La fonction de partition NPT de cet ensemble
est donne par lexpression :
( (
))
(II-28)
Notons que la quantit qui apparat dans lexponentielle de cette fonction, lorsquon value la moyenne,
donne lenthalpie thermodynamique H = E + PV . Dans cet ensemble, le volume V rejoint la liste des
quantits microscopiques (les coordonnes et les vitesses) dcrivant les tats quantiques du systme. Le
potentiel thermodynamique associ lensemble NPT est lnergie libre G. En effet, on peut montrer
facilement que :
G = kBT ln NPT = H T S
iv.
(II-29)
Lensemble grand canonique VT
Pour quelques applications, il est commode de travailler dans un ensemble o les systmes peuvent
changer non seulement de lnergie mais aussi de la matire (exemple au chapitre IV). Un ensemble
similaire est appel ensemble grand canonique qui correspond un nombre important de systmes
ouverts de volume constant, en quilibre interne et pouvant changer de lnergie et de la matire (les
molcules) avec le milieu extrieur.
Pour dcrire ltat thermodynamique dun tel ensemble, il est pratique dutiliser comme variables
caractristiques indpendantes la temprature T, le volume V et les potentiels chimiques 1, 2, des
constituants 1, 2, du systme. La fonction de partition est donne alors par la relation :
( (
))
(II-30)
Le potentiel thermodynamique associ cet ensemble est le produit -PV tant donn que :
PV = kBT ln Q VT
(II-31)
Il est possible de construire beaucoup dautres ensembles, mais ils ne sont pas tous intressants pour la
simulation molculaire. En effet, considrons lexemple de lensemble statistique PT o les paramtres
constants sont toutes des variables intensives ; dans ce cas, les quantits extensives correspondantes sont
sujettes des fluctuations, ceci signifie que la taille des systmes peut augmenter sans limite. Dautre part,
le vecteur potentiel chimique , la pression P et la temprature T du systme sont relis par une quation
dtat, et donc, mme si lexpression de cette quation tait connue, ces trois variables ne sont pas
II-10
indpendantes. Pour toutes ces raisons, la simulation dans un ensemble PT ou dautres ensembles
prsentant les mmes caractristiques nest pas pratique. Cest pourquoi, la majorit des ensembles usuels
en simulation molculaire ncessite la spcification, dau moins, une variable extensive (gnralement le
volume V ou le nombre de particule du systme N) afin de pouvoir contrler la taille du systme tudi.
Notons enfin que les couples typiques de variables intensives et leurs variables extensives conjugus sont :
(T, E), (P, V) et (, N) ce qui nous permet de retrouver, par combinaison directe, le choix des variables
caractristiques des ensembles classiques.
c.
La base thorique sur laquelle repose toute la thermodynamique est la mcanique statistique. Par
consquent, la mcanique statistique a t utilise comme point de dpart pour le dveloppement de
plusieurs quations dtat succinctement rsumes comme tant de la forme f(P,V,T) = 0.
En effet, comme nous lavons dcrit un peu plus loin, cest dans la fonction de partition dun ensemble
donn que se concentrent toutes les informations concernant le systme tudi. De ce fait, lorsquon fixe
la temprature T, le volume du systme V et le nombre de particules N comme variables caractristiques
indpendantes (c'est--dire lorsquon travaille dans lensemble canonique), toutes les informations
dcrivant le systme sont contenues dans la fonction de partition canonique ayant pour expression :
Q NVT = e
E j ( N ,V )
kB T
(II-19)
Le potentiel thermodynamique caractristique de cet ensemble est lnergie libre dHelmholtz, dont on
rappelle lexpression :
F (N , V , T ) = k B T ln Q NVT
(II-26)
Ainsi, une fois la fonction de partition QNVT ou, similairement, lnergie libre du systme F connue en
fonction des variables caractristiques N, V et T, toutes les autres proprits thermodynamiques seront
accessibles [ Sandler, 1994 ] ; en particulier, lquation dtat sera donne par la relation :
ln Q NVT
F
=
P (N ,V , T ) = k B T
V T , N
T ,N
V
(II-32)
Dautre part, lexpression de la fonction de partition peut tre prcise en sparant les contributions
cintique et potentielle lnergie du systme et en remplaant les sommations discrtes par des intgrales,
ce qui permet dcrire la fonction de partition sous la forme :
Q NVT =
q i (T )N
Z NVT
N!
(II-33)
II-11
Z NVT = L e
U (r1 , r2 ,Kr N )
kT
dr1dr 2 L dr N
(II-34)
est le facteur de Boltzmann de lnergie dinteraction U intgre sur toutes les positions des particules.
Si nous considrons un systme idal ou les particules ninteragissent pas entre elles, U prendra la valeur 0
et on obtiendra immdiatement Z = V N , ce qui, introduit dans lquation (II-32), donne lquation dtat
dun gaz parfait :
qi (T )N N
ln V
ln Q
ln
V = N
=
=
= NV
k BT V
!
V
N
V
T , N
T ,N
T , N
(II-35)
En multipliant cette quation par le nombre dAvogadro NA, on obtient lexpression classique de
lquation dtat dun gaz parfait :
PV = nRT ou
P
=
k BT
(II-36)
O est la densit molculaire. La constante des gaz parfaits R est gale au produit de la constante de
Boltzmann kB=1,38066.10-23 J/K et du nombre dAvogadro NA=6,023.10+23 entites. Il sagit dune
illustration remarquable de la faon dont une hypothse microscopique trs simple (pas dnergie
dinteraction) permet de driver une quation macroscopique importante.
Ce comportement est approch par tous les gaz rels aux faibles densits. A de plus fortes densits, les
particules se rapprochent et interagissent de faon non ngligeable. Dans ce cas, lintgrale de
configuration ZNVT nest plus seulement gale VN. Deux approches peuvent alors tre employes pour
caractriser le comportement non idal des fluides loin de la zone des faibles densits :
-
Lquation dtat du viriel dcrit les dviations du comportement idal par un dveloppement
limit en srie de Taylor de la densit molculaire :
P
= + B(T ) 2 + C (T ) 3 + ...
k BT
quation du viriel
(II-37)
Dans cette quation, B(T) est le second coefficient du viriel et est associ aux interactions entre
paire de molcules, C(T) est le troisime coefficient du viriel et est associ aux interactions entre
triplets de molcules, etc.
-
A lexception des cas de gaz rare ou prsentant une infime dviation par rapport au
comportement idal, lintgrale de configuration des fluides ne peut tre value exactement.
II-12
Du fait que la mcanique statistique ne fournie pas de solution exacte pour les fluides rels,
plusieurs mthodes dapproximations et des modles dinteraction sont employs dont la
thorie des tats correspondants, la thorie de perturbation et dautres mthodes
dapproximation de la fonction de partition comme la thorie de Van der Waals gnralise. Ils
permettent alors de driver des quations dtats spcifiques, par exemple lquation de Van der
Waals.
B.
Introduction
Introduction
Une fois que la modlisation relative un systme physique donn a t choisie, la deuxime tape
du travail consiste dterminer les proprits statistiques du modle en effectuant une simulation. Si l'on
s'intresse aux proprits statiques du modle, nous avons vu au sous-chapitre prcdent que le calcul de
la fonction de partition du systme se ramne au calcul d'une intgrale ou d'une somme discrte de
configurations de la forme :
Z = exp( U i )
(II-38)
II-13
Lutilisation de la mthode Monte Carlo pour modliser les problmes de physique nous permet dtudier
des systmes complexes en gnrant alatoirement un grand nombre de configurations parmi linfinit des
configurations que peut occuper un systme.
b.
Les gnrateurs de nombres alatoires sont au cur de toutes les simulations de Monte Carlo, non
seulement pour gnrer de nouvelles configurations mais aussi pour dcider de lacceptation ou du rejet de
nimporte quel mouvement effectu. On les retrouve galement dans linitialisation des simulations de
Dynamique Molculaire prsentes aprs pour attribuer les vitesses initiales aux diffrentes molcules du
systme.
Le prototype d'un gnrateur de nombres alatoires doit fournir des nombres issus d'une distribution
uniforme. Cependant, pour que ces nombres soient vraiment alatoires, il est ncessaire de satisfaire une
infinit de critres : la moyenne, la variance, mais aussi tous les moments de la distribution doivent tre
ceux d'une distribution uniforme. De plus, les suites de nombres doivent tre sans corrlation entre elles.
Comme les nombres sont reprsents par un nombre fini d'octets en informatique, les gnrateurs sont
forcment priodiques. Un critre ncessaire mais non suffisant est d'avoir une priode trs leve. Dans
les premiers temps de l'informatique, les gnrateurs utilisaient une reprsentation des nombres sur 8 bits,
impliquant des priodes trs courtes et des rsultats systmatiquement faux pour les simulations Monte
Carlo.
Cest pourquoi, lors des simulations ncessitant un grand nombre de tirages, il est indispensable de
s'assurer que la priode du gnrateur reste trs suprieure au nombre de tirages. Mais d'autres qualits
sont ncessaires, comme l'absence de corrlation entre les squences de nombres ; l'initialisation correcte
du gnrateur reste un point encore trop souvent nglig.
Il existe deux grands types d'algorithme pour obtenir des gnrateurs de nombres alatoires. Le premier
type est bas sur la congruence linaire :
x n +1 = (ax n + c ) mod(m )
(II-39)
mod tant la fonction modulo. Ce type de relation gnre une suite pseudo-alatoire de nombres
entiers compris entre 0 et (m-1). m donne la priode du gnrateur. Par consquent, il suffit de diviser xn+1
par m pour obtenir des nombres alatoires compris entre 0 et 1, puisque xn+1 ne peut excder la valeur de
m. Parmi les gnrateurs utilisant cette relation, on trouve les fonctions randu d'IBM, ranf du Cray,
drand48 sur les machines Unix, ran de Numerical Recipes (Knuth), etc. Les priodes de ces suites de
nombres vont de 229 (randu IBM) 248 (ranf).
La deuxime classe de gnrateurs est base sur le dplacement de registre travers l'opration ou
exclusif . Un exemple est donn par le gnrateur de Kirkpatrick et Stoll.
II-14
x n = x n 103 x n 250
(II-40)
Sa priode est grande, 2250, mais il ncessite de stocker 250 mots. Le gnrateur avec la priode la plus
grande est sans doute celui de Matsumoto et Nishimura connu sous le nom MT19937 (Mersenne Twister
generator). Sa priode est de 106000 ! Il utilise 624 mots par gnrateur et il est quidistribu dans 623
dimensions !
c.
Exemple illustratif
Pour illustrer lapplication dune mthode dchantillonnage base sur des nombres alatoires dans
le but de calculer une proprit macroscopique, prenons lexemple du calcul de la valeur de .
La Figure II-1.a reprsente un cercle unitaire inscrit sur un carr. On pourrait examiner ce problme en
terme du cercle complet et du carr, mais il est plus facile de traiter le problme en considrant
uniquement le quart, comme le montre la Figure II-1.b.
y
1
1
0
(II-41)
Lapplication dune technique de Monte Carlo faisant appel un gnrateur de nombres alatoires pour
cet exemple consiste la gnration de 2N nombres alatoires compris entre 0 et 1 lintrieur dune
distribution uniforme qui seront affects au couple de coordonnes (x,y).
II-15
Cet exemple est facile comprendre ; en simulation molculaire, par contre, nous sommes face des
problmes dvaluation de moyennes et dintgrales beaucoup plus complexes pour un nombre trs lev
de configurations. Cest pourquoi lapplication des mthodes de Monte Carlo pour la simulation
molculaire a fait lobjet dun travail particulier aboutissant des techniques spcifiques (lalgorithme de
Mtropolis, les biais statistiques, ) dont nous prsentons, ci-dessous, les principes thoriques.
d.
Echantillonnage uniforme
Tout d'abord considrons l'exemple le plus simple, celui d'une intgrale unidimensionnelle :
b
I = f (x )dx
(II-42)
I = (b a ) f (x )
O
(II-43)
(b a )
N
f (x )
(II-44)
i =1
Echantillonnage
uniforme
tat
b
Echantillonnage
alatoire
tat
tat
Echantillonnage
prfrentiel
f ( x) = 1 x 2
1/ 2
(II-45)
Une estimation de la valeur de par cette fonction donne la valeur de 3,14169, pour 107 tentatives, soit
quatre chiffres significatifs.
II-16
(II-46)
Pour approcher cette intgrale par cette mthode laide dune distribution uniforme, on procde selon les
tapes suivantes :
-
Pris par triplet, ces nombres spcifient les coordonnes des N molcules lintrieur du cube.
On calcule lnergie du systme grce un champ de force qui est fonction des positions des
diffrentes particules dans la configuration obtenue.
Cette procdure est rpte Nmax fois (Nmax >> 106), ce qui permet destimer lintgrale de
configuration par la relation :
Z NVT =
V
N max
N max
exp( U ( n ))
(II-47)
n =1
Hlas, cette valuation de lintgrale de configuration nest pas envisageable pour le calcul des proprits
thermodynamiques dun systme moins idal vu le trs grand nombre de configurations dun systme qui
ont une probabilit dexistence ou encore un facteur de Boltzmann quasi nul cause dune nergie trop
leve du fait du recouvrement nergtique entre les particules. Ces configurations peu probables sont
innombrables et donc longues calculer et pourtant ninterviennent pas ou peu dans le calcul des
grandeurs moyennes.
La convergence de cette mthode peut tre estime en calculant la variance 2 de la somme I. Or, on a :
2 =
1
N2
( f (x ) f (x ) ) ( f (x ) f (x )
N
i =1 j =1
(II-48)
Les points tant choisis indpendamment, les termes croiss qui qualifient la corrlation s'annulent, et on
obtient
2 =
1
N
( f (x )
N
i =1
f (x i )
(II-49)
La dpendance en 1/N donne une convergence a priori assez lente, mais il n'y a pas de modifications
simples pour obtenir une convergence plus rapide. On peut l'inverse modifier de faon importante
l'cart type. Pour cela, il parat clair que si la fonction f(x) ne prend des valeurs significatives que sur des
petites rgions de lintervalle [a,b], il est inutile de calculer la fonction en des points o sa valeur est trs
faible.
II-17
e.
Si par contre, on travaille avec une distribution alatoire non uniforme (Figure II-2) avec un poids
I =
a
f (x )
(x )dx
(x )
(II-50)
I =
a
f (x (u ))
du
(x (u ))
(b a )
N
(II-51)
f (x (u i ))
(x (u ))
i =1
(II-52)
Avec la distribution de poids (x), de manire similaire, la variance de l'estimation devient alors
1
=
N
2
i =1
f (x (u i ))
(x (u i ))
f (x (u i ))
(x (u i ))
(II-53)
En choisissant la fonction de poids (x) proportionnelle la fonction f(x), la variance s'annule. Cette
merveilleuse astuce n'est possible que dans un problme une dimension. En dimension suprieure, le
changement de variables dans une intgrale multiple fait intervenir la valeur absolue d'un jacobien et on ne
peut donc pas trouver de manire intuitive le changement de variable effectuer pour obtenir une
fonction de poids satisfaisante.
Pour des intgrations unidimensionnelles simples, la technique de Monte Carlo nest pas comptitive avec
les mthodes numriques de type Simpson par exemple bien connue des mathmaticiens pour calculer des
intgrales simples. Cependant, pour les intgrales de thermodynamique statistique, de dimension 6N pour
un systme N particules, la mthode sample mean integration avec un choix appropri de la fonction
de distribution de densit (x) est la seule fournir une solution acceptable.
Une solution au problme dvaluation de lintgrale de configuration ZNVT, introduite dans le paragraphe
prcdent, consiste ne gnrer que les configurations les plus stables et qui participent donc le plus
lvaluation des grandeurs thermodynamiques moyennes. Cest la stratgie adopte par la mthode
chantillonnage prfrentiel ( importance simpling ) introduite par Metropolis [ Metropolis et al., 1953 ]
(Figure II-2) et qui est employe pour la prdiction des quilibres liquide - vapeur.
II-18
f.
Reprenons notre problme de mcanique statistique : nous sommes intresss le plus souvent par le
calcul de la moyenne d'une grandeur A et non pas directement par la fonction de partition. Cette moyenne
scrit dans le cas de lensemble canonique sous la forme :
A =
A
i
exp( U i )
(II-54)
Z NVT
pi =
exp( U i )
Z NVT
(II-55)
p
i
= 1 ). Si l'on
tait capable de gnrer des configurations avec un mme poids peq, la moyenne de A serait estime par :
1
N
(II-56)
O N est le nombre de points calculs ; on serait alors ramen au calcul de la section prcdente.
L'astuce imagine par Metropolis, Rosenbluth et Teller en 1953 consiste gnrer une dynamique
stochastique Markovienne stationnaire, entre configurations successives, qui converge vers la distribution
d'quilibre peq.
Avant d'expliciter ce point, nous allons introduire quelques dfinitions. Considrant l'ensemble des
configurations i, on introduit un temps t prenant les valeurs discrtes associes au comptage des itrations
dans la simulation. Ce temps n'a pas de relation directe avec le temps rel du systme. On appelle p (i, t )
la probabilit du systme d'tre dans la configuration i au temps t.
Reprenons maintenant les termes de la dynamique choisie : dynamique stochastique signifie que le passage
d'une configuration une autre est le choix d'une procdure alatoire et Markovien signifie que la
probabilit d'aller vers une configuration j l'instant t+1, sachant que le systme tait dans la configuration
i l'instant t, ne dpend pas des configurations du systme pour des instants antrieurs (mmoire limite
l'instant t); cette probabilit conditionnelle est note w(i j ) . L'quation d'volution du systme est
alors donne par l'quation matresse suivante :
p(i, t + 1) = p(i, t ) +
(II-57)
II-19
Cette quation traduit le bilan suivant : l'instant t+1, la probabilit du systme d'tre dans l'tat i est gale
celle de l'instant prcdent, augmente par la possibilit que le systme qui se trouve dans n'importe
quelle autre configuration puisse aller dans l'tat i et diminue par la possibilit que le systme qui se
trouvait dans l'tat i puisse aller vers n'importe quelle autre configuration.
A linstant initial t = 0, le systme est plac dans une configuration initiale i0, ce qui signifie que, quel que
soit le choix de cette configuration, le systme ne satisfait pas P(i ) = i qui est la condition recherche.
w( j i )
j
= i
w(i j )
(II-58)
w( j i ) j = w(i j )i
micro-rversibilit
(II-59)
Cette relation, (II-59), est connue sous le nom de micro-rversibilit ou de bilan dtaill. Elle exprime le
fait que, dans l'tat stationnaire, (ou tat d'quilibre si le processus n'a pas engendr une brisure
d'ergodicit), la probabilit que le systme puisse aller d'un tat d'quilibre i vers un tat j est la mme que
celle d'aller d'un tat d'quilibre j vers un tat i. Ajoutons que cette condition n'est qu'une condition
suffisante, car nous n'avons pas prouv simultanment que la solution du systme d'quations (II-58) est
unique et que l'quation (II-59) est la meilleure solution. Pour des raisons pratiques, la quasi-totalit des
algorithmes de Monte Carlo reposent sur cette solution.
L'quation (II-59) peut tre aisment rcrite sous la forme :
w(i j ) j
=
= exp U j U i
w( j i ) i
( (
))
(II-60)
Cela implique que les inconnues w(i j ) que l'on cherche dterminer ne dpendent pas de lintgrale
de configuration ZNVT si difficile valuer, mais uniquement du facteur de Boltzmann reli lnergie de
chaque tat qui peut tre calcule.
g.
Algorithme de Metropolis
Le choix du processus Markovien stationnaire pour satisfaire l'ensemble des quations (II-59) est
l'une des solutions. Afin d'obtenir des solutions des quations (II-59) ou, en d'autres termes d'obtenir la
matrice de transition (w(i j )) , notons que le processus stochastique lmentaire dans un algorithme de
Monte Carlo est la succession de deux tapes :
II-20
A partir d'une configuration i, on tire au hasard une configuration j, avec une probabilit a
priori dessayer daller de ltat i ltat j (i j ) .
Cette nouvelle configuration est accepte avec une probabilit dacceptation (i j ) . Ainsi,
on a
w(i j ) = (i j ) (i j )
(II-61)
Dans l'algorithme original de Metropolis (et dans la plupart des algorithmes Monte Carlo), on choisit une
fonction symmtrique (i j ) = ( j i ) . Dans ce cas, les quations (II-60) sexpriment comme :
(i j )
= exp( (U j U i ))
( j i)
(II-62)
(i j ) = 1
si
U j Ui
(i j ) = exp( (U j U i ))
si
U j > Ui
(II-63)
Ainsi, dans la mthode de Metropolis la nouvelle configuration j est toujours accepte si son nergie est
plus faible que celle de ltat original i. Sinon, cest dire si lnergie de cette nouvelle configuration est
plus grande, le dplacement est accept suivant lquation (II-63) en comparant le facteur de Boltzmann
(II-64)
exp( U )
Rejeter
Toujours
accepter
Accepter
2
0
U (i j )
II-21
La Figure II-3 illustre le critre dacceptation dune transition de ltat i ltat j par la mthode de
Metropolis o i est un nombre alatoire compris entre 0 et 1.
Cependant, le calcul d'une moyenne ne peut commencer que lorsque le systme a atteint l'quilibre, c'est
dire quand p peq. Ainsi dans une simulation de Monte Carlo, il y a gnralement deux priodes : la
premire, o partant d'une configuration initiale, on ralise une dynamique afin d'amener le systme prs
de l'quilibre ; la seconde priode, o le systme volue au voisinage de l'quilibre, et o le calcul des
moyennes est ralis. En l'absence de critre prcis, la dure de la premire priode n'est pas facilement
prvisible. Une des techniques consisterait suivre l'volution de lnergie instantane du systme et
considrer que l'quilibre est atteint lorsque l'nergie se stabilise autour d'une valeur quasi-stationnaire.
3.
Dynamique molculaire
a.
Introduction
La limitation intrinsque de la simulation de Monte Carlo provient du fait qu'on ne considre pas la
dynamique ``relle'' du systme. Pour des systmes continus dfinis partir d'un Hamiltonien classique, il
est possible de rsoudre les quations du mouvement de la mcanique newtonienne pour l'ensemble des
particules. Cette mthode offre le moyen de calculer prcisment partir de corrlations temporelles les
proprits dynamiques du systme l'quilibre, grandeurs qui sont accessibles exprimentalement par
diffusion de la lumire ou des neutrons. Cette mthode permet aussi de calculer partir de corrlations
spatiales les grandeurs statiques d'quilibre comme dans une simulation Monte Carlo afin d'tre compares
directement l'exprience.
Pour des systmes atomiques en dehors de la rgion critique, cela est suffisant, mais a se rvle limitant
pour des difices molculaires plus complexes (molcules biologiques). En se basant sur un potentiel de
Lennard-Jones pour dcrire les interactions nergtiques, le temps obtenu partir d'une analyse
dimensionnelle des paramtres microscopiques du systme est :
(II-65)
m tant la masse d'une particule, son diamtre, et l'chelle d'nergie du potentiel d'interaction.
Le temps reprsente le temps pour un atome de se dplacer sur une distance gale sa taille avec une
vitesse gale la vitesse moyenne dans le fluide. Par exemple, pour l'argon, on a =3 , m=6.6310-23 kg
et =1.6410-20J, ce qui donne =2.810-14s. Pour une intgration numrique des quations du
mouvement, le pas d'intgration du temps doit rester une fraction du temps , typiquement t = 10-15 s,
voire plus petite. Le nombre de pas que l'on peut raliser en dynamique est typiquement de l'ordre de 105
107, ce qui nous amne pouvoir suivre un phnomne au maximum sur un intervalle de temps qui va
II-22
jusqu' 10-8 s. Pour de nombreux systmes atomiques, les temps de relaxation des phnomnes sont trs
infrieurs 10-8 s et la Dynamique Molculaire est un bon outil dtude.
b.
quations du mouvement
mouvement. La valeur instantane X(t) de Xobs peut tre exprime en fonction de ces 6N coordonnes. Le
principe de la Dynamique Molculaire est de rsoudre numriquement les quations de mouvement du
systme tudi. On obtient alors une trajectoire dans lespace des phases. Pour chaque point de cette
trajectoire, il est possible de calculer la valeur de X(t) et la grandeur Xobs mesurable exprimentalement
peut tre assimile la moyenne temporelle de X(t) :
X obs = X (t ) = lim
X (t )dt
(II-1)
= mi
dv i
dt
(II-66)
Les forces fi agissant sur la particule de masse mi sont gales la drive du potentiel U(r) agissant sur la
particule. Dans la direction x :
f x (r ) =
U ( r )
x dU (r )
=
x
r dr
(II-67)
12 6
, la
r
12
6
48 x
1
=
r
2 r
r2
Pour simuler un milieu infini, on utilise dans une simulation des conditions aux limites priodiques. Le
calcul de la force interagissant entre deux particules i et j se fait souvent entre l'image de j la plus proche de
i (voir le paragraphe sur les techniques de calcul).
II-23
Comme dans le cas dune simulation de Monte Carlo, le calcul de la force agissant sur la particule i
ncessite a priori le calcul de (N-1) forces lmentaires provenant des autres particules. L'utilisation d'un
potentiel tronqu dans lespace permet alors de limiter le calcul de la force aux particules entourant la
particule i dans une sphre dont le rayon est celui de la troncature du potentiel.
c.
Pour intgrer numriquement des quations diffrentielles, il est ncessaire de les discrtiser en
temps. Une grande varit de choix est possible a priori, mais comme nous allons le voir par la suite, il est
important que l'nergie du systme soit conserve au cours du temps (l'ensemble statistique est ici
lensemble microcanonique, NVE). L'algorithme propos par Verlet est historiquement l'un des premiers
introduit et il reste encore l'un des plus utiliss actuellement.
Pour des raisons de simplicit, nous considrons un systme constitu de N particules identiques et nous
appelons r, un vecteur 3N composantes: r = (r1, r2, , rN), o ri dsigne le vecteur position de la
particule i. L'quation d'volution du systme peut s'crire formellement comme :
d 2r
= f (r (t ))
dt 2
(II-68)
r (t + t ) = r (t ) + v(t )t +
3
f (r (t ))
(t )2 + d 3r (t )3 + O (t )4
2m
dt
(II-69)
3
f (r (t ))
(t )2 d 3r (t )3 + O (t )4
2m
dt
(II-70)
et de manire similaire,
r (t t ) = r (t ) v(t )t +
f (r (t ))
(t )2 + O (t )4
2m
(II-71)
Le calcul de la nouvelle position est donc effectu avec une prcision de l'ordre de (t)4. Cet algorithme a
lavantage de ne pas utiliser les vitesses des particules pour calculer les nouvelles positions, vitesses qui
requiert de calculer la drive des trajectoires r(t). On peut toutefois dterminer les vitesses a posteriori de la
manire suivante :
v(t ) =
r (t + t ) r (t t )
2
+ O (t )
2 t
(II-72)
La qualit d'une simulation de Dynamique Molculaire est lie la qualit de l'algorithme utilis et ses
proprits. La rapidit de l'excution du programme peut aussi tre dterminante. Notons que l'essentiel
II-24
du temps de calcul dans une Dynamique Molculaire est consomm dans le calcul des forces, ce qui
signifie que le cot du calcul des nouvelles positions est marginal.
La prcision du calcul pour l'algorithme de Verlet est approximativement donne par
t 4 N t
(II-73)
O Nt est le nombre de pas dintgration de la simulation. Le temps maximal coul dans la simulation est
donn par t N t . Il pourrait sembler intressant d'utiliser un algorithme faisant intervenir des drives
des coordonnes des ordres plus levs mais dans ce cas, la quantit d'information garder en mmoire
augmente rapidement. Ainsi, lalgorithme de Verlet est particulirement compact, ne ncessitant que 9
vecteurs de N lments (r x, r y, r z, r& x, r& y, r& z, et les nouvelles forces Fx, Fy, Fz), l o un algorithme type
Gear bas sur des quations analogues lquation (II-69) en ncessite 15 (r x, r y, r z, r& x, r& y, r& z, &r& x, &r& y, &r& z,
&r&& x, &r&& y, &r&& z et les nouvelles acclrations (3N) et nouvelles forces (3N)).
Les algorithmes d'ordres plus levs, comme celui de Gear, ont tendance fournir une dynamique aux
temps courts de meilleure qualit, mais l'nergie totale du systme tend ne pas rester constante aux
temps longs. L'algorithme de Verlet possde au contraire la vertu de conduire une faible drive de
lnergie totale aux temps longs. Or tudier le plus longtemps la dynamique du systme est prcisment ce
que lon cherche avec la Dynamique Molculaire.
Une symtrie particulirement importante, contenue dans les quations de Newton du systme, est la
symtrie par renversement du temps. Il est important de noter que l'algorithme de Verlet satisfait cette
symtrie. En effet, en changeant t t , l'quation (II-71) reste inchange. La consquence de cette
proprit est que si, un instant t de la simulation, on inverse le cours du temps, la trajectoire de la
Dynamique Molculaire revient sur ses pas. En pratique, les erreurs d'arrondi accumuls dans la simulation
limitent la rversibilit quand le nombre de pas dintgration devient important. En revanche, on peut, en
utilisant cette proprit de l'algorithme, tester l'importance des erreurs d'arrondi, en inversant le temps
dans la simulation pour des temps de plus en plus grands.
Signalons un algorithme, apparent celui de Verlet, connu sous le nom d'algorithme Leap frog
(algorithme saute-mouton) bas sur le principe suivant : les vitesses sont calcules pour des intervalles de
temps demi-entiers et les positions sont obtenues pour les intervalles de temps entiers. Si on dfinit les
vitesses pour les temps t + t 2 et t t 2 :
v(t + t 2) =
r (t + t ) r (t )
t
(II-74)
v(t t 2 ) =
r (t ) r (t t )
t
(II-75)
II-25
On obtient :
r (t + t ) = r (t ) + v(t + t 2 ) t
(II-76)
r (t t ) = r (t ) v(t + t 2 ) t
(II-77)
v(t + t 2) = v(t t 2 ) +
f (t )
t + O (t )3
2m
(II-78)
Puisqu'il est aussi bas sur l'quation (II-71), cet algorithme est identique l'algorithme de Verlet en ce qui
concerne le calcul des trajectoires (les valeurs des vitesses aux temps demi-entiers n'apparaissent que
comme des intermdiaires de calcul). Il peut tre toutefois diffrent pour le calcul des grandeurs
thermodynamiques car la moyenne de l'nergie potentielle peut tre calcule aux temps entiers (elle fait
intervenir les positions), tandis que la moyenne de l'nergie cintique est calcule aux temps demi-entiers
(elle fait intervenir les vitesses).
Linitialisation de la simulation en Dynamique Molculaire se fait en fixant la position de dpart des
molcules et en leur donnant une vitesse initiale. Les quations de Newton impliquent ensuite une
conservation de lnergie au cours du temps. C'est donc naturellement dans lensemble micocanonique
(c'est--dire nombre de molcule N, volume V et nergie E constants) que seffectue une simulation
en Dynamique Molculaire.
d.
Thermostats T et P
Toutefois, lensemble microcanonique nest pas toujours reprsentatif des systmes que lon
souhaite modliser. On peut avoir besoin de travailler non plus volume constant mais pression
constante. De plus exprimentalement, les systmes sont plutt en contact avec un thermostat et donc
temprature fixe plutt qu nergie fixe.
En effet, si lon part dun systme qui nest pas lquilibre et que lon fait une simulation dans lensemble
microcanonique, la stabilisation du systme va saccompagner dun changement important de la
temprature. Une faon simple de rguler la temprature du systme est de multiplier, aprs chaque pas de
calcul, les vitesses par un pas correctif
Tc (t ) =
i =1
mi vi (t )
kB N f
2
(II-79)
Nf tant le nombre de degr de libert, gal (3N-3) pour un systme de N particules avec une nergie
totale constante.
II-26
Cependant, il existe diffrentes mthodes qui permettent de travailler dans dautres ensembles que
lensemble microcanonique, mais lune des difficults de la Dynamique Molculaire dans un ensemble
statistique temprature ou pression constante est de maintenir T et ou P constant en moyenne, mme
si les valeurs instantanes de la pression et de la temprature peuvent fluctuer. En effet, dans un ensemble
canonique NVT, comment imposer une temprature prcise un systme dont les particules sagitent et
crent de ce fait un chauffement ? La solution trouve est de considrer que le systme est plong dans
un large bain thermostat.
Le thermostat dAndersen consiste choisir priodiquement une particule du systme choisie au hasard et
lui attribuer une vlocit rpondant une loi de distribution de Boltzmann. Cela revient considrer que
priodiquement, une particule du systme choisie au hasard entre en collision avec une particule
imaginaire du bain thermostat. Les inconvnients de cette mthode sont que lnergie totale nest plus
conserve et que la fonction dauto corrlation de la vitesse qui sert calculer le coefficient de diffusion
nest plus exacte puisque certaines vitesses instantanes sont modifies.
Le thermostat de Nos-Hoover propose un autre formalisme dans lequel les simulations sont ralises sur
la runion de lensemble statistique et du bain thermostat. Lnergie schange librement entre lensemble
statistique et le bain qui dispose dune certaine inertie thermique permettant de contrler la temprature au
sein du systme modlis. Une mthode similaire permet dimposer la pression du systme en couplant le
systme modle avec un large volume tampon.
4.
Ergodicit
Dans les diffrents ensembles statistiques que nous avons prsents au cours de la partie A, les
moyennes ont un sens purement statique, et sont dfinies sur l'ensemble des tats accessibles au systme.
Ceci ne correspond pas la moyenne que l'on effectue dans une exprience, dans une simulation de
Dynamique Molculaire, ou dans une simulation de Monte Carlo. Dans ces trois situations, on ralise une
srie de mesures ou de simulations durant des intervalles de temps finis ou sur un nombre de
configurations finis et on dtermine la moyenne sur ces mesures ou simulations.
Pour tre plus concret, considrons, dans un liquide simple, la densit locale moyenne obtenue par
Dynamique Molculaire. La densit instantane dpend des coordonnes r de toutes les particules du
systme. Lorsque le temps volue, les coordonnes atomiques changent (selon lquation de mouvement
de Newton), et donc la densit autour dun atome donn va changer aussi.
1
(r ) =
N
1
lim
dt (r, t )
t (t t )
i t
i =1
(II-80)
O i est un indice qui parcourt un ensemble de conditions initiales et ti, le temps initial correspondant.
Dans la mesure o l'ensemble des conditions initiales contient des lments compatibles avec l'ensemble
statistique appropri (dans le cas de la Dynamique Molculaire, il s'agit de l'ensemble microcanonique
NVE), on obtient :
II-27
1
dt (r, t )
t (t t )
1 t
(r ) = lim
(II-81)
NVE
Comme la moyenne d'ensemble ne dpend pas de l'instant initial, il y a une quivalence entre la moyenne
sur l'instant initial et la moyenne temporelle dfinie par lquation (II-3) de la partie A.
(r ) = (r )
ergodicit
NVE
(II-82)
Cette dernire quation implique que si le systme est capable d'explorer, durant le temps de l'exprience,
l'espace des phases de manire reproduire un ensemble reprsentatif de conditions initiales de la
moyenne d'ensemble, on a quivalence entre moyenne temporelle (obtenue par Dynamique Molculaire)
et moyenne d'ensemble (obtenu par les mthodes Monte Carlo). Cette galit est connue sous le nom du
thorme dergodicit . En pratique, conjointement sur le plan exprimental et sur le plan de la
simulation numrique (dans ce dernier cas, la dynamique est une dynamique fictive), l'mergence d'tats
mtastables conduit une brisure d'ergodicit, et il devient trs difficile, voire impossible d'chantillonner
correctement l'espace des phases, cest le cas aussi pour les systme solides qui ne sont pas ergodiques par
principe.
5.
r
r
dv
F = m dt
Configurations - trajectoire
t
?
Dynamique Molculaire
?
Monte Carlo
Configurations alatoires
II-28
transports) souvent values partir du dplacement moyen des molcules. La mthode de Monte Carlo
est moins coteuse en temps de calcul car elle ne ncessite pas le calcul des forces qui agissent sur les
molcules et reste particulirement adapte ltude des systmes en quilibre (phase isole ou quilibre de
phase). On choisira donc lune ou lautre des deux mthodes selon la nature des informations que lon
souhaite obtenir, et la nature des expriences avec lesquelles on comparera les rsultats des simulations.
Cest aussi lensemble statistique dtude qui peut guider le choix de la mthode dchantillonnage. La
simulation numrique dans un ensemble statistique donn se prte souvent mieux lune des deux
mthodes. Cest le cas par exemple des ensembles dans lesquels le nombre de molcules N fluctue,
comme lensemble grand canonique, qui sapplique mieux un algorithme de Monte Carlo qu une
simulation de Dynamique Molculaire.
6.
a.
La formule usuelle pour calculer une grandeur moyenne sur un ensemble de total configurations :
ensemble
total
X ( )
total = 1
1
moyenne
(II-83)
Un grand nombre dautres grandeurs peuvent tre calcules partir de la variance qui exprime les
fluctuations autour de la valeur moyenne :
total
1
2
2
2
=
X ( ) X ensemble = X 2 X
X
total = 1
variance
(II-84)
Les fluctuations dpendent de lensemble dans lequel elles sont calcules. En effet, la fluctuation de
lnergie na pas de sens dans un ensemble NVE. Ainsi, la capacit calorifique CV est lie E dans
lensemble NVT ; la compressibilit isotherme est lie N dans lensemble VT et V dans
lensemble NPT. On peut aussi calculer des fonctions structurales comme la fonction de distribution g(r))
indicatrice de la structure du systme et lie la probabilit de trouver une paire datomes spars par une
distance r.
Les coefficients de corrlations permettent daccder des grandeurs dcrivant ltat dynamique du
systme. Pour une proprit X (ex. la vlocit), la forme non normalise du coefficient de corrlation est la
suivante :
correl
XX
( ) = < X ( ) X ( ) > =
0
total
X ( ) X (
0
+ )
corrlation
(II-85)
total =1
II-29
Lintgration des coefficients de corrlations non normaliss permet de calculer les coefficients de
transfert macroscopique (coefficient de diffusion, viscosit, conductivit thermique). Leurs transformes
de Fourrier peuvent tre compares des spectres exprimentaux. La deuxime galit de lquation (II85) illustre comment on calcule en pratique le coefficient de corrlation pour un ensemble de total
configurations. Les formules sont dtailles dans les livres cits en tte de ce chapitre.
b.
Erreur statistique
La simulation molculaire est une exprience numrique. Par consquent, les rsultats sont sujets
des erreurs systmatiques et statistiques. Les erreurs systmatiques doivent tre values puis limines.
Elles ont pour causes des effets de taille, une mauvaise gnration alatoire, un mauvais quilibrage (voir
plus loin). Les erreurs statistiques interviennent puisquune majorit des proprits calcules sont des
moyennes sur le nombre de configurations certes important mais fini.
Dans lhypothse o une grandeur suit une loi de Gauss, lerreur statistique associe la valeur moyenne
est lcart type, racine de la variance, multipli par un facteur dpendant de l'indice de confiance attendu.
Cependant, lchantillonnage grand mais fini des configurations induit une corrlation entre les total
configurations stockes qui persiste pendant un certain nombre de configurations successives. Une faon
de tenir compte de ce problme est dintroduire un facteur dinefficacit statistique s qui value le nombre
de configurations successives qui sont corrles. Pour cela, on dcoupe les total configurations en nb blocs
de b configurations : total = nb.b.
On calcule sur lensemble des total configurations la moyenne de la proprit Xtotal (quation II-86) et sa
variance 2(X)total (quation II-87). Sur chaque bloc sont calcules la moyenne de la proprit Xb et sa
variance 2(Xb). En slectionnant plusieurs valeurs de b de plus en plus grandes, on peut calculer
linefficacit statistique s :
2 X b
b
s = lim
2 ( X )total
b
(II-86)
Alors, lerreur statistique 2( Xtotal ) associe la moyenne Xtotal sur lensemble des total configurations
est calcule :
X total =
s
total
c.
(x )total
(II-87)
Les cellules, ou boites, utilises dans des simulations de Monte Carlo ou de Dynamique Molculaire
sont de taille finie. Pour limiter les effets de bords, on applique les conditions priodiques aux limites.
Chaque cellule est rplique dans les trois dimensions, en gardant la mme position relative pour chaque
II-30
particule (Figure II-5). Au cours de la simulation, lorsquune particule se dplace dans la boite originale,
ses images priodiques dans chaque boite avoisinante se dplace exactement de la mme manire. Ainsi, si
lun des mouvements aura pour effet de dplacer lune des particules en dehors de la boite centrale, celuici aura pour effet secondaire de dplacer lune de ses images lintrieure de la cellule centrale ; le nombre
totale de particules reste ainsi inchang. Il nest pas ncessaire, non plus, de stocker les coordonnes de
toutes les particules images dans une simulation (un nombre infini !), mais uniquement celles de la boite
centrale.
surface ex. cristal
image
systme simul
Lnergie de la particule noire de la Figure II-5 doit tre calcule partir des interactions avec
toutes les particules, y compris les particules images.
Le systme priodique sera reprsentatif de la phase modle si la taille de la boite centrale est
suffisante et que les effets de bords sont attnus.
La contribution nergtique des interactions avec des particules images sera dautant plus importante que
le potentiel a une porte longue. Cest le cas pour les interactions coulombiennes et diffrentes techniques
existent, simples comme un rayon de coupure ; distance au-del de laquelle linteraction est suppose
nulle ; ou complexes (sommation dEwald). Dans le cas dun rayon de coupure, il faut inclure des
corrections longue distance.
Dans lexemple de la Figure II-5 o lon a fait volontairement figur une interface cristalline, les effets de
bord peuvent tre importants : linteraction de la particule noire avec la surface hachure est la somme sur
toutes les boites images de linteraction avec la surface par-dessous mais aussi par-dessus ! En outre, il
faudra prendre en considration le fait que larrangement des particules au voisinage de la surface sera
certainement diffrent de celui au sein du fluide homogne.
II-31
d.
Les interactions entre centres de force ne sont pas calcules pour toute la boite de simulation. Dans
la plupart des simulations molculaires, on impose un rayon de coupure rc pour le calcul explicite des
interactions ; ceci signifie que lorsquon calcule lnergie dinteraction dun centre de force i, on ne
considre que les particule dont la distance par rapport ce centre est infrieure rc, les interactions entre
centres de force situs une distance suprieure rc ne sont pas prises en compte. On introduit alors une
correction longue distance pour corriger cette approximation. Pour les interactions autres que les
interactions lectrostatiques et dinduction, on suppose quau-del de rc, les centres de force sont rpartis
uniformment, ce qui permet dvaluer cette correction par lexpression :
Ecorr = 2N r 2U (r )dr
(II-88)
rc
cran
i
exp( 2 r 2 )
= q i
( )
(II-89)
i (r ) = [ i (r ) + icran (r )] + [ icran (r )]
(II-90)
2 ,i (r ) = icran (r ) . La contribution lnergie lectrostatique nest pas calcule de la mme faon pour
ces deux termes. Le terme 1, i (r ) correspondant aux charges crantes converge rapidement dans
lespace direct, alors que la convergence de la contribution due 2 ,i (r ) et lente dans lespace direct mais
peut scrire comme une srie de Fourrier rapidement convergente dans lespace rciproque.
II-32
Plus de dtails techniques sur cette mthode et sa mise en uvre dans le code de calcul utilis dans ce
travail sont disponibles dans la thse de Bourasseau [ 2003 ] et le rapport de stage de Nicolas [ 2000 ] du
Laboratoire de Chimie Physique de lUniversit de Paris-sud.
e.
Le biais configurationnel
Le fait de simuler des molcules flexibles implique la prise en compte des interactions
intramolculaires qui requiert des mouvements de Monte Carlo permettant la relaxation interne des
molcules. Lorsque le mouvement consiste insrer une molcule dans une phase dense, la probabilit
dun tel mouvement ralis de faon alatoire est quasi nulle du fait du peu despace libre. Une parade
consiste introduire un biais statistique en considrant que la probabilit a priori dessayer daller de
ltat i ltat j (i j ) nest plus symtrique, mais biaise. Lquation (II-60) devient :
w(i j ) (i j ) j
=
w( j i ) ( j i ) i
(II-91)
(i j ) j
j
i
(
)
i
(i j ) = min 1,
(II-92)
Particulirement utile pour assurer une insertion de molcule dans une phase dense avec une bonne
probabilit, le biais configurationnel sera dtaill dans le paragraphe relatif la mthode de Monte Carlo
dans lensemble de Gibbs.
f.
Aprs avoir choisi un ensemble statistique, la configuration initiale consiste dfinir une boite, en
gnral cubique, dans laquelle sont places les particules, non pas alatoirement afin dviter les
chevauchements probables, mais selon un rseau priodique. Ensuite, on associe chaque particule les
paramtres du champ de force qui lui correspond pour pouvoir valuer lnergie du systme. Enfin, dans
le cas de la Dynamique Molculaire, on attribue (par exemple selon une distribution statistique de
Boltzmann) des vlocits initiales chaque particule afin de rsoudre lquation de Newton.
Puis la simulation est lance. Elle consiste en une phase dquilibrage et une phase de production.
La phase dquilibrage a pour but damener le systme dune configuration initiale une configuration
reprsentative du systme partir de laquelle on peut effectuer un chantillonnage respectant la condition
de micro-rversibilit : distribution alatoire des molcules et des vlocits au sein dun systme avec des
conditions thermodynamiques imposes (celle de lensemble statistique dans lequel sont ralises les
simulations). En Dynamique Molculaire dans un ensemble temprature T fixe, cela peut consister
chauffer graduellement ou par paliers le systme jusqu la temprature T puis attendre que lnergie
potentielle, que la pression, que la densit soient stabiliss, comme en Monte Carlo.
II-33
quilibrage
production
quilibrage
0,6
0,8
0,7
0,4
0,6
0,5
liquide 1
0,4
liquide 2
0,3
gaz
0,2
0,1
pression rduite
densit rduite
production
0,2
liquide 1
0
-0,2
liquide 2
0
500000
1000000
1500000
2000000
gaz
-0,4
0
0
500000
1000000
1500000
2000000
-0,6
nombre de configuration
nombre de configuration
0,8
fraction molaire
0,7
0,6
Ya gaz
Yb gaz
Xa liquide 1
Xb liquide 2
Xa liquide 2
Xb liquide 1
0,5
0,4
0,3
0,2
L1
L2
0,1
0
0
500000
nombre de configuration
Figure II-6 : Evolution de quelques grandeurs caractristiques lors de la simulation dun quilibre liquide
liquide vapeur pour un fluide modle de Lennard-Jones
La phase de production va consister enregistrer les caractristiques du systme pendant un grand
nombre de configurations afin de calculer in fine les proprits macroscopiques partir de moyennes, de
fluctuations ou de coefficients de corrlation. Typiquement en Dynamique Molculaire, une phase de
production dure quelques dizaines de nanosecondes, correspondant quelques dizaines de millions de pas
dintgration. En Monte Carlo, une phase de production correspond quelques dizaines de millions de
configurations.
La Figure II-6 prsente lvolution de la densit, de la pression et des fractions molaires lors de la
simulation dun quilibre liquide liquide vapeur dans lensemble de Gibbs (voir description aprs) pour
des particules modles de Lennard-Jones (les seules interactions prises en compte sont celles de Van der
Waals par un potentiel de Lennard-Jones). Le systme comprend trois boites de configuration initiale
identique (mme taille, nombre de molcules et positions). On constate que la pression lie au volume des
boites sest quilibre avant la densit et les fractions molaires. Les valeurs de la densit montrent bien
lexistence de deux boites denses correspondant aux liquides et dune boite de faible densit
correspondant la vapeur. Cette simulation est cependant inexploitable car on constate une drive dans
lvolution des fractions molaires qui peut tre li au choix de la configuration initiale, en loccurrence une
taille de boite et un nombre de molcules insuffisants , coupl une mauvaise valuation des corrections
pour les forces longue porte. Or un critre de succs dune simulation dont il faut sassurer est
II-34
C.
1.
La mcanique quantique
a.
Diffrentes approximations
Les principes de la mcanique quantique ne peuvent pas tre mis en uvre pour des systmes
reprsentatifs dune phase car le temps calcul est proportionnel au mieux Nlectron1,5 [ Leach, 1996 ].
Pourtant, leur utilit est relle en mcanique molculaire pour le calcul de lnergie dinteraction :
pour lobtention des charges atomiques ponctuelles, intervenant dans le calcul des interaction
colombiennes.
pour lobtention des diples et autres multiples permanents intervenant dans le calcul des
autres interactions lectrostatiques.
pour la rgression des constantes de rappels associes aux diffrents modes vibratoires des
molcules (longation, pliage, torsion).
Tout systme est considr sous la forme de noyaux autour desquels gravitent des lectrons. Les calculs de
mcanique quantique fournissent les proprits nuclaires et lectroniques du systme et lnergie totale
vraie du systme. Lnergie totale est relie la fonction donde gnrale (r,t) dpendante des
coordonnes et du temps au moyen de lquation de Schrdinger gnralise :
H = ih
d
dt
(II-93)
avec : h = h 2
O H est lHamiltonien, un oprateur mathmatique dcrivant les interactions du systme sous forme de
contributions cintiques et potentielles. (r,t) est interprte comme une amplitude de probabilit de
II-35
prsence et son carr|(r,t)|2 est donc la densit de probabilit de prsence de la particule, gale lunit
lorsque somme sur tout lespace.
Si lon suppose que (r,t) scrit comme le produit dune fonction (r) et f(t), on peut montrer que
lquation gnrale de Schrdinger se dcompose en deux quations, dont lquation de Schrdinger
indpendante du temps :
H ( r ) = E ( r )
(II-94)
C'est l'quation aux valeurs propres de l'oprateur H. Les nergies possibles sont donc les valeurs propres
de H et les conditions aux limites ne permettent E, lnergie totale du systme, de ne prendre que
certaines valeurs (quantification).
Cette quation a une solution analytique pour un systme un seul lectron, cest--dire pour le seul atome
dhydrogne H ou dhlium He+ ! Tout autre systme ncessite des techniques approximatives pour tre
tudi avec un compromis invitable entre le souhait dobtenir des rsultats rapidement et des rsultats
avec une grande prcision.
On distingue trois grands niveaux dapproximation selon que lon value les interactions entre les
lectrons plus ou moins finement qui conduisent aux mthodes suivantes dans lordre dapproximation
croissant :
semi-empirique
Parmi les mthodes ab-initio, les mthodes CI (interaction de configurations) sont celles qui obtiennent les
meilleurs rsultats mais sont extrmement lentes. Les valeurs de lnergie quelles calculent ont une
prcision comparable celles gnralement mesures exprimentalement. Les solutions des mthodes CI
sont obtenues en minimisant une combinaison linaire des fonctions dondes associes ltat
fondamental du systme et tous les tats excits. Cependant, leur application au-del de quelques
dizaines datomes devient impossible mme avec un super calculateur.
Pour des systmes de quelques centaines datomes, dautres approximations ont t dveloppes, parmi
lesquelles le modle appel Self Consistent Field Molecular Orbital (SCF-MO). Il utilise le concept
dorbitale atomique qui reprsente la fonction donde dun lectron se dplaant au sein dun potentiel
gnr par les noyaux et dun potentiel effectif moyen gnr par les autres lectrons du systme. Les
meilleures fonctions donde de ce type sont calcules sans aucun paramtre empirique. Elles sont appeles
des fonctions donde Hartree-Fock et rsolvent lquation de Schrdinger pour une configuration
lectronique donne (ex. ltat fondamental). Elles peuvent alors servir pour les calculs CI.
II-36
A lautre extrmit de lchelle des approximations des mthodes SCF, on trouve les modles semiempiriques dlectron . La thorie des orbitales molculaires de Huckel est la plus simple : les calculs
peuvent se faire sur une feuille de papier. Des modles semi-empiriques plus fins existent et permettent
dobtenir avec une bonne prcision des nergies dionisation, des gomtries optimales avec des valeurs
prcises des longueurs de liaisons et des angles de valence, des surfaces dnergie potentielle. Cependant,
elles prsentent le dsavantage de ne pas calculer directement les fonctions dondes mais essayent de
remplacer les diffrentes intgrales par des paramtres empiriques.
Entre les deux niveaux dapproximation, on trouve les mthodes du champ moyen dont la populaire DFT
thorie de la fonction de la densit.
b.
Lobjectif de toutes les mthodes de mcanique quantique numrique est de dterminer la structure
lectronique du systme. Pour cela, il leur faut rsoudre lquation de Schrdinger (II-94) dans laquelle
loprateur Hamiltonien H est gal (les termes noyau noyau sont omis) :
H =
1
2
i2
Z
+
ri R
j >i
12
ri r j
(II-95)
o Z est le nombre atomique et aussi la charge sur le noyau , ri donne la position de llectron par
rapport lorigine, R est la position du noyau par rapport lorigine, 2 est le laplacien. La sommation
sur i et j rfre aux lectrons et la sommation sur rfre aux noyaux. Le j > i vite de compter deux fois
les interactions entre les lectrons. Le terme associ au laplacien est celui de lnergie cintique du nuage
lectronique, les deux autres termes noyau lectron et lectron lectron tant lis lnergie potentielle.
La mthode Hartree Fock a t la premire proposer une faon de rsoudre lquation (II-94) avec
lHamiltonien de lquation (II-95). Pour cela, elle utilise lapproximation de Born-Oppenheimer qui
postule que les noyaux tant beaucoup plus lourds que les lectrons, ceux-ci sajustent instantanment
tout mouvement des noyaux ou inversement, que les noyaux paraissent figs. Ainsi, on peut sparer les
mouvements lectroniques et nuclaires. La fonction donde lectronique ne dpend que des positions des
noyaux et pas de leur nergie cintique.
De plus, chaque fonction donde multilectronique peut tre crite comme le produit de fonctions
donde mono lectroniques i. Cela implique que chaque lectron se dplace dans un champ moyen d
aux noyaux et aux autres lectrons. Cela caractrise ainsi les mthodes SCF.
Afin de satisfaire le principe dexclusion de Pauli, il faut que soit antisymtrique et on lcrit alors sous la
forme dun dterminant, =
N! 1 2
nergies associes i sont calcules en rsolvant N quations de Schrdinger associes chaque lectron :
II-37
quation de Fock
(II-96)
o VC est le potentiel de Coulomb et VXi est le terme dchange dinteraction discut aprs et qui prend
en compte tous les termes dchanges entre les lectrons et te linteraction entre un lectron et lui-mme.
Alors, un lectron est soumis un potentiel gnr par tous les autre lectrons et uniquement eux.
Les fonctions donde solutions de cette quation peuvent tre exprimes sous forme de fonctions
analytiques approximatives appeles les orbitales de Slater mais cest une approximation gaussienne
analytique de ces orbitales qui est le plus souvent employe. Selon le degr dapproximation gaussienne, on
obtient diffrents basis set. Ainsi la base STO-3G signifie quune combinaison de 3 fonctions
gaussiennes est utilise pour dcrire chaque orbitale atomique. Des noms comme 6-311++G(3df,3pd) ou
6-311++G** dcrivent des bases encore plus dtailles [ voir Leach, 1996 ].
c.
Les mthodes semi-empiriques drivent de la thorie de Huckel. La fonction donde i est calcule
comme une combinaison linaire dorbitales atomiques :
i = civ v
(II-97)
c (h
iv
i,effectif
i v = 0
(II-98)
c (h
iv
i S v ) = 0
(II-99)
O :
hv = hi,effectif v dv
(II-100)
Et :
S v = v dv
(II-101)
S est appel lintgrale doverlap (de chevauchement). On dtermine alors les nergies i en annulant le
dterminant de lquation (II-99) :
hv i S v = 0
(II-102)
Ensuite, chaque nergie i est substitue dans lquation (II-99) pour dterminer les coefficients ci
appropris utiliss pour reconstruire la fonction donde mono lectronique i par lquation (II-97).
II-38
Pour les mthodes semi-empiriques, les termes h et S sont remplacs par des paramtres empiriques
afin de faciliter la rsolution du dterminant de lquation (II-102). Elles donnent leur nom aux mthodes
empiriques : ZDO signifie Zero Differentiel Overlap (S =0), MINDO Modified Intermediate Neglect of
Differential Overlap, Parmi les plus populaires, citons AM1, MINDO/3, MNDO, PM3 et MM2.
Toutes ces mthodes sont disponibles au travers de trs nombreux logiciels de calculs de mcanique
quantique (ex. Hyperchem, GAUSSIAN, ).
d.
VXi 1 / 3
(II-103)
Analyse de population
Parmi les proprits que la mcanique quantique numrique permet de calculer, la densit
lectronique occupe une place de choix. Dcrivant la rpartition statistique des lectrons du systme, elle
se reprsente sous la forme de lignes disodensit. Elle permet aussi de raliser une analyse de population
qui vise partager la densit lectronique en charges ponctuelles places sur des centres atomiques rels
mais aussi parfois virtuels. Rappelons cependant, que la charge atomique na aucune ralit physique.
Pourtant bien utile dans les modles de mcanique molculaire, elle ne peut tre estime quen faisant des
approximations.
II-39
i.
Analyse de Mulliken
Pour construire un modle simple pour calculer les charges que nous voulons attribuer chaque
atome de la molcule, considrons une molcule diatomique AB. Les orbitales molculaires sont la
combinaison linaire de deux orbitales atomiques A et B centres sur A et B respectivement.
Considrons une orbitale molculaire i contenant ni lectrons :
i = c iA A + c iB B
(II-104)
(II-105)
SAB tant le recouvrement entre les deux orbitales atomiques. On peut considrer que ciA2 reprsente la
probabilit de trouver un lectron, situ dans i, localis sur latome A et ciB2 de le trouver localis sur
latome B. Le terme de recouvrement restant reprsente plus ou moins une probabilit de dlocalisation
sur A et B.
Lhypothse de Mulliken, est de considrer que lon peut partager cette partie restante par moiti sur
chacun des atomes A et B. Les probabilits de prsence de llectron sur chacun des atomes scrivent
alors :
(II-106)
Si on a ni lectrons sur lorbitale molculaire numro i, et si lon a N orbitales molculaires contenant des
lectrons, on peut calculer la population lectronique donc la charge lectronique se trouvant sur chacun
des atomes ; par exemple sur A :
q A = e
n p
i
(II-107)
iA
q A = e
n (c
2
iA
+ ciA ciB S AB
(II-108)
ou en charge nette :
QA = e Z A e
n (c
i
B A
2
iA
+ ciA ciB S AB
avec :
=0
(II-109)
Par dfinition, les charges de Mulliken dpendent naturellement de la base dorbitale atomique utilise, les
bases les plus tendues donnant parfois les valeurs les plus irralistes. Pourtant elle reste trs utilise du fait
de la simplicit de lanalyse qui la rend disponible dans tous les logiciels.
II-40
ii.
Cest une analyse visant sommer les degrs doccupation des NAO (Natural Atomic Orbitals). Les
orbitales molculaires sont dlocalises sur toute la molcule et nont gnralement aucune ressemblance
avec les liaisons covalentes localises ou si employes dans les raisonnements chimiques habituels.
Pourtant, le chevauchement entre des orbitales hybrides invent par L. Pauling pour dcrire les liaisons
localises a un fondement physique fort en la personne de la densit lectronique.
Lanalyse NBO consiste transformer les N orbitales atomiques en N Natural Atomic Orbitals, puis
combiner les NAO en NHO (Natural Hybrid Orbitals) de faon dcrire limplication des atomes dans la
densit lectronique de la molcule. Enfin, les NHO donnent les NBO en se recouvrant. Orbitales
molculaires ou NBO forment chacune un ensemble dorbitales solution de lquation de Schrdinger.
Une diffrence essentielle est que les orbitales molculaires possde 0, 1 ou 2 lectrons tandis que les
NBO peuvent possder un degr doccupation fractionnel entre 0 et 2 inclus. Les NBO fractionnelles
dcrivent gnralement des orbitales anti-liantes responsables des phnomnes de donneur , de
conjugaison, Cette similarit avec les raisonnements chimiques (liaisons covalentes, , conjugaison,
donneur, accepteur) sont lorigine de la popularit grandissante des orbitales NBO. Elles sont calcules
notamment par le logiciel GAUSSIAN.
Les charges NBO sont trs peu sensibles la base dorbitale atomique utilise, ce qui en fait une analyse
allchante pour comparer des travaux effectus avec des logiciels diffrents et des bases diffrentes.
In fine, les valeurs de charges ponctuelles varient normment dune mthode danalyse lautre. Toutes
sont des reprsentations approximatives de la densit lectronique. Par consquent, le choix dune
mthode est impossible a priori. Dautres critres doivent tre utiliss ; reproduction dune grandeur
mesurable exprimentalement comme le moment dipolaire ou bonne prdiction de proprits
thermodynamiques par simulation molculaire comme nous proposons de le faire.
II-41
2.
La mcanique molculaire
Le but de la mcanique molculaire est de dcrire les interactions nergtiques au sein dun systme
molculaire. Elle se diffrencie ainsi de la mcanique quantique o lnergie totale du systme satisfait
lquation de Schrdinger qui comporte un hamiltonien rassemblant lensemble des interactions (noyau
noyau, noyau lectron, lectron lectron) dans le systme et une fonction donde permettant de dcrire
le systme dans son ensemble. Une formulation alternative plus simple de cette diffrence dapproche est
quen mcanique quantique on sintresse au mouvement des lectrons, alors quen mcanique molculaire
les seuls mouvements pris en compte sont ceux des noyaux. Les lectrons sont simplement supposs
distribus de faon optimale autour des noyaux atomiques et peuvent alors tre assimils des charges
ponctuelles. Comme consquence, lnergie totale en mcanique quantique est lnergie vraie, mesurable
physiquement tandis que lnergie en mcanique molculaire est assimilable lnergie potentielle du
systme et est relative un tat de rfrence : ltat dquilibre du systme au repos.
Grce ces hypothses, la mcanique molculaire permet dtudier des systmes de plusieurs milliers
datomes avec un temps de calcul proportionnel M2, M tant le nombre de centres considrs (atomes
ou groupes datomes, par exemple 1 carbone), sans commune mesure avec le temps calcul ncessaire en
mcanique quantique qui dpend principalement du nombre dlectron N2,5 N6 (par exemple 12
lectrons sur un atome de carbone => 122,5498 1262 985 984 !). De plus lnergie potentielle du
systme permet daccder par les relations thermodynamiques de nombreuses autres grandeurs.
Les utilisations de la mcanique molculaire se sont accentues avec lavnement des ordinateurs et des
simulations numriques. On peut distinguer :
En mcanique molculaire une structure molculaire est donc dcrite comme un assemblage de billes
relies par des ressorts qui reprsentent explicitement les liaisons. La ncessit de formaliser les liaisons est
imprative alors quen mcanique quantique ce formalisme est inutile et la liaison entre deux atomes
dcoule de la probabilit leve de trouver des lectrons entre les deux atomes en question.
Les interactions nergtiques sont dcrites par des fonctions de potentiel mcanistiques (do le nom de
mcanique molculaire) et sont uniquement fonction de la gomtrie de la structure. Alors, les variables
II-42
mises en jeu dans l'expression mathmatique du potentiel dinteraction nergtique sont les seules
coordonnes internes de la structure ; ces coordonnes tant elles-mmes dfinies par la seule position des
noyaux entre eux :
Le modle mathmatique reprsentant l'nergie d'une molcule en mcanique molculaire est appel
champ de force. Sous ce terme sont en fait regroups deux lments : d'une part l'expression des
diffrentes fonctions contribuant au calcul nergtique et d'autre part les valeurs des diffrentes constantes
paramtrant ces fonctions. Ces paramtres sont identifis partir de donnes exprimentales ou valus
thoriquement. Lnergie totale du systme est donc une nergie potentielle multidimensionnelle dcrivant
les interactions intermolculaires ou interactions liantes (impliquant des atomes relis par des liaisons
explicites) et les interactions intramolculaires ou interactions non-liantes (impliquant des atomes non lis
directement par des liaisons explicites) :
(II-110)
+
-
longation
Interaction lectrostatique
Flexion
Torsion
Interaction liante
II-43
Par consquence, le champ de force est un modle du systme rel que lon veut tudier. De sa qualit
descriptive dpend la qualit du rsultat des simulations numriques par comparaison avec les valeurs
exprimentales. Lamlioration des champs de forces, modles de systme rel, reprsente un enjeu majeur
pour le succs des techniques de la simulation molculaire prdire avec prcision les proprits et le
comportement de systmes toujours plus nombreux et plus complexes.
La porte dune interaction intermolculaire qui concerne surtout les interactions non-liantes est
dtermine par lexposant de la distance entre les particules concernes : un potentiel en 1/r comme le
potentiel de coulomb gnrera une force en 1/r2 qui sera ressentie beaucoup plus loin (porte jusqu
25) quun potentiel en 1/r6 comme le terme dispersif du potentiel de Lennard-Jones (porte de 1 5).
a.
Les interactions entre atomes lis correspondent des nergies de dformation des liaisons, des
angles de valence et aux nergies de torsion :
(II-111)
1
2
k long (l l0 )
2 liaisons
(II-112)
Avec : klong la constante dlongation pour la liaison considre ; l la longueur relle de la liaison considre
et l0 la longueur de rfrence pour la liaison considre.
De mme, l'nergie de pliage est calcule en considrant l'ensemble des angles de valences de la structure
molculaire selon la formule :
E pliage =
1
2
k pliage ( 0 )
2 angle
(II-113)
Avec : kpliage la constante de pliage pour la liaison considre ; la valeur relle de langle de valence
considr et 0 la valeur de rfrence pour langle de valence considr.
II-44
Ce potentiel est souvent aussi exprim sous une autre forme harmonique quivalente :
E pliage =
1
2
k pliage (cos( ) cos( 0 ))
2 angle
(II-114)
En effet, un dveloppement limit lordre 2 au voisinage de 0 montre que ces deux dernires quations
sont quivalentes [ Ungerer, 1999 ] si :
(II-115)
Enfin, le potentiel de torsion permet de rendre compte de la dformation due aux variations de langle
didre entre deux liaisons spares par une troisime. Il est gnralement exprim par un dveloppement
en srie de la fonction cosinus selon les deux formes suivantes :
ETorsion =
Vn
2 (1 + cos(n ))
(II-116)
n=0
ETorsion =
cos( )
(II-117)
n=0
Dans ces expressions, Vn et Cn sont des constantes paramtriques, reprsente langle de torsion, est le
facteur de phase qui reprsente langle o lnergie de torsion passe par un minimum et n la multiplicit
qui reprsente le nombre de points minimum de la fonction lorsque la liaison est tourne de 360.
La dtermination des valeurs de rfrence, longueur et angle est ralise soit exprimentalement soit par
des calculs en mcanique quantique en calculant la variation de lnergie en fonction des diffrentes
longueurs et angles de faon dterminer les valeurs minimales qui deviennent les valeurs de rfrence et
aussi les constantes harmoniques permettant de lisser la forme quadratique du diagramme dnergie alors
obtenu.
La modlisation des interactions entres les atomes spars par plus de trois liaisons, E1-4, est effectue par
un potentiel de Lennard-Jones, du mme type que celui utilis pour modliser les interactions
intermolculaires de rpulsion-dispersion (prsentes juste aprs). Ce dernier terme permet dviter que les
atomes situs en bout de chane viennent se superposer avec ceux de lautre bout.
b.
Interaction lectrostatique
Linteraction lectrostatique entre diffrents centres de force dune mme molcule ou de deux
molcules, peut tre calcule par la loi de Coulomb comme une somme dinteractions entre des paires de
charges ponctuelles qi habituellement localises sur les atomes :
N A NB
ECoulomb =
i =1 j =1
qi q j
4 0rij
(II-118)
II-45
O 0 reprsente la permittivit du vide, et NA, NB le nombre de charges ponctuelles dans les deux
molcules. Ces charges peuvent tre dtermines par lissage des surfaces de potentiel lectrostatique
obtenues par des calculs ab-initio, ou dtermines partir de llectrongativit des atomes, ou ajustes de
faon reproduire un certain nombre de proprits thermodynamiques. Dans ce dernier cas, il nest pas
certain que les charges obtenues aient une signification physique. Par exemple, si on nimpose pas aux
charges de reproduire le moment dipolaire de la molcule isole, les interactions lectrostatiques ne
pourront pas prendre en compte la partie induction (polarisation) de lnergie dinteraction de manire
effective. Des exemples connus utilisant ce type de potentiel sont le modle TIP4P de Jorgensen
[Jorgensen et al., 1983] et le potentiel SPC de Berendsen pour la molcule deau [Berendsen et al., 1981]
pour laquelle la valeur des charges partielles et leur position sont dtermines par optimisation de lnergie
et de la structure en phase liquide. Ceci donne un diple effectif permanent environ 20% suprieur sa
valeur mesure exprimentalement en phase gaz. Cette procdure permet de donner des potentiels peu
coteux en temps de calcul mais qui risquent dtre peu transfrables.
ii.
Interaction de polarisation
Le phnomne de polarisation apparat lorsque les molcules du systme sont soumises un champ
lectrique extrieur. Dans ce cas, les nuages lectroniques des molcules se dforment. Cette dformation
du nuage lectronique se traduit par lapparition de moments induits. Lnergie de polarisation peut se
dcomposer en deux contributions : la premire, positive, reprsente lnergie quil faut fournir au systme
pour quil se polarise (car le nuage lectronique dform ne correspond pas une situation dquilibre de
la molcule isole) ; la deuxime, ngative, traduit la stabilisation due linteraction des moment induits
avec les autres charges et moments multipolaires prsents dans le milieu. Le systme se polarise de telle
sorte que la somme des deux (soit lnergie de polarisation) est minimale. Cette nergie est reprsente par
une fonction de potentiel en 1/r3 semblable la loi de Coulomb et a donc une porte moindre que
linteraction lectrostatique en 1/r. En pratique, ces interactions ne sont prises en compte que pour des
cas particuliers. En effet, il est souvent remarqu que les forces de polarisation contribuent de manire peu
importante lnergie totale dinteraction.
c.
Les interactions lectrostatiques et les nergies de polarisation elles seules ne peuvent pas tenir
compte de toutes les interactions non liantes du systme. Un exemple vident est celui des atomes dun
gaz rare o tout les moments multipolaires et les charges ponctuelles sont nuls. Il est cependant clair quil
y a des interactions entre les atomes, sinon comment peut-il apparatre sous forme liquide, solide, et
montrer un cart du comportement du gaz idal ? Ces dviations du comportement dun gaz idal ont t
quantifies par Van der Waals au 19me sicle.
II-46
Les interactions de Van der Waals sont la combinaison de deux termes : un terme rpulsif traduisant le
recouvrement des nuages lectroniques de mme polarit courtes distances et un terme attractif d aux
forces dispersives qui apparaissent instantanment durant les fluctuations des nuages lectroniques.
En effet, les forces rpulsives sont les plus locales et les plus intenses, leur porte est trs courte et elles
dlimitent une sphre quasi impntrable autour du noyau des atomes. Lnergie de rpulsion est une
fonction rapidement croissante lorsque la distance r, sparant les deux atomes, diminue. Elle est assez bien
reprsente par une fonction exponentielle :
VdW
Erep
(r ) = exp( r )
(II-119)
(r ) =
r
(II-120)
(II-121)
La plus connue des fonctions de potentiel de type Van der Waals est la fonction de Lennard-Jones 12-6,
dans laquelle les nergies de dispersion et de rpulsion sont dfinie par une seule et mme expression qui
scrit pour les deux molcules, neutres et non polaires, i et j sous la forme :
EijVdW LJ 12 6
= Erep + Edisp =
ij
4 ij
rij
12
ij
rij
(II-122)
O lindice i reprsente les atomes de la premire molcules et lindice j ceux de la deuxime molcules.
Cette quation contient deux paramtres ajustables : le diamtre de collision (ij) qui est la distance
minimale dapproche entre les deux atomes i et j pour laquelle lnergie entre deux atomes est nulle et la
II-47
profondeur du puits (ij) qui reprsente le minimum de lnergie potentielle, ce qui correspond
linteraction la plus stable. La Figure II-8 donne la forme de ce potentiel en fonction de r.
1/r12
U(r)
U(r)
U(r)
1/r12
r r
-1/r6
ij =
ii + jj
2
ij = ii jj
(II-123)
(II-124)
Les paramtres (ii), (jj), (ii) et (jj) sont les paramtres nergtiques du potentiel Lennard-Jones
dinteraction entre deux atomes de type identique. Une partie de notre travail a consist dterminer ces
paramtres pour les atomes de carbone et dazote qui forme la liaison -CN des nitriles. Nous discuterons
de la mthode utilise un peu plus loin. En effet, malgr le caractre physique illustr sur la Figure II-8, il
nexiste aucune mthode de calcul direct de ces paramtres. Leur obtention passe par une rgression de
faon reproduire des donnes exprimentales.
Dautres rgles de mlanges plus complexes existent mais cest la rgle de Lorentz-Berthelot qui reste la
plus utilise, prsentant lavantage dtre simple et fournissant des rsultats satisfaisants.
Ces rgles de mlange simples mettent en vidence que ltude dun systme multiconstituant comportant
M centres de forces diffrents ne ncessite que la connaissance des 2M paramtres n , n correspondant
et des ventuelles charges ponctuelles.
La diffrence est notable pour des systmes multiconstituants avec lapproche thermodynamique de type
quation dtat qui ncessite M (M + 1) 2 paramtres dinteractions binaires pour simuler un systme de
M composs chimiques. Par ailleurs, contrairement aux paramtres dinteraction binaire, les paramtres n,
II-48
capacits dextrapolation de la simulation molculaire dans une trs large plage de conditions de
temprature et de pression pour un systme donn au del de la zone pression - temprature pour laquelle
on dispose des donnes exprimentales.
Enfin, notons ici que bien que le potentiel de Lennard-Jones soit le plus utiliser pour dcrire les
interaction de Van der Waals, vu sa simplicit et le peu de paramtres quil fait intervenir ; il existe
plusieurs autres formulations dans lesquelles le terme en r12, traduisant la rpulsion, est remplac par un
terme en exponentielle. Parmi les potentiels adoptant ce formalisme, on retrouve le potentiel de
Buckingham, appel aussi potentiel exp-6, o les nergies de dispersion-rpulsion sont exprimes par
lquation suivante
Eexp 6 (r ) =
6
6
r rm
exp 1
1 6
rm r
pour
r > rmax
pour
r < rmax
(II-125)
Le rayon de coupure rmax reprsente la plus courte distance pour laquelle dEexp 6 (r ) dr = 0 obtenue par
rsolution itrative de lquation (II-125).
Ce potentiel semble donner des rsultats lgrement meilleurs celui de Lennard-Jones [ Errington and
Panagiotopoulos, 1999 ]. En effet, tout en tant utilis avec un modle UA (United atoms, voir plus loin),
les dviations par rapport aux donnes exprimentales sont comparables celle que prsentent les
modles AA (All Atoms, voir plus loin), tenant compte explicitement des atomes dhydrogne ; mais il
ncessite un nombre plus important de paramtres (les trois paramtres , et rm).
ii.
La liaison dhydrogne
La liaison hydrogne est une interaction entre un hydrogne globalement dficient en lectron et un
atome de forte densit lectronique portant un doublet dlectrons libres. Le modle lectrostatique,
appliqu lorsque la distance intermolculaire est grande, nest pas suffisant pour dcrire ces interactions
particulires. A plus courte distance, les phnomnes de rpulsion et de dlocalisation lectronique
interviennent.
Plusieurs champs de force tiennent compte explicitement de la liaison dhydrogne en ajoutant au terme
de Lennard-Jones 12-6 un nouveau terme souvent dcrit par une expression de la forme Lennard-Jones
12-10 entre latome dhydrogne donneur et lhtro-atome accepteur.
E (r ) =
A C
r12 r10
(II-126)
La majorit des champs de force considre que les liaisons hydrognes sont implicitement dcrites par le
potentiel de Van der Waals Lennard-Jones 12-6 au travers des valeurs rgresses des paramtres ii, ii.
II-49
d.
Il existe de nombreux modles permettant de dcrire les molcules [ Panagiotopoulos, 2001 ]. Tous
ces modles peuvent tre regroups en trois grandes familles selon que tous les atomes sont explicitement
pris en compte ou pas et selon la disposition des centres de force au sein de la molcule. Lagrgation
datomes en groupe a naturellement pour consquence la rduction du nombre de centres de force et
donc du temps de calcul.
i.
La solution la plus raliste du point de vue physique consiste utiliser un potentiel de type All
Atoms (AA) qui tient compte de tous les atomes. Ainsi, en ce qui concerne les hydrocarbures, les atomes
dhydrogne et les atomes de carbone sont deux types de centres de force diffrents, et tous les atomes
dune molcule sont reprsents par un centre de force. Plusieurs modles ont t dvelopps selon ce
schma pour prdire les proprits des alcanes saturs et insaturs. Les modles les plus utiliss lheure
actuelle sont les suivants : OPLS [ Jorgensen et al., 1984 ], TRAPPE-EH [ Chen et Siepmann, 1999 ],
AMBER [ Cornell et al., 1995 ], CHARMM [ Mackerell et al., 1995 ] et COMPASS [ Sun, 1998 ]. Tous ces
modles utilisent le potentiel de Lennard-Jones sous sa forme 12-6, mise part le champ de force
COMPASS qui utilise le potentiel 9-6. Linconvnient de cette reprsentation est quelle implique un
temps de calcul plus important que pour les modles agrgeant les atomes en groupes et demande un
effort didentification des paramtres aussi plus important.
ii.
Une deuxime approche, couramment utilise, consiste dcrire la molcule par des centres de
force regroupant plusieurs atomes. Pour le cas des hydrocarbures par exemple, chaque groupement CH3,
CH2 et CH est considr comme un pseudo atome, et un seul centre de force situ sur latome de carbone
reprsente le groupe. Il sagit du modle united atoms (UA). Cette reprsentation est plus conomique que
le modle AA puisque le nombre de centre de force est divis par deux, trois ou quatre selon les
molcules. Ce type de potentiel a donn de bons rsultats pour la description des hydrocarbures. Les
principaux potentiels de type Lennard-Jones utilisant ce modle sont le potentiel TRAPPE [ Martin et
Siepmann, 1998 ], et le potentiel NERD [ Nath et al., 1998 ] alors que le modle exp-6 [ Errington et
Panagiotopoulos, 1999 ] est bas sur la forme exponentielle dveloppe plus haut. Cependant, le modle
UA ne prend en compte que les distances entre atomes de carbone, et ne tient pas compte de la
disposition des atomes hydrogne. Le modle UA apparat moins apte reprsenter correctement les
interactions entre molcules que les modles AA, notamment haute pression lorsque la distance entre les
molcules est rduite.
iii.
Finalement, le modle utilis dans notre travail est celui propos par Toxvaerd [ 1990, 1997 ] et dj
utilis par le laboratoire de Chimie Physique dans des travaux antrieurs. Il consiste placer les centres de
II-50
force des pseudo atomes prs du centre gomtrique des groupements CH, CH2 et CH3. Ce potentiel
nomm Anisotropic United Atom (AUA) permet dune part de considrablement diminuer le temps de
calcul par rapport au potentiel de type AA, et dautre part de mieux diffrencier les groupements de type
diffrent que ne le fait le potentiel UA, en prenant en compte linfluence des atomes dhydrogne. Le
modle AUA donne notamment de bons rsultats pour les alcanes purs haute pression.
Molcule relle
Modle AA
Modle UA
Modle AUA
Figure II-9. Diffrents type de modles de champ de force [ Bourasseau, 2003 ].
La Figure II-9 reprsente les trois types de modle molculaire prsents au-dessus. La Table 1 prsente
les valeurs des paramtres de Lennard-Jones actuellement publis pour le champ de force AUA4.
e.
II-51
()
groupe
/kB (K)
()
Rfrence
CH3
3,6072
120,15
0,21584
CH2
3,4612
86,291
0,38405
CH
3,3625
50,980
0,64599
CH2 cycliques
3,4610
90,090
0,336
CH2 olfines
3,48
111,1
0,295
CH olfines
3,32
90,6
0,414
C olfines
3,02
61,9
>CH aromatiques
3,2464
89,415
0,407
Contreras [ 2002 ]
>C aromatiques
3,0648
42,079
Contreras [ 2002 ]
Cependant, confrer un pouvoir prdictif est plus que souhaitable pour un jeu de paramtres. Ainsi, les
paramtres optimiss sur une molcule seront transfrables pour dautres de la mme famille. Cest
lintrt espr du potentiel AUA4 que nous avons utilis dans ce travail. En effet, Ungerer [ Ungerer et
al., 2000 ] a montr quil tait possible doptimiser les paramtres des centres de force en se basant
uniquement sur quelques donnes exprimentales de trois molcules alcanes, lthane, le pentane et le
dodcane, pour les appliquer avec succs ensuite, tels quels, la prdiction des proprits dquilibre de
tous les alcanes linaires de C2 C20 et de nombreux alcanes ramifis. Ainsi, il a pu amliorer les rsultats
obtenus par le potentiel prcdent AUA3 pour les alcanes linaires (groupements CH3 et CH2) en se
basant uniquement sur trois grandeurs dquilibre (soit ln(Psat), Hvap et liq pour les trois groupes
datomes CH, CH2, CH3) deux tempratures diffrentes.
Pour loptimisation des paramtres Lennard-Jones des nitriles, nous avons adopt la mme procdure que
propose par Ungerer et applique aux alcanes. Il sagit dune mthode de minimisation de gradient qui
consiste optimiser simultanment les paramtres sigma et epsilon du potentiel de Lennard-Jones
ainsi que le paramtre delta spcifique au potentiel AUA.
La mthode consiste minimiser la fonction derreur adimensionnelle F, par rapport aux trois
paramtres tels que les drives partielles de F soient nulles :
F=
1
n
i =1
(f
cal
f i exp
si2
(II-127)
O si est lincertitude statistique estime sur la variable calcule fcal, tandis que fexp est la proprit associe
mesure dans les mmes conditions de tempratures et de pressions.
II-52
Deux mthodes de calcul de ces drives ont t proposes, elles ncessitent toutes les deux de partir dun
point de dpart doptimisation, c'est--dire un jeu de paramtres pour lequel tous les paramtres du
potentiel sont dfinis.
i.
Le calcul des drives par diffrence finie, utilise durant ce travail, a t la premire mthode
propose et utilise par Ungerer pour loptimisation des paramtres des alcanes, point de dpart trivial
pour ltablissement nimporte quel nouveau champ de force vu le vaste champ dapplication de ce type de
molcules.
Au cours de cette procdure, Le calcul des drives partielles seffectue par diffrence finie :
f i
x j
=
xk j
( )
f i x j + x j f i x j
(II-128)
x j
O fi est la variable thermodynamique simule (soit ln(Psat), Hvap et liq), xj le paramtre optimis (sigma
, epsilon ou delta ), fi(xj+xj) reprsente la valeur de la grandeur fi pour une variation de xj applique
au paramtre initial xj en gardant les deux autres constant (xk avec k j).
Cette technique prsente quelques inconvnients bien quelle a permis de donner de bons rsultats pour
les alcanes saturs. En effet, la prcision de lapproximation par diffrences finies dpend de xj et
ncessite dvaluer deux fois la fonction fi. Ainsi, pour optimiser simultanment les trois paramtres AUA
(sigma , epsilon et delta ) dun centre de force donn en utilisant deux molcules de rfrence et
deux tempratures diffrentes ; il faudra alors effectuer 16 simulations compltes (les quatre simulations
initiales pour chaque temprature et pour chaque molcule, en plus des trois autres simulations par couple
de temprature et de molcule ncessaire lvaluation des drives partielles) ; chaque simulation
ncessitant de gnrer plusieurs millions de configurations pour permettre une moyenne statistique
correcte.
ii.
Afin damliorer la qualit des rsultats et de rduire le temps doptimisation, une deuxime
technique pour lvaluation des drives partielles a t dveloppe par Bourasseau et al. [ 2003 ], base sur
la mthode des fluctuations qui permet dobtenir les drives partielles par lexpression gnrale suivante :
fi
x j
O
U
f i
fi
fi
x j
x j
fi
reprsente la moyenne de la grandeur fi optimiser (cette fois Psat, Hvap et liq) prise dans
U
x j
(II-129)
II-53
Cette fois, les drives sont considres comme des grandeurs statistiques de lensemble de travail, ce qui
donne plus de cohrence aux rsultats obtenus. De plus, cette technique permet de rduire
considrablement le nombre de simulations, puisquune seule et unique simulation permet lvaluation de
toutes les drives correspondantes une molcule donne pour des conditions de temprature et de
pression donnes.
Plus de dtails techniques concernant lvaluation des drives dans diffrents ensembles sont disponibles
dans la thse de Bourasseau [ 2003 ] qui a appliqu cette mthode avec succs pour loptimisation des
paramtres Lennard-Jones des olfines.
f.
Mme si les principales formes fonctionnelles des potentiels (longation, pliage, torsion, Van der
Waals, Coulomb) sont prsentent dans tous les champs de force cits, le choix doit tre fait en connaissant
le type de donnes exprimentales qui ont t utilises pour rgresser les paramtres de Lennard-Jones et
la faon dont sont calcules et attribues les charges ponctuelles.
Un autre dbat est la prise en compte explicite de la liaison hydrogne ou implicite en lincluant dans
lnergie de Van der Waals. Cela a des consquences dans le choix dun champ de force. Celui-ci
comportera forcment des paramtres semi-empiriques tels que les paramres de Lennard-Jones et . Il
faut viter de mlanger les paramtres de Van der Waals de champs de force qui diffreraient dans la prise
en compte ou non de certains termes tels que la liaison hydrogne ou une extension du terme rpulsif.
Les plus anciens des champs de force (OPLS ou AMBER par exemple) ont t dvelopps en ne tenant
que des donnes relatives la phase liquide dans des conditions proches de la temprature ambiante, et
non lquilibre liquide vapeur sur une large gamme de temprature. Actuellement, un effort
considrable est en cours au niveau mondial pour tendre les domaines dapplications des champs de
force isols afin de les harmoniser.
II-54
D.
Cette partie met en lumire quelques principes fondamentaux des quilibres liquide vapeur
(paragraphe 1) et de leur modlisation au moyen de modles macroscopiques (paragraphe 2) et par
simulation molculaire (paragraphe 3). Les rfrences bibliographiques des paragraphes 1 et 2 ne sont pas
cites dans le texte mais se trouvent dans les excellents livres de Sandler [ Sandler, 1993 ] et de Vidal
[ Vidal, 1997 ].
1.
Degr de libert
Il est ncessaire, pour dfinir compltement un systme, dindiquer les valeurs numriques dun
certain nombre de variables, que lon choisit par exemple parmi la pression, la temprature, lnergie
interne, lentropie et les concentrations des diffrentes phases. Parmi les variables, on distingue les
variables extensives (qui dpendent de la masse) et les variables intensives dites variables dtat. Un
quilibre entre phase ne dpendant pas des quantits relles des diverses phases en prsence, on ne tiendra
compte que des variables intensives, comme la temprature, la pression et la concentration pour le dcrire.
On peut faire varier un certain nombre de ces variables de manire indpendante, mais les valeurs des
autres sont dtermines par les valeurs choisies et par les conditions thermodynamiques dquilibre. Alors,
le degr de libert ou variance du systme est dfini par le nombre de variables dtat quon peut faire
varier de manire indpendante sans modifier le nombre de phase.
b.
La rgle de phase
La rgle de phase propose par Gibbs fournit une relation gnrale entre le nombre de degr de
libert dun systme, le nombre de phase et le nombre de constituants C. Elle exprime quon a toujours
pour un systme non ractionnel :
= C + 2
(II-130)
Par exemple, pour un systme binaire (C=2) en quilibre liquide - vapeur (=2), la rgle de phase indique
que la variance est gale : =2. Ceci signifie que lon peut faire varier de manire indpendante deux
II-55
variables dtat sans modifier le nombre de phases. Ces variables sont en gnral choisies parmi la
temprature, la pression et la composition.
c.
dP sat
H
=
dT
T V
(II-131)
d ln P sat H vap
=
dT
RT 2
(II-132)
Lune des relations les plus importantes fournissant cette variation pour des tempratures nettement
infrieures la temprature critique est lquation empirique dAntoine, de la forme :
log P sat = A
B
T +C
(II-133)
Pour chaque constituant i du mlange, la condition dquilibre thermodynamique entre phase est
donne par lgalit des potentiels chimiques des phases qui se rduit, par dfinition de la fugacit,
lgalit des fugacits dans chaque phase pour chaque constituant. Dans le cas dun quilibre liquide
vapeur :
f iV = f i L
(II-134)
Principe de Le Chtelier : toute modification dun des facteurs dterminant ltat dun systme en
quilibre provoque un dplacement de lquilibre dans un sens qui tend sopposer la variation du facteur
considr .
II-56
Le problme fondamental est de relier ces fugacits aux compositions de chacune des phases. Dans ce qui
suit, nous ngligeons les effets dus aux forces de surface, gravit, champs lectriques ou magntiques et
membranes semi-permables. Alors, la fugacit dun constituant dans un mlange ne dpend que de la
temprature, de la pression et de la composition de ce mlange. On exprimera les compositions de la
vapeur et du liquide par les fractions molaires, yi et xi respectivement.
Si la pression est faible et le mlange idal (peu dinteractions entre les molcules prsentes dans le
mlange), lquation (II-134) se simplifie en la loi de Raoult :
yi P = xi Pi 0
Loi de Raoult
(II-135)
o P est la pression totale, yi P reprsente la pression partielle du constituant i dans la phase vapeur et Pi 0
la pression de vapeur saturante du constituant i ; une fonction de la temprature T que lon peut calculer
par la loi dAntoine.
2.
Comme le montre la Figure II-10 pour un certain nombre de mlanges binaires, lhypothse dun
mlange idal nest valable que dans un nombre trs limit de cas (constituants similaires, par exemple : nbutane et isobutane) et peut introduire des comportements non idaux induisant lapparition dazotrope
htrogne, dazotrope homogne temprature dbullition minimale ou maximale ou minimale avec
dmixtion dans une certaine plage de composition. Sans apparition dazotrope, dautres mlanges
peuvent prsenter des pincements sur le diagramme composition y-x. Des phnomnes de double
azotropie ont aussi t observs exprimentalement. Tous les mlanges non idaux dcrits dans la Figure
II-10 sont difficiles ou impossible sparer par un procd de distillation classique et requirent des
procds de distillation complexes (Gerbaud, HDR 2003).
De faon pratique, lcart lidalit est introduit dans la loi de Raoult par des coefficients correctifs
qui tendent vers 1 pour un mlange idal. On distingue deux approches :
1. Le modle commun dans lequel les deux phases sont dcrites par la mme quation dtat
iV yi P = iL xi P
(II-136)
2. Le modle combin suivant lequel la phase vapeur est dcrite par une quation dtat et la phase
liquide par un modle thermodynamique.
iV yi P = i xi f i 0 L
(II-137)
O i est le coefficient de fugacit et tend vers 1 lorsque P tend vers 0 ; et i est un coefficient dactivit de
la phase liquide gal 1 pour un mlange idal.
II-57
Min T htroazotrope
dviation/Raoult [+]
eau hexane
Mlange idal
n-butane n-heptane
T
Vapeur
P = cst.
Vapeur
P = cst.
Courbe de rose
L2+V
L1+V
Tazeo
Courbe de bulle
a
L1
b
Liquide
xazeo
xA, yA
eau thanol
T
Vapeur
chloroforme - actone
T
P = cst.
Vapeur
P = cst.
xA, yA
Min T homoazotrope
avec rgion LLV
dviation/Raoult [+]
eau 2-butanone
Max T homoazotrope
dviation/Raoult [-]
Min T homoazotrope
dviation/Raoult [+]
L2
L1 + L 2
Vapeur
P = cst.
Tazeo
L1
L2+V
L1+V
Tazeo
Tazeo
Liquide
xazeo xA, yA
Mlange zotropique
avec pincement
propane propylne
Vapeur
Liquide
T
P = cst.
f
xA, yA
Vapeur
Liquide
L1 + L 2
L2
xazeo
xA, yA
xazeo
Mlange zotropique
temperature dbullition proche
T
Liquide
xA, yA
Mlange zotropique
avec pincement
eau actone
T
P = cst.
g
xA, yA
Vapeur
P = cst.
Liquide
h
xA, yA
Figure II-10. Equilibres liquide vapeur pour diffrents systmes binaires non idaux
Pour les systmes de constituants peu polaires (hydrocarbures, ) dans lesquels les interactions
intermolculaires sont faibles, on privilgiera le modle commun (approche par quation dtat).
En gnral, les quations dtat ne dcrivent pas correctement le comportement des systmes de
constituants polaires ou avec de fortes interactions (alcools, eau, nitriles, ). On privilgiera alors
lapproche combine.
II-58
A pression faible (jusqu 3-5 atm), considrant que les molcules dans la phase vapeur interagissent peu,
on peut faire lhypothse dun gaz parfait qui implique que le coefficient de fugacit de la phase vapeur iV
est gal 1 et disparat des quations (II-136) et (II-137).
a.
2P
P
= 2
=0
V T =Tc V T =Tc
(II-138)
Ceci est d au fait que lisotherme critique possde une pente horizontale et un point dinflexion au point
critique. La Figure II-11 prsente le diagramme de Clapeyron (P-V) de lisobutane en coordonnes
logarithmiques.
II-59
10 000 P, psi
Point critique
Enveloppe du domaine exprimental
liquide-vapeur
Isotherme exprimentale pour 190F
Equation de Redlich Kwong pour 190F
isobutane
Tc
T = 190F
1 000
Pc
0
P (190 F)
100
vapeur
sature
Liquide
satur
10
Vexp Vcal
10
100
1 000
V, ft3 /lb mole
La premire quation dtat reliant les proprits macroscopiques dun systme est la loi des gaz
parfait dcoulant de la loi de Boyle-Mariotte :
Pv = RT
(II-139)
II-60
Pv = ZRT
(II-140)
Van der Waals a constat que le facteur de compressibilit pour des fluides simples avec des molcules
sphriques comme les gaz rares et le mthane ne dpend que des paramtres rduits Tr et Pr dfinis de la
manire suivante :
Tr =
T
Tc
et Pr =
P
Pc
(II-141)
avec : Prsat =
P sat
Pc
= log Prsat
Tr = 0,7
(II-142)
(II-143)
Le premier terme de cette galit (Z(0)) reprsente le facteur de compressibilit pour une molcule
sphrique, alors que Z(1) exprime la dviation par rapport la sphricit idale de la molcule. Pitzer et ses
collaborateurs ont rparti dans un tableau les valeurs de Z(0) et Z(1) en fonction des coordonnes rduites
Tr et Pr, alors que Edmister propose les mmes valeurs sous forme de diagrammes.
ii.
Plusieurs quations dtat empiriques ont t proposes, mais la seule qui reste fonde sur des bases
thoriques solides est lquation du viriel. Cette quation est un dveloppement en srie de puissance du
facteur de compressibilit autour du point o fluide rel et gaz parfait se confondent, ce qui correspond
une masse volumique nulle. On crira alors :
Z=
Pv
B C
=1+ + 2 +L
RT
v v
(II-144)
II-61
O les coefficients B, C, sont appels respectivement second, troisime, coefficient du viriel et sont
associs respectivement aux interactions deux deux, trois trois, . Pour un fluide pur, les coefficients
du viriel ne dpendent que de la temprature.
On peut galement dfinir un dveloppement en pression sous la forme :
Z = 1 + BP + C P 2 + L
(II-145)
dont les paramtres sont lis aux coefficient du viriel par des relations dont nous ne citons que les deux
premires : B =
B
RT
C =
et
C B2
R 2T 2
Ces dveloppements en srie ne sauraient tre considrs comme une quation dtat, puisquils
comportent un nombre infini de termes dont la valeur est inconnue et la sommation impossible. Par
contre, leur troncature a donne naissance aux quations dtat du viriel que lon caractrise par lordre
de la troncature et par la nature du dveloppement dont elles sont drives, selon quil soit en volume ou
en pression.
Au-del de son application pratique restreinte quelques corps purs (NH3, CO2, H2O, ) cette quation
permet de dmontrer la justification des rgles de mlange de Van des Waals qui restent universellement
utilises pour dcrire le comportement de mlanges de fluides partir dquations dtats de fluides purs.
iii.
Lquation des gaz parfaits nest utilisable que dans un domaine de pression restreint (celui de
basses pressions), pour cette raison cette quation a t modifie sous la forme propose par Hirn :
(P + )(v b ) = RT
(II-146)
o b appel covolume reprsente la part du volume molaire qui nest pas rellement disponible cause de la
prsence dautres molcules, et est un terme correctif de pression d lattraction mutuelle des
molcules, appel pression interne.
La premire quation dtat fournir une description quantitative des deux phases vapeur et liquide et une
prdiction de la transition entre phases est celle de Van der Waals en 1873 qui a tabli thoriquement que
la pression interne est inversement proportionnelle au carr du volume molaire, do son quation dtat :
P + 2 (v b ) = RT
v
Ou bien, P =
RT
a
2
(v b ) v
(II-147)
(II-148)
Dans cette dernire quation, le premier terme associ au covolume b caractrise la rpulsion et le second
associ la pression interne lattraction. Si cette formule ne reprsente pas parfaitement les isothermes
II-62
dun fluide (voir Figure II-11), elle ne conduit pas des absurdits physiques, do sa validit et son
application pour des calculs approximatifs.
Pour la dtermination des paramtres de cette quation (et ceux des autres quations dtat du mme
type), il existe deux alternatives :
1. la premire consiste choisir des paramtres qui permettent de reproduire des donnes
exprimentales, gnralement la pression vapeur et la densit liquide ou vapeur. Avec lquation
de Van der Waals dont les coefficients ne dpendent pas de la temprature, une seule temprature
suffit, pour les autres quations, par contre, o ces paramtres seront fonction de la temprature,
il faudra rgresser ces paramtres sur un intervalle de temprature.
2. la deuxime procdure, qui reste la plus recommande, utilise la condition du point critique,
donne par lquation (II-138), pour obtenir les valeurs de ces paramtres aux points critiques ac et
bc ; et ensuite applique une correction qui dpendra de la temprature (c'est--dire a(T)=ac(T)),
gale lunit au point critique et ajust pour permettre de mieux reproduire tout lintervalle de
temprature.
Hormis son caractre historique, lquation de Van der Waals nest pas quantitativement prcise. Par
exemple, elle prdit que le facteur de compressibilit critique :
Zc =
Pc vc
RTc
(II-149)
est gal 3/8 = 0,375 pour tous les fluides ; alors que sa valeur pour diffrents hydrocarbures varie de
0,24 0,29, et cet intervalle est plus grand si dautres fluides non hydrocarbures sont considrs. De
mme, les pressions de saturation prdites sont incorrectes. A cause de ces insuffisances, de nombreux
chercheurs ont tent de modifier la forme de lquation de Van der Waals afin damliorer sa prcision.
Dans le but de mieux reprsenter les proprits volumtriques des fluides, Redlich et Kwong ont propos
une modification empirique du terme dattraction de Van der Waals en introduisant galement une
variation du paramtre a avec la temprature. Ceci permet damliorer les calculs des masses volumiques,
en particulier ceux de la phase vapeur et de fournir une meilleure estimation du facteur de compressibilit
critique (Zc = 0,333). Par contre, les prdictions de la masse volumique liquide et les tensions de vapeurs
restent mdiocres.
P=
RT
a
(v b ) v(v + b )
(avec = T
12
(II-150)
Exprime en fonction du facteur acentrique, lquation de Redlich-Kwong a une forme cubique en posant
A=
aP
(RT )
et B =
bP
:
RT
II-63
Z 3 Z 2 + A B B 2 Z AB = 0
Redlich-Kwong
(II-151)
La reprsentation de cette quation sur la Figure II-11 illustre bien les trois racines de cette quation sur
lisotherme dquilibre liquide vapeur : deux racines correspondent aux points de bulle et de rose, tandis
que la troisime situe dans la rgion liquide vapeur na pas de signification physique. Comme le montre
galement la Figure II-11, lquation de Redlich-Kwong ne sapplique bien qu la phase vapeur.
Pour la phase liquide, des relations telles que celle de Rackett [ Vidal, 1997 ] dcrivent le volume liquide
saturation.
Cependant, on a remarqu que les quations dtat cubique deux paramtres ne peuvent reprsenter avec
une prcision satisfaisante la fois le comportement volumtrique et les tensions de vapeur, et quil
convient de choisir une priorit. Cest ce qua compris Soave en 1972 qui a appliqu lquation dtat de
Redlich-Kwong au calcul des tensions de vapeur des hydrocarbures en modifiant lexpression du
paramtre devenu fonction de la temprature et du facteur acentrique de Pitzer :
[ (
) (
)]
Soave-Redlich-Kwong
(II-152)
Soave-Redlich-Kwong
(II-153)
Z 3 Z 2 + A B B 2 Z AB = 0
P=
RT
a
(v b ) v(v + b ) + b(v b )
Peng-Robinson
(II-154)
Peng-Robinson
(II-155)
Peng-Robinson
(II-156)
Avec :
[ (
) (
)]
) (
Z 3 (1 B )Z 2 + A 2 B 3B 2 Z AB B 2 B 3 = 0
II-64
Les quations de Peng-Robinson (PR) et de Soave-Redlich-Kwong (SRK) sont largement utilises dans
lindustrie, particulirement pour le raffinage et la simulation des rservoirs. Leurs avantages rsident dans
le fait quelles ncessitent peu de donnes exprimentales (les coordonnes critiques et le facteur
acentrique), un temps de simulation relativement court et surtout quelles conduisent une bonne
estimation des quilibres liquide vapeur pour les hydrocarbures qui sont particulirement importants
pour la conception des procds. Par ailleurs, ces quations prsentent dimportantes limitations. Par
exemple, les densits liquides ne sont pas correctement values, les paramtres ne sont pas adquats pour
les autres fluides non hydrocarbures, notamment les fluides polaires. De plus, quelles que soient les
modifications apportes, il apparat impossible de reprsenter de manire satisfaisante la rgion critique.
Abbott (1979), dans une tude comparative des quations dtat cubique, a conclu que la forme la plus
adquate doit faire intervenir des paramtres supplmentaires. De nouvelles quations ont t
dveloppes, toutes convergeant vers les quations SRK ou PR lorsque des valeurs particulires sont
attribues ces paramtres additionnels. Ce changement permet damliorer simultanment la qualit de
prdiction de la pression de saturation et du volume molaire liquide. Par contre, le progrs au voisinage du
point critique est moins marquant (par rapport celui des quations de SRK et PR), car, comme signal
par Abbott [ 1979 ], toutes les quations dtat cubiques doivent surestimer la compressibilit critique pour
quelles soient prcises dans les autres rgions.
Avant de finir cette courte synthse sur les quations dtat cubiques, citons la mthode base sur le
concept de translation de volume (Peneloux et al., 1982). Cette translation, opre sur les quations de Van
der Waals, SRK et PR, a pour effet de niveler les diffrences quelles prsentaient pour le calcul des
volumes molaires. Par exemple, la translation v v + c et b b + c applique lquation SRK donne :
P=
a
RT
(v b ) (v + c )(v + b + 2c )
(II-157)
O le paramtre de translation de volume c est choisi pour prdire la valeur correcte du volume molaire
liquide une temprature donne. On prend gnralement une temprature proche de la temprature
dbullition, ce qui correspond une temprature rduite au voisinage de 0,70. Sur la Figure II-11, le
paramtre c pour lquation de Redlich-Kwong correspondrait la diffrence entre VLexp et VLcal si la
temprature de 190F correspondait Tr=0,70. Cette procdure permet damliorer la densit liquide tout
en ayant une incidence ngligeable sur celle de la phase vapeur. Dautre part, les diffrences entre les
proprits thermodynamiques des phases liquide et vapeur en quilibre ne sont pas modifies par la
translation ; ainsi lenthalpie de vaporisation reste inchange.
En ce qui concerne leur application, les quations dtats reprsentent bien le comportement des
constituants non ou faiblement polaires dans une large gamme de pression. Elles ne sont pas adaptes
pour les constituants polaires dans lesquelles les interactions en phase liquide sont importantes et ncessite
un autre formalisme comme celui des modles de coefficients dactivit prsents aprs.
II-65
Enfin, pour un certain nombre de composs dimportance industrielle, tels que leau, lammoniac, les
hydrocarbures lgers, on dispose de donnes exprimentales abondantes et varies (mesures PVT ,
tension de vapeur, capacits calorifiques). En sappuyant sur ces donnes, pralablement soumises une
valuation rigoureuse, on a propos des quations dtat susceptibles de les restituer avec une excellente
prcision ; on peut ainsi considrer le rsultat du calcul comme une donne quasi exprimentale . Ces
quation dtat sont souvent spcifiques, sinon dun constituant, du moins dun petit groupe de composs,
et se prtent mal la reprsentation des mlanges, en particulier au calcul des quilibres liquide vapeur.
b.
Lors de la prsentation des quations dtat, nous avons insist sur le fait que ltat du gaz parfait ne
pouvait prsenter une approximation acceptable du gaz rel qu basse pression. De mme, en ce qui
concerne les mlanges, le concept de solution idale ne constitue quune tape dans lvaluation des
proprits dun mlange. Les termes qui dcrivent la diffrence entre les proprits dun mlange rel et
celles dun mlange idal sont appels, quand il sagit de proprits extensives, grandeurs dexcs .
On dfinit alors respectivement par :
V E = V V id ; S E = S S id ; G E = G G id
(II-158)
Le volume dexcs, lentropie dexcs et lenthalpie libre dexcs. Cette dernire proprit est
particulirement employe lors de ltude des quilibres liquide vapeur du fait quelle est lie aux
coefficients dactivit i par la relation pour un mlange de m constituants comportant ni moles de
constituant i :
G E = RT
(n ln )
i
(II-159)
i =1
La notion de non idalit a impos dintroduire un coefficient dactivit i pour la phase liquide. Ce
coefficient apparat aussi dans la relation thermodynamique de Gibbs-Duhem qui scrit pour un mlange
de deux constituants :
ln 1
ln 2
= x2
x1
x1 T , P
x2 T , P
(II-160)
II-66
G E
RT ln i =
ni
T , P,n j ,i j
(II-161)
Les quations (II-159) et (II-161) sont utiles dans la pratique puisquelles permettent dinterpoler ou
dextrapoler des donnes incompltes par rapport la composition, pourvu que lon connaisse une
expression mathmatique de GE en fonction de la composition, appel modle de coefficient dactivit.
iii.
les modles de Van Laar et de Margules dont lquation de lenthalpie libre dexcs est exprime
laide de deux paramtres qui sont fonction de la temprature et qui peuvent tre obtenus partir
des donnes exprimentales (plus prcisment dilution infinie). Ces modles sont limits en
gnral aux systmes binaires.
Les modles fonds sur le concept de la composition locale introduit par G. M. Wilson en 1964.
Ce concept est lorigine des modles qui se sont rvls les meilleurs pour la corrlation et la
prdiction des dviations lidalit de la phase liquide et sont intensivement utiliss dans
lindustrie chimique. Il est bas sur ltablissement dun bilan nergtique du mlange pour lequel
on doit introduire des approximations et proposer une valuation des compositions locales. Il
conviendra ensuite de passer du bilan nergtique lenthalpie libre dexcs. Chacune de ces
tapes dterminera un modle empirique dont la valeur se jugera par son sens physique et son
pouvoir de reprsentation des valeurs exprimentales. On retrouve principalement le modle de
Wilson propos en 1963, le modle NRTL (Non Random Two Liquids, propos par Renon et
Prausnitz en 1968) et le modle UNIQUAC (UNIversal QUAsi Chemical, propos par Abrams et
Prausnitz en 1975 et amlior par Maurer et Prausnitz en 1978, ce modle a t dvelopp sur une
base thorique de thermodynamique statistique)
Le modle UNIQUAC tablit une quation en deux parties, une partie combinatoire qui dcrit la
contribution entropique dominante, une partie rsiduelle qui est due principalement aux forces
intermolculaires qui sont responsables de lenthalpie de mlange, mesure de lnergie dgage la suite du
mlange de constituants purs. La part combinatoire est dtermine par la composition xi et par la forme et
la taille des molcules. Elle requiert uniquement des donnes relatives aux constituants purs. La part
rsiduelle dpend en plus des forces entre les molcules, prises en compte par lintermdiaire de deux
paramtres dinteraction binaire ij et ji.
II-67
GE GE
=
RT RT
GE
combinatoire RT
rsiduelle
(II-162)
G
= f combinatoire ( xi , r, q ) + f rsiduelle ( xi , q' , ij , ji )
RT
Les paramtres r, q et q sont des constantes molculaires, pour chaque constituant pur, relatives
respectivement sa taille, sa surface gomtrique externe et sa surface dinteraction. q peut tre diffrent
de q, notamment pour les molcules polaires : ainsi, pour les alcools, linteraction molculaire est domine
par le groupe hydroxyle OH susceptible de former des liaisons hydrogne. q (0,92 pour lthanol) est
alors plus petit que q (1,97 pour lthanol) puisque les autres groupes chimiques des alcools (CH3 et CH2
principalement) interviennent trs peu dans linteraction intermolculaire mais comptent nanmoins pour
valuer la forme de la molcule.
Les deux paramtres dinteraction binaire ij et ji sexpriment en fonction des nergies dinteraction
Uij entre une molcule i et une molcule j, Uii entre deux molcules i et Ujj entre deux molcules j.
ij = exp
1
2
(U
1
2
(U
ji = exp
ij
U jj N A
RT
U ii N A
RT
ji
(II-163)
O NA est le nombre dAvogadro. Faute de pouvoir calculer systmatiquement les nergies dinteractions
Uij, Uii et Ujj, par exemple par des calculs de chimie quantique, ij et ji deviennent des paramtres
ajustables quil faut identifier partir de la rgression de donnes exprimentales.
Linconvnient des modles semi-prdictifs rside dans la ncessit de disposer de paramtres
dinteraction binaire pour prdire le comportement de mlanges de fluides. Des sources importantes
existent (DECHEMA Chemical series) mais restent trs insuffisantes : on trouve moins de 50 000
systmes binaires recenss lorsque la combinatoire des plus de 15 millions de constituants chimiques
existant avec dinnombrables conditions opratoires (T,P) possibles en laisse prsager une infinit. En
outre, il est avr quune infinit de valeurs de paramtres dinteraction binaire peuvent reprsenter le
mme ensemble de donnes exprimentales. Il est donc conseill de nutiliser les modles semi-prdictifs
que pour interpoler des donnes existantes, et jamais pour extrapoler.
En termes dapplication, les modles de coefficient dactivit en gnral donnent des rsultats mitigs
forte pression mais sont adapts pour la prdiction des proprits des constituants polaires non
lectrolytiques faible pression. La reprsentation de systmes lectrolytiques est cependant possible, par
exemple en ajoutant une contribution lectrolytique aux termes combinatoire et rsiduels de lexpression
du modle UNIQUAC.
II-68
iv.
ln i = ln i, combinatoire + ln i, rsiduel
(II-164)
Les diffrences entre ces deux mthodes, ASOG et UNIFAC, dcoulent des modles qui ont t choisis
pour exprimer le terme combinatoire et rsiduel. Chacune proposent des tables de paramtres pour une
liste de groupes chimiques. En particulier, la mthode UNIFAC a fait lobjet de publication de tables de
paramtres restreintes et des tables tendues sont disponibles pour les membres du consortium UNIFAC
(http://134.106.215.86/UNIFAC/).
Parmi les inconvnients de ces mthodes, citons :
des incohrences structurelles qui ont conduit publier des tables de paramtres UNIFAC
diffrents pour la reprsentation des quilibres liquide liquide et liquide vapeur , mais aussi
ont conduit dfinir des groupes principaux et des groupes secondaires, et aussi des paramtres
fonction de la temprature,
dans les mthodes initiales seules les contributions additives sont prises en compte sans
considrer les contributions lies aux interactions entre les groupes chimiques. Une
consquence est limpossibilit de distinguer les isomres. Les mthodes avec des groupes
secondaires ou dcrivant la connectivit des groupes sont ncessaires.
Constantinou et Gani et leurs collgues [ Constantinou et Gani, 1994 ; Constantinou et al., 1996 ] ont fait
une revue critique des mthodes de contribution de groupe et prsent une nouvelle mthode incluant
groupes primaires et secondaires et son application pour la prdiction de proprits dquilibre.
II-69
Parmi les avantages, ces mthodes peuvent directement tre incorpores aux quations dtats travers la
dfinition de rgles de mlange complexes telles que celles de Huron-Vidal et ses modifications MHV1 et
MHV2, de Wong et Sandler et de Gmehling (PSRK [ Chen et al., 2002 ]) afin de combiner les avantages
des quations dtats (utilisation aux faibles et hautes pressions) et des modles de coefficient dactivit
(description de la non idalit de la phase liquide).
Enfin, signalons que dans le cadre de la conception de nouvelles molcules pour lesquelles bien souvent
aucune donne nest disponible, les mthodes de contribution de groupes sont employes dans les phases
de recherche prliminaire et quensuite, la modlisation molculaire ou les expriences simposent afin
dobtenir des valeurs des proprits plus prcises pour affiner le choix des molcules [ Harper et al.,
1999 ].
c.
Comme nous lavons suggr lors de la prsentation du modle UNIQUAC, il est possible de relier
les paramtres dinteraction binaire aux nergies dinteraction entre les molcules semblables et
dissemblables [ Vidal, 1997 ]. Par consquent, il est thoriquement possible de calculer ces nergies par
simulation molculaire. Plusieurs tentatives ont t ralises avec un rsultat convaincant par des
mthodes ab-initio et moins convaincant par des mthodes de mcanique molculaire.
Dautres approches existent comme les modles de type COSMO [ Klamt and Eckert, 2000 ] o lon
value le changement denthalpie libre lors du transfert dune molcule en phase gaz idale dans un solvant
trait comme un continuum polarisable. Une revue rcente met en vidence les avantages et les
inconvnients de ces approches avec des mthodes ab-initio pour des systmes dintrt industriel [
Gerbaud, 2003 ; Sandler, 2003 ]. Sandler dclare que avec les futurs dveloppements esprs, les
modles COSMO remplaceront les modles prdictifs comme UNIFAC bass sur le concept de
contribution de groupe. .
Enfin, les quation de type SAFT (Softely Associated Fluid Thory) connaissent un dveloppement
important depuis 15 ans [ Mller and Gubbins, 2001 ]. Elles valuent lnergie libre de Helmoltz comme la
somme de plusieurs contributions formalisant une reprsentation des molcules et de leur interaction :
contribution de monomre, contribution de liaison de ces monomres pour former une molcule ;
contribution dassociation intermolculaire entre sites dinteraction. Ce formalisme est proche de celui
imagin en mcanique molculaire et lun des atouts de lapproche SAFT est le fait que les prdictions
dune telle quation peuvent tre confrontes des simulations molculaires utilisant le mme formalisme.
d.
Comme nous lavons signal au dbut de ce chapitre, les quations dtat ont t dveloppes
initialement pour des fluides purs puis tendues pour les mlanges. Cette extension ncessite
lintroduction de rgles de mlange permettant le calcul des paramtres de mlange en fonction de ceux
des corps purs.
II-70
La rgle de mlange la plus rpandue est sans doute la rgle de mlange conventionnelle o rgle de Van
der Waals qui permet de calculer les coefficients de mlange am et bm par les deux relations :
am =
x x
i
aii a jj 1 kij
bm =
1
2
x x (b
i
ii
)(
+ b jj 1 lij
(II-165)
kij et lij sont les paramtres dinteraction binaire, ils expriment les forces dinteraction exerces entre les
diffrentes types de particules i et j qui constituent le mlange. Ils doivent tre rgresss partir de
donnes exprimentales du binaire considr. Gnralement, lij est pris gal zro, ce qui permet dcrire
le second coefficient sous la forme : bm =
xb
i ii
La justification thorique de ces rgles de mlange provient de lquation dtat du viriel donne par
lquation (II-142). En effet, la mcanique statistique montre que le second B et le troisime coefficient C
du viriel ne dpendent que de la temprature. Dans le cas des mlanges ils scrivent en fonction de la
composition du mlange sous la forme :
B=
x x B (T )
C=
x x x C (T )
(II-166)
ij
j k
ijk
Dautre part, en faisant un dveloppement limit, par exemple de lquation de Van der Waals en
puissance de b
( )
a
1
pv
b
=
=1+
v
vRT
RT 1 + b
n =0
v
a
vRT
(II-167)
B(xi , T ) =
x x B (T ) = b RT
i
ij
(II-168)
Notons que ces rgles de mlange ne sont applicables quaux mlanges prsentant une dviation
relativement faible par rapport lidalit ; car comme nous venons de le voir le dveloppement limit
nest valable qu faibles densits. Cependant, jusqu prsent, ce sont ces rgles de mlange qui sont les
plus utilises nimporte quelle densit, du moment o la mcanique statistique ne fournie pas encore
suffisamment dinformation valide densits leves pour pouvoir vrifier la crdibilit des rsultats
obtenus.
Lie aux quations dtat de type Van der Waals, cette rgle de mlange classique ne doit sappliquer
quaux systmes non polaires (mme domaine dapplication que les quations dtat dont elle drive) ou
dans le cas des solutions prsentant une lgre dviation par rapport lidalit et qui peuvent tre dcrits
par la thorie de la solution rgulire. Cependant, la plupart des solutions intrt industriel sont polaires
II-71
ou prsentent une importante dviation par rapport lidalit et sont gnralement dcrites par les
modles de coefficient dactivit. Vidal en 1978 et plus tard Huron et Vidal en 1979 ont introduit une
nouvelle gnration de rgles de mlange combinant une quation dtat avec un modle de coefficient
dactivit. Cette combinaison est ralise par le biais de lexpression de la fonction de lenthalpie libre
dexcs calcule partir dune quation dtat pour dvelopper une rgle de mlange qui reproduit le
mme rsultat :
G E = RT ln x i ln i
i
(II-169)
II-72
kij = K ij K ij K ji xi
(II-170)
et de Schwartzentruber et Renon :
kij = K ij + lij
mij xij m ji x j
mij xij + m ji x j
(x + x )
i
(II-171)
Les mthodes de simulation molculaire pour le calcul des quilibres liquide vapeur peuvent tre
divis en deux types : les mthode de calcul direct et les mthodes de calcul indirecte o le potentiel
chimique est calcul pour un intervalle de conditions et lquilibre de phase est dduit de ces valeurs. Nous
prsentons successivement les mthodes de calcul direct et indirect sachant que parmi celles-ci il en existe
quelques unes qui sappliquent pour la phase solide aussi.
a.
La mthode la plus triviale pour tudier les quilibres entre deux phases est daccomplir la
simulation avec une seule boite, les deux phases tant spares par une interface. Dans la procdure
classique, N molcules sont places dans une boite rectangulaire de dimension avec des conditions
priodiques limites appliques aux directions x et y. la phase liquide est confine au centre de la cellule
avec une phase vapeur de part et dautre par rapport la direction z [ Goujon et al., 2000 ; 2002 ].
Cette mthode prsente la caractristique attractive dtre naturelle, dans le sens que systme simul
reproduit le systme diphasique exprimental. De plus, dans le principe, cette mthode peut sappliquer
aux phases liquides de nimporte quelle densit, ce que beaucoup dautres mthodes prsentes par la suite
nautorisent pas. Cependant, elle prsente un nombre dinconvnients significatifs. En effet, les
simulations sont lentes car les systmes doivent tre larges et ncessitent de trs longue dure de
simulation pour quilibrer linterface prsente. De plus, il devient difficile de raliser un systme
diphasique stable lorsque les densits des deux phases sont similaires, au voisinage du point critique.
Enfin, le confinement des phases fluides entre des plaques parallles peut influencer les valeurs de la
densit et la pression dquilibre, de telle sorte que lon simule un fluide confin plutt que deux phases en
quilibre.
II-73
De nos jours, cette mthode parat plus approprie pour ltude de la rgion interfaciale et ltude du
comportement des fluides dans des systmes poreux, mais peut devenir, dans un future proche, plus
attractive pour le calcul des quilibres entre phase avec les progrs continus en puissance de calcul.
ii.
La temprature du systme est fixe ds le dpart ; les deux conditions dquilibre restantes sont
respectivement satisfaites par deux types de dplacement alatoires :
Transfert des particules entre les deux rgions pour assurer lgalit des potentiels chimiques.
Cependant, lnergie interne de chaque phase est stabilise par des dplacements des particules par
translation, rotation et changement de conformation des molcules dans chacune des deux rgions.
Smit et al. [ 1989 ] ont prouv que lensemble de Gibbs est un ensemble statistique formellement
quivalent lensemble canonique NVT, et que sa densit de probabilit est proportionnelle :
N!
+ N1 ln V1 + N 2 ln V2 U1 U 2
exp ln
N1! N 2!
(II-172)
(II-173)
II-74
2. Pour ltape de changement de volume durant laquelle, volume total constant, le volume de la
rgion 1 augmente de V, avec une diminution correspondante du volume de la rgion 2, le
critre dacceptation est donn par :
V + V
V V
+ N 2 ln 2
pvolume = min1, exp U1 U 2 + N1 ln 1
V1
V2
(II-174)
N 2V1
exp( U1 U 2 )
ptransfert = min1,
(N1 + 1)V2
(II-175)
Cette dernire quation peut tre gnralise pour des systmes plusieurs constituants, en remplaant N1
et N2, par respectivement N1j et N2j, pour lespce j.
La Figure II-12 montre ces diffrentes tapes.
tat
initial
Phase I
Phase II
Type de
mouvement
dplacements
internes
V1 = -V2
transfert de
particule
Condition
dquilibre
Ei minimale
P1 = P2
i,1 = i,2
Figure II-12. Principes de la simulation dans lensemble de Gibbs pour le calcul dquilibre liquide vapeur
ractif ou non.
Ces quations sont valables uniquement pour une simulation temprature, volume et nombre de
particule constants (Gibbs NVT). Pour des systmes un seul constituant pur cest la seule possibilit, car
la rgle de phase exige quune seule variable intensive (ici la temprature) soit fixe indpendamment
lorsque deux phases coexistent : la pression est alors un rsultat de la simulation.
Pour des systmes multiconstituants, en tenant compte toujours de la rgle de phase, la pression peut aussi
tre spcifie en avance, ce qui permet de raliser des simulations temprature, pression et nombre de
particule constants (Gibbs NPT). Dans ce cas, la densit de probabilit du systme est proportionnelle :
II-75
N!
exp ln
+ N1 ln V1 + N 2 ln V2 U1 U 2 PV1 PV2
N1! N 2!
(II-176)
Pour Gibbs NPT, la seule modification de lalgorithme intervient au niveau de ltape du changement de
volume, qui se fait maintenant indpendamment dans les deux phases : en effet, si le volume de la rgion 1
change de V1, et celui de la rgion 2 de V2, la probabilit dacceptation de cette tape serait donne
par :
V + V1
V + V2
+ N 2 ln 2
P (V1 + V2 )
min1, exp U1 U 2 + N1 ln 1
V1
V2
(II-177)
Dans la pratique, il faut pondrer limportance relative de chaque dplacement. Le moindre changement
de volume par exemple doit tre suivi de nombreuses translations - rotations pour rtablir la stabilit de
lnergie interne. En plus V est born en gnrant un nombre alatoire uniforme entre 0 et 1 et Vmax
un changement de volume maximal ajust pour obtenir un pourcentage dacceptation fix :
V = Vmax min(V1,V2 )
(II-178)
La mme procdure est applique pour choisir les amplitudes maximales de translation et de rotation.
Une des principales difficults de cette mthode rside dans le transfert de particules pour obtenir lgalit
des potentiels chimiques. Cest particulirement difficile dinsrer une particule et fortiori une molcule
polyatomique dans une phase dense comme la phase liquide. Une alternative est de casser le caractre
alatoire de cette insertion en cherchant les endroits prfrentiels (des espaces libres) o on peut insrer la
particule. Dtruire le caractre alatoire signifie lintroduction dun biais statistique. Un exemple trs utile
est celui du biais configurationnel introduit par Smit et al. [ 1995 ] qui consiste insrer segment par
segment la molcule considre.
iii.
Biais configurationnel
La mthode du biais configurationnel consiste insrer segment par segment une molcule dans
une phase. Supposons que toutes les molcules ont l segments et que lon veut transfrer la molcule de la
rgion 2 vers la rgion 1. Linsertion du premier segment l1 se fait alatoirement et on dfinit une fonction
de pondration w1(l1) :
E1,l1
w1 (l1) = exp
k BT
(II-179)
o E1,l1 est lnergie dinteraction du segment l1 avec le reste des molcules de la rgion 1.
La deuxime insertion, le segment l2 est ralise en slectionnant m directions possibles et en choisissant
une en fonction de lnergie de pliage et de torsion, par exemple la kime avec la probabilit :
II-76
pacc,k
E1k,l 2
exp
k T
B
=
w1 (l 2)
(II-180)
o E1k, l 2 est lnergie dinteraction du segment l2 insr dans la direction k avec le reste des molcules de
la rgion 1 et la fonction de pondration w1(l2) :
w1 (l 2) =
1
m
E1i,l 2
exp
k T
B
i =1
(II-181)
Le mme scnario est rpt pour les autres segments. Le facteur de pondration total W1 qui est
simplement lnergie dinteraction totale de la molcule insre avec le reste de la rgion est calcul par :
W1 =
w ( j)
(II-182)
j =1
N 2 V1 W1
pacc,transfert = min1,
(N1 + 1) V2 W2
(II-183)
Notons que cette quation nest valable que pour un champ de force dcrit par un potentiel datome unifi
(UA). Cependant, une formule un peu plus complique sapplique avec les atomes unifis anisotropes.
iv.
Nous avons vu que dans lensemble de Gibbs, les deux boites de simulations convergent vers deux
phases diffrentes. Cette convergence permet en effet de minimiser lnergie du systme, alors quil
apparat plus couteux de crer des interfaces. Cependant, lapproche du point critique, lnergie
ncessaire la cration de ces interfaces, c'est--dire lnergie ncessaire pour initier lquilibre de phase,
est considrablement rduite. Par consquent, il est frquent, dans ce domaine de temprature et de
pression, de voir apparatre des interfaces au sein des deux boites et mme dobserver des changes
rapides didentit entre les deux boites. Ainsi, au cours dune simulation dans lensemble de Gibbs une
temprature proche du point critique, une boite de simulation voit sa densit osciller entre la densit de la
phase liquide et celle de la phase vapeur.
II-77
Tr = 0.80
Tr = 0.94
Tr = 0.98
Figure II-13. volution de la densit de deux phases en quilibre au voisinage du point critique
La Figure II-13 illustre ce cas pour des simulations dquilibre liquide vapeur ralises dans lensemble
de Gibbs Monte Carlo pour un fluide de Lennard-Jones. Il est difficile dans ces conditions deffectuer une
prise de moyenne correcte sur un domaine stable de taille suffisante. Ce phnomne dindistinguabilit des
phases est dailleurs parfaitement observ exprimentalement.
Les diffrents auteurs utilisant lensemble de Gibbs en pris lhabitude de limiter la temprature de leurs
simulations une temprature rduite infrieure 0,95.
Il est cependant intressant dobtenir les proprits critiques des fluides tudis ; dune part parce que
pour les molcules usuelles, les donnes exprimentales correspondantes sont souvent disponibles, et nous
pourrons donc nous y comparer ; dautre part, parce que pour les molcules moins bien renseignes, ou
instables haute temprature, ces donnes sont souvent inexistantes, alors quelles sont largement utilises
en thermodynamique classique, notamment pour les quation dtat. La simulation molculaire prsente
donc un moyen pour les obtenir.
Pratiquement, on accde aux proprits critiques par extrapolation des densits de coexistence obtenues
grce lensemble de Gibbs pour des tempratures infrieures. La temprature critique est obtenue en
ajustant la loi dchelle suivante :
liq vap = A (T Tc )
(II-184)
O est une constante appele facteur dIsing qui a t value gale 0,32 pour les alcanes. Mme si une
procdure a t dcrite [ Frenkel et Smit, 1996 ], aucune tude dtaille na t ralise pour identifier le
facteur dIsing pour dautres familles chimiques, notamment pour les nitriles. La valeur de 0,32 est donc
frquemment prise titre indicatif pour valuer les coordonnes critiques.
Puis, pour dterminer la densit critique (c), nous ajustons les paramtres de la loi des diamtres
rectilignes :
II-78
liq + vap
2
= c + B (T Tc )
(II-185)
Dans les mthodes de simulation indirecte, le potentiel chimique ou lnergie libre des constituants
purs est calcul dans la simulation pour une srie de points ; lquilibre de phase est ensuite dtermin
dans lespace dtat o les deux phases possdent le mme potentiel chimique, la mme temprature et la
mme pression. La majorit des applications des quilibres entre phases sont alors plus lentes calculer
que par les mthodes de calcul directes.
Par ailleurs, les mthodes conventionnelles de Monte Carlo et de Dynamique Molculaire ne fournissent
pas directement les proprits thermiques comme lnergie libre et lentropie, mais donnent uniquement
les valeurs des proprits mcaniques comme la pression et lnergie. Dans le cas des mthodes de Monte
Carlo, ceci est d au fait que lalgorithme dchantillonnage de Metropolis utilis est conu pour permettre
dobtenir une moyenne sur lnergie de configuration du systme et non pas sur la fonction de partition ou
sur lnergie libre des configurations. Des problmes similaires sont rencontrs dans les simulations de
Dynamique Molculaire. Cest pourquoi plusieurs techniques dchantillonnage ont t proposes pour
calculer de telles proprits thermiques.
i.
Le potentiel chimique est gnralement calcul avec la technique propose par Widom [ Widom,
1963 ], dans laquelle on calcule la variation de lnergie potentielle du systme due linsertion dune
particule test. Cette approche est valable pour les deux mthodes de simulation molculaire savoir : la
Dynamique Molculaire et les simulations de type Monte Carlo, notamment la mthode de lensemble de
Gibbs (paragraphe D.3.a.ii de ce mme chapitre).
Dmontrons la relation de Widom dans lensemble canonique NVT, sachant que le raisonnement est
similaire pour les autres ensembles.
Soit un systme reprsentant une phase de volume V, temprature donne T et contenant N particules
monoatomiques de types diffrents. Nous avons dj dmontr que le potentiel thermodynamique
quivalent lensemble canonique est lnergie libre de Helmholtz F, relie au potentiel chimique du
constituant i (dans lune des deux phases liquide ou vapeur) par la relation thermodynamique classique :
i =
N i T ,V
(II-186)
l
,V , N j , j i
II-79
Dautre part, nous avons montr que le potentiel thermodynamique caractristique de lensemble
canonique est reli sa fonction de partition QNVT par lquation :
F = kBT ln Q NVT
(II-26)
i =
N i T ,V
,V v , N j , j i
ln Q NVT
= kBT
N i
Q
kBT ln N VT
Q NVT
T ,V l ,V v , N j , j i
(II-187)
Avec : N = N + 1
En sparant les contributions cintique et potentielle lnergie du systme dans lexpression de la
fonction de partition ; en remplaant les sommations discrtes par des intgrales et en introduisant la
variation de lnergie potentielle totale Ut associe linsertion dune molcule test , le potentiel
chimique se calcule nouveau par lexpression :
q (N i +1)
i
exp( (U + U ))d 3 N i dV
ext
t
(N i + 1)! V
V =0 i
i = kBT ln
N
V =V
i
q
3 Ni
i Ni i ! exp( U ext )d dV
V =0
V
V =V
= kBT ln
qi
(N i + 1) exp( U t )
(II-188)
qi = V 3 o : i = h
reprsente la longueur donde de De Broglie pour une
(2 mi kT )12
i
i = kBT ln
V
exp( U t )
(N i + 1)3i
(II-189)
En choisissant comme tat de rfrence ltat du gaz parfait correspondant une nergie potentielle totale
nulle ( U t = 0 ), le potentiel chimique rsiduel r ( = id , avec id reprsentant le potentiel chimique
du gaz idal la mme temprature et densit) dun constituant scrira sous la forme :
II-80
r = kBT ln exp( U t )
(II-190)
Cette quation a t largement employe par les mthodes de Monte Carlo et de Dynamique Molculaire
pour calculer le potentiel chimique des fluides de constituants purs et des mlanges.
Il est plus avantageux de travailler dans lensemble isobare isotherme (NPT) au lieu de lensemble
canonique. Ceci est d au fait que cet ensemble permet des fluctuations de la densit au travers du volume
du systme, alors que lensemble canonique ne permet pas de telles fluctuations. Ceci est particulirement
important lorsquon travaille dans des conditions proche de la rgion critique. Lexpression donne par
lquation scrit alors sous la forme :
V exp( U t )
r = k BT ln
ii.
(II-191)
Lorsquon dispose dun point de la courbe dquilibre entre deux phases, on peut dterminer le
reste de cette courbe sans calculs supplmentaires de lnergie libre du systme. Kofke a propos en 1993
une mthode numrique qui permet de raliser cette intgration [Kofke, 1993a, 1993b]. Dans sa forme la
plus simple, cette technique est quivalente lintgration numrique de lquation de Clausius-Clapeyron
(bien que Kofke a choisi dappeler son approche lintgration de Gibbs-Duhem).
Rappelons la drivation de lquation de Clausius-Clapeyron. Lorsque deux phases et coexistent
temprature T et pression P donnes, leurs potentiel chimique doit tre le mme. Si on change et la
temprature et la pression du systme dune quantit infinitsimale dT et dP respectivement ; la diffrence
entre le potentiel chimique des deux phases devient alors :
= (S S )dT + (V V )dP
(II-192)
S S
dP
H
=
=
dT V V TV
(II-193)
II-81
H
d ln P
=
d1 T
PV T
(II-194)
Kofke et al. ont appliqu cette mthode pour dterminer les courbes dquilibre liquide vapeur [Kofke,
1993a, 1993b] et dquilibre solide liquide pour les fluides de Lennard-Jones [Agrawal et Kofke, 1995].
Prcisons que cette mthode nest en aucun cas limite au calcul des diagrammes dquilibre pression
temprature. Cependant, bien quelle soit trs efficace pour tracer les courbes de coexistence, elle nest
pas trs robuste numriquement. En effet, les erreurs numriques dans lintgration de lquation (II-193)
peuvent tre lorigine dune dviation significative entre les valeurs calcules et les donnes
exprimentales. De plus, cest une mthode qui dpend de la qualit du point initial : une mauvaise
localisation de ce dernier conduira une estimation incorrecte de toute la courbe dquilibre. Cest pour
cette raison, quil est important de raliser quelques calculs prliminaires additionnels de lnergie libre du
systme pour fixer deux ou trois points supplmentaires sur la courbe dquilibre.
iii.
Pseudo-ensembles
II-82
iv.
Semigrand ensemble
Cette mthode est base sur lobservation quune fois que le potentiel chimique de lun des
constituants du mlange est fix ; le potentiel chimique de tous les autres constituants peut tre impos en
ralisant des mouvements alatoires de changement didentit des particules.
Rappelons le rsultat que nous avons dmontr pour lexpression du potentiel chimique dexcs en
fonction du facteur de Boltzmann correspondant au mouvement alatoire dinsertion dune particule test :
ex = kBT ln exp( U )
(II-195)
Supposons quon dsire simuler le comportement dun mlange binaire, nous devons calculer lnergie
libre de Gibbs par mole de mlange, donne en fonction de la composition par :
G (x A ) = x A A + x B B
(II-196)
G (x )
= A B
x T ,P,N
= ( A B )id + ( A B )ex
(II-197)
Connaissant les potentiels chimiques des constituants purs, il suffira de calculer uniquement la quantit
ex = ( A B )ex . On peut penser alors utiliser la mthode dinsertion dune particule test pour
calculer sparment le potentiel dexcs des deux particules A et B et faire la diffrence par la suite ; ce qui
est correct dans le principe mais ni prcis ni optimal en terme de temps de calcul. En effet, il possible de
calculer directement ex en valuant le facteur de Boltzmann correspondant un mouvement alatoire
virtuel qui consiste choisir une particule du type B du systme et de changer son identit pour quelle
devienne du type A. Le mme raisonnement que celui appliqu pour la technique dinsertion dune
particule test permet daboutir :
r = k B T ln
NB
exp U +
N A +1
(II-198)
II-83
E.
Conclusion
A travers ce chapitre, nous pouvons conclure que la simulation molculaire se prsente comme une
troisime voie dexploration du monde rel, intermdiaire entre exprience et thorie. La figure II-14
confronte ces trois approches possibles du rel [ Allen et Tildesley, 1987 ], elle met aussi en vidence les
deux utilisations principales de la simulation molculaire, savoir :
Tester la validit du modle dun systme rel en comparant les mesures exprimentales avec les
prdictions numriques.
Tester la validit dune thorie base sur le mme modle que celui qui sert la simulation
molculaire. Lexemple le plus courant est le test des quations dtat thermodynamique.
Systme
rel
MODELE
Systme
modle
EXPERIENCE
SIMULATION
THEORIE
Rsultats
exprimentaux
Rsultats
exacts pour le
modle
Prdictions
thoriques
comparaison
comparaison
Test des
modles
Test des
thories
Figure II-14 : La place de la simulation molculaire dans les diffrentes approches de la ralit.
Par ces deux techniques dexploration de lespace des phases et celui des configurations : la Dynamique
Molculaire et la mthode de Monte Carlo, la simulation molculaire est une technique dexprimentation
numrique permettant dobtenir les proprits physiqo-chimiques des systmes macroscopiques partir de
deux ingrdients indissociables :
II-84
Les interactions nergtiques lmentaires sont rassembles dans un champs de force qui prend en
compte tout ce qui contribue lnergie interne du systme : nergies liantes (vibration, longation,
torsion) et nergies non liantes (interaction de Van der Waals, interaction coulombienne, liaison
hydrogne).
La thermodynamique statistique permet de faire le lien entre les interactions lmentaires et les proprits
macroscopiques. Partout o les phnomnes nergtiques ont une place prdominante, la simulation
molculaire mrite dtre considre pour tudier et approfondir la connaissance des phnomnes au cur
des procds [ Gerbaud, 2003 ].
Un intrt particulier a t port dans ce chapitre lun de ces phnomnes : celui des quilibres liquide
vapeur, dimportance particulire pour la conception des procds. La liste des techniques de simulation
molculaire le traitant est loin dtre exhaustive, montrant la diversit et les potentialits de cette nouvelle
approche.
II-85
A.
Introduction
III-1
calcul des charges ponctuelles atomiques ont t raliss en chimie quantique avec une mthode de
fonctionnelle de la densit, une mthodologie doptimisation des paramtres de Lennard-Jones semblable
celle employe par Ungerer et al. [ 2000 ] est suivie.
In fine, la capacit prdictive de la mthode de l'ensemble de Gibbs Monte Carlo couple un modle
appropri d'interaction est dmontre pour les nitriles.
B.
III-2
consquence. Mais il est utile d'essayer de prvoir d'autres proprits qui n'ont pas t incluses dans le
procd de paramtrage. Cest mme un des avantages intrinsques de la simulation molculaire. Par
consquent, la transfrabilit de la forme fonctionnelle et des paramtres du champ de force est une
proprit importante d'un bon champ de force. La transfrabilit signifie que le mme ensemble de
paramtres peut tre employ pour modliser une srie de molcules homologues, plutt que de dfinir un
nouveau jeu de paramtres pour chaque molcule individuelle.
Par exemple, nous souhaiterions pouvoir employer pour tous les nitriles linaires le mme jeu de
paramtres de Lennard-Jones CH3, CH2 du champ de force AUA4 et dvelopper un couple de
paramtres de Lennard-Jones pour le groupe nitrile CN en essayant de reproduire les nombreuses
donnes exprimentales disponibles pour lactonitrile. Puis nous comptons prdire, sans autre
rajustement des paramtres CN, les proprits dautres nitriles linaires comme le propionitrile et le
n-butyronitrile.
La transfrabilit est clairement importante si nous voulons employer le champ de force pour faire de
vraies prdictions de donnes physico-chimiques. Nous verrons par ailleurs que ces prdictions atteignent
une prcision remarquable. Ce nest seulement que pour quelques petites molcules quun travail dtaill
de paramtrage est important (H2O, H2S ou de CO2 par exemple), faute de pouvoir bien les dcomposer
en groupes chimiques gnriques.
Remarque importante : quelle que soit leur source (donnes exprimentales ou autre champ de force
rput), les paramtres de champ de force, non dvelopps pour le champ de force qui nous intresse
(AUA4), doivent seulement tre employs comme le point de dpart dune optimisation. Ainsi, nous
avons employ comme point de dpart les paramtres des nitriles publis par Allen et Tildesley [ 1987 ]
ct des paramtres CH3, CH2 publis dans le champ de force AUA4.
2.
1
N
(X
mod
i
X iexp )
s i2
(II-127)
O si est une estimation de lincertitude associe la diffrence entre la valeur calcule et la valeur
exprimentale de la proprit (Xmod - Xexp), elle reflte la fois lincertitude exprimentale et lincertitude
III-3
statistique. Ainsi, F est considre comme une fonction des couples de paramtres de Lennard-Jones
(N,N) et (C,C) que lon veut optimiser.
Les valeurs exprimentales choisies sont la densit liquide liq, lenthalpie de vaporisation Hvap et le
logarithme de la pression de vapeur saturante ln(Psat). Les incertitudes associes aux proprits de
rfrence sont : si = 0,1 pour ln(Psat), 1 kJ.mol-1 pour Hvap et 10 kg.m-3 pour liq.
Leur choix est dict par le fait que nous souhaitons prdire des quilibres liquide vapeur :
-
La pression de vapeur saturante est une des caractristiques essentielles des proprits
dquilibre liquide vapeur puisquelle dcrit la courbe de rose qui dlimite la rgion de
coexistence des deux phases en quilibre. En outre, les mesures exprimentales sont aises et
nombreuses. Ce nest pas le cas pour la densit vapeur qui aurait pu tre une variable de
rfrence alternative.
Enfin, la densit molculaire liquide est fortement corrle aux autres proprits de la matire.
Pouvoir la prdire de faon prcise est un prrequis pour prtendre calculer prcisment dautre
proprits importantes comme la pression de vapeur saturante dans ce chapitre, mais aussi
comme dautres lon fait, des proprits mcaniques, des nergies de mlanges, viscosit ou
index de rfraction. Une densit molculaire correcte est aussi lillustration que la conformation
des molcules individuelles et leur empilement molculaire sont bien dcrits, signe dune bonne
modlisation des interactions inter et intramolculaires.
Les gradients sont calculs aprs simulation par diffrence finie, faute de disposer de leur valuation au
cours des simulations comme cela a t fait rcemment pour les olfines [ Bourrasseau et al., 2003 ].
Chaque cycle doptimisation consiste alors lancer 4 simulations pour chaque temprature pour
[(N+N,N) ; (C,C)], [(N,N+N) ; (C,C)], [(N,N) ; (C,C+C)], [(N,N) ; (C+C,C)], soit 16
simulations dans notre cas (avec = 0,05 et /kB = 10K).
3.
III-4
Un autre critre de premire importance tait la disponibilit de donnes exprimentales fiables, qui
reprsente le problme principal de n'importe quelle approche semblable. La recherche bibliographique a
montr que la banque de donnes chimiques des nitriles est pauvre. Seul l'actonitrile a t le sujet de
beaucoup d'expriences qui permettent une bonne dtermination de sa courbe d'quilibre liquide vapeur
et des proprits en phase liquide pour un large ventail de valeurs de temprature et de pression.
L'actonitrile est aussi typique d'un liquide aprotique dipolaire dont les proprits thermodynamiques,
structurales et de dynamiques des fluides ont t tudies. Son moment dipolaire lev et sa constante
dilectrique importante lui permettent de dissoudre les substances polaires et mme ioniques et font de lui
un solvant majeur en chimie et en gnie chimique.
En outre, pour les nitriles, un unique article, crit en russe de surcrot [Chakhmuradov et Guseinov, 1984],
fournit des donnes exprimentales de pression de vapeur saturante sur une plage de temprature et de
pression significative pour les molcules proponitrile et n-butyronitrile. Cest un des facteurs qui ont
motiv le choix de la molcule d'actonitrile comme composant de rfrence pour optimiser les
paramtres de Lennard-Jones. Leur transfrabilit qui constitue lobjectif ultime sera examine en
prdisant plusieurs proprits du proponitrile et du n-butyronitrile.
Les valeurs exprimentales doivent tre aussi prcises que possible. Cela nous a amen les extraire de
plusieurs publications, en allemand [ Francesconi et al ., 1975 ], anglais [ Mousa, 1981 ; Kratzke et Mller,
1985 ; Warowny, 1994 ] et russe [ Chakhmuradov et Guseinov, 1984 ].
Les grandeurs exprimentales de rfrence utilises pour loptimisation sont :
-
Bohm et al. ont utilis un modle six sites de type All Atoms (les 3 H, C de CH3, C et N de
CN) pour l'actonitrile dans des simulations de Dynamique Molculaire [ Bohm et al., 1983 ;
1984 ]. Ces auteurs ont prouv que leur modle donne une excellente description de la phase
liquide.
Jorgensen et al. ont employ un modle trois sites de centres de force de type United Atoms
(CH3, C et N) avec un potentiel intermolculaire comprenant un terme de Lennard-Jones et un
potentiel de coulomb [ Jorgensen et al., 1988 ]. Ils ont effectu des simulations de Monte Carlo
dans lensemble NPT deux tempratures (298,15K et 343,15K) la pression atmosphrique.
Les proprits liquides simules (densit et enthalpie de vaporisation) taient en accord avec les
valeurs exprimentales avec des erreurs moyennes denviron 1 2%. En outre, leur tude met
en vidence limportance des interactions lectrostatiques sur la structure de la phase liquide.
Par rapport Bohm et al. la rduction du nombre de site a permis un gain de temps calcul
valu comme permettant des calculs 3-4 fois plus rapides.
III-5
Rcemment, Hloucha et al. ont ralis une tude complte afin d'essayer d'tablir le lien entre
linformation contenue dans lnergie dinteraction obtenue par des calculs quantiques ab initio
et la prdiction de proprits macroscopiques [ Hloucha et al., 1997 ; 2000 ]. Dans la dernire
publication, ils ont ralis des calculs de simulation molculaire en utilisant un potentiel bas
sur les calculs ab initio de Bukowski et al. [ Bukowski et al., 1999 ] qui a appliqu la Symmetry
Adapted Perturbation Theory (SAPT) afin de calculer les nergies dinteraction pour des
centaines de configurations de paires de molcules dont la paire actonitrile actonitrile.
Comme dcrit dans leur papier original, Hloucha et al. [ 2000 ] ont raliss trois rgressions
diffrentes des nergies dinteractions ab initio de Bukowski et al. [ 1999 ] pour obtenir deux
potentiels de Lennard-Jones (nomms ACN1 et ACN2) et un potentiel de Buckingham modifi
(nomm ACN3). Ils ont ajout chaque fois un potentiel coulombien multi-site. Hloucha et al.
utilisent des charges ponctuelles rgresses partir de surfaces de potentiel lectrostatique
calcules au niveau ab-initio (QCISD/aug-cc-pVDZ). Les potentiels ACN1, ACN2 et ACN3
ont t utiliss dans des simulations de Monte Carlo dans l'ensemble de Gibbs Monte Carlo
pour calculer le diagramme dquilibre liquide vapeur de l'actonitrile. Les rsultats varient
fortement selon les trois modles, le modle ACN1 donnant les meilleurs rsultats. Cest
pourquoi nous comparerons plus loin nos rsultats ceux obtenus avec le modle ACN1.
4.
2. Prparation du modle de mcanique molculaire pour une molcule de nitrile. Les paramtres
de ltape 1 sont utiliss. Pour les paramtres du potentiel de Lennard-Jones des groupes CH3
ou CH2 des molcules dactonitrile (CH3CN), de propionitrile (CH3 CH2 CN) et de
butyronitrile (CH3 CH2 CH2 CN), les valeurs gnriques AUA4 sont utilises. Pour le groupe
nitrile -CN, deux couples (N,N) et (C,C) sont initialiss partir de valeurs trouves dans la
littrature de champs de force plus anciens. Les molcules sont considres comme flexibles.
La forme gnrale du potentiel est la somme des contributions intramolculaires de pliage et de
torsion et des contributions intermolculaires de Van der Waals et lectrostatiques. Les
longueurs des liaisons tant fixes, le potentiel intramolculaire dlongation tant nglig, ce
qui permet dexprimer lnergie totale par la sommation :
III-6
( ) k2 ( )
V rN =
angles
WN
(1 + cos(n ))
torsions 2
12 6 M M q q
+ 4 ij ij ij + k l
r
r
i =1 j = i +1
ij k l = k +1 4 0 rkl
ij
N
(III-1)
3. Simulations dquilibres liquide vapeur dans lensemble de Gibbs Monte Carlo NVT pour
plusieurs conditions opratoires de temprature (pour un corps pur, il nest pas possible de
fixer la pression, seules des simulations dans lensemble Gibbs Monte Carlo NVT sont
autorises par la rgle de phase).
4. Rptition de ltape 2 4 afin de minimiser lcart avec des donnes exprimentales de
rfrences reprsentatives de lquilibre liquide vapeur : la pression de vapeur saturante
ln(Psat), lenthalpie de vaporisation Hvap et la densit liquide liq. On minimise la somme des
moindres carrs pondrs pour optimiser les valeurs des deux couples (N,N) et (C,C).
5. Lorsque les deux couples (N,N) et (C,C) sont dtermines, on teste leur caractre gnrique
en les utilisant tels quels sans aucun autre ajustement pour prdire lquilibre liquide vapeur et
les proprits liquide et vapeur dautres nitriles (propionitrile, n-butyronitrile).
C.
Optimisation de la gomtrie
Loptimisation de gomtrie est faite avec le code DeMon [ St-Amant et Salahub, 1990 ; St-Amant,
1991 ] et une mthode de Fonctionnelle de la Densit (DFT) de type PD86 avec une base dorbitales
atomiques TZVP. Elle vise trouver lnergie minimale correspondant la conformation la plus stable de
la molcule. La recherche du minimum nergtique utilise une mthode de gradient conjugu. Les rsultats
correspondent la conformation stable dans le vide. Ce nest pas un obstacle significatif pour la
simulation de la molcule en phase condense puisque dans les simulations, cette conformation est
autorise changer comme le dcrit lnergie totale (quation (III-1)).
Les rsultats sont prsents dans la Figure III-1, la Table III-1 et la Table III-2 pour lactonitrile, dans la
Figure III-2, la Table III-3 et la Table III-4 pour le propionitrile et dans la Figure III-3 et la Table III-5
pour le n-butyronitrile.
III-7
a.
Actonitrile
H1
180
1,098
C1
C2
1,457
1,168
N1
H2
H3
Figure III-1. Conformation stable de lactonitrile calcule par DFT (PD86 base TZVP).
CH3CN
qMEPa
C1
0,000
0,000
0,000
-0,62
-0,40
-0,464
C2
0,000
+1,457
+0,001
+0,47
+0,15
+0,441
H1
+1,030
- 0,382
0,000
+0,21
+0,16
+0,169
H2
- 0,514
- 0,384
+0,891
+0,21
+0,16
+0,169
H3
- 0,515
- 0,377
- 0,894
+0,21
+0,16
+0,169
N1
0,000
-0,48
-0,23
-0,485
3,98
4,01
3,95
+2,626
- 0,001
Moment dipolaire (Debye) exp. 3,92
qMullikena
qab initiob.
analyse de population MEP et Mulliken avec une base dorbitale atomique TZVP
b Hloucha et al. [2000]
a
Distance ()
Exp
MM3 a
Ce travail b
C1C2
C1H1,2,3
C2 N1
1,468
1,095
1,159
1,470
1,108
1,158
1,457 ( - 0,75% )
1,098 ( +0,27% )
1,168 ( +0,77% )
Angles (deg)
Exp
MM3 a
C2C1H
109,7
110,0
HC1H
108,9
108,9
C1C2 N
-----180,0
a Goldstein et al. [1996]
b calcul DFT PD 86 avec une base dorbitale atomique TZVP
Ce travail b
110,30 ( +0,55 % )
108,63 ( - 0,25 % )
179,83
Table III-2. Distances et angles exprimentaux, optimiss par mcanique molculaire et par DFT de
lactonitrile.
III-8
b.
Propionitrile
H4
H1
1,100
1,096
C1
1,547
C2
180
C3
1,463
1,169
N1
H2
H3
H5
Figure III-2. Conformation stable du propionitrile calcule par DFT (PD86 base TZVP).
CH3CH2CN
qMEP a
qMulliken a
C1
- 1,395
- 0,313
- 0,592
- 0,23
-0,31
C2
0,000
0,000
0,000
- 0,05
-0,25
C3
0,000
0,000
+1,463
+0,08
+0,13
H1
- 2,130
+0,439
- 0,281
+0,08
+0,12
H2
- 1,750
- 1,296
- 0,258
+0,08
+0,12
H3
- 1,345
- 0,316
- 1,687
+0,08
+0,12
H4
+0,741
- 0,740
- 0,340
+0,04
+0,15
H5
+0,359
+0,981
- 0,346
+0,04
+0,15
N1
- 0,019
- 0,010
+2,632
Moment dipolaire calcul (Debye) (exp : - )
-0,12
-0,23
4,09
analyse de population MEP et Mulliken avec une base dorbitale atomique TZVP
Table III-3. Coordonnes cartsiennes (en ) et charges partielles du propionitrile dans la conformation
stable (calcul DFT).
Distance ()
Exp
MM3 a
Ce travail b
C1C2
1,548
1,533
1,547 ( - 0,06 % )
C2C3
1,474
1,473
1,463 ( - 0,74 % )
C1H
-------
1,110
1,096
C2H
1,091
1,113
1,096 (+ 0,46 % )
C3 N1
1,157
1,158
1,169 (+ 0,77 % )
Angles (deg)
Exp
MM3 a
Ce travail b
C1C2C3
110,3
110,5
112,48 (+ 1,97 % )
C3C2H4,5 ___
109,3
108,15
C1C2H4,5
-------
110,0
110,70
C2C1H1,2,3
-------
111,5
110,58
H4C2H5
109,2
107,3
106,38 ( - 2,58 % )
C2C3 N
-------
180,0
178,93
a
b
Table III-4. Distances et angles exprimentaux, optimiss par mcanique molculaire et par DFT du
propionitrile.
III-9
c.
n-butyronitrile
H4
H1
H6
1,099
1,097
C1
1,535
C2
1,551
1,102
C3
1,460
180
C4
1,169
N1
H2
H3
H5
H7
Figure III-3. Conformation stable du n-butyronitrile calcule par DFT (PD86 base TZVP).
CH3(CH2)2CN
C1
+0,189
- 1,392
- 2,128
-0,33
C2
+0,150
- 1,426
- 0,594
-0,16
C3
0,000
0,000
0,000
-0,25
C4
- 0,037
0,000
+1,460
+0,12
H1
+0,286
- 2,402
- 2,546
+0,11
H2
- 0,729
- 0,948
- 2,538
+0,11
H3
+1,039
- 0,798
- 2,493
+0,11
H4
- 0,686
- 2,046
- 0,242
+0,12
H5
+1,065
- 1,883
- 0,192
+0,12
H6
+0,835
+0,641
- 0,326
+0,14
H7
- 0,976
+0,472
- 0,373
+0,14
N1
- 0,063
- 0,031
+2,628
Moment dipolaire calcul (Debye) (exp : - )
-0,23
qMulliken a
4,21
Toutes les conformations calcules par DFT sont en excellent accord avec les conformations optimales
calcules en mcanique molculaire avec le champ de force MM3 et les dterminations exprimentales
cites dans Goldstein et al. [ 1996 ]. Lerreur maximale est gale 0,77% pour les distances et 2,58% pour
les angles.
Pour ces diffrentes molcules appartenant la mme srie homologue des nitriles linaires, on constate
sur la Figure III-1, Figure III-2 et Figure III-3 lhomognit des distances interatomiques et des angles.
Ainsi la distance CN optimale est 1,168 ou 1,169 (1,157 1,159 exprimentalement) et langle C
CN proche de 180, soit un alignement parfait des trois atomes.
III-10
2.
Charges ponctuelles
Le champ de force AUA4 dcrit les interactions lectrostatiques par un potentiel de Coulomb et
donc ncessite des charges ponctuelles. Les calculs de mcanique quantique raliss par DFT (thorie de la
fonction de densit) dterminent des surfaces de potentiel lectrostatiques partir desquelles une analyse
de population peut tre effectue (chapitre II). Cette analyse de population consiste attribuer des
centres de charges (dans notre cas chaque noyau atomique) la fraction dlectron ayant une probabilit de
rsidence importante sur chacun dentre eux. Deux analyses de population de charges ponctuelles ont t
ralises aprs des calculs quantiques de DFT employant une base dorbitale atomique TZVP et une
fonctionnelle PD86 (code de calcul DeMon) :
-
lanalyse de Mulliken qui est une analyse simple ne tenant compte que du diamtre de Van der
Waals de chaque atome est fortement dpendante de la base dorbitales atomique utilise.
Pour lactonitrile, lanalyse de population MEP (Molecular Electrostatic Potential) donne des charges
suprieures en valeur absolue celle de lanalyse de Mulliken. Du fait de la proportionnalit du potentiel
coulombien vis--vis des charges partielles, linteraction lectrostatique calcule avec les charges
ponctuelles MEP sera donc suprieure celle pour les charges ponctuelles Mulliken distance fixe. Dans
la Table III-1, les charges ab-initio calcules par Hloucha et al. [ 2000 ] se situent entre lanalyse de Mulliken
et lanalyse MEP. Lazote, lectrongatif, se voit attribuer des charges partielles ngatives, de mme que le
carbone mthylique. Par ailleurs, le carbone du nitrile, situ entre les deux, porte des charges positives. Le
groupe mthyle CH3 est globalement neutre pour les distributions MEP et Hloucha et al. [ 2000 ] alors
quil est lgrement positif pour la distribution Mlliken.
Pour le propionitrile, les charges MEP sont cette fois infrieures aux charges de Mulliken. En outre, elles
sont trs diffrentes des charges MEP pour lactonitrile sur les mmes atomes. De mme CH3 est
globalement neutre, CH2 lgrement positif et CN lgrement ngatif avec un azote lectrongatif. Cest
le mme constat pour la molcule n-butyronitrile sauf que le groupe CN est encore plus ngatif.
Discussion
Le choix dune distribution ou dune autre est toujours sujet dbat et reste important puisquil fixe
de facto limportance des interactions lectrostatiques mais aussi parce que les paramtres de
Lennard-Jones qui vont tre rgresss dpendront de ce choix. Ces analyses de population ne sont
que des modles imparfaits de linformation contenue dans une surface de potentiel lectrostatique
qui dpend elle aussi des approximations ralises lchelle quantique au travers du choix de la
base dorbitale atomique. Plusieurs critres peuvent tre regards pour faire un choix partir des
valeurs fournies dans les Table III-1, Table III-3 et Table III-5 :
III-11
Pour lactonitrile, cest le cas pour les trois analyses de population MEP et Mulliken et celle de
Hloucha et al. [ 2000 ]. En consquence, dautres arguments, comme lutilisation dans un champ de
force gnrique, doivent alors tre utiliss afin de prfrer lune ou lautre des analyses. Les calculs
prdictifs des proprits des nitriles prsents ci-dessous contribueront prfrer lanalyse de
population de Mulliken, pourtant la plus simple et la plus criticable (voir chapitre II).
Pour complter la discussion, nous avons trs rcemment ralis dautres calculs de conformation
optimale avec dautres bases des mthodes de DFT en utilisant le logiciel Gaussian 03. Des analyses
de population de Mulliken et NBO ont t faites et sont prsentes en annexe. Les rsultats des
analyses de population de Mulliken sont trs diffrents de ceux prsents dans la Table III-1,
confirmant lextrme dpendance de cette analyse au choix de la base et de la fonctionnelle, mais
ajoutant notre indcision quant au choix dun jeu de charge. En revanche, lanalyse de population
NBO savre insensible au choix de la taille de la base dorbitale atomique et donne des valeurs
voisines mais diffrentes de celles des analyses de population MEP et de Hloucha et al. [ 2000 ] pour
lactonitrile. En outre les charges NBO sont cohrentes entre les diffrents nitriles et vis--vis de
llectrongativit des atomes. Cest un rel progrs pour lutilisation raisonne en mcanique
molculaire de charges ponctuelles. Ces rsultats trop rcents nont pas t exploit dans ce travail
de thse mais constituent une piste explorer.
Remarquons que chaque atome, notamment les atomes hydrognes, possde une charge ponctuelle
associe et interviendra donc dans le calcul du potentiel de Coulomb alors que pour le potentiel de
Lennard-Jones, des groupes datomes CH3 et CH2 sont considrs comme des centres de force
uniques.
Remarquons aussi que dans les conditions des simulations, plusieurs facteurs vont influencer la
surface de potentiel lectrostatique et par consquent la distribution des charges ponctuelles :
-
Des phnomnes de polarisation induite peuvent survenir, surtout en phase condense o les
molcules sont proches. De nouveau, un calcul ractualis des charges serait ncessaire, non
plus sur une molcule mais sur un large ensemble de molcules.
La mthode quantique que nous utilisons est incompatible avec ces ractualisations vu son
temps calcul dj important pour le calcul dune molcule. Dautres mthodes destimation
des charges permettent une ractualisation [ Leach, 1996 ] mais, bass sur des concepts
dlectrongativit des atomes, elles ne fournissent pas une distribution de charge reposant
III-12
sur des bases thoriques aussi solides que les analyses de population des surfaces de potentiel
lectrostatiques calcules en mcanique quantique.
In fine, il nexiste pas de mthode universelle dans la littrature pour dterminer les charges ponctuelles.
3.
Les valeurs des constantes intramolculaires de pliage et de torsion utilises sont prsentes dans la
Table III-6. Ces valeurs issues de la littrature ont t rgresses sur des surfaces de potentiel de pliage ou
de torsion dtermines partir de calculs de mcanique quantique. Nous avons galement calcul la
constante de pliage XCN par le champ de force MM1 en utilisant le logiciel Hyperchem et nous avons
trouv Kbend = 27688 K.rad-1 et 0 = 178,8, en accord avec les valeurs de la Table III-6.
Prcisons que le potentiel de torsion nest pris en compte que pour la chane CH3 CH2 CH2 C du
n-butyronitrile tant donn que les trois atomes CCN sont parfaitement aligns.
Pliage
CH3
15,030
CH2
14,030
C (C N)
12,011
N (C N)
14,006
C C C
CCC
CN
Torsion
C CH2 CH2 C
0 (deg)
112,0 a
Kbend (K)
62500 a
0 (deg)
108,8 b
Kbend (K)
69500 b
0 (deg)
180,0 b
Kbend (K)
24300 b
a0 (K)
1001,35 c
a1 (K)
2129,52 c
a2 (K)
-303,06 c
a3 (K)
-3612,27 c
a4 (K)
2226,71 c
a5 (K)
1965,93 c
a6 (K)
-4489,34 c
a7 (K)
-1736,22 c
a8 (K)
2817,37 c
Toxvaerd [ 1997 ]
Table III-6. Masses molaires des groupements chimiques et les constantes intramolculaires de pliage
et de torsion du champ de force AUA4.
III-13
En conclusion partielle, les calculs prliminaires sur les molcules, raliss en mcanique
quantique par la thorie de la fonction de la densit (DFT) sont en accord avec les donnes
conformationnelles exprimentales et fournissent des valeurs pertinentes des charges
ponctuelles et des constantes intramolculaires qui sont requises par le champ de force de
mcanique molculaire utilis dans les simulations molculaires.
Le choix du jeu de charge ponctuelle est impossible a priori. Les simulations molculaires
seront prcieuses pour faire ce choix.
D.
Simulations molculaires
1.
Mthodes dchantillonnage
la mthode de lensemble de Gibbs Monte Carlo NVT intgrant un biais configurationnel pour
calculer les quilibres entre phase,
la mthode de Monte Carlo dans lensemble isobare isotherme NPT pour calculer les
proprits monophasiques, notamment en phase liquide.
Comme indiqu prcdemment, la mthode de lensemble de Gibbs Monte Carlo NPT nest pas utilisable
pour un constituant seul comme cest notre cas puisquelle viole la rgle de phase. Ces mthodes sont
programmes dans le code de calcul GIBBS dvelopp par lInstitut Franais du Ptrole (IFP) et le
Laboratoire de Chimie Physique (LCP) dORSAY (voir annexe II).
Lincertitude statistique est rapporte dans toutes les tables sous forme de la notation suivante : une
densit vapeur gale 12,25 81 kg.m-3 (table III-8) indique que la valeur moyenne est gale 12,25 kg.m-3 et
que lerreur statistique est gale 0,81 kg.m-3. Cette erreur est dtermine partir de la variance par blocs
non corrls selon la procdure dcrite dans Allen et Tildesley [ 1987 ] et rsume dans le chapitre II.
b.
Nombre de molcules
Les simulations ont t effectues avec un nombre total de 200 molcules. Il reprsente un bon
compromis entre la taille du systme et le temps calcul tout en permettant dutiliser des tailles de boites
suffisantes. En effet, en considrant 50, 200 et 400 molcules dans une boite cubique darrte de 28
pour des simulations NPT en phase liquide condense T = 325,809 K et P = 4,1317 MPa avec les
charges MEP, nous obtenons les rsultats prsents dans la Table III-7. Le nombre de configurations
pour les trois simulations varie de 2.106 8.106 selon le nombre total de molcules pour permettre une
meilleure exploration de lespace de configurations. Cependant, la moyenne est prise partir de 1.106
configurations aprs stre assur de la stabilisation de lnergie du systme.
III-14
On constate quavec 50 molcules, les valeurs moyennes de densit et de pression sont loin des valeurs
exprimentales. En outre, lcart type est important, compar aux autres simulations.
Taux de transfert
Prise de
moyenne
Dure/10+6 configurations
(min)
liq (kg.m-3)
Pmoy (MPa)
2)10+6
750,41
762,52 3,18
4,1317
2,042 1,927
17
200
(1 - 4)10+6
754,48 0,20
3,872 579
124
400
(1 - 8)10+6
753,20 0,88
4,282 342
368
valeur exprimentale
50
(1 -
La Table III-8 illustre linfluence de la distribution statistique des mouvements de Monte Carlo dans
lensemble de Gibbs Monte Carlo sur les rsultats.
Distribution (%)
transl/rot/int/vol/tranfert
liq (kg.m-3)
vap (kg.m-3)
valeur exprimentale
613,37
11,92
10/10/10/0,5/69,5
612,90 1,49
12,88 78
20/20/10/1/49
612,18 3,14
30/30/29/1/10
32/32/30/1/5
Pliq (kPa)
-
Taux de
transfert
Pvap (kPa)
Dure
729,42
3331,8 1082,2
1016,5 58,1
2%
10h
12,25 81
219,7 936,2
968,8 61,8
2%
9h
606,86 2,20
11,12 82
259,6 1214,7
888,1 50,9
2%
8h
608,70 2,12
11,35 49
1126,7 674,7
864,4 36,8
2%
7h15
Table III-8. Influence de la distribution statistique des mouvements de Monte Carlo dans lensemble de
Gibbs Monte Carlo NVT 433,15 K avec les charges Mlliken.
On constate que :
-
le taux dacceptation des mouvements de transfert de particules entre les boites est voisin de
2% quelle que soit la distribution choisie. Cette valeur est similaire ce qui est rapport en
gnral dans la littrature pour les simulations dans lensemble de Gibbs Monte Carlo.
La dure des simulations est rduite lorsque le pourcentage de tentatives de transfert diminue,
ce qui est normal puisque les tentatives de transfert utilisent un biais configurationnel et sont
donc coteuses en temps de calcul.
III-15
avec 49% ou 10% de tentatives de transfert, les pressions des deux phases sont gales aux
incertitudes prs, mais les fluctuations de la pression liquide sont trs leves.
Avec 5% de tentatives de transfert, les pressions des deux phases sont gales aux incertitudes
prs. Les fluctuations sont aussi rduites, notamment du fait du plus grand nombre de
tentatives de translation, rotation et relaxation interne qui rendent la phase plus stable en
rduisant son nergie propre.
La meilleure distribution entre les mouvements est la dernire, imposant environ 5% de tentatives de
transfert. Lavantage est net sur la dure des simulations tout en garantissant un taux de transfert russi
convenable. Quant aux valeurs moyennes des densits et pressions, une distribution accordant trop
dimportance aux tentatives de transfert conduit des fluctuations importantes de la pression de la phase
liquide, voire des valeurs errones. Une bonne procdure de choix de la distribution des mouvements est
de dmarrer avec un taux de transfert dau moins 50% puis de le diminuer tant que le taux dacceptation
effectif est satisfaisant (> 2%). En plus, et notamment pour de faibles valeurs de taux de transfert (5 ou 10
%), on sassurera galement que le nombre de transfert accept est suffisant.
iii.
Le potentiel de Van der Waals et le potentiel lectrostatique ayant une longue porte, nous avons
employ des corrections longues distances classiques, savoir avec un rayon de coupure gal la moiti
de la largueur de la bote de simulation [ Allen et Tidldesley, 1987 ]. Lutilisation de tels cut-off est un
moyen simple de ne pas considrer la contribution dune interaction longue porte au-del dun certain
rayon autour du centre de force considr. Cependant il existe dautres corrections moins drastiques,
thoriquement plus correctes mais qui rallongent de beaucoup les temps de calcul (sommation Ewald, voir
chapitre II). Leur impact est prsent dans les rsultats.
iv.
La sommation Ewald
III-16
Pvap (kPa)
Dure
exp
613,37
11,92
729,42
Sans Ewald
607,90 1,75
12,00 60
Avec Ewald
609,58 1,40
12,25 95
3049,6 512.8
exp
583,19
17,59
Sans Ewald
589,17 1,42
17,64 45
21h
Avec Ewald
590,07 1,81
18,26 43
3668,4 582.0
31h
Une simple comparaison permet dobserver que lintroduction de cette technique de calcul ici namliore
pas la qualit des rsultats, notamment la pression dquilibre de la phase liquide, mais comme on sy
attendait augmenter considrablement le temps de calcul.
Ceci a motiv notre dcision de travailler avec un rayon de coupure classique (gale la demi taille de le
boite) pour les deux potentiels intermolculaires : de Lennard-Jones et de Coulomb.
v.
Le nombre de configurations chantillonnes lors des simulations a fait lobjet de plusieurs essais
prliminaires. In fine, pour lactonitrile, il est en moyenne de 5.106 configurations pour les simulations
dquilibre liquide vapeur ralises dans lensemble de Gibbs Monte Carlo NVT (moyenne entre 1 et
5.106) et de 2.106 configurations pour les simulations de la phase liquide par lensemble isotherme - isobare
(NPT) (moyenne entre 0,2 et 2.106). Pour le propionitrile, il est en moyenne de 7.106 configurations dans
les simulations NVT de lensemble de Gibbs Monte Carlo (moyenne entre 2 et 7.106) et de 2.106
configurations pour les simulations NPT de la phase liquide (moyenne entre 1 et 2.106), et enfin pour le nbutyronitrile il est en moyenne de 9.106 configurations dans les simulations de lensemble de Gibbs Monte
Carlo NVT (moyenne entre 2 et 9.106). La longueur des simulations a t augmente pour le propionitrile
et le n-butyronitrile afin de maintenir le mme degr de convergence et de prcision que celui de
lactonitrile vu que ces deux molcules possdent plus de degr de libert.
vi.
La conformation des molcules est diffrente selon les simulations. En effet, sappuyant sur la
valeur de la constante de pliage de langle XCN gale 180 (Table III-6) qui indique un alignement
parfait, les simulations avec les charges ponctuelles MEP ont t ralises sans cette contribution de pliage.
Ainsi, la molcule dactonitrile devient rigide et les autres nitriles semi-flexibles, nayant comme degr de
libert de pliage que celui concernant les atomes de carbone CCC.
III-17
Avec les charges Mulliken, tous les degrs de libert ont t pris en compte et les trois molcules sont
donc flexibles.
c.
Les simulations ont t ralises avec deux jeux de charges ponctuelles, issus des analyses de
population MEP et Mulliken. Les simulations avec les charges de Mulliken ont t ralises aprs celles
avec les charges MEP pour des raisons expliques aprs.
Avec les charges ponctuelles MEP ; pas de pliage de langle XCN :
-
simulations dquilibre liquide vapeur dans lensemble de Gibbs Monte Carlo NVT deux
tempratures (433,15 K et 453,15 K) pour loptimisation puis cinq autres tempratures (Table
III-12 et Table III-13),
simulations de la phase liquide dans lensemble de Monte Carlo NPT deux tempratures
(273,150 K / 3,700 kPa sur la courbe de saturation et 298,405 K / 24,714 MPa - phase liquide
condense) (Table III-12).
simulations de la phase liquide condense dans lensemble de Monte Carlo NPT sept autres
couples de temprature et pression (Table III-14).
simulations dquilibre liquide vapeur dans lensemble de Gibbs Monte Carlo trois
tempratures (433,15 K, 443,15 K 453,15 K) pour loptimisation puis cinq autres
tempratures (Table III-18 et Table III-19).
simulations de la phase liquide condense dans lensemble de Monte Carlo NPT pour trois
couples de tempratures et pressions (Table III-20).
d.
trois centres de force de Lennard-Jones : le groupe CH3, latome de carbone du groupe nitrile
CN et latome dazote du groupe nitrile CN.
Pour les autres nitriles (propionitrile et n-butyronitrile), le groupe CH2 est rajout ainsi que trois centres de
charge sur chaque atome de CH2. Les paramtres des groupes CH3 et CH2 pour le champ de force AUA4
sont rappels dans la Table III-10.
III-18
()
/kB (K)
()
CH3
3.6072
120.15
0.21584
CH2
3.4612
86.29
0.38405
Table III-10. Les paramtres nergtiques de Lennard-Jones des groupes CH2 et CH3 dans le champ de
force AUA4 [ Ungerer et al., 2000 ]
Les charges ponctuelles sont calcules au pralable par analyse de population issue de calculs de
mcanique quantique. Deux jeux de charges issus des analyses de populations MEP et Mulliken sont
compars.
Les constantes harmoniques du potentiel intramolculaire sont fournies dans la Table III-6, en prenant en
compte le pliage de langle XCN uniquement pour les simulations avec les charges Mulliken.
Pour le groupe CH3, les paramtres de Lennard-Jones sont pris gaux ceux du champ de force AUA4.
Dtermins pour des alcanes, ils sont ici considrs comme suffisamment gnriques pour tre employs
avec dautres classes de molcules. Les rsultats prsents ci-dessous conforteront leur gnricit.
Pour le groupe nitrile CN et avec les charges ponctuelles MEP, la valeur initiale des paramtres de
Lennard-Jones est fixe gale celle de Allen et Tildesley (C/kB = 51.20 K, N/kB = 37.30 K, C = 3.35
, N = 3.31 ). Les simulations avec ces valeurs ont conduit une valeur de la fonction objectif F =
4,76. Elle reprsente la borne suprieure de F avant optimisation des paramtres. Les rsultats ci-dessous
montre quavec les valeurs optimises des paramtres, on atteint une fonction objectif F = 1,35.
Pour le groupe nitrile CN et avec les charges ponctuelles Mulliken, la valeur initiale des paramtres de
Lennard-Jones est diffrente (C/kB = 100.0 K, N/kB = 140.0 K, C = 3.30 , N = 3.35 ) et choisie en
fonction de lexprience acquise lors des simulations avec les charges ponctuelles MEP, simulations
ralises avant celles avec les charges ponctuelles Mulliken.
e.
Dans les simulations dquilibre entre phase, la pression moyenne a t prise gale celle de la
pression de la phase vapeur. En effet, la pression de la phase liquide fluctue significativement du fait de la
densit plus importante. Dautre part, un ajustement prcis de loccurrence des transferts entre boites doit
tre fait afin de limiter les fluctuations de la pression de la phase liquide et davoir concordance entre la
valeur de pression moyenne de la phase liquide et de la phase vapeur.
L'enthalpie molaire de vaporisation a t simplement calcule comme diffrence entre les enthalpies
molaires moyennes des phases liquide et vapeur. La densit liquide moyenne a t dtermine directement
comme le rapport de la masse moyenne de la bote liquide et de son volume moyen.
III-19
Lincertitude statistique calcule selon la mthode de variance par bloc est fournie dans les tables de
rsultats. En moyenne elle est gale 5%, 2% et 0.5% pour la pression, lenthalpie de vaporisation et la
densit liquide respectivement, dans la gamme favorable des tempratures rduites pour la mthode
d'ensemble de Gibbs Monte Carlo (0.60 < T/Tc < 0.95), c.--d., loin des tempratures rduites leves
(proche de la temprature critique) o il est difficile de distinguer les deux phases, et galement loin des
basses tempratures o la mthode d'ensemble de Gibbs Monte Carlo procde un chantillonnage des
configurations de mauvaise qualit en raison dun trop faible taux de transfert entre la boite liquide et la
bote vapeur.
2.
Simulations molculaires
a.
Comme indiqu prcdemment, les simulations avec les charges ponctuelles MEP ont t ralises
sans tenir compte du pliage de langle XCN qui est donc considr comme linaire et fixe. Ainsi, la
molcule dactonitrile devient rigide et les autres nitriles nont comme degr de libert de pliage que celui
concernant les atomes de carbone CCC.
i.
La Table III-11 montre lvolution de la fonction objectif et des valeurs des paramtres de LennardJones du groupe nitrile au cours de leur optimisation avec les charges ponctuelles MEP.
N (C N)
C (C N)
Initialisation
1re optimisation
2nde optimisation
/kB (K)
()
/kB (K)
()
51.200
51.200
50.677
3.350
3.350
3.504
37.300
65.470
65.470
3.310
3.307
3.307
F
4,76
2,33
1,35
Table III-11. Evolution des paramtres de Lennard-Jones du groupe nitrile. Charges ponctuelles MEP.
Lordre de grandeur de ces valeurs est critiqu plus loin lors de lvaluation de la transfrabilit des
paramtres. On notera la rduction significative de la fonction objectif au fur et mesure des cycles
doptimisation, ce qui signifie que la procdure doptimisation remplit son rle.
Les rsultats des simulations dquilibre liquide vapeur dans lensemble de Gibbs Monte Carlo
(433,15 K et 453,15 K) et des simulations de la phase liquide dans lensemble de Monte Carlo NPT
(273,15 K et 298,00 K) ayant servies lors de loptimisation des paramtres sont prsents dans la Table
III-12.
III-20
T(K)
Conditions
273,15
NPT
(3,7 kPa)
433,15
298,40
453,15
GEMC
NPT
(24,7 MPa)
GEMC
Valeurs
exprimentales
Point initial
1re optimisation
2nde
optimisation
803,00a
800,62 0,296
811,40 1,046
799,25 0,467
Hvap
34176,25b
32219,97 5,724
34541,97 1,070
34838,97 1,939
613,3691c
573,54 6,493
619,02 0,921
619,15 0,942
Hvap
Ln (Psat)
25120,8b
19807,6 21,151
23647,55 5,865
24360,64 3,026
13,5007c
14,2866 5,821
13,8663 2,708
13,7751 2,032
795,97d
791,02 0,622
808,03 1,515
795,03 0,118
Hvap
33039,35b
31586,94 4,396
34126,94 3,292
34303,94 3,828
583,1970c
527,42 9,564
581,55 0,282
586,77 0,613
Hvap
23465,86b
16456,1 29,872
21363,76 8,958
22547,03 3,916
Ln (Psat)
13,8934c
14,6636 5,544
14,2888 2,846
14,2149 2,314
4,76
2,33
1,35
Type de
donnes
Table III-12. Evolution des principales grandeurs physico-chimiques de lactonitrile avec les charges
MEP en fonction des cycles doptimisation des paramtres de Lennard-Jones du groupe nitrile.
Lerreur relative moyenne finale calcule sur les grandeurs de la Table III-12 est infrieure 2%, une
excellente valeur confirme par les trs faibles carts relatifs des grandeurs prises individuellement. Les
carts les plus importants sont constats pour la pression de la phase vapeur.
ii.
Le cycle doptimisation achev, le champ de force ainsi complt est utilis pour prdire dautres
proprits de lactonitrile. Les Table III-13 et Table III-14 prsentent des simulations complmentaires
dquilibre liquide vapeur dans lensemble de Gibbs Monte Carlo et de la phase liquide condense dans
lensemble de Monte Carlo isobare isotherme (NPT).
III-21
Vap (kg.m-3)
T (K)
exp.a
348,15
350,00
423,15
433,15
443,15
448,15
453,15
463,15
473,15
500,00
520,00
523,15
----9,78
11,92
14,45
--17,59
20,31
---------
exp b
1,25
--9,88
----15,22
----22,5
68,0
Liq (kg.m-3)
Simule
--2,118 347
--14,609 977
17,458 5541
--21,956 1,494
23,101 837
--50,804 1,601
74,315 2,408
---
exp.a
----627,98
613,37
599,90
--583,19
567,15
---------
exp b
727,5
--635,0
----598,0
----556,2
----436,0
Pression(kPa)
Simule
Pliq
Psat exp.a
Psat
Table III-13. Donnes dquilibre liquide vapeur exprimentales et simules dans lensemble de Gibbs
Monte Carlo. Actonitrile avec charges ponctuelles MEP.
T (K)
287,987
298,405
325,809
327,577
398,150
410,140
435,196
464,279
a
P (MPa)
45,836
24,714
4,1317
60,236
36,087
24,826
7,5078
10,414
liqsim (kg.m-3)
818,97 2,03
795,03 1,34
749,30 2,79
787,82 2,13
715,25 3,08
679,26 2,81
609,65 3,23
595,25 2,93
liqexp (kg.m-3) a
819,57
795,97
750,41
794,06
710,59
685,00
625,75
593,80
Erreur (%)
0,07
0,12
0,15
0,79
0,66
0,84
2,57
0,24
Table III-14. Densits du liquide simules et exprimentales ( Kratzke et Mller [1985] ). Actonitrile avec
charges ponctuelles MEP.
La Figure III-4, la Figure III-5 et la Figure III-6 prsentent ces valeurs respectivement sous formes de
courbes de coexistence liquide vapeur, de diagramme de Clapeyron pression densit et de courbe de
pression de vapeur saturante. Les erreurs statistiques associes aux valeurs calcules ne sont pas
distinguables des points reports.
III-22
Francesconi 75
Warowny 94
simulations (charges MEP)
point critique estim (charges MEP)
simulations (charges Mulliken)
point critique estim (charges Mulliken)
600
Temprature (K)
550
500
450
400
350
300
250
0
200
400
600
800
Densit (kg.m-3)
Figure III-4. Courbe de coexistence liquide vapeur de lactonitrile. Valeurs calcules avec deux jeux de
charges ponctuelles (MEP & Mulliken) et donnes exprimentales de Francesconi et al. [ 1975 ] et de
Warowny [ 1994 ].
La Table III-13 et la Figure III-4 montrent lcart entre les densits de la phase vapeur simules et
exprimentales lorsque la temprature augmente. En revanche, les densits liquides sont en parfait accord
avec les densits exprimentales sur le diagramme densit - temprature.
Pourtant, sur la Figure III-5, ces mmes points se dtachent de la courbe de coexistence pression densit
exprimentale. Cela traduit le fait que les pressions dquilibre calcules (qui sont prises gales la pression
dans la phase vapeur) sont suprieures aux pressions exprimentales. Lcart est dailleurs systmatique et
suggre que son origine proviendrait dune surestimation systmatique de la pression de vapeur saturante
plutt que dune survaluation de lampleur des interactions dans la phase liquide. Reproduire les pressions
simules de la phase liquide sur le graphe ne suggre rien du fait des fluctuations trs importantes de la
pression dans la phase liquide.
Lcart constat sur la pression vapeur reste cependant raisonnable comme le montre la Figure III-6 :
le modle que nous proposons se compare favorablement avec les rsultats simuls par Hloucha et al.
[ 2000 ] avec son meilleur modle ACN1. Rappelons que ce modle ACN1 est issu de la rgression des
paramtres de Lennard-Jones partir dnergies dinteraction calcules ab-initio. Leur dmarche trs
rigoureuse dans son approche thorique donne pourtant de moins bons rsultats que notre modle.
III-23
100000
pression [kPa]
10000
1000
10
0
100
200
300
400
500
600
700
800
-3
densit [kg.m ]
Figure III-5. Diagramme de Clapeyron de lactonitrile. Valeurs calcules avec charges ponctuelles MEP,
valeurs exprimentales [ Francesconi, 1984 ; Kratzke et Mller, 1985 ; Warovny, 1994 ].
10000
Warowny 94
Chakhmuradov 84
Simulations (charges Mulliken)
ACN1 model of Hloucha
1000
0,003
0,0028
0,0022
1/500
0,002
Pression (kPa)
100
0,0018
1/Temprature [K]-1
Figure III-6. Courbe de pression de vapeur saturante de lactonitrile. Valeurs calcules avec deux jeux de
charges ponctuelles (MEP & Mulliken), valeurs exprimentales [ Chakhmuradov et Guseinov, 1984 ;
Warovny, 1994 ] et valeurs calcules par Hloucha et al. [ 2000 ] avec le modle ACN1.
III-24
Les enthalpies de vaporisation sont reportes sur la Figure III-7. On remarque lexcellent accord entre les
valeurs calcules avec les charges MEP et la corrlation DIPPR base sur des donnes exprimentales.
39000
36000
Hvaporisation (J/mol)
33000
30000
27000
24000
21000
18000
15000
12000
200
250
300
350
400
450
Temprature (K)
500
550
600
Figure III-7. Enthalpie de vaporisation des nitriles. Valeurs calcules avec deux jeux de charges
ponctuelles (MEP & Mulliken), valeurs exprimentales [ NIST, 2004 ] et corrlation.
iii.
Lestimation du point critique selon la mthode prsente dans le chapitre II est reporte sur la
Figure III-4 et dans la Table III-15. La pertinence dutilisation de la mthode dIsing est vrifie en
estimant les coordonnes critiques partir des valeurs exprimentales de Francesconi et al. [ 1975 ] et de
Warowni [ 1994 ]. On note une surestimation systmatique de la densit critique avec la mthode dIsing.
Cette constatation nous permet daffirmer que lestimation laide de nos rsultats simuls est en excellent
accord avec le point critique exprimental dtermin par Francesconi et al. [ 1975 ].
a
b
exp.a
Estimation calculs
avec MEP
TC (K)
547,85
547,50
548,42
543,28
C (kg.m-3)
237,1
242,78
241,61
247,56
Francesconi et al [ 1975 ]
Warowni [ 1994 ]
Table III-15. Estimation des coordonnes critiques de lactonitrile avec charges MEP.
III-25
iv.
La transfrabilit des paramtres de Lennard-Jones bass sur les charge MEP est value en les
employant tels quels pour calculer des proprits du propionitrile :
a Chakhmuradov
T (K)
P (MPa)
289,96
290
289,96
289,96
300
314,02
314,02
314,02
320
342,62
342,62
342,62
342,62
350
368,19
368,19
390,7
390,7
400
417,33
417,33
0,1
1,0
5,0
50,0
50,0
50,0
30,0
40,0
35,0
30,0
40,0
10,0
20,0
15,0
10,0
20,0
20,0
30,0
25,0
20,0
30,0
liqsim (kg.m-3)
756,74 2,14
786,9
815,4
802,07 2,35
796,2
783,0
789,6
757,48 1,41
758,4
766,0
740,9
750,2
688,38 1,87
715,2
725,9
704,0
715,1
640,21 1,90
676,3
689,9
et Guseinov [ 1984 ]
Table III-16. Comparaison des densits liquides obtenues par simulation pour le propionitrile, corrles et
exprimentales. Paramtres de Lennard-Jones optimiss avec les charges ponctuelles MEP.
Equilibre liquide vapeur par simulation dans lensemble de Gibbs Monte Carlo.
En dpit des efforts importants consentis, aucune simulation na donn de rsultat. Plus prcisment,
nous navons pas trouv de conditions initiales des boites conduisant la sparation des deux phases.
Les densits liquides simules dans lensemble isobare isotherme NPT du propionitrile sont compares
avec les donnes exprimentales de Chakhmuradov et Guseinov [ 1984 ] dautres tempratures et
pressions (Table III-16).
Un examen attentif de ces rsultats sur le diagramme de Clapeyron pression densit de la Figure III-8
indique un dcalage important par rapport aux isothermes exprimentales notamment pour les
III-26
tempratures 320 K, 350 K et 400 K. Dautre part, une comparaison avec les donnes exprimentales de la
phase liquide de lactonitrile montre une autre incohrence. En effet, le rallongement de la chane alkyle
(passage de lactonitrile au propionitrile) devrait se traduire par une augmentation de la densit liquide
dans les mmes conditions de temprature et de pression alors que cest le contraire qui est constat : la
table III-14, donne une densit liquide exprimentale de lactonitrile gale 685 kg.m-3 pour les conditions
thermodynamiques de 410,14 K et 24,826 MPa alors que les simulations du propionitrile dans des
conditions presque similaires, le dernier point de la Table III-16, pour une temprature de 400 K et une
pression de 25 MPa donne une densit liquide infrieure gale 640,21 kg.m-3.
100000
300 K
320 K
400 K
1000
252,14 K
289,96 K
314,02 K
342,62 K
368,19 K
390,70 K
417,33 K
pression [kPa]
443,17 K
10000
443,17 K
350 K
290 K
simulation molculaire charges MEP
densit liquide bouillant (DIPPR)
isotherme exprimentale (Chakhmuradov, 1984)
100
400
450
500
550
600
650
700
750
800
850
-3
densit [kg.m ]
Figure III-8. Courbe de coexistence du propionitrile. Corrlation DIPPR, valeurs prdites avec les charges
ponctuelles MEP et valeurs exprimentales [Chakhmuradov et Guseinov, 1984].
III-27
v.
La Figure III-9 et la Figure III-10 comparent les valeurs des paramtres de Lennard-Jones du
groupe nitrile optimiss avec les charges ponctuelles MEP avec les valeurs de ces paramtres pour dautres
groupes chimiques du champ de force AUA4 issu des alcanes.
3,80
3,60
sigma ()
3,40
3,20
3,00
2,80
2,60
CH4
CH3
CH2
CH
C (C=C)
C (C N) N (C N) C (C N) N (C N)
MEP
MEP
Mlliken Mlliken
Figure III-9. Paramtres de Lennard-Jones sigma pour le carbone et lazote du groupe nitrile optimiss
avec les charges MEP ou Mulliken et pour diffrents groupes chimiques du champ de force AUA4.
180,00
epsilon/kB (K)
150,00
120,00
90,00
60,00
30,00
0,00
CH4
CH3
CH2
CH
C (C=C)
C (C N) N (C N) C (C N) N (C N)
MEP
MEP
Mlliken Mlliken
Figure III-10. Paramtres de Lennard-Jones epsilon pour le carbone et lazote du groupe nitrile
optimiss avec les charges MEP ou Mulliken et pour diffrents groupes chimiques du champ de force
AUA4.
Discussion
Le paramtre sigma est gnralement associ dans le potentiel de Lennard-Jones au rayon
dinteraction du centre de force.
III-28
Il est anormalement lev pour latome de carbone isol du groupe nitrile, se situant entre les
valeurs pour les groupes de trois atomes CH2 et de quatre atomes CH3. Mme si ces
paramtres sont rgresss, une logique physique voudrait que la valeur de C dans -CN soit
plutt voisine de petits groupes chimiques (C dans C=C ou CH). La valeur finale optimise
est dailleurs largement suprieure aux valeurs rapportes dans la littrature pour latome de
carbone (Table III-11 o les valeurs dinitialisation sont prises dans la littrature).
Il est comparable pour lazote terminal du groupe nitrile aux valeurs rapportes dans la
littrature pour latome dazote comme le montre la Table III-11. Nanmoins, les valeurs de la
littrature correspondent une autre rgression que celle entreprise ici.
La contribution nergtique de Van der Waals lui est directement proportionnelle. Par
consquent, il semblerait logique que, pour une molcule polaire, comme lactonitrile, la
valeur de ce paramtre soit plus importante que dans des molcules peu polaires comme les
alcanes. Ce nest pas le cas avec le jeu de paramtre bas sur des charges MEP. En effet, les
valeurs du groupe nitrile obtenues sont faibles et comparables aux valeurs de petits groupes
chimiques (C dans C=C ou CH) du champ de force AUA4 et celles rapportes dans la
littrature (Table III-11).
vi.
Les charges MEP permettent de reproduire avec une prcision comparable aux mesures
exprimentales les proprits de lactonitrile en phase liquide condens ou en quilibre
liquide vapeur.
Cependant, la transfrabilit des paramtres de Lennard-Jones trouvs est mise en dfaut,
dune part par limpossibilit de simuler lquilibre liquide vapeur du propionitrile,
dautre part par des valeurs des paramtres, certes rgresss, qui, compars aux valeurs
dautres groupes chimiques, ne prsentent pas la logique physique escompte.
Aprs constat de la faible transfrabilit des paramtres de Lennard-Jones obtenus avec les
charges ponctuelles MEP, de nombreuses tentatives ont t ralises pour amliorer ces
paramtres en conservant les charges ponctuelles MEP. Elles ont choues. Nous sommes
arrivs la conclusion que la solution ne rsidait pas dans lvaluation de linteraction de
Van der Waals au travers des paramtres de Lennard-Jones, mais plutt dans lvaluation
de linteraction lectrostatique :
Pour lactonitrile, les charges ponctuelles MEP conduisent une valuation de la
contribution lectrostatique de lnergie dinteraction intermolculaire trop importante (de
III-29
b.
La Table III-17 montre lvolution de la fonction objectif et des valeurs des paramtres de LennardJones du groupe nitrile au cours de leur optimisation avec les charges ponctuelles de Mulliken.
C (C N)
N (C N)
/kB (K)
()
/kB (K)
()
Initialisation
100,00
3,300
140,00
3,350
2,91
1re optimisation
98,28
3,383
152,72
3,542
2,04
2nde optimisation
95,52
3,218
162,41
3,564
1,63
Table III-17. Evolution des paramtres de Lennard-Jones du groupe nitrile. Charges ponctuelles Mulliken.
La Figure III-9 et la Figure III-10 comparent les valeurs des paramtres de Lennard-Jones du groupe
nitrile optimiss avec les charges ponctuelles Mulliken et ceux optimiss avec les charges MEP, avec les
valeurs de ces paramtres pour dautres groupes chimiques du champ de force AUA4 issu des alcanes.
Ayant ralis cette optimisation aprs avoir constat lchec de la transfrabilit des paramtres obtenus
avec les charges MEP, la transfrabilit des paramtres de la Table III-17 a t value la fin de chaque
III-30
cycle doptimisation en vrifiant la possibilit de dmarrer une simulation dquilibre liquide vapeur pour
le propionitrile et le n-butyronitrile.
Discussion
Le paramtre sigma :
-
Contrairement aux valeurs optimises avec les charges ponctuelles MEP, le paramtre sigma
du carbone C (-CN) est du mme ordre de grandeur que les valeurs pour de petits
groupes chimiques (C dans C=C ou CH).
le paramtre sigma de lazote N (-CN) est plus grand que le paramtre sigma du
carbone C (-CN). Cest un rsultat logique dans la mesure o le nuage lectronique de
lazote est plus gros que le nuage lectronique du carbone en terme de rayon de Van der
Waals.
Le paramtre epsilon :
-
Les deux valeurs du groupe nitrile sont cette fois-ci plus leves que les valeurs optimises
avec les charges ponctuelles MEP. En parallle avec des charges ponctuelles plus faibles pour
Mulliken que pour MEP, les valeurs du paramtre epsilon optimises avec les charges
ponctuelles Mulliken compensent probablement la rduction de linteraction lectrostatique.
i.
Les rsultats des simulations dquilibre liquide vapeur dans lensemble de Gibbs Monte Carlo
(433,15 K, 443,15 K et 453,15 K) sont prsentes dans la Table III-18 ainsi que les valeurs exprimentales.
Ils montrent que la valeur de la fonction objectif diminue au fur et mesure des cycles doptimisation.
Lerreur relative moyenne finale et calcule sur les grandeurs de la Table III-12 est infrieure 3%. Cest
un excellent rsultat mais pour les grandeurs prises individuellement, on constate que les rsultats sont
moins bons que ceux obtenus avec les charges MEP. De nouveau, les carts les plus importants sont
constats pour la pression de la phase vapeur.
III-31
T(K)
Type de
donnes
Conditions
433,15
443,15
453,15
GEMC
GEMC
Valeurs
exprimentales
Point initial
2nde
optimisation
1re
optimisation
613,37a
601,22 1,981
618,45 0,828
608,47 0,799
Hvap
Ln (Psat)
25120,8b
24634,9 1,934
27779,2 10,582
26371,1 4,977
13,5007a
13,9780 3,535
13,5269 0,194
13,7570 1,898
599,90a
592,11 1,299
607,24 1,224
600,50 0,100
Hvap
Ln (Psat)
24319,7b
24205,9 0,468
27226,2 11,951
25816,7 6,156
13,6963a
14,1955 3,645
13,7745 0,571
13,9803 2,074
583,197a
579,09 0,704
598,82 2,679
588,19 0,856
Hvap
23465,86b
23436,2 0,126
26663,8 13,628
25015,5 6,604
Ln (Psat)
13,8934a
14,3659 3,401
13,9976 0,750
14,1280 1,689
GEMC
1,63
Table III-18. Evolution des principales grandeurs physico-chimiques en fonction des cycles
doptimisation des paramtres de Lennard-Jones du groupe nitrile. Actonitrile avec charges
ponctuelles Mulliken.
ii.
Le cycle doptimisation achev, le champ de force ainsi complt est utilis pour prdire dautres
proprits de lactonitrile. La Table III-19 et la Table III-20 prsentent des simulations complmentaires
de la phase liquide condense dans lensemble de Monte Carlo NPT et dquilibre liquide vapeur dans
lensemble de Gibbs Monte Carlo.
T (K)
Vap (kg.m-3)
exp.a
exp b
Liq (kg.m-3)
Simule
exp.a
exp b
Pression(kPa)
Simule
Pliq
Psat exp.a
Psat
413,15
7,88
---
7,945 288
641,49
---
629,36 1,56
423,15
9,78
9,88
9,423 357
627,98
635,0
617,91 1,73
433,15
11,92
613,37
---
608,47 2,03
443,15
14,45
599,90
---
600,50 1,75
448,15
---
---
598,0
453,15
17,59
583,19
---
588,19 1,15
463,15
20,31
567,15
---
580,52 1,04
473,15
---
---
---
556,2
500,00
---
35,132 900
---
520,00
---
49,700 1,301
523,15
---
---
15,22
22,5
68,0
---
---
---
---
---
---
---
---
---
526,09 1,94
2736,8 493,1
---
---
497,86 3,35
3572,3 403,7
---
436,0
---
---
---
---
---
Table III-19. Donnes dquilibre liquide vapeur exprimentales et simules dans lensemble de Gibbs
Monte Carlo. Actonitrile avec charges ponctuelles Mulliken.
III-32
T (K)
P (MPa)
Psim (MPa)
liqsim (kg.m-3)
liqexp (kg.m-3)
Erreur (%)
298,405
24,714
24,938 495
741,47 1,48
795,97
6,9
325,809
4,1317
3,861 780
709,10 93
750,41
5,5
435,196
7,5078
7,448 688
614,64 1,99
625,75
1,8
Table III-20. Densits du liquide simule et exprimentales [ Kratzke , 1985 ]. Actonitrile avec les charges
ponctuelles Mulliken
La Figure III-4 dcrivant la courbe de coexistence liquide vapeur montre que globalement, les
simulations avec les charges MEP sont en meilleur accord avec les densits liquides que les simulations
avec les charges Mulliken. Linverse est constat pour les densits de la phase vapeur. En particulier, les
densits liquides calcules avec les charges Mulliken sont trop leves haute temprature et trop faibles
basse temprature. Ces constatations sont encore plus videntes sur la Figure III-11 qui prsente les
rsultats des simulations dquilibre liquide vapeur sur le diagramme pression densit. Suivant les
conventions usuelles, les pressions de la phase liquide sont prises gales celle de la phase vapeur.
100000
pression [kPa]
10000
1000
10
0
100
200
300
400
500
600
700
800
-3
densit [kg.m ]
Figure III-11. Diagramme de Clapeyron de lactonitrile. Valeurs calcules avec charges ponctuelles
Mulliken, valeurs exprimentales [ Francesconi, 1975 ; Kratzke et Mller, 1985 ; Warovny, 1994 ].
III-33
Les enthalpies de vaporisation calcules avec les charges ponctuelles Mulliken sont reportes sur la Figure
III-7. Elles scartent significativement de la corrlation, confirmant la tendance observe sur la courbe de
coexistence liquide vapeur de la Figure III-4 et confirmant le fait que le champ de force avec les charges
de Mulliken dcrit trop approximativement le comportement thermodynamique de lactonitrile.
iii.
Lestimation du point critique selon la mthode prsente dans le chapitre II est reporte sur la
Figure III-4 et dans la Table III-21. Lestimation de la temprature critique est mdiocre (survalue de
33K) compare au point critique exprimental dtermin par Francesconi et al. [ 1975 ] mais aussi
lestimation ralise avec les charges MEP . On pouvait sy attendre cause des densits liquides simules.
Lcart sur la densit critique est faible au regard des remarques faites sur la mthode dIsing qui semble
surestimer cette grandeur.
exp.a
Estimation calculs
avec Mulliken
TC (K)
547,85
580,06
C (kg.m-3)
237,1
248,40
Francesconi et al [ 1975 ]
Table III-21. Estimation des coordonnes critiques de lactonitrile avec charges Mulliken
iv.
Discussion
On constate que le modle avec les charges ponctuelles Mulliken reprsente moins bien les
proprits de lactonitrile, notamment en phase liquide (estimation des coordonnes critiques
mdiocre, cart sur les densits liquides, enthalpies de vaporisation dcales) que le modle avec les
charges ponctuelles MEP. Or, la faiblesse relative de ses paramtres de Lennard-Jones indiquait que
le modle avec les charges MEP semblait exagrer linteraction lectrostatique. Nous pouvons
proposer comme critique et comme hypothse que le modle avec les charges ponctuelles Mulliken
ne prend pas assez en compte les interactions lectrostatiques en phase liquide. Six corrections
peuvent tre imagines :
1. valuer les interactions lectrostatiques en nimposant pas un seuil de coupure mais en
employant des mthodes rputes moins approximatives telle que la sommation Ewald. La
nouvelle version logicielle du code GIBBS propose cette option. La mise en uvre de la
sommation Ewald rallonge considrablement le temps de calcul comme nous lavons
indiqu auparavant.
III-34
De nouveau, la transfrabilit des paramtres de Lennard-Jones du groupe nitrile bass sur les
charges ponctuelles de Mulliken est value en simulant dans lensemble de Gibbs Monte Carlo NVT les
proprits dquilibre liquide vapeur du propionitrile et du n-butyronitrile. Ces paramtres sont pris
gaux la valeur trouve pour les simulations de lactonitrile. Les autres paramtres sont dtermins
comme auparavant (distribution de charges ponctuelles issue dune analyse de population de Mulliken
partir de surface de potentiel lectrostatique calcul en DFT, constantes harmoniques et conformation
dquilibre calcule par DFT, autres paramtres de Lennard-Jones pris dans le champ de force AUA4).
La Table III-22 et la Table III-24 prsentent les rsultats des prdictions des quilibres liquide vapeur
pour le propionitrile et le n-butyronitrile respectivement. Les valeurs sont compares aux donnes
exprimentales dquilibre liquide - vapeur sur les Figure III-12 et Figure III-15, de pression de vapeur
III-35
saturante sur les Figure III-13 et Figure III-16. La Table III-23 prsente les rsultats des prdictions de la
densit du propionitrile liquide condens et les donnes sont reproduites sur la Figure III-14. Les
enthalpies de vaporisation sont reportes sur la Figure III-7.
T (K)
400
425
450
460
480
500
510
520
Vap (kg.m-3)
Liq (kg.m-3)
Psat (kPa)
Pliq (kPa)
5,729 191
9,347 321
14,552 273
18,636 488
24,786 559
38,294 1,372
44,273 1,531
50,539 1,536
650,64 81
623,15 1,42
596,12 1,52
588,09 1,33
559,76 1,74
541,75 2,64
526,77 2,55
509,30 2,44
315,67 9,21
545,25 22,11
871,35 16,55
1110,30 23,29
1483,90 28,03
2110,30 75,29
2423,20 63,20
2778,60 91,28
-90,7 294,2
481,0 483,0
143,5 425,0
1613,9 328,8
1699,2 302,9
2968,3 454,2
3246,2 372,6
3165,6 305,6
Hvap (J,mol-1)
29 913,0
28 360,0
26 634,3
25 840,8
24 055,4
21 983,8
20 814,3
19 802,1
Table III-22. Prdictions dquilibre liquide vapeur du propionitrile avec les paramtres de LennardJones optimiss pour lactonitrile avec les charges Mulliken. Moyenne prise entre 2 et 7 millions de
configurations.
10000
pression [kPa]
1000
100
0
100
200
300
400
500
600
700
800
Figure III-12. Courbe de coexistence du propionitrile. Corrlation DIPPR et valeurs prdites avec les
charges ponctuelles de Mulliken.
III-36
10000
Chakhmuradov 84
1000
Pression (kPa)
100
0,0027
0,0026
0,0025
1/400
0,0024
0,0023
0,0022
0,0021
0,002
1/500
0,0019
0,0018
0,0017
1/Temprature [K]-1
Figure III-13. Pression de vapeur saturante du propionitrile exprimentale et prdite laide des
paramtres de Lennard-Jones optimiss pour lactonitrile avec les charges Mulliken.
T (K)
289,96
290
289,96
314,02
314,02
320
342,62
342,62
342,62
342,62
350
368,19
368,19
390,7
390,7
400
417,33
417,33
a Chakhmuradov
P (MPa)
0,1
1,0
5,0
30,0
40,0
35,0
30,0
40,0
10,0
20,0
15,0
10,0
20,0
20,0
30,0
25,0
20,0
30,0
liqsim (kg.m-3)
liqexp(a) (kg.m-3)
784,0
744,07 99
786,9
783,0
789,6
741,49 84
758,4
766,0
740,9
750,2
709,34 1,13
715,2
725,9
704,0
715,1
671,94 1,28
676,3
689,9
et Guseinov [ 1984 ]
Table III-23. Comparaison des densits liquide obtenues par simulation pour le propionitrile et corrles.
Paramtres de Lennard-Jones optimiss pour lactonitrile avec les charges Mulliken.
III-37
100000
320 K
400 K
1000
252,14 K
289,96 K
314,02 K
342,62 K
390,70 K
417,33 K
443,17 K
443,17 K
pression [kPa]
368,19 K
350 K
10000
290 K
simulation molculaire charges Mulliken
densit liquide bouillant (DIPPR)
isotherme exprimentale (Chakhmouradov, 1984)
100
400
450
500
550
600
650
700
750
800
850
-3
densit [kg.m ]
Figure III-14. Densit du liquide condens du propionitrile avec des paramtres de Lennard-Jones
optimiss pour lactonitrile avec les charges Mulliken.
T (K)
400
425
450
465
482
500
520
540
Vap (kg,m-3)
3,361 151
5,743 130
9,885 221
12,977 397
18,101 344
24,082 471
33,149 536
48,340 1,408
Liq (kg,m-3)
671,21 1,70
649,32 1,13
624,28 1,00
604,09 1,86
587,62 1,26
568,05 1,40
541,52 1,32
512,46 2,83
Psat (kPa)
155,73 7,16
276,76 5,74
488,97 8,29
642,51 15,84
905,29 14,59
1193,70 21,76
1603,00 25,97
2222,10 42,37
Pliq (kPa)
Hvap (J.mol-1)
-86,6 277,2
451,0 411,9
972,5 326,3
116,1 357,3
645,2 259,7
1270,9 247,3
1438,7 304,1
2568,9 278,2
34 063
32 544
30 669
29 394
28 118
26 528
24 427
21 661
Table III-24. Prdictions dquilibre liquide vapeur du n-butyronitrile laide des paramtres de
Lennard-Jones optimiss pour lactonitrile avec les charges Mulliken. Moyenne prise entre 2 et 9 millions
de configurations.
Une remarque concerne les fluctuations de la pression en phase liquide. Nous avons constat que ces
fluctuations sont toujours trs importantes compares ceux de la phase vapeur. Cela est d la proximit
des molcules entres elles qui est responsable de fortes contraintes : dans certains cas, ceci sest traduit par
des valeurs moyennes ngatives (la Table III-24 pour T = 400 K) ; dans dautres cas, ltat dquilibre nest
pas atteint puisque les deux pressions liquide et vapeur sont diffrentes mme en tenant compte des
incertitudes, cest le cas pour les trois molcules mme en augmentant les longueur des simulations.
Un remde possible serait de jouer sur le taux de transfert, ici pris gal 69,5 % pour les trois molcules.
III-38
10000
pression [kPa]
1000
100
0
100
200
300
400
500
600
700
800
-3
densit [kg.m ]
Figure III-15. Courbe de coexistence du n-butyronitrile. Corrlation DIPPR et valeurs prdites laide des
paramtres de Lennard-Jones optimiss pour lactonitrile avec les charges Mulliken.
10000
Chakhmuradov 84
1000
Pression (kPa)
100
0,0026
0,0025
1/400
0,0024
0,0023
0,0022
0,0021
0,002
1/500
0,0019
0,0018
0,0017
1/Temprature [K]-1
Figure III-16. Pression de vapeur saturante du n-butyronitrile exprimentale et prdite laide des
paramtres de Lennard-Jones optimiss pour lactonitrile avec les charges Mulliken.
III-39
Les valeurs sont de vraies prdictions. Lexamen des tables et des graphes montre que laccord est
remarquable en ce qui concerne la pression de vapeur saturante et moins bon en ce qui concerne la densit
du liquide bouillant. Pour le propionitrile liquide condens 290 K, la valeur de la densit est sousestime. Le mme phnomne, pour la mme temprature avait t constat avec les charges MEP lors de
lvaluation de la transfrabilit des paramtres MEP. En revanche, les autres carts avec les isothermes
exprimentales sont plus faibles quavec les charges MEP et on remarque que la densit liquide du
propionitrile pour le point 400 K et 25 MPa est plus leve que celle obtenue avec les charges MEP,
mais reste toujours plus faible que celle de lactonitrile dans des conditions presque semblable (cf.
commentaire page III-27 pour les charges MEP).
Cependant, les enthalpies de vaporisation sont en bon accord avec les donnes de la littrature, un accord
meilleur que pour lactonitrile sans que nous lexpliquions.
La mauvaise prdiction des densits du liquide bouillant se traduit comme pour lactonitrile par une
mauvaise estimation du point critique reporte dans la Table III-25 pour le propionitrile et le nbutyronitrile, selon la mme mthode (chapitre II). Comme pour lactonitrile, lestimation est mdiocre
en ce qui concerne la temprature critique, la surestimant denviron 20K. Lcart sur la densit critique est
faible au regard des remarques faites sur la mthode dIsing qui semble surestimer cette grandeur.
propionitrile
n-butyronitrile
exp.a
exp.b
Estimation
calculs avec
Mulliken
exp.a
exp.b
Estimation
calculs avec
Mulliken
TC (K)
564,4
561,20,2
584,49
582,85
585,40,2
606,50
C (kg.m-3)
240,52
255,05
248,80
a
b
252,05
Les charges Mulliken reproduisent avec satisfaction les proprits de lactonitrile liquide
condens ou en quilibre liquide vapeur mais avec une prcision significativement
moindre de celle obtenue avec les charges MEP comme le montrent principalement les
valeurs de la densit du liquide bouillant et de lenthalpie de vaporisation.
Cependant, la transfrabilit des paramtres de Lennard-Jones avec les charges Mulliken
est avre : les prdictions vraies (sans nouvel ajustement de paramtres) des proprits
III-40
E.
Excellent rsultat pour lactonitrile mais incohrence pour les autres nitriles . Telle pourrait tre
lapprciation des rsultats obtenus en optimisant des paramtres de Lennard-Jones avec les charges
ponctuelles MEP. La procdure doptimisation des paramtres donne effectivement des valeurs de
paramtres epsilon et sigma de Lennard-Jones pour le carbone et lazote du groupe nitrile qui
permettent de reproduire avec prcision les donnes exprimentales disponibles pour lactonitrile
lorsquon les combine avec des paramtres de Lennard-Jones gnriques du champ de force AUA4 et des
constantes intramolculaires issues de calculs quantiques.
Cependant, la transfrabilit des paramtres associs aux charges MEP pour prdire les proprits du
propionitrile est mise en dfaut. Cela nous a conduit raliser une autre optimisation avec des charges
ponctuelles de Mulliken issues dune analyse de population trs simple.
Bien mais peut mieux faire . Telle pourrait tre lapprciation des rsultats obtenus en optimisant des
paramtres de Lennard-Jones avec les charges ponctuelles Mulliken. Certes, les paramtres sont
transfrables et donnent des prdictions remarquablement en accord avec les rares donnes
exprimentales disponibles pour le propionitrile et le n-butyronitrile. Mais, les prdictions des proprits
de lactonitrile avec les charges Mulliken sont globalement moins bonnes quavec les charges MEP,
notamment en ce qui concerne la phase liquide et cela est flagrant puisque la temprature critique est
systmatiquement surestime avec les charges Mulliken denviron 20K pour les trois nitriles tudis.
En termes de valeurs numriques des paramtres epsilon et sigma de Lennard-Jones obtenus pour le
carbone et lazote du groupe nitrile, les valeurs avec les charges de Mulliken sont plus cohrentes que
celles obtenues avec les charges MEP lorsquon les compare aux autres valeurs obtenues pour des groupes
du champ de force AUA4.
III-41
A posteriori, nous avons constat que la distribution de charges ponctuelles MEP pour le propionitrile nest
pas cohrente avec celle de lactonitrile, ce qui pourrait expliquer la non transfrabilit du jeu de
paramtres MEP malgr son excellente reprsentation des proprits de lactonitrile.
Par consquent, nous pouvons proposer deux critres pour lobtention de paramtres de Lennard-Jones
gnriques :
1. homognit des distributions de charges sur les sries chimiques homologues lors de leur
obtention par analyse de population des surfaces de potentiel lectrostatique. Lanalyse de
Mulliken semble tre trop simpliste et se traduit par des rsultats bons mais pas excellent. Les
distributions de charges NBO mriteraient dtre testes.
2. vrification de la pertinence physique des paramtres optimiss par comparaison dautres
valeurs obtenues pour le mme champ de force.
En perspective, il semble opportun de raliser une nouvelle optimisation de paramtres avec les charges
NBO prsentes en annexe. Les charges NBO sont plus fortes que les charges MEP sur les groupes
mthyle (latome C1 et les hydrognes) mais infrieure pour le groupe nitrile, tout en tant suprieure aux
charges de Mulliken (voir Table III-1 et Table 1 en annexe III). Mettre en uvre dautres corrections
voques dans ce chapitre comme lutilisation de la sommation Ewald ne semble pas une priorit court
terme.
III-42
A.
Introduction
L'importance des quilibres entre phases pour valuer la faisabilit et pour concevoir des procds
de distillation est bien connue. En particulier, essentielle est la capacit de prvoir si un mlange donn
forme un ou plusieurs azotropes et de calculer les conditions appropries de composition, temprature et
pression qui y sont associes [ Widagdo et Seider, 1996, Kiva et al., 2003 ]. La Figure IV-1 rappelle les
diffrents types dazotropes binaires :
-
Mlange idal
n-butane n-heptane
T
Vapeur
P = cst.
Vapeur
P = cst.
Courbe de rose
L2+V
L1+V
Tazeo
Courbe de bulle
a
L1
b
Liquide
xazeo
xA, yA
eau thanol
T
Vapeur
chloroforme - actone
T
P = cst.
Vapeur
P = cst.
xA, yA
Min T homoazotrope
avec rgion LLV
dviation/Raoult [+]
eau 2-butanone
Max T homoazotrope
dviation/Raoult [-]
Min T homoazotrope
dviation/Raoult [+]
L2
L1 + L 2
Vapeur
P = cst.
Tazeo
L1
L2+V
L1+V
Tazeo
Tazeo
Liquide
c
xazeo xA, yA
Liquide
xazeo
xA, yA
L1 + L 2
xazeo
L2
xA, yA
Figure IV-1. Equilibres liquide vapeur pour diffrents systmes binaires azotropiques
IV-1
B.
Dveloppement de la mthodologie.
1.
Travaux antrieurs
La seule utilisation dune technique dchantillonnage base sur des simulations dans lensemble de
Gibbs Monte Carlo nest pas adquate pour dterminer un azotrope. En premier lieu, dans lhypothse
IV-2
o il existe un azotrope (ce qui nest pas connu a priori), il y a une probabilit infime que la simulation
dans lensemble de Gibbs Monte Carlo converge exactement au point azotropique. En deuxime lieu, les
simulations dans lensemble de Gibbs Monte Carlo se font dans les conditions NVT pour les corps purs
et NVT ou NPT pour les mlanges de faon satisfaire la rgle de phase (chapitre II) indiquant que pour
un systme diphasique le nombre de proprits intensives que lon peut fixer est gal 1 pour un systme
monoconstituant et 2 pour un systme binaire.
Or, nous pouvons dmontrer que raliser directement la simulation du point azotropique ncessite de se
placer dans les conditions PT, ce qui viole la rgle de phase puisqualors quatre proprits intensives sont
fixes (1, 2, P et T) : en effet, une simulation dans lensemble de Gibbs Monte Carlo NPT permet
dobtenir a posteriori (1, 2) mais sans tre forcment au point azotropique. Pour cela, il faudrait ajouter
aux traditionnels mouvements dans lensemble de Gibbs des mouvements de cration/destruction de
particules de faon changer la composition totale N=NAI+NAII+NBI+NBII pour quau final elle satisfasse
les conditions azotropiques (XAi=XBi) avec Xji=Nji/Ni. Ne pas fixer Ntotal a priori revient imposer a priori
(1, 2) et donc faire des simulations diphasiques binaires dans un ensemble de Gibbs PT, en
contradiction avec la rgle de phase.
Une dernire possibilit de simulation directe dans lensemble de Gibbs serait de fixer Ntotal a priori et de ne
modifier que lidentit des molcules au sein des boites pour modifier la composition. Cependant, cela
ncessite de partir dun nombre de molcules de chaque constituant identique dans chaque boite et qui
satisfassent la composition de lazotrope, point que nous ne connaissons pas et cherchons prcisment
dterminer.
Le calcul direct des azotropes par la simulation molculaire a pourtant t ralis avec succs en utilisant
l'intgration de Gibbs Duhem qui fournit des informations de coexistence semblables l'ensemble Monte
Carlo de Gibbs (voir chapitre II) mais ne rend pas ncessaire l'change de particules pour atteindre
l'quilibre potentiel chimique puisque celui-ci est assur par l'intermdiaire de lintgration de l'quation de
Clapeyron [ Pandit et Kofke, 1999 ] suivant la ligne de coexistence vers l'azotrope.
Ces deux constatations nous ont conduit imaginer une approche alternative base sur une adaptation de
la mthode de simulation dans lensemble de Gibbs Monte Carlo volume impos couple avec des
mouvements de changements didentit appliqus chacune des deux phases sparment. Cette ide est
galement analogue celle imagine par Teja et Rowlinson [ 1973 ] qui ont calcul les azotropes
homognes de mlanges binaires en utilisant une quation d'tat en tant que modle thermodynamique :
ils fixaient la temprature et changeaient la composition et le volume jusqu' ce que les conditions
dazotropie soient satisfaites.
La mthodologie peut tre dcrite succinctement par :
1. valuation de lexistence dun azotrope par calcul de la valeur des constantes dquilibre aux
bornes du diagramme binaire par simulation dans lensemble de Gibbs Monte Carlo NVT,
IV-3
2. calcul prcis des caractristiques de lazotrope par une succession de simulation dans
lensemble de Gibbs Monte Carlo NVT sur deux boites et des mouvements de changement
didentit raliss par boite indpendamment afin de progresser le long de la courbe de
coexistence liquide vapeur vers le point azotropique.
2.
Nouvelle
particule
i VT = Q N i VT exp i N i
Ni
i
O Q N i VT =
(IV-1)
Ej
S = kB ln i VT + E i N i
i
(IV-2)
IV-4
kBT .
qi N i
i N ! exp U ext + i i N i
i
O qi = V
(IV-3)
(2mi kT )12
(n )
p (o n ) = min 1,
(o )
(IV-4)
(o )
q1 (o )N 1 (o ) q 2 (o )N 2 (o )
exp( (U ext (o ) (1 N 1 (o ) + 2 N 2 (o ))))
(
)
(
)
N
o
!
N
o
!
1
2
q (n )N 1 (n ) q (n )N 2 (n )
(IV-5)
(IV-6)
De plus,
(i)
q 2 (o ) = q 2 (n ) = q 2 .
(ii)
N 2 (o ) = N 2 (n ) + 1 , donc N 1 (o ) + N 2 (o ) = N 1 (n ) + N 2 (n ) .
En dfinitive, le mouvement de changement didentit doit tre accept avec une probabilit :
IV-5
N (o) q1
N1 (o ) + 1 q2
(IV-7)
Auparavant, Kofke et Glandt [ 1988 ] ont propos la mme ide dans un ensemble semi-grand canonique,
un ensemble qui incorpore des caractristiques des ensembles Canonique et Grand Canonique.
Remarques :
1. Un mouvement dinsertion dune particule de type 1 dans lensemble Grand Canonique a une
probabilit dacceptation :
(IV-8)
tandis que la destruction dune particule de type 2 est accepte avec une probabilit :
(IV-9)
Il est vident de ces deux dernires quations que la probabilit du mouvement prsent est le
produit des probabilits classiques de l'insertion et de la destruction :
(IV-10)
2. Comme nous venons de le dmontrer, le mouvement de changement didentit ralis est bas sur
la fonction de partition de lensemble grand canonique et sa probabilit dacceptation est gale au
produit des deux probabilits dacceptation des deux mouvements dinsertion et de destruction
spcifiques cet ensemble statistique. Nanmoins, ce mouvement ne suffit pas simuler lensemble
grand canonique VT car, comme le montre lquation (IV-7), il nest pas ncessaire de fournir les
deux potentiels 1 et 2 pour appliquer ce mouvement mais seulement leur diffrence (2-1). Dans
ce travail, ces deux potentiels chimiques sont calculs partir des simulations effectues dans
lensemble de Gibbs Monte Carlo NVT selon la mme procdure que celle propose par
Panagiotopoulos et al. [ 1988 ]. Par consquence, ceci revient effectuer ce mouvement dans un
pseudo-ensemble N(2-1)VT sous condition que le principe de micro rversibilit soit respect,
c'est--dire raliser autant de tentative de changement didentit dans un sens que dans lautre. Ce
nest pas le cas dans ce travail et cest pour cette raison quon parlera uniquement de mouvements
de changement didentit permettant une rinitialisation des deux phases sans prciser lensemble
statistique dans lequel ces mouvements sont oprs.
3. Nous pouvons faire lhypothse dune sparation complte entre les contributions cintiques et
m
exp( U int )dr j et se
3i r mj
IV-6
3i
Equilibre chimique entre phases (nergie minimale, galit des tempratures, pressions et
potentiels chimiques de chaque constituant dans chaque phase),
De lgalit des potentiels chimique dcoule lexpression macroscopique de lquilibre entre deux phases
avec lintroduction dune constante dquilibre Ki :
yi = K i . x i
(IV-11)
En supposant que le constituant A est plus volatil que le constituant B, la Figure IV-1 montre que :
-
Au point azotropique KA = 1.
Ainsi, nous pouvons dfinir une procdure dvaluation rapide de lexistence dun azotrope dans un
mlange binaire :
-
C.
Description de lalgorithme
IV-7
au voisinage des composants purs, par exemple ZH = 0,975 et ZL = 0,025, les indices H et L signifiant
respectivement haut et bas.
ZL = 0.025
ZH = 0.975
Gibbs NVT
Gibbs NVT
KL = (x/y)L
stocke la configuration
finale (XYZ)L
KH = (x/y)H
stocke la configuration
finale (XYZ)H
OUI
KL = 1
Coordonnes azotropiques
OUI
KH = 1
NON
pas dazotrope
STOP
NON
KL 1
<0
KH 1
1RE TAPE
I=I+1
Changement
didentit
sur (XYZ)L
OUI
ABS(ZI-ZL) <
ABS(ZI-ZH)
Changement
didentit
sur (XYZ)H
NON
Gibbs NVT
KH = KI
KL = KI
KI = 1
2nde TAPE
ZL = ZI
NON
Coordonnes azotropiques
NON
KI 1
<0
KH 1
OUI
ZH = ZI
Figure IV-3. Algorithme pour trouver lventuel azotrope dun mlange binaire A-B
Naturellement, la densit initiale des botes doit tre choisie pour provoquer un quilibre liquide vapeur.
Aprs avoir fait la moyenne des compositions rsultantes x et y en quilibre, les constantes d'quilibre KH
et KL sont calcules :
-
IV-8
Sinon, on compare les deux constantes dquilibre KH et KL selon le critre dcrit auparavant
pour dterminer lexistence dun azotrope et son type.
Si lazotrope existe, la deuxime tape permet de calculer ses caractristiques composition, pression et
temprature (fixe a priori) en progressant le long de la courbe de coexistence liquide vapeur :
-
Puis, la dernire configuration obtenue par les mouvements de changement didentit est
utilise comme point de dpart dune simulation dans lensemble de Gibbs Monte Carlo NVT
qui a pour but de satisfaire la condition dquilibre liquide vapeur et de rvaluer la constante
dquilibre K :
o
Sinon, une nouvelle tape 2 est ritre sur une gamme de composition plus troite parce
qualors la configuration prcdente ayant une valeur de constante dquilibre place
comme KH ou KL par rapport lunit est remplace par la configuration finale obtenue de
la simulation dans lensemble Gibbs.
La progression le long de la courbe de coexistence liquide vapeur procde ainsi comme une mthode par
dichotomie sauf que le nouveau point nest pas situ systmatiquement au milieu des deux points
prcdents mais est choisi alatoirement.
D.
Rsultats
1.
Mlange 11
12
22
11
12
22
1,000
0,750
1,000
1,000
1,000
1,000
Oui
II
1,000
1,000
1,000
1,000
0,885
0,769
Non
III
1,000
0,773
0,597
1,000
0,884
0,768
Non
IV
1,000
0,900
1,000
1,000
1,050
1,000
Oui
1,000
0,900
1,000
1,000
1,050
1,150
Oui
Azotrope
IV-9
Les rsultats de simulations sur les fluides de Lennard-Jones sont toujours exprims en coordonnes
rduites par rapport aux paramtres et [ Frenkel et Smit, 1996 ].
Une simulation reprsentative concerne N = 800 particules. Pour les deux premiers points initiaux de la
premire tape (ZH = 0,975 et ZL = 0,025), nous dmarrons la simulation dans lensemble de Gibbs
Monte Carlo volume impos avec une configuration identique dans chacune des deux boites. Ainsi, avec
ZH = 0,975, chaque boite contient 390 particules de type A et 10 particules de type B places selon un
rseau cubique face centr respectant la densit initiale choisie. Le choix de cette densit initiale dans les
simulations de lensemble de Gibbs NVT est important pour assurer lchantillonnage correct dun
systme de deux phases en coexistence. Lorsque le diagramme de phase du fluide nest pas connu, cela
peut tre ralis par plusieurs simulations exploratoires de la coexistence des phases. Les diagrammes de
phase des cinq mlanges de la Table IV-1 ont t publis dans la littrature [ Panagiotopoulos et al., 1988 ;
Pandit et Kofke, 1999 ].
Le second paramtre dune simulation dans lensemble de Gibbs dlicat initialiser est le nombre de
tentatives de transfert de particule entre les boites. En effet, Panagiotopoulos [ 1987 ] a observ qu'
tempratures rduites leves, o la probabilit d'accepter un change est relativement grande, des
rsultats incorrects sont trouvs pour la phase liquide avec des pressions liquides moyennes parfois
ngatives et une trs grande incertitude. Ce problme est attnu en ajustant le nombre d'changes essays
(voir une illustration plus dtaille dans le chapitre III).
Une volution typique des densits des deux rgions en fonction du nombre de configurations produites
est prsente sur la Figure IV-4 pour le mlange I pendant une itration de ltape 2. Aprs une priode
d'quilibrage denviron 106 configurations, les proprits des deux rgions convergent vers des valeurs
stables tout fait diffrentes des conditions initiales.
0,8
reduced density
0,7
0,6
0,5
0,4
0,3
0,2
0
Figure IV-4. Densit rduite pour le mlange I. Simulation dans lensemble de Gibbs Monte Carlo NVT.
Plus de 3.106 configurations ont t produites pour chaque simulation, et des graphes tels que celui
reprsent sur la Figure IV-4 sont tracs pour vrifier que la convergence est satisfaisante.
IV-10
2.
La Table IV-2 prsente les rsultats de la premire tape de la mthodologie qui value l'existence
d'azotrope. Les nombres suivent la notation pour des valeurs d'incertitude utilises dans la simulation
molculaire (0,607823 = 0,60780,0023). La notation * indique des units rduites [ Frenkel et Smit, 1996 ].
1st Step
liq
vap
P*liq
P*vap
xA
yA
KA
Mixture I
ZL
0.607823
0.07678
0.056743
0.06195
0.00993
0.04405
>1
T* = 1.15
ZH
0.605636
0.310929
0.064754
0.059026
0.98375
0.96488
<1
Mixture II
ZL
1.302047
0.150371
0.128577
0.123035
0.02914
0.01418
<1
T* = 1.15
ZH
0.614828
0.07537
0.057635
0.06075
0.97832
0.97073
<1
Mixture III
ZL
1.146674
0.202549
0.119841
0.103916
0.05088
0.00401
<1
T* = 0.75
ZH
0.603833
0.065214
0.058647
0.05598
0.98481
0.935111
<1
Mixture IV
ZL
0.635625
0.048510
0.041430
0.04157
0.98004
0.961111
<1
T* = 1.10
ZH
0.629821
0.32387
0.020643
0.013125
0.02056
0.02895
>1
Mixture V
ZL
0.428021
0.03595
0.038224
0.03142
0.02032
0.04738
>1
T* = 1.10
ZH
0.641119
0.058316
0.074448
0.06782
0.97884
0.963612
<1
Table IV-2. Rsultats de la premire tape pour examiner l'existence des azotropes homognes.
En accord avec la Table IV-1, trois parmi les cinq mlanges prsentent un azotrope homogne (mlanges
I, IV et V). Comme indiqu par leurs paramtres de Lennard-Jones, les mlanges I et IV sont symtriques
et leurs coordonnes azotropiques devraient tre trouves pendant les itrations de la deuxime tape
x = y = 0,50. C'est pourquoi nous choisissons seulement les mlanges I et V pour illustrer la deuxime
tape de la mthodologie. En outre, les tudes prcdentes ont montr le point d'azotrope du mlange V
se situe autour de xA = yA = 0,70 [ Pandit et Kofke, 1999 ].
Comme dj observ dans la littrature, la pression de la phase vapeur fluctue moins que celle de la bote
liquide o les densits leves soumettent des particules de plus grandes contraintes lors du transfert
dune particule partir de la bote de faible densit vers la bote forte densit. Les compositions xA et yA
en quilibre encadrent la composition initiale (ZH ou ZL) comme dans n'importe quel calcul de flash
macroscopique qui conduit un quilibre entre phases.
3.
Les Table IV-3 et Table IV-4 prsentent les rsultats des simulations pour toutes les itrations
ralises au cours de la recherche des azotropes selon la mthodologie de ltape 2.
IV-11
Mlange I
*liq
*vap
P*liq
P*vap
xA
yA
KA
ZL = 0,025
0,607823
0,07678
0,056743
0,06195
0,00993
0,04405
>1
ZH = 0,975
0,605636
0,310929
0,064754
0,059026
0,98375
0,96488
<1
0,5124107
0,365161
0,107444
0,102337
0,705583
0,631053
<1
0,602331
0,07634
0,069254
0,06332
0,02337
0,095710
>1
0,4512176
0,406885
0,116675
0,115060
0,4963175
0,497195
0,468564
0,395631
0,115031
0,114040
0,620184
0,583054
<1
0,4370201
0,418251
0,119459
0,115354
0,547965
0,547161
0,428357
0,4085180
0,112536
0,114740
0,528536
0,5214116
Z1=0,666
Z2=0,053
Z3=0,497
Z4=0,599
Z5=0,530
Z6=0,505
ZL = 0,025
ZH = 0,975
ZL = 0,025
ZH = 0,666
ZL = 0,053
ZH = 0,666
ZL = 0,497
ZH = 0,666
ZL = 0,497
ZH = 0,599
ZL = 0,497
ZH = 0,530
Mlange V
*liq
*vap
P*liq
P*vap
xA
yA
KA
ZL = 0,025
0,428021
0,03595
0,038224
0,03142
0,02032
0,04738
>1
ZH = 0,975
0,641119
0,058316
0,074448
0,06782
0,97884
0,963612
<1
0,547122
0,054411
0,045527
0,04566
0,706921
0,702867
0,489441
0,056523
0,043333
0,046010
0,471310
0,550643
>1
0,504913
0,06269
0,051933
0,04895
0,519319
0,595457
>1
0,540117
0,057415
0,046625
0,04718
0,685317
0,700660
>1
0,551330
0,068214
0,054927
0,05226
0,695818
0,698249
Z1=0,7068
Z2=0,4929
Z3=0,5343
Z4=0,6898
Z5=0,6991
ZL = 0,025
ZH = 0,975
ZL = 0,025
ZH = 0,707
ZL = 0,493
ZH = 0,707
ZL = 0,534
ZH = 0,707
ZL = 0,689
ZH = 0,707
IV-12
0,12
Pressure
0,1
0,08
0,06
ZL
Z2
Z3 Z5
Z4
Z1
ZH
0,04
0
0,1
0,2
0,3
0,4
0,5
0,6
0,7
0,8
0,9
Figure IV-5. Courbe de coexistence liquide vapeur pour le mlange I T*=1.,15. Nos simulations
(symboles en croix et points pour les phases liquide et vapeur respectivement) compares aux simulations
obtenues par Panagiotopoulos et al. [ 1988 ] (symbole triangle).
Comme le montre la Figure IV-6, la mthodologie propose permet datteindre la rgion de lazotrope en
quelques itrations, partir de laquelle une simulation plus longue permettra dchantillonner prcisment le
point azotropique.
Dans sa publication comme le montre la Figure IV-5, Panagiotopoulos et al. [ 1988 ] n'a pas effectu de
simulation prs de l'azotrope, traant seulement une ligne de guidage vers une pression d'azotrope
lgrement infrieure P* = 0,11. Nos calculs prouvent que la pression azotropique est voisine de P* = 0,115.
KL > 1
ZL
0,1
0
KL > 1
ZL
0
KL > 1
0,2
0,3
0,4
0,1
0,2
0,3
0,4
0,3
0,4
0,1
0,2
0,3
Z3
KL > 1
0,1
0,2
0,3
0,4
ZL
0,1
0,2
0,3
0,4
K4 < 1
ZL
K5 < 1
Z5 Z H
Z L ZH
0,5
0,7
Z4 ZH
0,7
0,6
0,5
0,5
KL > 1
0,8
0,9
1,0
0,8
0,9
1,0
0,8
0,9
1,0
0,9
1,0
KH < 1 ZH
0,6
0,5
0,4
0,7
0,6
0,5
K3 > 1
0,2
KH < 1 ZH
KH < 1 ZH
ZL
0,1
0,7
0,6
0,5
Z2 K2 > 1
KL > 1
Z1
K1 < 1
KH < 1
0,8
KH < 1
0,6
0,7
0,8
0,9
1,0
0,7
0,8
0,9
1,0
KH < 1
0,6
Figure IV-6. Exploration de lespace des compositions au cours des itrations pour le mlange I.
IV-13
La Figure IV-6 montre aussi que l'azotrope est approch la troisime itration et que les itrations
suivantes permettent seulement de rtrcir la fentre de composition. Ce manque d'exactitude au point
azotropique est partiellement d la symtrie particulire du mlange I qui provoquent un aplatissement
plus large qu'habituel de la courbe de bulle l'azotrope.
La prcision avec laquelle est dtermine le point azotropique a pu galement tre amliore en
augmentant le nombre de particules. En effet, dans notre cas, avec 800 particules au commencement, un
changement d'identit concernant 1 seule particule change la composition dau moins 1/800, soit 0,00125.
Cette prcision minimale est analogue lincertitude intrinsque dun appareil de mesure et sajoute donc
lerreur de mesure qui est rapporte dans les Table IV-3 et calcule comme une variance par bloc (voir
chapitre II).
Pour le mlange V, les rsultats sont prsents dans la Table IV-4. La composition azotropique est
calcule 0,695 pour le liquide et 0,6982 pour la vapeur, en accord avec les valeurs rapportes par Pandit
et Kofke [ 1999 ] savoir 0,6982 pour le liquide et 0,7094 pour la vapeur.
E.
Discussion et conclusions
Ce chapitre prsente une mthodologie pour trouver l'existence et la composition d'azotrope pour
des fluides de Lennard-Jones. Il combine efficacement des proprits des azotropes macroscopiques et
des principes de modlisation molculaire.
En effet, la constante dquilibre macroscopique change de signe par rapport lunit de part et dautre de
lazotrope. Cette proprit est dj employe dans des calculs classiques de gnie chimique avec les
modles thermodynamiques macroscopiques comme les modles de coefficient d'activit ou d'quation
d'tat. Ici, les modles macroscopiques sont remplacs par des simulations molculaires dans l'ensemble de
Gibbs pour valuer la valeur de la constante dquilibre au voisinage des corps purs. Cette recherche
grossire de la prsence dazotrope constitue la premire tape de la mthodologie.
La deuxime tape de la mthodologie cherche dterminer avec prcision la composition et la pression
de lazotrope. Elle met en uvre un processus itratif gouvern par un algorithme de type recherche par
dichotomie qui progresse le long de la courbe de coexistence liquide vapeur en ralisant successivement
des simulations dans lensemble de Gibbs Monte Carlo afin dvaluer la constante dquilibre, et des
mouvements de changement didentit afin de pouvoir modifier la composition au sein de chacune des
deux boites de simulation. La temprature azotropique est maintenue constante tout au long des
simulations.
Bien que couronne de succs, la mthode propose souffre de plusieurs insuffisances :
IV-14
1. La temprature de lazotrope est fixe et la pression est cherche, or dans les installations de
distillation rencontres en gnie chimique, cest plus souvent la pression qui est fixe et la
temprature que lon cherche.
2. L'tape d'initialisation, propose ici, ne permet pas la dtection d'un point azotropique dont les
coordonnes sont infrieures 0,025 fraction molaire en composant pur. Relaxer cette
hypothse exigerait dutiliser des nombres de molcules plus grands pour assurer un nombre
suffisant de chaque composant dans les boites lquilibre liquide vapeur.
3. Le mouvement de changement d'identit est appliqu la configuration finale obtenue par les
calculs dans l'ensemble de Gibbs Monte Carlo qui ne sont pas obligatoirement les plus stables.
Cest un choix arbitraire : on pourrait chercher raliser cette simulation en se basant sur la
configuration la plus stable ou sur une configuration choisie alatoirement parmi lchantillon
des configurations obtenues dans lensemble de Gibbs. Lincidence sur le rsultat final est
probablement ngligeable.
4. Au cours de ltape de changement didentit o le potentiel chimique est postul comme tant
constant, il na pas t vrifi a posteriori que le potentiel chimique reste constant. Ceci peut
induire une petite erreur qui est aisment compense parce que le calcul dans lensemble de
Gibbs va conduire une nouvelle galit des potentiels chimiques entre les phases.
Les insuffisances 3 et 4 montrent que nous ne progressons pas exactement le long des courbes de
coexistence liquide vapeur. Mais le point azotropique final qui est atteint est dtermin exactement par
le dernier calcul dquilibre liquide vapeur dans lensemble de Gibbs Monte Carlo.
Lextension de cette mthodologie pour traiter des fluides rels sera possible avec quelques modifications,
notamment le remplacement du mouvement de changement didentit des molcules A en B par la
destruction de A suivie de la cration de B. En effet, pour des molcules polyatomiques le changement
didentit est quasi impossible. La destruction et la cration devront aussi tre ralise avec un biais
configurationnel (chapitre II) pour avoir des chances de succs raisonnables.
IV-15
Chapitre V.
Le chapitre III contribue llaboration dun champ de force gnrique en optimisant des paramtres de
Lennard Jones de la contribution de Van der Waals pour le groupe nitrile. Des simulations dquilibre
liquide vapeur dans lensemble de Gibbs Monte Carlo volume impos (NVT) et des simulations de la
phase liquide condense dans lensemble isobare isotherme (NPT) sont ralises avec le code Gibbs du
laboratoire de Chimie Physique dOrsay (LCP) et lInstitut Franais du Ptrole (IFP). Le caractre polaire
des nitriles ncessite de prendre en compte la contribution lectrostatique au travers de charges
ponctuelles issues, dans ce travail, dune analyse de population de type MEP ou Mulliken.
Excellent rsultat pour lactonitrile mais incohrence pour les autres nitriles . Telle pourrait tre
lapprciation des rsultats obtenus en optimisant des paramtres de Lennard Jones avec les charges
ponctuelles MEP. La procdure doptimisation des paramtres donne effectivement des valeurs de
paramtres epsilon et sigma du potentiel de Lennard Jones pour le carbone et lazote du groupe nitrile
(CN) qui permettent de reproduire avec prcision les donnes exprimentales disponibles pour
lactonitrile lorsquon les combine avec des paramtres de Lennard Jones gnriques du champ de force
AUA4 et des constantes intramolculaires issues de calculs quantiques. Cependant, la transfrabilit des
paramtres associs aux charges MEP pour prdire les proprits du propionitrile est mise en dfaut. Cela
nous a conduit raliser une autre optimisation avec des charges ponctuelles de Mulliken.
Bien mais peut mieux faire . Telle pourrait tre lapprciation des rsultats obtenus en optimisant des
paramtres de Lennard Jones avec les charges ponctuelles Mulliken. Certes, les paramtres sont
transfrables et donnent des prdictions en bon accord avec les maigres donnes exprimentales
disponibles pour le propionitrile et le n-butyronitrile. Mais, les prdictions des proprits de lactonitrile
V-1
avec les charges Mulliken sont globalement moins bonnes quavec les charges MEP, notamment en ce qui
concerne la phase liquide et cela est flagrant puisque la temprature critique est systmatiquement
surestime avec les charges Mulliken denviron 20K pour les trois nitriles tudis. En termes de valeurs
numriques des paramtres epsilon et sigma de Lennard Jones obtenus pour le carbone et lazote du
groupe nitrile, les valeurs avec les charges de Mulliken sont plus cohrentes que celles obtenues avec les
charges MEP lorsquon les compare aux autres groupes du champ de force AUA4.
A posteriori, nous avons constat que la distribution de charge ponctuelle MEP pour le propionitrile nest
pas cohrente avec celle de lactonitrile, ce qui pourrait expliquer la non transfrabilit du jeu de
paramtre MEP malgr son excellente reprsentation des proprits de lactonitrile.
Par consquent, nous pouvons proposer deux critres pour lobtention de paramtres gnriques :
1. homognit des distributions de charges sur les sries chimiques homologues lors de leur
obtention par analyse de population des surfaces de potentiel lectrostatique. Lanalyse de
Mulliken semble tre trop simpliste et se traduit par des rsultats bons mais pas excellents.
2. vrification de la pertinence physique des paramtres optimiss par comparaison dautres
valeurs obtenues pour le mme champ de force.
Le chapitre IV propose une mthodologie pour la prdiction des azotropes de mlanges binaires. Au-del
de leur signification de non idalit du mlange, limportance des azotropes est reconnue puisquelle
conditionne directement la faisabilit et la conception des procds de distillation. En particulier,
essentielle est la capacit de prvoir si un mlange donn forme un ou plusieurs azotropes et de calculer
les conditions appropries de composition, temprature et pression qui y sont associes. La mthodologie
est base sur une combinaison de proprits macroscopiques des azotropes frquemment utilises en
Gnie des Procds et de proprits de thermodynamique statistique permettant dchantillonner
correctement un systme molculaire possdant un azotrope en progressant le long de la courbe de
coexistence liquide vapeur.
La premire tape exploite le fait que la constante dquilibre macroscopique change de signe par rapport
lunit de part et dautre de lazotrope par des simulations molculaires dans l'ensemble de Gibbs Monte
Carlo volume impos (NVT) pour valuer la valeur de la constante dquilibre au voisinage des corps
purs. Cette recherche grossire de la prsence dazotrope constitue la premire tape de la mthodologie.
La deuxime tape de la mthodologie cherche dterminer avec prcision la composition et la pression
de lazotrope. Elle met en uvre un processus itratif gouvern par un algorithme de type recherche par
dichotomie qui progresse le long de la courbe de coexistence liquide vapeur en ralisant successivement
des simulations dans lensemble de Gibbs Monte Carlo afin dvaluer la constante dquilibre, et des
mouvements de changement didentit afin de pouvoir modifier la composition au sein de chacune des
boites de simulation. Nous avons dmontr que la probabilit dacceptation de ce mouvement est gale au
V-2
produit des deux mouvements classiques dinsertion et de destruction spcifiques lensemble grand
canonique (VT). La temprature azotropique est maintenue constante tout au long des simulations.
Bien que couronne de succs, la mthode propose souffre de plusieurs insuffisances :
1. La temprature de lazotrope est fixe et la pression est cherche, or dans les installations de
distillation rencontres en gnie chimique, cest plus souvent la pression qui est fixe et la
temprature que lon cherche.
2. Les bornes de ltape d'initialisation devront dans certains cas tre rapproches des corps purs,
ce qui peut exiger dutiliser des nombres de molcules plus grands pour assurer un nombre
suffisant de chaque composant dans les boites lquilibre liquide vapeur.
3. Le choix des configurations employes entre chaque changement de technique de simulation
est un choix arbitraire : on pourrait chercher raliser les mouvements de changement
didentit sur la configuration la plus stable ou sur une configuration choisie alatoirement
parmi lchantillon des configurations obtenues dans lensemble de Gibbs Monte Carlo.
Lincidence sur le rsultat final est probablement ngligeable.
4. Au cours de ltape de changement didentit o le potentiel chimique est postul comme tant
constant, il na pas t vrifi a posteriori que le potentiel chimique reste constant. Ceci peut
induire une erreur qui est aisment compense parce que le calcul dans lensemble de Gibbs
Monte Carlo va conduire une nouvelle galit des potentiels chimiques entre les phases.
Les insuffisances 3 et 4 montrent que nous ne progressons pas exactement le long des courbes de
coexistence liquide vapeur. Mais de toute faon, le point azotropique final qui est atteint est exactement
dtermin par le dernier calcul dquilibre liquide vapeur dans lensemble de Gibbs Monte Carlo.
En perspective, nous pouvons citer pour chacune des deux contributions respectivement :
-
Vrifier les charges MEP du propionitrile avant de raliser une nouvelle optimisation de
paramtres avec les charges NBO prsentes en annexe. Les charges NBO sont plus fortes que
les charges MEP sur les groupes mthyle (latome C1 et ses hydrognes) mais infrieures pour
le groupe nitrile, tout en tant suprieures aux charges de Mulliken (voir table 1 du chapitre III
et table 1 en annexe IV). Mettre en uvre dautres corrections voques dans ce chapitre
comme lutilisation de la sommation Ewald ne semble pas une priorit court terme.
Lextension de la mthodologie de recherche des azotropes pour traiter des fluides rels. Elle
sera possible avec quelques modifications, notamment le remplacement du mouvement de
changement didentit par une combinaison des deux mouvements classiques de lensemble
grand canonique : insertion destruction raliss avec un biais configurationnel pour avoir des
chances de succs raisonnables.
V-3
Rfrences
Agrawal R. and Kofke D.A., 1995, thermodynamic and structural properties of model systems at solidfluid coexistence : II. Melting and sublimation of the Lennard-Jones system. Mol. Phys., 85, 43-59.
Allen M.P. and D.J. Tildesley, 1987, Computer Simulation of Liquids, Oxford University Publications, ISBN
0-19-855645-4, New York, Royaume-Uni. [nouvelle dition en 2000].
Berendsen H.J.C., Postma J.P.M., Van Gunsteren W.F. and Hermans H.J., 1983, In intermolecular forces,
Pullman, B., Ed., Dordrecht.
Bhm H.J., McDonald I.R. and Madden P.A., 1983, An effective pair potential for liquid acetonitrile, Mol.
Phys., 49(2), 347-360.
Bhm H.J., Lynden-Bell R.M., Madden P.A. and McDonald I.R., 1984, Molecular motion in a model for
liquid acetonitrile, Mol. Phys., 51(3), 761-777.
Bourasseau E., Ungerer P., Boutin A., Fuchs A.H., 2002a, Monte Carlo simulation of branched alkanes
and long chain n-alkanes with ansotropic united atoms intermolecular potential, Mol. Sim., 28(4),
317-336.
Bourasseau E., Ungerer P., Boutin A., 2002b, Prediction of Equilibrium Properties of Cyclic Alkanes by
Monte Carlo Simulations New Anisotropic United Atoms Intermolecular Potential New
Transfer Bias Method, J. Phys. Chem. B, 106, 5483-5491.
Bourasseau E., Haboudou, M., Boutin A., Fuchs A.H., Ungerer P., 2003, New optimization method for
intermolecular potentials: Optimization of a new anisotropic united atoms potential for olefins:
Prediction of equilibrium properties, J. Chem. Phys., 118(7), 3020-3034.
Bourasseau E., 2003, Prdiction des proprits thermodynamiques dquilibres de fluides par simulation
molculaire, rapport de thse, Universit de Paris-sud.
Bukowsky R., Szalewicz K., Chabalaowski C.F., Ab initio interaction potentials for simulations of
dimethylnitramine solutions in supercritical carbon dioxide with cosolvents, 1999, J. Phys. Chem. A.,
103, 7322-7340.
Case F., Chaka A., Friend D.G., Frurip D., Golab J., Johnson R., Moore J., Mountain R.D., Olson J.,
Schiller M. and Storer J., 2004, The first industrial fluid properties simulation challenge, Fluid Phase
Equilibria, 217(1), 1-10. et http://www.cstl.nist.gov/FluidSimulationChallenge/
Chakhmuradov C.G. and Guseinov S.O., 1984, Iz. Vys. Uc. Zav., 65-69. (en russe)
Chapman R.G. and Goodwin S.P., 1993, A general algorithm for the calculation of azeotropes in fluid
mixtures, Fluid Phase Equilibria, 85, 55-69.
Chen B. and Siepmann J.I., 1999, Transferable potentials for phase equilibria. 3. Exliplcit-hydrogen
description of normal alkanes, J. Phys. Chem. B., 103, 5370-5379.
Chen J., Fisher K., Gmehling J., 2002, Modification of PRSK mixing rules and results for vapor-liquid
equilibria, enthalpy of mixing and activity coefficients at infinite dillution, Fluid Phase Equilibria, 200,
411-429.
Chen C.C. and Mathias P.M., 2002, Applied Thermodynamics for Process Modelling, AIChE Journal,
48(2), 194-200.
COMSEF. Computational Molecular Science and Engineering Forum. Cours en ligne relatifs la
simulation molculaire. http://www.ecs.umass.edu/che/am3/AIChE.html
Constantinou L. and Gani R., 1994, New Group Contribution Method for Estimating Properties of Pure
Compounds, AIChE Journal, 40(10), 1687-1710.
Constantinou L., Bagherpour K., Gani R., Klein J.A., Wu D.T., 1996, Computer Aided Product Design:
Problem Formulations, Methodology and Applications, Comp. Chem. Eng., 20(6/7), 685-702.
Contreras O., 2002, Determinacion del equilibrio liquido-vapor de agua, aromaticos y sus mezclas
mediante simulacion molecular, thse de lUniversit de Tarragone, Espagne.
Cornell W.D., Cieplak P., Bayly C.L., Gould I.R., Merz K.M., Ferguson D.M, Spellmeyer D.C, Fox T.,
Caldwell J.W. and Kollman P.A, 1995, A second generation force fields for the simulation of
proteins, nucleic acids and organic molecules, J. Am. Chem. Soc., 117, 5179-5197.
Delhommelle J., Granucci G., Brenner V., Millie P., Boutin A., Fuchs A.H., 1999, A new method for
deriving atomic charges and dipoles for alkanes: investigation of transferability and geometry
dependence, Mol. Phys., 97(10), 1117-1128.
Delhommelle J., Tschirwitz C., Ungerer P., Granucci G., Millie P., Pattou D., Fuchs A.H., 2000,
Derivation of an Optimized Potential Model for Phase Equilibria (OPPE) for Sulfides and Thiols. J.
Phys. Chem. B, 104, 4745-4753.
De Pablo J.J. and Escobedo F.A., 2002, Perspective: Molecular Simulations in Chemical Engineering:
Present and Future, AIChE Journal, 48 (12), 2716-2721.
Diu B., Guthmann C., Lederer D., Roulet B., 1989, Physique statistique, Hermann, ISBN 2 7056 6065 8,
Paris.
Doherty M.F. and Perkins J.D., 1979, On the Dynamics of Distillation Process - III. The Topological
Structure of Ternary Residue Curve Maps, Chem. Eng. Sci., 34, 1401-1414.
Dugas H., 2000, Principes de base en modlisation molculaire. Aspects thoriques et pratiques, 5ime Edition,
disponible sur le WEB ladresse : http://www.centrcn.umontreal.cal~dugas.
Eckert E. and Kubicek M., 1996, Computing Heterogeneous Azeotropes in Multicomponent Mixtures.
Comp. Chem. Eng., 21, 347-350.
Errington J.R. and Panagiotopoulos A.Z., 1999, A New Potential Model for the n-Alkanes Homologous
Series, J. Phys. Chem. B, 103, 6314-6322.
Escobedo F.A. and de Pablo J.J., 1997, Pseudo-ensemble simulations and Gibbs-Duhem integrations for
polymers, J. Chem. Phys., 106, 2911-2923.
Escobedo F.A., 1998, Novel Pseudo-ensembles for Simulation of Multicomponent Phase Equilibria, J.
Chem. Phys., 108, 8761-8772.
Fidkowski, M.F., Malone, M.F. and Doherty M.F., 1993, Computing Azeotropes in Multicomponent
Mixtures, Comp. Chem. Eng., 17, 1141-1155.
Francesconi A.Z., Franck E.U., Lentz H., 1975, Pressure volume temperature data for acetonitrile up
to 450 C and 2500 bar, Ber. Bunsen-Ges. Phys. Chem., 79, 897-901. (en allemand)
Frenkel D. and Smit B., 1996, Understanding Molecular Simulation. From Algorithms to Applications, Academic
Press, ISBN 0-12-267370-0, San Diego. [nouvelle dition en 2002].
Fuchs A., Boutin A., Rousseau B., 1997, Simulation molculaire et Gnie des Procds, Entropie, 208, 5-12.
Gerbaud V., 2003, Apport de la simulation molculaire en gnie des procds. Habilitation Diriger les
Recherches, Institut National Polytechnique de Toulouse, France.
Gmehling J., Menke J., Krafczyk J. and Fischer K., 1994, Azeotropic Data; VCH Editor: Weinheim.
Goldstein E., Buyong M., Lii J.H. and Allinger N.L., 1996, Molecular Mechanics Calculations (MM3) on
Nitriles and Alkynes, J. Phys. Org. Chem., 9, 191-202.
Goujon F., Malfreyt P., Boutin A. and Fuchs A.H., 2001, Vapour liquid phase equilibria of n-alkanes by
direct Monte Carlo simulations, Mol. Simul., 27, 99-114.
Goujon F., Malfreyt P., Boutin A. and Fuchs A.H., 2002, Direct Monte Carlo simulations of the
equilibrium properties of n-pentane liquid vapour interface, J. Chem. Phys., 116 (18), 8106-8117.
Harding S.T., Maranas C.D., McDonald C.M. and Floudas C.A., 1997, Locating all homogenous
Azeotropes in Multicomponent Mixtures, Ind. Eng. Chem. Res., 36, 160-178.
Harper P.M., Gani R., Kolar P., Ishikawa T., 1999, Computer Aided Molecular Design with combined
molecular modeling and group contribution, Fluid Phase Equilibria, 158-160, 337-347.
Hloucha M., Deiters U.K., 1997, Monte Carlo simulations of acetonitrile with an anisotropic polarizable
molecular model, Mol. Phys., 90(4), 593-597.
Hloucha M., Sum A.K. and Sandler S.I., 2000, Computer Simulation of Acetonitrile and Methanol with abinitio Based Pair Potential, J. Chem. Phys., 113(13), 5401-5406.
Jorgensen W.L., Chandrasekhar J., Madura J.D., Impey R. W. and Klein M.L., 1983, Comparison of
simple potential functions for simulating liquid water, J. Chem. Phys., 79, 926-935.
Jorgensen W.L., Madura J.D. and Swenson C.J., 1984, Optimized Intermolecular Potential Functions for
Liquid Hydrocarbons, J. Am. Chem. Soc., 106, 6638-6646.
Jorgensen W.L., Briggs J.M., 1988, Monte Carlo simulations of liquid acetonitrile with a three-site model,
Mol. Phys., 63(4), 547-558.
Kiva V.N., Hilmen E.K. and Skogestad S., 2003, Azeotropic phase equilibrium diagrams: a survey, Chem.
Eng. Sc., 58, 1903-1953.
Klamt A. and Eckert F., 2000, COSMO-RS: a novel and efficient method for the a priori prediction of
thermophysical data of liquids, Fluid Phase Equilibria, 172, 43-72.
Kofke D.A., Glandt E.D., 1988, Monte Carlo simulation of multicomponent equilibria in a semigrand
canonical ensemble, Mol. Phys., 64(6), 1105-1131.
Kofke D.A., 1993a, Gibbs-Duhem integration : a new method for direct evaluation of phase coexistence
by molecular simulation. Mol. Phys., 78, 1331-1336.
Kofke D.A., 1993b, Direct evaluation of phase coexistence by molecular simulation via integration a long
the coexistence line. J. Chem. Phys., 98, 4149-4162.
Kratzke H. and Mller S., 1985, Thermodynamic properties of acetonitrile 2. (P, , T) of saturated and
compressed liquid acetonitrile, J. Chem. Thermodynamics, 17, 151-158.
Lagache M., 2003, La simulation de Mont Carlo et lindustrie ptrolire : dveloppement de potentiel
pour les composs organomercurs, Calcul de grandeurs thermodynamiques drives de gaz
condensat, rapport de thse, Universit de Paris-sud.
Leach A.R., 1996, Molecular Modelling. Principles and Applications, Longmann, Harlow, ISBN 0-582-23933-8,
Royaume-Uni.
Mackerell A.D., Wiorkiewiczkuczera J. and Karplus M., 1995, All-atom empirical energy function for the
simulation of nucleic acids, J. Am. Chem. Soc., 117, 11946-11975.
Martin M.G. and Siepmann J.I., 1998, Transferable Potentials for Phase Equilibria. 1. United Atom
Description of n-Alkanes, J. Phys. Chem. B, 102(14), 2569-2577.
Mc Quarry D.A., 1976, Statistical Mechanics, Harper and Collins Publishers, New York, ISBN 06-044366-9,
Etats-Unis.
Metropolis N., Rosenbluth A.W., Rosenbluth M.N., Teller A.H. and Teller E., 1953, Equation of state
calculations by fast computing machines, J. Chem. Phys., 21, 1087-92.
Mota J.P.B., 2003, Towards the atomistic description of equilibrium-based Separation Processes. 1.
Isothermal Stirred-Tank Adsorber., ESCAPE 13 Proceedings, Computer-Aided Chemical
Engineering, 14, 791-796, Eds. Kraslawski A., Turunen I., Elsevier Science B.V., Amsterdam.
Mousa A.H.N., 1981, Vapour pressure and saturated-vapour volume of acetonitrile, J. Chem. Therms.,
13, 201-202.
Mller E.A. and Gubbins K.E., 2001, Molecular-Based Equations of State for Associating Fluids: A
Review of SAFT and Related Approaches, Ind. Eng. Chem. Res., 40, 2193-2211.
Nath S.K., Escobedo F.A. and De Pablo J.J., 1998, On the Simulation of Vapor-Liquid Equilibrium for
Alkanes, J. Chem. Phys., 108, 9905-9911.
Nicolas C, 2000, Modlisation et simulation molculaire de leau : prise en compte explicite de lnergie
dinduction, Rapport de stage DEA, Universit de Paris-sud.
NIST chemistry webbook, 2004, http://webbook.nist.gov/chemistry/
Panagiotopoulos A.Z., 1987, Direct determination of phase coexistence properties of fluids by Monte
Carlo in a new ensemble, Mol. Phys., 61, 813-826.
Panagiotopoulos A.Z., Quirke N., Stapleton M., Tildesley D.J., 1988, Phase equilibria by simulation in the
Gibbs Ensemble Alternative derivation, generalization and application to mixture and membrane
equilibria, Mol. Phys., 63(4), 527-545.
Panagiotopoulos A.Z., 2000, Monte Carlo Methods for Phase Equilibria, J. Phys. Cond. Matter, 12, R25R52.
Panagiotopoulos A.Z., 2001, Force Field Development for Simulations of Condensed Phases. In First
International Conference on Foundations of Molecular Modeling and Simulation, AIChE Symp. Ser.,
97, 61-70.
Pandit S.P. and Kofke D.A., 1999, Evaluation of a Locus of Azeotropes by Molecular Simulation, AIChE
Journal, 45(10), 2237-2244.
Peneloux A., Rauzy E. and Freze R, 1982, A consistent correction for Redlich-Kwong-Soave volumes,
Fluid Phase Equilibria, 8, 7-23.
Prausnitz J.M., Lichtenthaler R.N. and de Azevedo E.G., 1999, Molecular Thermodynamics of Fluid Phase
Equilibrium. 3rd edition. Prentice Hall International, Upper Saddle River. ISBN 0-13-977745-8.
Sandler S.I., 1994, Models for thermodynamic and Phase Equilibria Calculations, Marcel Dekker Inc., New York,
ISBN 0-8247-9130-4, Etats-Unis.
Sandler S.I., 2003, Quantum mechanics: a new tool for engineering thermodynamics, Fluid Phase Equilibria,
201, 147-160.
Smit B., de Smedt Ph. and Frenkel D., 1989, Computer simulations in the Gibbs ensemble., Mol. Phys., 68,
931-950.
Smit B., Karaboni S., Siepmann J.I., 1995, Computer simulation of vapor liquid phase equilibria of nalkanes, J. Chem. Phys., 102, 2126-2140.
Shulgin I., Fischer K., Noll O. and Gmehling J., 2001, Classification of Homogeneous Azeotropes, Ind.
Eng. Chem. Res., 40, 2742-2747.
St-Amant A., Salahub D.R., 1990, New algorithm for the optimization of geometries in local density
functional theory, Chem. Phys. Letters, 169, 387-392.
St-Amant A., 1991, rapport de thse, Universit de Montral.
Sun H., 1998, COMPASS: an ab initio force field optimized for condensed-phase applications-overview
with details on alkane and benzene compounds, J. Phys. Chem. B., 102, 7338-7364.
Teja A.S. and Rowlinson J.S., 1973, The prediction of the thermodynamic properties of fluids and fluid
mixtures IV. Critical and azeotropic states, Chem. Eng. Sci., 28, 529-538.
Toxvaerd S., 1990, Equations of state of alkanes I. J. Chem. Phys., 93, 4290-4295.
Toxvaerd S., 1997, Equations of state of alkanes II. J. Chem. Phys., 107, 5197-5204.
Ungerer P., 1999, Calcul des proprits dquilibre et de transport de mlanges dhydrocarbures par
simulation molculaire, rapport de thse, Universit de Paris-sud.
Ungerer P., Boutin A. and Fuchs A. H., 1999, Direct calculation of bubble points by Monte Carlo
simulation, Mol. Phys., 97(4), 523-539.
Ungerer P., Beauvais C., Delhommelle J., Boutin A., Rousseau B. and Fuchs A.H., 2000, Optimisation of
the anisotropic united atoms intermolecular potentiel for n-alkanes. J. Chem. Phys., 112(12), 54995510.
Ungerer P., Boutin A. and Fuchs A. H., 2001, Direct calculation of bubble points for alkane mixtures by
molecular simulation, Mol. Phys., 99(17), 1423-1434.
Vidal J., 1997, Thermodynamique. Application au Gnie Chimique et lindustrie ptrolire, dition Technip, ISBN
2-7108-0715-7.
Warowny W., 1994, Volumetric and phase behavior of acetonitrile at temperatures from 363K to 463K, J.
Chem. Eng. Data, 39, 275-280.
Wasylkiewicz S.K., Doherty M.F. and Malone M.F. 1999, Computing All Homogeneous and
Heterogeneous Azeotropes in Multicomponent Mixtures. Ind. Chem. Eng. Res., 38, 4901-4912.
Widagdo S. and Seider W.D., 1996, Azeotropic Distillation, AIChE J., 42, 96-146.
Widom B., 1963, Some topics in the theory of fluids. J. Chem. Phys., 39, 2802-2812.
Lapproximation de Stirling
ln(m) .
1
m =1
ln ( N !) = ln (m ) ln ( x )dx = N ln ( N ) N
Qui est lapproximation de Stirling de ln ( N !) .
Remarques
Ici la borne infrieure de lintgrale a tait prise nulle puisque N est trs grand, et on sait que x. ln ( x ) tend
vers zro quand x tend zro.
1
Un calcul plus affin de cette approximation donne : ln ( N !) N ln ( N ) N + 1 ln (2 N ) 2 mais ce
AI.2
n1 + n2 + n3 + L = N
Nous prsentons ici une dmonstration de la formule utilise dans notre dmonstration. Tout dabord
calculons le nombre de permutations possibles de N objets distincts, ce nombre est donn par :
Ensuite, calculons le nombre de possibilits de partager les N objets en deux groupes diffrents, lun
contenant N1 objets et lautre le reste, c'est--dire N 2 = N N 1 .
Annexe I-1
Il y a [N .(N 1)(
. N 2 )L (N (N1 1))] manires de former le premier groupe, et N 2!= ( N N1 )!
manires de former le second. Le nombre total est donc le produit :
N .(N 1)(
. N 2 )L ( N ( N1 1)) ( N N1 )!= N !
( N N1 )!= N !
( N N1 )!
Mais de cette manire, nous avons dnombr tous les cas radicalement en faisant mme tenir compte de
lordre des diffrents lments dans les diffrents groupes. Mais puisque ce dtail est immatriel, il suffit
de diviser le nombre total (N !) par N 1! N 2 ! (o N 1! reprsente lordre dans le premier groupe et
N!
N1!( N N1 )!
= N!
N1!.N 2 !
N!
=
N 1! N 2 !L N r !
N!
r
N !
i
i =1
Annexe I-2
=
liq
M N liq
V liq N A
) NN
A
liq
Annexe II-3
et
int
ext
H vap = U vap
+ U vap
+ Pvap V vap
) NN
A
vap
Uint sont les nergies intramolculaire et Uext les nergies intermolculaires des phases liquide et vapeur
selon lindice.
En supposant que le volume de la phase liquide ngligeable devant celui de la phase vapeur
( V liq <<< V vap ), en ngligeant lnergie intermolculaire en phase vapeur (condition dun gaz parfait,
Uext=0 qui peut sappliquer aux gaz rels faible pression) et en considrant lgalit de la contribution
int
= U liqint , lexpression de lenthalpie de vaporisation se
intramolculaire en phase liquide et vapeur U vap
rduit :
H vap U liq + RT
Ceci signifie quon peut dterminer, dans ces conditions, les valeurs de la masse volumique liquide et de
lenthalpie de vaporisation par une seule simulation NPT en phase liquide.
3. Pression de vapeur saturante
La pression est value au cours de la simulation laide du thorme du viriel :
P=
1
3V
dU ij
3 NkBT ri r j
dr
i
j
ij
Annexe II-4
Le code Gibbs utilis dans sa version 2.5.1, conu par le laboratoire de chimie physique de
lUniversit de Paris XI dOrsay, est destin calculer les proprits l'quilibre de fluides par simulation
molculaire en utilisant une mthode de Monte Carlo. Ses principales possibilits sont les suivantes :
calcul de proprits monophasiques par une simulation de Monte Carlo dans les
ensembles statistiques NVT, NPT.
Les fluides qu'il est possible de prendre en compte sont des corps purs ou des mlanges dans lesquels les
molcules peuvent tre rigides ou flexibles.
Les formes suivantes d'nergie peuvent tre prises en compte :
nergie de polarisation (ou induction), c'est--dire l'action du champ lectrique sur les
centres de Lennard-Jones, considrs comme polarisables.
Dans le cas des molcules rigides, les positions des centres de Lennard-Jones et des charges
lectrostatiques peuvent tre choisies librement.
Dans le cas des molcules flexibles, construites sur le modle des alcanes, les centres de Lennard-Jones
peuvent tre placs sur les atomes de carbone (potentiel d'Atomes Unifis ou UA) ou sur la bissectrice
externe des angles C-C-C (potentiel d'Atomes Unifis Anisotropes ou AUA). La distance entre deux
centres de force voisins est constante pour une mme molcule. Les charges lectrostatiques et les diples,
Annexe III-5
si l'on veut les prendre en compte, sont localises sur les carbones et les hydrognes ; la distance C-H tant
impose.
Afin de gnrer les configurations, les mouvements suivants sont implants :
rotation d'un centre autour de ses plus proches voisins pour les molcules flexibles ;
changements de volume ;
Les mouvements de recroissance partielle et de transfert mettent en uvre des techniques de biais
configurationnel pour les molcules flexibles. L'amplitude des mouvements de translation, de rotation et
de changement de volume est optimise au cours de la simulation pour obtenir le taux moyen
d'acceptation dsir.
Le code Gibbs, crit en Fortran 90, fonctionne en mode batch : il lit les donnes dans trois fichiers
d'entre (control.dat, potparam.dat et molecule.dat dont nous prsenterons des exemples ci-aprs) et crit
des fichiers de rsultats. Pour visualiser et traiter les rsultats, nous avons eu recours aux utilitaires
comme :
plotmtv et xmgrace pour afficher l'volution des grandeurs telles que la masse
volumique, la pression, lnergie, etc.
AIII.2
Organigramme simplifi
Nous avons prsent ci-dessous un organigramme simplifi qui met en vidence les diffrentes
tapes de calcul suivies par le code, et ladaptation de lalgorithme de Metropolis lensemble de Gibbs.
Annexe III-6
Initialisation
Translation et rotation
des centres de forces
Rejet/acceptation de la
configuration selon le
critre de Metropolis.
Changement de
volume
Moyenne statistique,
Stockage des informations
Graphiques
Non
N
itrations
Oui
Traitement
statistique final
Visualisation
des boites.
AIII.3
Fichiers dentre
AIII.3.a. molecule.dat
Le fichier aceto.dat avec les charge Mulliken, lorsque la molcule est considre flexible.
3
1 1 2
2 2 1
3 1 2
6
-0.40
0.16
0.16
0.16
0.15
-0.23
1
1
8
8
3
0
2
2
2
0
0
Annexe III-7
Le fichier aceto.dat avec les charge MEP, lorsque la molcule est considre rigide.
3
! Nombre de sites Lennard-Jones
1
-1.192 0.151 -0.472
1
2
-1.192 1.619 -0.473
8
3
-1.192 2.796 -0.472
8
6
! Nombre de charges ponctuelles
-0.62
-1.192 0.151 -0.472
0.47
-1.192 1.619 -0.473
0.21
-0.157 -0.222 -0.471
0.21
-1.708 -0.224 0.424
0.21
-1.709 -0.224 -1.367
-0.48
-1.192 2.796 -0.472
1
AIII.3.b. control.dat
.true.
1
1
.false.
!
!
!
!
Annexe III-8
AIII.3.c. potparam.dat
1
3.6072
15.03
120.15
2.22
2
3.2183
12.011
95.519
3
3.5638
14.0067
162.41
! nbre de dist
1 2
1.457
! entre i et j
2 3
1.168
! entre i et j
! nbre de ben
1 2 3
1
24300.0
180.0
0
8
1
2
3
4
5
6
7
8
! nbr de torsion
0.21584
0.38400
0.336
0.64599
0.000
0.000
0.000
0.000
!
!
!
!
!
!
!
!
!
number of
AUA tupe,
AUA type,
AUA type,
AUA type,
AUA type,
AUA type,
AUA type,
AUA type,
Annexe III-9
Annexe IV.
Annexe IV-10
des diffrences significatives pour lanalyse de Mulliken qui nest pas explique (mme nom de
la fonctionnelle, mme base dorbitale atomique) si ce nest que les logiciels sont diffrents.
une homognit entre les diffrents nitriles pour les charges sur chaque atome. Ainsi pour la
fonctionnelle BP86, la charge NBO sur latome dazote varie entre 0,3 et 0,31 pour les trois
nitriles linaires.
linsensibilit des charges NBO en fonction de la fonctionnelle, tandis que les charges Mulliken
y sont sensibles.
Les charges NBO sont plus fortes que les charges MEP sur les groupe mthyle (latome C1 et
les hydrognes) mais infrieure pour le groupe nitrile, tout en tant suprieure aux charges de
Mulliken (voir table 1 chapitre II).
In fine, les charges NBO reprsentent un ensemble cohrent, qui satisfait donc un des critres de
transfrabilit noncs dans le chapitre III et, par sa nature moins simpliste que lanalyse Mulliken
mriterait dtre teste pour ralise des simulations molculaires des nitriles.
Annexe IV-11