Perméabilité des gaz dans les polymères
Perméabilité des gaz dans les polymères
3 Méthodes retenues 37
3.1 Mécanique statistique et modélisation moléculaire . . . . . . . . . . . 37
3.2 Méthode Monte Carlo . . . . . . . . . . . . . . . . . . . . . . . . . . 38
3.2.1 L’algorithme de Metropolis . . . . . . . . . . . . . . . . . . . . 38
3.2.2 Mouvements Monte Carlo . . . . . . . . . . . . . . . . . . . . 39
3.2.3 Les ensembles statistiques . . . . . . . . . . . . . . . . . . . . 42
3.2.4 Le potentiel chimique . . . . . . . . . . . . . . . . . . . . . . . 45
3.3 Dynamique moléculaire . . . . . . . . . . . . . . . . . . . . . . . . . . 45
3.3.1 Équations du mouvement . . . . . . . . . . . . . . . . . . . . 46
3.3.2 Thermostat de Nosé-Hoover . . . . . . . . . . . . . . . . . . . 47
3.3.3 Intégration des équations du mouvement . . . . . . . . . . . . 48
3.4 Modèles moléculaires . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
3.4.1 Potentiels inter-moléculaires . . . . . . . . . . . . . . . . . . . 50
3.4.2 Potentiels intra-moléculaires . . . . . . . . . . . . . . . . . . . 54
3.4.3 Paramètres de potentiel . . . . . . . . . . . . . . . . . . . . . 55
3.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
4 De l’amorphe au semi-cristallin 59
4.1 Données expérimentales de solubilité . . . . . . . . . . . . . . . . . . 62
4.2 Calcul de la solubilité dans l’ensemble osmotique . . . . . . . . . . . . 63
4.3 Gonflement du polymère en présence du gaz . . . . . . . . . . . . . . 66
4.4 Calcul de la solubilité dans l’ensemble Grand Canonique . . . . . . . 68
Table des matières 3
Bibliographie 131
1 Introduction
1.1 Motivation industrielle
La perméabilité des matériaux polymères aux gaz et aux liquides est une pro-
priété qui est mise à profit dans de nombreux domaines industriels. Grâce à leur
caractère imperméable, les polymères rentrent dans la fabrication de la vaisselle, des
réservoirs de liquides, des vêtements, etc. Des polymères relativement imperméables
comme le polyéthylène, le polypropylène ou le polyamide sont par exemple utilisés
dans l’emballage alimentaire (cf. figure 1). Dans le domaine médical les polymères
sont utilisés pour isoler le médicament avant l’absorption finale dans le corps. Dans
le secteur du transport d’énergie par exemple, l’étanchéité des tuyaux est assurée
par des polymères comme le polyéthylène haute densité, le polyamide 11 ou le po-
lyfluorure de vinylidène.
Dans tous les exemples précédemment cités, moins les polymères sont perméables
plus ils sont sollicités comme barrière étanche. Cependant, il y a des applications
dans lesquelles la perméabilité des polymères vis-à-vis de certains gaz est favorisée.
C’est le cas de la séparation des gaz par les membranes polymères. Dans ce cas les po-
lymères sont utilisés sous forme de membranes fines perméables pour certains gaz et
étanches pour d’autres gaz. La séparation des mélanges CO2 /H2 par les membranes
de polyéthylène terephthalate en est un exemple. Les membranes polymères sont
également utilisées pour la désalinisation de l’eau. L’utilisation de polymère dans
toutes ces applications demande une connaissance approfondie de la perméabilité des
gaz dans ce type de matériau. Cela peut se faire par le biais d’études expérimentales
et théoriques sur les systèmes gaz/polymère.
Fig. 1 – L’application des polymères dans la vie quotidienne pour leur étanchéité a)
à gauche : l’emballage des produits alimentaires par le polyéthylène b) au centre :
les vaisselles à base de polypropylène c) à droite : isolation du fromage par deux
couches de polymère. La couche extérieure est du polyamide pour empêcher l’air
d’entrer dans le sac. Le polyéthylène est utilisé dans la couche intérieure pour être
en contact avec la nourriture.
1.2 Contexte de la thèse 2
c
Fig. 2 – Les différents types de plateformes sous-marines BOEMRE
La pression autour d’une conduite flexible est très grande, d’une part à cause de
la pression hydrostatique de l’eau profonde à l’extérieure de la conduite et d’autre
part à cause de la pression des gaz contenus dans le brut à l’intérieur. La pression
peut valoir jusqu’à 100 MPa, et la température peut aller de 4˚C en eaux profondes
jusqu’à 180 ˚C pour certains puits. La problématique est l’étanchéité des flexibles
dans ces conditions de pression et de température à la fois vis-à-vis de l’eau externe
et vis-à-vis du gaz interne, ainsi que de garder la conduite toujours flexible. Pour
cela différents polymères sont employés dans la fabrication des conduites flexibles
1.2 Contexte de la thèse 3
c
Fig. 3 – Structure classique d’une conduite flexible Technip
J = −D∇C (2)
où C est la concentration locale de gaz dissous dans la membrane polymère et D est
le coefficient de diffusion fickien. La concentration de gaz dissous dans le polymère
est une propriété d’équilibre et est reliée à la pression de la phase gazeuse en contact
par une relation linéaire :
C = S(P ) P (3)
où S(P ) est une quantité d’équilibre appelée solubilité. A haute pression, plu-
sieurs modes d’absorption interviennent ; la solubilité peut alors varier avec la pres-
sion. Mais, à basse pression, conformément à la loi de Henry, la concentration dépend
linéairement de la pression, et la solubilité est une constante indépendante de la pres-
sion. Si S est indépendante de la pression, les relations (2) et (3) donnent le flux de
gaz à travers une membrane polymère comme étant directement lié à la différence
de pression des deux cotés de la membrane :
J = −DS∇P (4)
Le facteur de proportionnalité dans la relation 4 s’appelle perméabilité P e. Il est
définit comme :
Pe = S × D (5)
De manière courante, la concentration de gaz dissous est exprimée par le rap-
port du volume de gaz dans les conditions standards (cm3 (STP)) sur le volume de
polymère cm3 (polymère) à 0˚C sous pression atmosphérique. Elle peut être aussi
exprimée par la proportion massique du gaz dans 100 grammes de polymère. Par
1.4 Intérêt d’une approche par modélisation 5
aléatoirement (régions amorphes) et des régions contenant des chaı̂nes orientées sur
un réseau (régions cristallines).
Dans ce travail de thèse, on s’intéressera à la perméabilité des polymères au gaz
dans leur état semi-cristallin. Un problème majeur introduit par l’étude de polymère
semi-cristallin est sa nature extrêmement hétérogène, aussi bien d’un point de vue
structural que d’un point de vue dynamique.
La figure 4 est une image d’un polypropylène semi-cristallin prise par un mi-
croscope optique. On observe à l’échelle microscopique des sphérolites, constituées
de lamelles cristallines de quelques dizaines de nanomètres d’épaisseur, séparées de
régions amorphes de dimensions variables. La structure sphérolitique n’est pas la
seule structure des polymères semi-cristallins. Suivant les types de polymères et les
conditions de préparation, il peut y avoir d’autres structures plus complexes. Les
régions amorphes et les régions cristallines s’alternent pour former des agrégats de
géométries différents. Cependant, les sphérolites sont les formes les plus souvent
rencontrées. Ces hétérogénéités structurales d’un polymère semi-cristallin sur des
grandes échelles d’espace ne sont pas accessibles par la simulation moléculaire dans
un temps raisonnable.
L’enjeu de cette thèse est d’obtenir la perméabilité aux gaz d’un polymère semi-
cristallin sans passer pas la modélisation de toutes les régions. Les paramètres de
transport du gaz seront uniquement calculés dans la phase amorphe. La présence
1.6 Plan du manuscrit 7
des régions cristallines sera prise en compte implicitement selon des hypothèses
adéquates, c’est ce qui constitue l’originalité de ce travail.
peuvent pas être éliminées au cours de la cristallisation. Ces structures sont rejetées
par les cristallites en cours de croissance et forment une proportion importante de
régions non-cristallines.
expliquée par le fait que les chaı̂nes gaussiennes, pendant la période de cristalli-
sation, ont une certaine probabilité de se replier et de rentrer dans la surface du
cristal [13–15]. Cette probabilité est une fonction décroissante de la rigidité de la
chaı̂ne. Ce modèle montre que l’épaisseur des lamelles augmente avec la rigidité des
chaı̂nes. L’augmentation de Lc avec la rigidité des chaı̂nes a été observée pour des
polymères de natures différentes [12, 16].
Fig. 7 – Comparaison des lamelles ayant des interfaces enchevêtrées (à gauche)
et non-enchevêtrées (à droite) et températures de fusion correspondantes à chaque
structure obtenue par calorimétrie [17]
des lamelles épaisses, l’angle d’inclinaison est de l’ordre de 19-20˚. Cet angle s’accroı̂t
graduellement jusqu’à ∼45˚ avec la baisse de la température de fusion [19, 20].
Les polymères semi-cristallins sont toujours caractérisés par leurs taux de cris-
tallinité. La courbe d’évolution de densité en fonction de température (figure 5) a
montré qu’en général le polymère ne se cristallise pas entièrement. Le taux de cris-
tallinité (Xc) est défini comme la fraction massique ou volumique du système à l’état
cristallin. Dans un polymère ne contenant que les phases amorphe et cristal, le taux
d’amorphe sera 1-Xc. La subtilité est qu’il y a souvent une fraction non négligeable
d’une troisième phase avec une structure intermédiaire entre l’amorphe et cristal.
Cette phase s’appelle «interface». Alors, pour accentuer la nature multi phase des
polymères semi-cristallins, le terme «composition de phase» est plus approprié que
celui de la simple «cristallinité». L’existence de l’interface ainsi que ses propriétés
seront discutées dans les parties suivantes.
La cristallinité diminue avec la masse moléculaire, ce qui est observé pour le cas
des polyéthylènes [21, 22]. L’effet de branchements dans les chaı̂nes est également
étudié sur le taux de cristallisation. Il a été déjà montré que des branchements
comme des groupements méthyls, et ethyls dans une certaine mesure, pouvaient
être incorporés dans la phase cristal. Néanmoins, les branchements contenant plus
de deux carbones sont essentiellement exclus de cette phase [23–26]. La présence de
branchements, en général, va diminuer le taux de cristallisation et l’épaisseur des
lamelles [27].
La structure microscopique d’un polymère semi-cristallin est souvent assez com-
plexe et elle dépend de plusieurs facteurs comme : le mode de cristallisation, l’his-
toire thermique, la présence de branchements, la masse moléculaire, la nature du
polymère, etc. Dans la partie suivante nous introduirons les différentes techniques
expérimentales pour étudier la structure des polymères semi-cristallins.
thermique du cristal est différente selon les différents axes de la maille. Dorset [29]
applique la technique de diffraction d’électrons (ED)7 pour étudier la structure d’un
cristal. Contrairement à la diffusion de rayons X, le facteur de diffusion des électrons
par les atomes d’hydrogène reste pratiquement constant dans tous les angles, ce qui
rend l’ED plus favorable pour détecter les atomes d’hydrogène. Tashiro et al. [30]
ont effectué une expérience de diffraction de neutrons à grand angle (WAND)8 sur
le cristal du polyéthylène dans lequel les atomes d’hydrogène sont remplacés par des
atomes de deutérium. Ils obtiennent les paramètres de maille suivants : a=7.43 Å,
b=4.93 Ået c=2.545 Å. Ils présentent également les coordonnées réduites de chaque
atome dans la maille et attribuent un groupe d’espace «Pnam» au cristal. Leurs
résultats sont en bon accord avec les résultats obtenus par diffraction des électrons.
La figure 8 montre un schéma de la maille élémentaire du cristal de polyéthylène.
suite, l’épaisseur des lamelles est obtenue en multipliant L par le taux de cristallinité
mesuré par d’autres techniques. C’est ainsi que Robelin et al. [12] se servent de la
diffusion de rayons X à petite angle (SAXS) pour mesurer la périodicité à longue
portée du poly(tetrahydrofuran), du poly(éthylène téréphtalate), du polypropylène
et du polyéthylène en fonction de leur masse molaire. Puis en utilisant le taux de
cristallinité de chaque échantillon, ils estiment l’épaisseur des lamelles cristallines et
celle de la partie amorphe. Ils établissent une corrélation linéaire entre ces épaisseurs
et la racine carrée des masses molaires des polymères. Ils montrent également que
la polydispersité augmente la périodicité à longue portée. Pour les polyéthylènes
de masse molaire comprise entre 10000 et 250000 gr/mol, Robelin et al. obtiennent
où Ll est l’épaisseur de la phase amorphe non perturbée.
2.1 Morphologie 17
des périodicités à longue portée (L) qui varient de 200 Å à 450 Å ainsi que des
épaisseurs de lamelles qui varient de 120 Å à 190 Å. Cela montre que le changement
de la périodicité à longue portée (L) dû à la masse molaire provient majoritairement
du changement de l’épaisseur de la phase amorphe. Voigt-Martin et al. [19] utilisent
la microscopie électronique (TEM) pour obtenir la distribution des épaisseurs des
lamelles pour des polyéthylènes monodisperses avec différentes masses molaires. La
figure 9 montre une image du polyéthylène semi-cristallin obtenu par microscopie
électronique. Elle montre une distribution relativement large des épaisseurs des la-
melles. Les épaisseurs obtenues par Voigt-Martin et al. sont en bon accord avec
les résultats de Robelin et al. L’accès à l’épaisseur moyenne de lamelles cristallines
peut se faire également par calorimétrie. Pour cette méthode l’enthalpie de fusion
du cristal et l’énergie libre de l’interface doivent être connues. Il reste finalement à
mesurer le décalage de température de fusion par rapport à la température de fusion
idéale du cristal et ensuite de le relier à l’épaisseur des lamelles par la relation de
Gibbs-Thomson. Hedenqviste et al. [31] utilisent cette technique pour obtenir Lc
dans le polyéthylène. Ils mesurent des épaisseurs de lamelles allant de 80 Å à 200
Å pour différentes masses moléculaires.
La variété des techniques expérimentales disponibles pour étudier la structure
des polymères semi-cristallins montre l’intérêt d’une bonne connaissance de l’hétéro-
généité du polymère à l’échelle microscopique. Le comportement du polymère à cette
échelle a aussi des conséquences importantes sur les propriétés macroscopiques. Les
phases amorphes non perturbées et cristallines sont assez connues. A l’opposée,
l’étude des zones interfaciales est relativement récente. Dans la partie suivante, nous
discuterons de la structure de l’interface ainsi que des différents facteurs qui influent
sur la conformation des chaı̂nes de polymère dans cette région.
Pour conclure, la dissipation d’ordre lors du passage d’une région cristalline à une
région amorphe a lieu sur une certaine longueur. Comme la figure 12 le montre, le
passage d’un cristal à la région amorphe n’est pas brutal et s’effectue dans une région
de transition, l’interface. Les différentes propriétés intensives du polymère, y compris
la densité, changent graduellement de la phase cristalline vers la phase amorphe non
perturbée. En plus de la phase cristalline et de la phase amorphe non perturbée,
2.2 Morphologie et perméabilité 22
la prise en compte de l’interface comme troisième phase dans les polymères semi-
cristallins est essentielle. Dans la suite, l’ensemble de la phase amorphe non perturbée
plus l’interface sera appelé «phase amorphe».
S = αSa (7)
KD = αKD,a (11)
où KD,a est la constante de solubilité de Henry dans le PET 100% amorphe.
Revenons à l’équation (7) qui est basée sur l’hypothèse d’imperméabilité de la
phase cristalline. Cette hypothèse n’est pas toujours valable. Certains polymères ont
des structures cristallines plus ouvertes laissant entrer les petites molécules. Comme
nous l’avons montré dans la partie 2.1.3, il existe des polymères pour lesquels la
densité de la phase cristalline (ρc ) est plus petite que celle de la phase amorphe
2.2 Morphologie et perméabilité 24
D = D∞ exp(γC) (12)
2.2 Morphologie et perméabilité 25
Dans cette relation C est la concentration de gaz, D∞ est la diffusion à dilution infinie
et γ est une constante qui caractérise la plastification. Selon Vittoria et al. [59], γ
augmente très rapidement avec le taux de cristallinité. En revanche, D∞ diminue
avec la cristallinité.
Le processus de diffusion dans les polymères semi-cristallins a été étudié par
plusieurs modèles. Il a été montré que la diffusion du gaz s’effectue par des sauts
successifs dans la matrice polymère. Il existe des modèles théoriques pour calculer la
diffusion en fonction de la longueur moyenne de l’amplitude des sauts [71–73]. Dans
le cas de la diffusion de l’éthylène dans un polyéthylène semi-cristallin ayant une
cristallinité de 50% ces modèles conduisent à une longueur de saut d’environ 12 nm.
Cette valeur diminue jusqu’à 5,5 nm pour une cristallinité de 70% [74]. L’effet de la
cristallinité est donc pris en compte par une diminution de la longueur moyenne de
l’amplitude de sauts. Ce modèle n’a pourtant aucun aspect prédictif car l’amplitude
moyenne des sauts n’est pas reliée d’une manière claire à la cristallinité ou à la
morphologie du polymère.
Comme on l’a vu dans la partie précédente, les cristallites sont souvent considérées
comme des barrières imperméables. Les particules de diffusant ne peuvent pas pénét-
rer dans les cristallites. Sur la figure 13 on compare la trajectoire du diffusant dans
le polymère purement amorphe (couleur rouge) avec sa trajectoire dans le polymère
semi-cristallin (couleur noir). L’imperméabilité des cristallites impose un chemin
plus long aux diffusants qui doivent éviter les zones cristallines. Cette allongement
du chemin conduit à la réduction du coefficient de diffusion effectif. Cette réduction
peut être représentée par un facteur géométrique. Klute [70] a introduit le facteur
2.2 Morphologie et perméabilité 26
où B est un facteur qui dépend du diamètre du gaz. Il s’appelle «facteur bloquant».
Ce facteur est associé à la distribution des régions amorphes entre les régions cris-
tallines. Les zones amorphes constituent les zones de passage pour le diffusant. Ces
zones vont laisser passer ou non les molécules qui diffusent selon la dimension des
ces molécules par rapport à la dimension de ces zones amorphes. La figure 14 illustre
un schéma de cette situation. En conséquence, le facteur bloquant augmente avec la
taille des molécules.
Malgré les interprétations différentes, les deux formulations (13) et (14) sont
équivalentes si :
1
τ= et β=B (15)
ψ
Étant donné une littérature riche sur la notion de tortuosité et de facteur d’immo-
bilisation, nous tenterons d’approfondir ces notions dans les paragraphes suivants.
Le facteur β a été identifié par Myers [78] mais il a été mesuré pour la première
fois par Michaels [54]. Nous décrirons ci-dessous les techniques expérimentales pour
mesurer les différents facteurs de l’équation (13).
Michaels et Parker [54] ont extrapolé la perméabilité de la phase amorphe à
partir de l’état fondu par une loi d’Arrhenius (cf. annexe C). Ensuite, ils ont divisé
la perméabilité extrapolée par la solubilité dans la phase amorphe pour obtenir Da .
L’autre méthode est celle de la mesure de la diffusion dans des polymères analogues
100% amorphes. Prenons comme exemple le cas de la diffusion dans le caoutchouc
qui est équivalente à la diffusion dans le polyéthylène amorphe [69]. Mogri [68] a
2.2 Morphologie et perméabilité 28
βHe = 1 (16)
Michaels et Parker [54] ont utilisé les équations suivantes pour estimer τ et β :
P eHe /P e
β= caoutchouc
(17)
P eHe /P ecaoutchouc
Da
τ= (18)
Dβ
où l’exposant caoutchouc désigne les grandeurs associées dans le caoutchouc. Dans ces
équations il est supposé que la solubilité des gaz dans le caoutchouc est égale à la so-
lubilité dans la phase amorphe du polyéthylène. Pourtant, cette dernière hypothèse
n’est pas justifiée par l’expérience et elle semble erronée. Quelques années plus tard,
Michaels [69] a amélioré sa méthode s’affranchissant ainsi de cette hypothèse incor-
recte. Il a calculé les nouveaux τ et β par :
caoutchouc
DHe
τ= (19)
DHe
Dcaoutchouc
β= (20)
Dτ
Ces équations peuvent être démontrées à partir des équations (13) et (16). Cette
méthode est généralement plus exacte et elle donne des tortuosités plus grandes
comparativement à l’ancienne méthode.
Après avoir déterminé les différents facteurs de la relation (13), il reste à les relier
à la morphologie du cristal. Il existe des modèles macroscopiques et rigoureux pour
calculer la tortuosité dans les matériaux poreux en fonction de la morphologie du
système [80–82]. Néanmoins, ces modèles ne marchent pas bien pour les polymères
semi-cristallins. Les tortuosités τ obtenues dans les différents types de polymères
n’ont pas montré une corrélation simple avec le taux de cristallinité. L’analogie avec
2.2 Morphologie et perméabilité 29
Fig. 15 – Illustration de l’effet des régions cristallines (en noir) sur le gonflement
du polymère semi-cristallin. D’après de Pablo [91]
Ces observations montrent que la solubilité dans la phase amorphe des po-
lymères semi-cristallins ne peut généralement pas être traitée comme la solubilité
dans la phase amorphe libre. Une description plus réaliste de la solubilité dans le
polymère semi-cristallin est alors nécessaire. Cette description a été réalisée soit par
des modèles élastiques soit par des modèles basés sur la présence de l’interface. Dans
les paragraphes suivants, nous présenterons d’abord les modèles élastiques. Ensuite,
nous introduirons le modèle basé sur les trois phases.
Le modèle élastique consiste à attribuer une élasticité au réseau des chaı̂nes
pontantes entre les cristallites et donc à ajouter une contribution énergétique au
gonflement du polymère. Dans ce modèle, l’enthalpie libre de mélange ∆Gm s’écrit
sous la forme suivante :
∆Gm = ∆Gfmree + ∆Gel m (24)
dans laquelle l’exposant f ree signifie dans la phase amorphe libre. Le terme ∆Gel m est
le coût énergétique supplémentaire dû au gonflement d’un réseau élastique. Ainsi, le
potentiel chimique et l’activité du soluté sont donnés par la définition :
a1 (C1 ) = aGaz
1 (P, T ) (26)
2.2 Morphologie et perméabilité 32
Ici, aGaz
1 (P, T ) est l’activité du soluté dans son état gazeux qui est une fonction
de la température T et la pression P . Selon la relation (25), le potentiel chimique
du gaz peut être écrit comme la somme d’une contribution dans l’état librement
amorphe et d’une contribution élastique. Cette relation associée à la relation (26)
montre que l’augmentation d’activité par le terme élastique diminue la concentration
d’équilibre de gaz dans le polymère.
L’activité du soluté dans un polymère librement amorphe ln(af1 ree ) peut être cal-
culée par différents modèles thermodynamique [95]. Pour le terme élastique ln(ael 1 ),
il n’existe, à ce jour, que deux modèles : celui de Flory et Rehner(F-R) [96, 97] et
celui de Michaels et Hausslein (M-H) [98].
Le deuxième modèle élastique est le modèle de M-H [98]. Dans ce modèle l’élasti-
cité du réseau possède une origine enthalpique. Les chaı̂nes reliant les cristallites dans
ce modèle sont considérées comme des ressorts «Hookiéens». La force de traction
sur ces chaı̂nes est en équilibre avec la force induite par l’énergie libre de cristallisa-
tion. Cette force dépend de l’écart de température par rapport à la température de
fusion. Elle est plus importante aux basses températures (cf. figure 17). Michaels et
2.2 Morphologie et perméabilité 34
où ∆H2f est la chaleur de fusion par unité de masse du cristal, φ1 est la fraction
volumique du gaz dans la phase amorphe, φ2 est la fraction volumique du polymère
amorphe (= 1 − φ1 ) et χ est le paramètre d’interaction de Flory-Huggins. Il reste
f qui est le paramètre ajustable du modèle et représente la fraction de chaı̂nes
élastiques entres les cristallites. Le paramètre f dépend uniquement de la morpho-
logie, de la cristallinité et de l’histoire thermique du polymère. Cette expression
montre que le terme élastique tend vers zéro quand la température s’approche de
la température de fusion. L’ajustement du modèle d’élasticité M-H sur la solubilité
de para xylène dans le polyéthylène a conduit à une fraction de chaı̂nes élastiques
f comprise entre 0,25 et 0,35 en fonction de la vitesse de refroidissement du po-
lymère [98]. Pour les polyéthylènes cristallisés par l’action simultanée de l’étirement
et du refroidissement, la cristallinité est très élevée et les chaı̂nes de la phase amorphe
sont très orientées. Dans ce type de polyéthylène le paramètre f est de l’ordre de
1 et le terme élastique devient plus important [98]. Banaszak et al. [99] ont montré
que la reproduction des solubilités d’éthylène et de 1-hexene dans le polyéthylène
en dessous de la température de fusion par la méthode PC-SAFT [100, 100, 101]
demande l’ajout du terme élastique de M-H. Serna et al. [102] ont calculé la solu-
bilité de différents gaz dans le polyéthylène avec l’équation d’état de GCLF [89] en
ayant recours à la méthode UNIFAC-VdW-FV [103]. Ils ont montré que la contribu-
tion élastique de M-H est obligatoire pour reproduire les solubilités expérimentales.
Yao et al. [104] ont étudié également les solubilités des gaz dans le polyéthylène
pour des températures allant de 60 à 90 ˚C. Étant donné que la cristallinité du
polyéthylène dans cette région change en fonction de la température, la fraction des
chaı̂nes élastiques ne sera plus constante. Ils ont amélioré l’expression élastique de
M-H en tenant compte de la dépendance de f en température.
Les modèles élastiques de M-H [98] et F-R [96, 97] sont des modèles concur-
rents. Certains auteurs ont comparé ces deux modèles élastiques. Nous citons par
exemple les travaux de Doong et al. [95] où les solubilités de nombreux gaz dans
le polyéthylène ont été ajustées par les deux modèles élastiques. Ils ont pu montrer
que la dépendance en température de la solubilité dans le polyéthylène est mieux
reproduite par le terme d’élasticité de M-H. En revanche, dans le polypropylène,
c’est le modèle F-R qui représente mieux la solubilité [105]. Il en résulte que le choix
de l’un ou l’autre de ces deux modèles dépend du polymère étudié.
Outre les modèles élastiques, l’autre manière d’interpréter la différence entre
la solubilité dans l’amorphe libre S f ree et la solubilité expérimentale corrigée Sa
est la présence de l’interface. Hedenqvist et al. [31] ont évalué les taux d’inter-
face pour différents polyéthylènes par spectroscopie Raman. Ils ont ensuite mesuré
2.3 Conséquence sur la modélisation 35
les solubilités d’oxygène et d’hexane. Ils montrent que la solubilité corrigée par le
taux d’amorphe (contenant la phase amorphe non perturbée et l’interface) est en-
viron deux fois plus faible que la solubilité corrigée par le taux d’amorphe non
perturbée. En d’autres termes, la solubilité dans la phase amorphe du polyéthylène
semi-cristallin (Sa ) est plus petite que la solubilité dans une phase librement amorphe
(S f ree ). A l’inverse des modèles élastiques, la réduction de solubilité est cette fois-ci
indépendante du taux de gonflement.
Jusqu’à présent l’impact de l’environnement cristallin a été étudié uniquement
sur la solubilité dans la phase amorphe. Cet environnement peut éventuellement
influencer la diffusion dans la phase amorphe. Différents calculs de coefficients de
diffusion à basse température par simulation moléculaire [106–108] n’ont pas pris
en compte cette influence. L’étude de la diffusion des gaz dans le polyéthylène
semi-cristallin par le modèle de Cohen-Turnbull-Fujita (CTF) [109] a montré que
la fraction de volume libre 13 de la phase amorphe du polyéthylène semi-cristallin
est environ la moitié de celle obtenue dans le caoutchouc [31]. Cela signifie que la
diffusion dans le caoutchouc qui est l’équivalent du polyéthylène 100% amorphe est
supérieure à celle dans la phase amorphe du polyéthylène in situ. Hedenqvist et
al. [31] ont estimé la fraction de volume libre de l’interface comme étant entre 20
et 50% de celle de liquide suivant le taux de branchement ou de la morphologie du
polyéthylène, ce qui conduit à une diffusion plus lente dans la phase amorphe du
polyéthylène semi-cristallin par rapport à l’amorphe libre.
Pour conclure, on a montré dans cette partie que la présence des cristallites
peut changer les différentes propriétés de l’amorphe. L’impact de l’environnement
cristallin sur la solubilité et la diffusion des gaz dans l’amorphe a été étudié dans
différents travaux. Cette influence doit être aussi prise en compte dans la simulation
moléculaire de la phase amorphe, ce qui sera évoqué dans la partie suivante.
La contrainte ainsi exercée sur la phase amorphe, peut être calculée théoriquement
dans le cadre des théories élastiques de F-R ou M-H [99]. Ces contraintes seront di-
rectement liées au taux de gonflement du polymère induit par l’absorption de gaz
et elles dépendront de la nature et de la pression du gaz. En revanche, l’effet de
l’interface sur la solubilité dans la phase amorphe est indépendant de gonflement.
Dans ce travail σ est un paramètre ajustable qui à priori contient les deux effets.
C’est l’interprétation des valeurs obtenues de σ qui clarifie son origine à posteriori.
La méthodologie et la technique permettant d’exercer une telle contrainte supplé-
mentaire au cours de la simulation moléculaire seront présentées dans le chapitre
suivant.
37
3 Méthodes retenues
3.1 Mécanique statistique et modélisation moléculaire
La modélisation moléculaire permet d’étudier des systèmes complexes pouvant
contenir plusieurs milliers de particules. Tout en s’appuyant sur l’expérience pour
étalonner ses modèles, elle apporte des informations sur des grandeurs parfois non
accessibles expérimentalement. L’expérience permet de mesurer des grandeurs ma-
croscopiques à des échelles de temps et d’espace macroscopiques. Comme il a été
dit en introduction, la modélisation moléculaire se place à l’échelle nanométrique,
en temps et en espace. La thermodynamique statistique permet de faire le lien
entre modélisation moléculaire et expérience. Il s’agit du fondement des méthodes
de modélisation que sont la méthode de Monte Carlo et la Dynamique Moléculaire.
La thermodynamique statistique décrit l’état d’équilibre thermodynamique com-
me un ensemble de micro-états dont la contribution est fonction de leur poids de
Boltzmann. L’objectif de la modélisation moléculaire est de construire une col-
lection de ces micro-états, représentative de la thermodynamique macroscopique.
L’ensemble des configurations accessibles à un système est appelé espace des confi-
gurations. Lorsque les vitesses des particules sont prises en considération, on parle
d’espace des phases.
Soit une observable thermodynamique O qui dépend de la configuration Λ dans
laquelle le système se trouve. Puisque la configuration des particules évolue au
cours du temps t, la grandeur O(Λ) évolue aussi. Ainsi, la grandeur mesurable
expérimentalement est la moyenne de O(Λ) dans le temps :
Z Tobs
1
hOitemps = lim O(Λ(t))dt (29)
Tobs →∞ Tobs 0
où Tobs est le temps d’observation qui est souvent très grand par rapport au temps de
relaxation du système à l’échelle moléculaire 14 . Dans l’approximation classique, les
équations de Newton régissent la dynamique des particules et permettent de prédire
leur évolution temporelle.
La thermodynamique statistique montre que si on laisse le système évoluer
indéfiniment, il passera par une multitude de configurations qui contribuent à la
grandeur moyenne hOi. Dès lors, on peut s’affranchir de la notion de temps et en-
visager d’évaluer hOi par une moyenne d’ensemble. Ainsi, avec ρ(Λ) la densité de
probabilité de la configuration Λ, on peut écrire :
X
hOiens = O(Λ)ρ(Λ) (30)
Λ
e−βu(Λ)
ρ(Λ) = (31)
Q
dans laquelle u(Λ) est l’énergie potentielle associée à la configuration Λ, β est
1/kB T et Q est la fonction de partition définie comme :
X
Q= e−βu(Λ) (32)
Λ
Dans l’algorithme de Metropolis original (cf. référence [110]), α est une matrice
symétrique telle que α(a → n) = α(n → a). Dans le cas où α est symétrique,
l’équation (34) s’écrit :
figure 23) est effacé, puis il est reconstruit entre les atomes i et j pour former une
3.2 Méthode Monte Carlo 42
nouvelle chaı̂ne (chaı̂ne foncée, figure 24). Il en est de même pour le trimère compris
entre les atomes j2 et j qui est reconstruit entre j2 et i2 pour former une nouvelle
chaı̂ne (chaı̂ne claire, figure 24).
La conformation de ces deux chaı̂nes est complètement modifiée par ce mouve-
ment. L’échantillonnage des configurations issu de l’utilisation de ce mouvement est
de meilleure qualité et permet d’atteindre une meilleure précision dans l’évaluation
des grandeurs étudiées. Il est important de noter que ce mouvement peut engendrer
une certaine polydispersité dans les longueurs de chaı̂ne. Toutefois, la polydisper-
sité est totalement contrôlable, en limitant la section de la chaı̂ne dans laquelle les
atomes i et j sont choisis. Au cours de ce travail, les chaı̂nes sont considérées comme
strictement monodisperses.
IDR est la variante intra-moléculaire de DB. Les deux trimères sont reconstruits
au sein de la même chaı̂ne. Les figures 25 et 26 illustrent ce mouvement.
pour chaque ensemble statistique une fonction de partition Q et une densité de pro-
babilité ρ(Λ). Trois ensembles statistiques sont présentés ci-après : l’ensemble de
Gibbs, l’ensemble Grand Canonique et l’ensemble osmotique. Ils ont été largement
utilisés pour étudier le phénomène d’absorption.
L’ensemble osmotique
Dans le cadre de la simulation d’un tel système, seule la phase contenant tous les
n composés est considérée (boı̂te de droite de la figure 27). Les potentiels chimiques
µi=1..p des p composés sont imposés. Etant donné qu’ils ne sont pas connus a priori,
une simulation NPT préalable de la phase vapeur permet de les calculer comme dans
l’ensemble Grand Canonique. Par ailleurs, la boı̂te de simulation est soumise à une
contrainte isotrope σ, ce qui permet d’observer des changements de volume.
Pour chacun des n − p constituants ne pouvant pas traverser la paroi poreuse,
le nombre de molécules Nj=p+1..n est constant i.e. dNp+1 = · · · = dNn = 0. Puisque
l’ensemble osmotique comprend des paramètres extensifs fixés, l’ensemble ne peut
3.3 Dynamique moléculaire 45
être dégénéré. Un schéma des ensembles Grand Canonique et osmotique est représen-
té sur la figure 27.
où ∆U est le changement d’énergie totale du système suite à l’ajout d’une particule
(= U (N + 1) − U (N )).
Le permier terme s’appelle potentiel chimique idéal (µid i ). Ce terme est calculable
analytiquement. Par contre, le second terme, appelé potentiel chimique d’excès (µex i ),
s’obtient par le test d’insertion.
En pratique, il est plus simple de manipuler la fugacité plutôt que le potentiel
chimique. La fugacité fi d’une espèce i est assimilable à une pression. Elle se déduit
par l’expression suivante :
fi = P0i eµi /kB T (40)
La fugacité peut être également exprimée en fonction du potentiel chimique d’excès :
Ni kB T µex N kB T µex
fi = e i /kB T = xi e i /kB T (41)
V V
où xi est la fraction molaire de l’espèce i dans le système. Il faut noter que pour un
gaz parfait, le potentiel chimique d’excès µexi vaut zero et la fugacité de l’espèce i se
réduit à sa pression partielle.
où mi est la masse de la particule i et U est le potentiel de pair. Avec cette définition,
l’équation 42 prend la forme suivante :
mi r¨i = Fi (46)
3.3 Dynamique moléculaire 47
Donc, à chaque pas, l’énergie potentielle et la force peuvent être évaluées pour chaque
particule. Il ne reste plus qu’à intégrer les équations du mouvement pour déterminer
le déplacement de la particule dans le temps ∆t.
mhvα2 i = kB T (48)
Dans ces équations, Fi est la force exercée sur la particule i. La forme générale de
cette force est la suivante :
où −∇ri U est la partie conservative de la force. Ficontrainte est la force exercée sur
l’atome i pour le maintenir à une distance constante de ses deux plus proches voisins
à l’intérieur de la chaı̂ne de polymère. La version RATTLE de SHAKE [125] a été
3.3 Dynamique moléculaire 48
utilisée pour prendre en compte ces contraintes. Enfin, L est le nombre de degrés de
liberté du système. La quantité d’énergie conservée du système s’écrit :
N
X p2i p2ξ
HN H = + U(rN ) + + LkB T ξ (54)
i=1
2mi 2Q
Il faut noter que HN H n’est pas un vrai Hamiltonien puisqu’il n’est pas possible d’ob-
tenir les équations du mouvement par les dérivés partiels de HN H . La constante Q
est intitulée «constante de couplage». Elle indique l’inertie du thermostat. Les fluc-
tuations de la température sont en relation inverse avec la valeur de cette constante :
plus sa valeur est faible, plus les fluctuations de la température sont rapides. À l’in-
verse, une haute valeur de Q rend les fluctuations de température très lentes dans
une période longue. Ainsi, la bonne valeur de Q est celle qui permet de recréer les
fluctuations du système dans l’ensemble microcanonique. On peut également relier
Q à un temps T caractéristique de la relaxation du système par :
Q = LkB T T 2 . (55)
iL = Λ̇ .∇Λ (57)
3.3 Dynamique moléculaire 49
où Λ={ri ,pi ,ξ,pξ }. Sur le plan analytique, l’application de l’opérateur de Liouville
n’est possible que sur un intervalle de temps δt très petit. Si P est le nombre de
pas du calcul et t le temps total simulé on aura δt = t/P . On peut appliquer P fois
l’opérateur de Liouville sur l’intervalle de temps δt comme :
P
Y
Λ(t) = exp (iLδt)Λ(0) (58)
i=1
Étant donné que les opérations iLr , iLv et iLN H ne commutent pas entre elles, le
développement direct de l’exponentielle n’est pas possible. Pourtant, pour δt très
petit, l’application de la formule de Trotter indique que :
δt δt
exp (Aδt + Bδt) = exp (A ) exp (Bδt) exp (A ) (63)
2 2
Par conséquent, pour l’opérateur de Liouville on aura :
δt δt δt δt
e(iLδt) = e(iLN H 2 ) e(iLv 2 ) e(iLr δt) e(iLv 2 ) e(iLN H 2 ) (64)
L’opérateur iLN H peut à son tour être factorisé par le développement de Trotter :
δt δt δt δt δt
e(iLN H 2 ) = e(iL2 4 ) e(iL1 2 ) e(iL3 2 ) e(iL2 4 ) (65)
dans lequel :
N
X
iL1 = − ξ˙r˙i .∇r˙i (66)
i=1
N
1 X p2i
iL2 = − LkB T .∇ξ̇ (67)
Q i=1
mi
En dépit de leur apparence compliquée, ces formules sont finalement assez simples
à implémenter dans un code. De même, elles sont peu coûteuses en temps de calcul.
Ainsi, établir un algorithme d’intégration s’effectue facilement. L’opérateur iL est
appliqué en différentes étapes. Les règles de transformations sont les suivantes :
δt δt
e 4 iL2 : ξ˙ → ξ˙ + ξ¨ (69)
4
δt δt ˙
e 2 iL1 : r˙i → r˙i + exp (− ξ) (70)
2
δt δt ˙
e 2 iL3 : ξ→ξ+ ξ (71)
2
δt δt Fi
e 2 iLv : r˙i → r˙i + (72)
2 mi
eδtiLr : ri → ri + δtr˙i (73)
L’intégration des équations du mouvement demande d’appliquer les opérateurs δt4 iL2 ,
δt
2
iL1 , δt2 iL3 et δt2 iLv selon l’ordre indiqué dans les équations (64) et (65). Ensuite,
l’application de l’opérateur de δtiLr donne les nouvelles positions des particules.
Dans cette étape, on calcule de nouvelles forces sur chaque particule et on continue
à appliquer les opérateurs restants pour obtenir l’état du système {ri ,pi , ξi et pξ } à
l’instant t + δt.
Dispersion-répulsion
Les termes de dispersion et de répulsion sont regroupés au sein d’un même po-
tentiel. Le potentiel de Lennard-Jones 6-12 est le plus utilisé. Selon ce potentiel,
pour deux particules i et j en interaction à la distance rij , l’expression du potentiel
VijLJ est :
" #
σ ij σij
VijLJ = 4ij ( )12 − ( )6 (74)
rij rij
où ij est la profondeur du puits et σij représente la distance à laquelle l’énergie
d’interaction est nulle.
La partie dispersive correspond ici au terme en 1/r6 alors que la partie répulsive
est représentée par le terme en 1/r12 . De plus, l’adoption de la puissance 12 pour
la partie répulsive s’adapte très bien au calcul numérique puisqu’il suffit d’élever le
terme en 1/r6 au carré pour l’obtenir, entraı̂nant ainsi un gain de temps considérable.
Deux constantes caractérisent le potentiel de dispersion-répulsion. Elles dépendent
du couple de centres de forces (i et j) mis en jeu : ij et σij . Pour des raisons de com-
modité et d’universalité des champs de force, les paramètres sont souvent déterminés
pour chaque centre de force i pris séparément. Dans ce cas, les paramètres d’interac-
tion entre deux centres de force de différents types sont évalués à travers l’applica-
tion de règles de combinaison. Signalons que ces règles, ou formules, sont empiriques.
Même s’il en existe un grand nombre, on ne recourt généralement qu’à quelques-unes
d’entre elles : Lorentz-Betrhelot, Kong [128] et les moyennes géométriques.
La règle de Lorentz-Berthelot est la suivante :
√
ij = i j (75)
1
σij = (σi + σj ) (76)
2
Celle de Kong [128] est :
12 13
ii σii12 jj σjj
1
ij σij12 = 13 1 + ( ) 13 (77)
2 ii σii12
3.4 Modèles moléculaires 52
1
ij σij6 = (ii σii6 jj σjj
6 2
) (78)
Nous avons opté dans ce travail pour les règles de mélange de Kong. Ce choix tient
du fait qu’elles offrent de meilleurs résultats que celles de Lorentz-Bertelot pour les
simulations d’équilibre liquide-vapeur de mélanges [129] et notamment pour CO2
dans les alcanes [130].
L’évaluation des interactions entre particules est l’opération la plus coûteuse.
Alors, afin de réduire le nombre d’interactions à calculer il est possible d’introduire
un rayon de coupure, noté rc , dans les simulations. En effet, on considère le potentiel
nul pour une distance supérieure à rc . Toutefois, cette méthode entraı̂ne une erreur
systématique dans l’évaluation de l’énergie potentielle du système. Il est possible de
supposer que, du point de vue d’une particule, le milieu est homogène au-delà de rc .
On a alors :
Nρ ∞
X Z
U= U(rij < rc ) + U(r)4πr2 dr (79)
j>i
2 rc
on pose :
ρi,1 = qi δi (r) + gi (r) (84)
ρi,2 = −gi (r) (85)
La densité de charge liée à une charge ponctuelle qi s’exprime comme la somme de
ces deux termes. La contribution de l’énergie électrostatique totale provenant du
premier terme ρi,1 converge rapidement dans l’espace directe, et nous pouvons donc
choisir d’appliquer un rayon de coupure au delà duquel la contribution est jugée
négligeable. La contribution à l’énergie électrostatique totale provenant du second
terme ρi,2 ne converge pas facilement dans l’espace direct, mais il est possible de
l’expliciter sous la forme d’une série de Fourier rapidement convergente dans l’espace
réciproque.
3.4 Modèles moléculaires 54
~ r) = −∇φ(~r) = −∇(φ1 + φ2 )
E(~ (86)
N
X qi
φ1 (~r) = [1 − erf (α | ~r − ~ri |)] (87)
i=1
4π0 | ~r − ~ri |
N
1 X X qj 2 2
φ2 (~r) = 2
exp [i~k.(~r − ~rj )e−k /4α ] (88)
0 V j=1
k
~k6=0
3.5 Conclusion
Le calcul de la perméabilité des polymères aux gaz par simulation moléculaire
peut être effectué en deux étapes. Dans la première étape, la méthode de Monte
Carlo permet de calculer la solubilité. Afin de contrôler précisément la composition
des mélanges de gaz en équilibre avec la phase polymère, ainsi que pour imposer une
contrainte isotrope sur la phase amorphe du polymère indépendante de la pression
de gaz, l’ensemble osmotique a été utilisé. Toutes les simulations de Monte Carlo
dans ce travail ont été effectuées à l’aide du code Gibbs. Ce code a été développé
conjointement par le Laboratoire de Chimie Physique de l’Université de Paris-Sud
et par IFP Energies nouvelles. Ce code a été modifié au cours de la thèse pour
calculer l’énergie électrostatique intra moléculaire avec la méthode de la sommation
d’Ewald et effectuer le mouvement de «stretching». Ces modifications initialement
programmés pour simuler un polymère fluoré, ne sont pas utilisées pour l’étude de
l’absorption de gaz dans le polyéthylène.
Dans une deuxième étape, la méthode de dynamique moléculaire permet d’accéder
au coefficient de diffusion de gaz dans le polymère. La densité du polymère ainsi que
la concentration de gaz à l’équilibre sont celles obtenues dans la simulation Monte
Carlo. Toutes les simulations de dynamique moléculaire ont été effectuées à l’aide
du code Newton. Ce code est développé dans le Laboratoire de Chimie Physique de
l’Université Paris-Sud.
Finalement la perméabilité est le produit de la solubilité et du coefficient de
diffusion. Les méthodes introduites dans ce chapitre ont été employées pour obtenir
la solubilité et la diffusion de gaz dans le polyéthylène. Les résultats obtenus sont
présentés dans les chapitres suivants.
3.5 Conclusion 58
59
4 De l’amorphe au semi-cristallin
Afin de calculer la concentration de gaz dissous dans le polymère amorphe à
l’équilibre, et donc la solubilité, nous avons opté pour la méthode de Monte Carlo
dans l’ensemble osmotique. Cet ensemble permet de fixer la température T , le
nombre de chaı̂nes de polymères Np , le potentiel chimique de chaque espèce de
gaz µgaz et la contrainte isotrope exercée sur la phase polymère σ. Cette contrainte
est égale à la pression de gaz avec lequel le polymère est en contact. La valeur
ainsi imposée pour la contrainte conduit à une solubilité assez proche de la valeur
expérimentale pour le polyéthylène à l’état fondu. Cet accord a été vérifié dans le
cadre de la thèse de F. Faure [94]. On a représenté sur la figure 29 la concentration
de CO2 et de CH4 dans le polyéthylène fondu, à 433 K, en fonction de la pression
de gaz. Les résultats de simulation dans l’ensemble osmotique sont comparés aux
résultats expérimentaux de Sato [141] pour CO2 et de Lundberg [142] pour CH4 . Les
droites représentent les valeurs obtenues par un calcul à dilution infinie, dans l’en-
semble NPT. Comme on peut le voir sur cette figure, un accord raisonnable avec les
résultats de la littérature est obtenu, ce qui valide à la fois l’aspect méthodologique
et le choix des potentiels [94].
Des difficultés apparaissent quand on veut reproduire la solubilité à basse tempé-
rature. La figure 30 montre les résultats de simulation pour les solubilités du méthane
et du dioxyde de carbone obtenues dans l’ensemble osmotique. Cela montre une
surestimation d’un facteur 5 de la solubilité du CO2 et d’un facteur 4 de la solubilité
du méthane comparativement aux données expérimentales.
La surestimation de solubilité par rapport à l’expérience à basse température
est un problème qui ne saurait être lié à la qualité du champ de force : les mêmes
potentiels ont permis d’obtenir des résultats très satisfaisants dans le polyéthylène
fondu [94]. L’origine de cette différence peut être attribuée à la présence de régions
cristallines à basse température. En d’autres termes, la solubilité dans la phase libre-
ment amorphe et dans la phase amorphe du polymère semi-cristallin sont différentes
(cf. chapitre 2 ; relation (23)). Pour résoudre ce problème, on propose d’exercer une
contrainte supplémentaire sur la phase amorphe. Par cette contrainte supplémentaire,
on souhaite prendre en compte implicitement l’influence des régions cristallines sur
la phase amorphe. La contrainte isotrope appliquée s’écrit comme :
σ = P + σP E (93)
où P est la pression du gaz et σP E est la contrainte supplémentaire. σP E est un
paramètre ad-hoc. Une contrainte σP E nulle correspond à la simulation d’une phase
amorphe libre.
En général, la contrainte supplémentaire change la solubilité. Afin d’observer
une modification significative de la solubilité, il faut appliquer une contrainte σP E
bien plus importante que la pression de gaz. Le rôle de la contrainte supplémentaire
σP E est de diminuer le volume accessible pour les gaz. Cela conduit à l’expulsion
de certaines molécules de la matrice polymère et la diminution de la solubilité (cf.
60
15
H2S osmotique
H2S/nC100 D∞
12 CO2 osmotique
CO2/nC100 D∞
0
0 1 2 3 4 5 6
Pression (MPa)
25
CO2 exp Ash
CO2 osmotique
CH4 osmotique
20 CH4 exp Naito
C(gr/100gr)
15
x5
10
x4
0
0 2 4 6 8
P(MPa)
figure 31). Par conséquent, plus σP E augment, plus la solubilité diminue. La valeur
de contrainte σP E peut être ajustée afin d’obtenir la solubilité expérimentale. Cette
valeur ajustée est appelée dans la suite la contrainte effective (σef f ).
La valeur de contrainte effective signifie l’ampleur de l’influence des régions cris-
tallines sur la phase amorphe. Cette influence dépend de la morphologie du polymère
semi-cristallin. Tous les facteurs qui sont censés affecter la morphologie du polymère
peuvent changer la valeur de σef f . Ettant donné que la morphologie du polymère
dépend de la température, on s’attend à une contrainte effective dépendant de la
température. L’effet de température sur la valeur de la contrainte effective sera étudié
à la fin de ce présent chapitre.
Selon la définition donnée pour la contrainte effective, sa valeur peut changer en
fonction des données expérimentales choisies. L’écart entre les différentes données
expérimentales peut induire une incertitude importante sur la valeur de contrainte
effective. Plus on dispose de points expérimentaux, plus la valeur estimée de σef f
sera précise. Par contre, s’il n’y a pas assez de données expérimentales fiables, la
valeur de la contrainte effective devient très incertaine. C’est pour cette raison qu’il
est nécessaire de faire un point sur les différentes données expérimentales de so-
lubilité dans la littérature avant d’estimer la valeur de contrainte effective. Ainsi,
on commence ce chapitre par la présentation de différentes mesures expérimentales
de solubilité. A partir de ces données, une solubilité moyenne a été calculée. Cette
moyenne sera utilisée pour obtenir la contrainte effective. Afin de prendre en compte
la dispersion entre les résultats des différents auteurs, on calcule la déviation stan-
dard de la valeur moyenne qui représente l’erreur sur la solubilité expérimentale.
4.1 Données expérimentales de solubilité 62
Sa Sa Sa Sa
CO2 (293K) CH4 (298K) N2 (298K) H2 (313K)
Li [146] - 0,233 0,260 -
Michaels [33] 0,932 0,153 0,054 -
Naito [144] 0,966 0,143 0,052 -
Ash [143] 1,180 - 0,061 0,0035
IFP [148, 150] 0,911 - - 0,0026
Togawa [147] 1,379 - - -
Solms [145] 0,772 0,096 - -
< Sa > 1,023 0,131 0,056 0,0031
δSa 0,081 0,014 0,002 0,0004
gr(gaz)
Tab. 4 – Solubilités expérimentales (en 100gr(polymère) MPa
) des gaz dans le
polyéthylène semi-cristallin. Ces valeurs ont toutes été obtenues à partir de la pente
de la courbe de concentration du gaz en fonction de la pression. Toutes les solu-
bilités sont corrigées par la fraction volumique de phase amorphe du polymère et
représentent la solubilité de la phase amorphe Sa . Les valeurs représentées par la
couleur rouge sont très loin des autres données expérimentales et elles ne sont pas
prises en compte dans le calcul de solubilité moyenne.
14 P = 0.8 MPa
P = 2.9 MPa
P = 5 MPa
12
10
C(gr/100gr)
8
0
0 20 40 60 80 100 120
σ PE(MPa)
5
P=2.9 MPa
P=5 MPa
4 P=8.4 MPa
C (gr/100gr)
3
0
0 20 40 60 80 100 120
σPE(MPa)
1,5
P = 2.9 MPa
P = 5 MPa
P = 8.4 MPa
1
C(gr/100gr)
0,5
0
0 20 40 60 80 100 120
σPE(MPa)
0,012
0,01
C (gr/100g)
0,008
0,006
0,004
0,002
0
0 20 40 60 80 100 120
σPE (MPa)
libre. Cette contrainte unique est nécessaire pour représenter correctement l’influence
de la phase cristalline sur la solubilité dans la phase amorphe.
Comme la valeur de contrainte effective (σef f ) est indépendante de la quantité
de gaz dissous et donc du taux de gonflement, au moins dans la gamme de pres-
sion étudiée, celle-ci ne peut pas être due à l’élasticité du réseau semi-cristallin. La
contrainte effective a donc pour effet d’augmenter la densité de la phase amorphe. En
effet, la densité moyenne de la région amorphe dans le polyéthylène semi-cristillin
est différente de celle d’une phase amorphe libre. L’ordre imposé par les chaı̂nes
appartenant aux régions cristallines doit persister sur des longueurs de l’ordre de
quelques segments de Kuhn16 . Cela conduit à une fraction importante d’interface
amorphe-cristal qui augmente la densité moyenne de la phase perméable.
par une augmentation du volume. Le taux de gonflement (Qa ) ainsi obtenu dans la
simulation est celui de la phase amorphe. Il est défini par la relation suivante :
a a
< Vswell > − < Vbulk >
Qa = a
(94)
< Vbulk >
a
dans laquelle < Vswell > est le volume moyen de la boı̂te de simulation en présence
a
du gaz et < Vbulk > est le volume moyen de la boı̂te contenant uniquement les
chaı̂nes de polymère sous la pression atmosphérique. Le volume moyen du système
de polyéthylène pur peut être obtenu par une simulation à part dans l’ensemble
NPT.
Afin de relier le gonflement de la phase amorphe au gonflement dans un polymère
semi-cristallin, il faut également prendre en compte le volume des régions cristallines.
Le volume total est la somme du volume de la phase amorphe Va , et de la phase
cristalline Vc . Le gonflement global (Q) qui est observable expérimentalement est :
a a
(< Vswell > +Vc ) − (< Vbulk > +Vc )
Q= a
(95)
< Vbulk > +Vc
Du fait que les régions cristallines du polyéthylène sont imperméables vis-à-vis
des gaz ; leurs volumes ne changent pas au cours de l’absorption. En conséquence,
la relation 95 se réduit à l’équation suivante :
Q = φa Qa (96)
où φa est la fraction volumique de la phase amorphe dans le matériau semi-cristallin
pur.
La figure 36 montre le volume moyen d’une boı̂te de simulation pour le système
CO2 -PE en fonction de la contrainte appliquée, pour différentes pressions de gaz. On
voit qu’à une contrainte donnée, le volume moyen augmente avec la pression de gaz.
Par ailleurs, pour une pression donnée de gaz, le volume diminue avec la contrainte
exercée. Comme l’absence de contrainte supplémentaire (σP E = 0) conduit à une
surestimation de quantité de gaz dissous, le taux de gonflement correspondant est
aussi très grand. Par contre, dès qu’une contrainte supplémentaire est appliquée,
un gonflement beaucoup plus petit est observé. Ainsi, pour le dioxyde de carbone
à la pression de 5 MPa on obtient Qa = 20 ± 1,6% en l’absence de contrainte
supplémentaire et Qa = 3 ± 0,8% pour une contrainte effective de 80 MPa.
Étant donné que l’existence d’une contrainte effective a été nécessaire pour ob-
tenir la bonne quantité de gaz dissous, on pense que le gonflement à contrainte
effective est proche de la valeur du gonflement expérimental. Pour le système CO2 -
PE, on a trouvé une donnée expérimentale du gonflement. Boyer et al. [154] ont
mesuré la solubilité du CO2 dans le polyéthylène semi-cristallin (φa = 0,5). Ensuite,
ils ont ajusté leurs données par les équations d’état de Sanchez-Lacombe (SL-EOS)
en réglant le paramètre d’interaction k12 entre le CO2 et le polyéthylène. L’équation
de Sanchez-Lacombe (SL-EOS) a été enfin utilisée pour prédire le degré de gonfle-
ment Q. Ce taux a été estimé à 2,1% pour une pression de 5 MPa à 333 K. Cette
4.4 Calcul de la solubilité dans l’ensemble Grand Canonique 68
valeur est très proche de celle obtenue dans ce travail sous contrainte effective. On
en déduit que l’utilisation d’une contrainte effective dans notre modèle donne une
valeur de gonflement réaliste.
36000
P = 0.8 MPa
P=2.9 MPa
P=5 MPa
PE pur
34000
32000
V(A3)
30000
28000
0 20 40 60 80 100 120
σPE(MPa)
Les volumes moyens des systèmes CH4 -PE sont également illustrés sur la fi-
gure 37. Le gonflement de la phase amorphe induit par du méthane à une pression
de 8,4 MPa est de 11,2 ± 2,6% à contrainte supplémentaire nulle et 2,4 ± 1,4% pour
une contrainte effective de 80 MPa. La figure 38 montre les résultats obtenus pour le
mélange N2 -PE. Pour une pression de 8,4 MPa de N2 , on obtient Qa = 5± 2,8% en
l’absence de contrainte supplémentaire et Qa = 1 ± 1% pour une contrainte effective
de 80 MPa. Le gonflement induit par l’azote est pratiquement négligeable dans cette
gamme de pression. Pour les systèmes méthane-polyéthylène et azote-polyéthylène,
on n’a pas trouvé des données expérimentales de gonflement dans la littérature.
Ces résultats sont alors prédictifs. Toutefois, il serait intéressant de disposer de ces
données expérimentales sur les taux de gonflement afin de les comparer avec les
simulations.
Les gonflements des systèmes étudiés sont donnés dans le tableau 5. Trois gran-
deurs sont présentées : Qa calculé en l’absence de contrainte supplémentaire (σP E = 0
MPa), Qa calculé sous contrainte effective (σef f = 80 MPa) et Q calculé pour une
fraction volumique d’amorphe φa = 0, 5 avec Qa obtenu sous contrainte effective.
33000
P=2.9 MPa
P=5 MPa
P=8.4 MPa
32000 PE pur
31000
V(A3)
30000
29000
28000
27000
0 20 40 60 80 100 120
σPE(MPa)
31000
P=2.9 MPa
P=5 MPa
P=8.4 MPa
PE pur
30000
V(A3)
29000
28000
27000
0 20 40 60 80 100 120
σPE(MPa)
dans la gamme de pressions étudiées et sous une contrainte effective de 80 MPa est
faible. Selon le tableau 5, ce gonflement n’excède pas un taux de 3%. Des taux de
gonflement aussi faibles suggèrent que les solubilités pourraient être obtenues par
des simulations à volume constant comme les simulations dans l’ensemble Grand
Canonique (GCMC).
4.5 Effet de longueur de chaı̂ne 70
N2 CH4 CO2
P (MPa) 8,4 8,4 5
T (K) 298 298 293
Qa , σP E = 0 5% ± 2,8 11,2% ± 2,6 20% ± 1,6
Qa , σP E = σef f 1±1 2,4% ± 1,4 3% ± 0,8
Q, σP E = σef f 0,5% ± 0,5 1,2% ± 0,7 1,5% ± 0,4
10
CO2
8
6
4
2
0
Concentration (gr/100gr)
CH4
1,5
0,5
0
0,6 N2
0,5
0,4
0,3
0,2
0,1
0
0 2 4 6 8
Pression (MPa)
entes ont été réalisées sur une boı̂te de simulation contenant des chaı̂nes de nC70 .
Donc, il nous semble important d’étudier l’effet de la longueur de chaı̂ne sur les
résultats obtenus. Selon Dee et al. [155], la densité d’un polyéthylène fondu croı̂t
avec la longueur de chaı̂ne, rapidement pour les chaı̂nes courtes jusqu’à un plateau
à partir de quelques centaines de carbones par chaı̂ne. Cette tendance n’a pas été
vérifiée expérimentalement à basse température. Cela vient de difficulté de la mesure
de densité de la phase amorphe du polyéthylène semi-cristallin.
En simulation moléculaire, le temps de calcul dépend directement de la longueur
des chaı̂nes (à nombre de chaı̂nes constant). Une boı̂te de simulation contenant des
chaı̂nes courtes sous-estime la densité expérimentale ; tandis qu’un modèle contenant
des chaı̂nes longues nécessite une très forte puissance de calcul. La longueur de 70
carbones par chaı̂ne qui a été choisie dans ce travail est une longueur optimale. La
densité correspondant à nC70 est seulement 2% inférieure à la densité des chaı̂nes
infiniment longues. En revanche, cette longueur peut être étudiée dans un temps
raisonnable, à savoir de l’ordre de deux à trois semaines pour une simulation dans
l’ensemble osmotique sur un processeur Intel actuel.
Revenons à la densité expérimentale du polyéthylène en fonction de la longueur
de chaı̂ne. Le polyéthylène à 298 K est à l’état semi-cristallin donc la mesure directe
de la densité de la phase amorphe à cette température n’est pas possible. Néanmoins,
en utilisant les données expérimentales de Dee et al. [155] pour le polyéthylène fondu,
la densité de la phase amorphe à basse température peut être estimée en utilisant
4.5 Effet de longueur de chaı̂ne 72
une équation d’état [156]. Les densités du polyéthylène amorphe à haute et basse
température sont présentées sur la figure 40. De même, cette figure montre l’évolution
de la densité en fonction de l’inverse du nombre de carbones par chaı̂ne Np . Comme
la figure 40 le montre, la densité de nC70 obtenue par simulation NPT correspond
parfaitement avec la densité extrapolée. Elle est toutefois inférieure à la densité d’un
polyéthylène réel pour lequel 1/Np tend vers zéro.
880
σeff > 0
860 σeff = 80 MPa
840
820
ρPE(kg/m3)
800
780
720
0 0,005 0,01 0,015 0,02 0,025
1/NP
Sa Sa Sa Sa
(CO2 ) (CH4 ) (N2 ) (H2 )
0, 906(a) ± 11% 0, 125(b) ± 20% 0, 074(a) ± 10% 0, 00338(c) ± 12%
333 K
0, 763(b) ± 20% 0, 045(b) ± 34% 0,00423(a) ± 15 %
0, 640(b) ± 30% 0, 131(b) ± 30% 0, 058(b) ± 30% 0, 00375(c) ± 12%
353 K
0, 700(c) ± 20% 0, 156(c) ± 20% 0, 0446(c) ± 20%
298 K
ln Sa(gr/100grMPa)
-2,8
-3
-3,2
-3,4
-3,6
-3,8
-4
2,8 2,9 3 3,1 3,2 3,3
-1
1000/T(K )
20
CO2 P = 5 MPa
Ash
15 Serpe
10
Concentration (gr/100gr)
0
CH4 P = 5 MPa
Serpe
1,5
0,5
0
0 20 40 60 80 100 120
σPE(MPa)
N2
P = 5 MPa
Ash
1
Concentration (gr/100gr)
0,5
0
H2 P = 0,6 MPa
0,0035
Flaconneche
Ash
0,003
0,0025
0,002
0,0015
0,001
0 20 40 60 80 100 120
σPE(MPa)
30
25
20
C (gr/100gr) 15
10
-5
0 20 40 60 80 100 120 140 160 180 200
σPE(MPa)
la solubilité de l’azote à 353 K, la valeur extrapolée a été utilisée (voir les paragraphes
précédents). D’après les figures 45 et 46, la valeur de contrainte effective est de
l’ordre de 40 ± 10 MPa à 353 K. Certaines déviations sont observées par rapport
à ces valeurs (par exemple pour CO2 à 353 K ou H2 à 333 K), mais cela peut être
relié à la quantité et à la qualité des données expérimentales. Il est important de
noter que peu de données expérimentales sont disponibles à ces températures, ce
qui rend la détermination de σef f incertaine. Cependant, la tendance générale est
une diminution nette de la contrainte effective, de 80 ± 10 MPa à 298 K jusqu’à
à peu près 40 ± 10 MPa à 353 K. Le tableau 7 résume la valeur de contrainte
effective pour chaque température. Ce comportement soutient notre interprétation
sur la signification de cette contrainte : en augmentant la température, la cristallinité
diminue ainsi que la densité effective des phases perméables ; par conséquent, σef f
diminue.
On peut aller encore plus loin en associant une énergie d’activation (∆Hef f ) à
la contrainte effective. En ajustant les valeurs de σef f par une loi d’Arrhenius :
on arrive à une énergie d’activation d’environ -10,4 ± 6,2 kJ/mol (cf. figure 47). Cette
énergie pourrait être considérée comme la barrière énergétique liée à la relaxation
de la densité effective de la phase perméable du polyéthylène semi-cristallin.
4.7 Conclusion 77
8 CO2 P = 4 MPa
Flaconneche
Serpe
6
Concentration (gr/100gr)
2
0
CH4 P = 4 MPa
Flaconneche
Serpe
1
0,5
0
0 20 40 60 80 100 120
σPE(MPa)
2
N2
P = 10 MPa
1,5 Extr
1
Concentration (gr/100gr)
0,5
0
H2
P = 2,1 MPa
Flaconneche
0,015
0,01
0,005
0
0 20 40 60 80 100 120
σPE(MPa)
4.7 Conclusion
Au cours de ce chapitre, nous avons exposé l’application des simulations de
Monte Carlo afin d’obtenir les solubilités de gaz simples dans le polyéthylène semi-
cristallin. Nous avons eu pour dessein d’étudier l’influence de la morphologie des
polymères semi-cristallins sur la solubilité des gaz. Il est connu à partir des observa-
4.7 Conclusion 78
4,6
4,4
4,2
298 K
ln σeff(MPa)
4 353K
333 K
3,8
3,6
3,4
3,2
3
2,8 2,9 3 3,1 3,2 3,3
-1
1000/T (K )
tions expérimentales [2, 9, 22, 31, 157] que les polymères semi-cristallins contiennent
des régions amorphes, des régions cristallines et des régions interfaciales. Les régions
perméables sont formées des régions amorphes et interfaciales.
Nos simulations, en accord avec des travaux antérieurs [90–94], montrent que les
solubilités calculées dans la phase amorphe libre sont plus grandes que les données
expérimentales ; alors que les résultats de simulation pour l’état fondu sur le même
système sont en bon accord avec l’expérience [137]. Il s’en suit que les différences
observées ne seraient être attribuées au champ de force ou aux mouvements de
Monte Carlo utilisés. Le problème relève certainement de la description de la phase
amorphe dans l’état semi-cristallin.
En utilisant l’ensemble osmotique d’une manière originale, nous avons démontré
que les solubilités expérimentales étaient reproduites quand une contrainte supplém-
entaire unique σef f est exercée sur la phase polymère. Le point crucial réside dans
le fait que cette contrainte ne dépend pas de la pression de gaz étudiés (dans la
gamme de 0-8 MPa) ou de la nature du gaz (N2 , CO2 , CH4 et H2 ). Elle est plutôt la
caractéristique du polymère à une température et une cristallinité donnée. De plus,
nous avons observé que cette contrainte devait être appliquée même pour les gaz très
peu solubles (comme N2 et H2 ) qui induisent pratiquement aucun gonflement dans
4.7 Conclusion 79
Ci = Si Pi = Si (xi P ) (99)
où, Si représente la solubilité partielle de l’espèce i. Lorsqu’il n’y a pas de
compétition entre les différentes espèces, la solubilité partielle d’un composant est
indépendante de la présence des autres composants. Dans ce cas, la solubilité par-
tielle de l’espèce i dans le mélange (Si ) est égale à la solubilité de cette espèce
lorsqu’elle est à l’état pur (Sipur ) :
Si = Sipur (100)
La relation (100) s’appelle «condition d’idéalité». La condition d’idéalité pour
différents mélanges sera étudiée au cours de ce chapitre.
Dans la deuxième approche, la concentration dissous d’une espèce est reliée à sa
fugacité. Ainsi, on a :
CH4 -CO2
f = XP
Simulation NPT
2,5
fCH4 (MPa)
1,5
0,5
0
0 0,2 0,4 0,6 0,8 1
G
X CH4
f = XP
Simulation NPT
2,5
2
fCO2 (MPa)
1,5
0,5
0
0 0,2 0,4 0,6 0,8 1
G
X CO2
obtenus dans l’ensemble NPT. Sur la figure 48 la fugacité du méthane est tracée en
5.3 Mélanges binaires : résultats des simulations dans l’ensemble
osmotique 85
0
(Sim) 80% CH4, 20% CO2
2,5
(Sim) 95% CH4, 5% CO2
2 (Exp) 80% CH4, 20% CO2
(Exp) 95% CH4, 5% CO2
1,5
0,5
0
0 20 40 60 80 100 120
σPE(MPa)
Les fugacités ainsi obtenues peuvent être utilisées comme paramètres d’entrée
des simulations dans l’ensemble osmotique. La phase polymère contient 15 chaı̂nes
de nC70 . La figure 50 illustre le résultat issu de la simulation dans l’ensemble os-
motique. Elle présente la concentration totale du gaz absorbé en fonction de la
contrainte supplémentaire σP E sur la phase amorphe. Chaque couleur représente une
composition différente et les symboles pleins représente les données expérimentales
correspondantes. Les résultats montrent qu’une seule contrainte effective suffit à re-
produire les solubilités expérimentales de différentes compositions de mélange. Une
contrainte de l’ordre de 80 MPa est en bon accord avec celle obtenue pour les gaz
purs dans le chapitre précédent. Cet accord confirme que la contrainte effective est
une caractéristique de la phase amorphe du polyéthylène et qu’elle ne dépend pas
de la nature de l’absorbant.
La simulation du mélange CH4 -CO2 permet également de déterminer la composi-
tion de l’absorbant dans le polymère. Les fractions de méthane et de carbone dioxyde
dans le polyéthylène sont calculées sous contrainte effective de 80 MPa. La figure
5.3 Mélanges binaires : résultats des simulations dans l’ensemble
osmotique 86
1
(Sim σPE = 80 MPa)
PE G
0.9 X =X
0.8
0.7
0.6
CH4
PE
0.5
X
0.4
0.3
0.2
0.1
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
G
X CH4
PE
Fig. 51 – Fraction molaire de méthane dans le polyéthylène (XCH4 ) en fonc-
G
tion de la fraction molaire du méthane dans la phase gazeuse (XCH4 ) pour un
mélange méthane-dioxyde de carbone à 3 MPa et 308 K. Les symboles ouverts
sont les résultats des simulations dans l’ensemble osmotique avec une contrainte
supplémentaire de 80 MPa. La ligne pleine à une composition identique dans les
deux phases. La ligne pointillée est une approximation polynomiale des résultats de
simulation.
0,4
0,3
CCH4 (gr/100gr)
0,2
0,1
0
0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1
G
X CH4
3
(Sim σPE = 80 MPa)
2,5
2
CCO2 (gr/100gr)
1,5
0,5
0
0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1
G
X CO2
0.5
20
0.5
0
0 1 2 3 4 5 6 7 8 9 10
r(A)
Fig. 54 – Fonction de distribution radiale g(r) CH4 -CH2 (en haut) et CO2 -CH2
(en bas) pour deux compositions différentes du mélange. La phase gazeuse est un
mélange binaire de méthane et de dioxyde de carbone à 3 MPa et 308 K. Les simu-
lations ont été effectuées dans l’ensemble osmotique avec une contrainte de 80 MPa.
Les résultats présentés jusqu’ici montrent que la présence d’une espèce dans le
polyéthylène n’affecte pas l’autre espèce. L’étude de la distribution des molécules par
le biais des fonctions de distribution radiale g(r), permet de vérifier cette hypothèse.
Afin d’étudier la perturbation éventuelle du système, les fonctions de distribution
radiale sont comparées pour les différentes compositions de mélange. Les figures 54 et
55 montrent les fonctions g(r) du CH4 et du CO2 par rapport aux groupements CH2
et CH3 du polyéthylène. Quand le nombre des molécules de gaz dans le polymère
est plus important, les courbes sont moins bruitées. Selon ces figures, le premier pic
de la fonction de distribution radiale g(r) CH4 -CH3 est plus grand que le premier
pic de la fonction g(r) CH4 -CH2 . Cela montre que les molécules de méthane ont
plutôt tendance à se regrouper autour des groupements CH3 , grâce à un effet de
volume libre plus important autour de ces groupements qu’autour de CH2 . En plus,
ces courbes montrent que les fonctions de distributions radiales ne sont pas affectées
par la composition du mélange. Autrement dit, les distributions du CH4 et du CO2
autour des chaı̂nes, ne sont pas perturbées par la quantité de l’autre espèce.
CH4 -H2
5.3 Mélanges binaires : résultats des simulations dans l’ensemble
osmotique 89
gCH4-CH3(r)
1.5
0.5
20
0.5
0
0 1 2 3 4 5 6 7 8 9 10
r(A)
Fig. 55 – Fonction de distribution radiale g(r) CH4 -CH3 (en haut) et CO2 -CH3
(en bas) pour deux compositions différentes du mélange. La phase gazeuse est un
mélange binaire de méthane et de dioxyde de carbone à 3 MPa et 308 K. Les simu-
lations ont été effectuées dans l’ensemble osmotique avec une contrainte de 80 MPa.
Simulation NPT
f = XP
1,5
fCH4 (MPa)
0,5
0
0 0,2 0,4 0,6 0,8 1
G
X CH4
Simulation NPT
f = XP
2,5
2
fH2 (MPa)
1,5
0,5
0
0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1
G
X H2
0,6
CCH4 (gr/100gr)
0,4
0,2
0
0 20 40 60 80 100 120
σPE (MPa)
0,002
CH2 (g/100g)
0,0015
0,001
0,0005
0
0 20 40 60 80 100 120
σPE (MPa)
fonction de distribution radiale g(r) H2 -CH3 est plus grand que le premier pic de la
fonction g(r) H2 -CH2 . Par conséquent, les molécules de gaz ont plutôt tendance à
se regrouper autour des groupements CH3 qu’autour de CH2 . Comme prévu, aucun
changement de distribution des gaz autour des groupements CH3 et CH2 des chaı̂nes
n’est observé. Cela tient du fait que la présence d’un gaz dans le polyéthylène n’af-
fecte pas la solubilité de l’autre.
CH4 -H2 S
0.7
0.6
0.5
H2
PE
X
0.4
0.3
0.2
0.1
-0.1
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
G
X H2
PE
Fig. 60 – Fraction molaire d’hydrogène dans le polyéthylène (XH2 ) en fonction de
G
la fraction molaire d’hydrogène dans la phase gazeuse (XH2 ) pour un mélange de
méthane-hydrogène à 2,1 MPa et 313 K. Les symboles ouverts sont les résultats
des simulations dans l’ensemble osmotique avec une contrainte supplémentaire de
80 MPa. La ligne pleine correspond à une composition identique du mélange dans
les deux phases. La ligne pointillée est une approximation polynomiale des résultats
de simulation.
0.2
CCH4 (g/100g)
0.15
0.1
0.05
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
G
X CH4
Flaconneche et al.
(Sim σPE = 80 MPa)
0,005
0,004
CH2 (g/100g)
0,003
0,002
0,001
0
0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1
G
X H2
CH4-H2 (20-80) / PE
1.5 CH4-H2 (90-10) / PE
gCH4-CH2(r)
0.5
20
H2 / PE
1.5 CH4-H2 (20-80) / PE
CH4-H2 (90-10) / PE
gH2-CH2(r)
0.5
0
0 1 2 3 4 5 6 7 8 9 10
r(A)
Fig. 63 – Fonction de distribution radiale g(r) CH4 -CH2 (en haut) et H2 -CH2 (en
bas) pour deux compositions différentes du mélange hydrogène-méthane et pour
l’hydrogène pur. La phase gazeuse est à 2,1 MPa et 313 K. Les simulations ont été
effectuées dans l’ensemble osmotique avec une contrainte de 80 MPa.
Après avoir déterminé les fugacités, les simulations d’absorption dans l’ensemble
osmotique peuvent être réalisées. Les simulations ont été effectuées sous une con-
5.3 Mélanges binaires : résultats des simulations dans l’ensemble
osmotique 95
gCH4-CH3(r)
1.5
0.5
40
H2 / PE
3 CH4-H2 (20-80) / PE
CH4-H2 (90-10) / PE
gH2-CH3(r)
0
0 1 2 3 4 5 6 7 8 9 10
r(A)
Fig. 64 – Fonction de distribution radiale g(r) CH4 -CH3 (en haut) et H2 -CH3 (en
bas) pour deux compositions différentes du mélange hydrogène-méthane et pour
l’hydrogène pur. La phase gazeuse est à 2,1 MPa et 313 K. Les simulations ont été
effectuées dans l’ensemble osmotique avec une contrainte de 80 MPa.
Simulation NPT
f = XP
2,5
2
fCH4 (MPa)
1,5
0,5
0
0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1
G
X CH4
Simulation NPT
f = XP
2,5
fH2S (MPa)
1,5
0,5
0
0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1
G
X H2S
PE G
X =X
(Sim σPE = 80 MPa)
0.8
0.6
H2S
PE
X
0.4
0.2
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
G
X H2S
PE
Fig. 67 – Fraction molaire de sulfure d’hydrogène dans le polyéthylène (XH2S ) en
G
fonction de la fraction molaire de sulfure d’hydrogène dans la phase gazeuse (XH2S )
pour un mélange méthane-sulfure d’hydrogène à 3 MPa et 308 K. Les symboles ou-
verts sont les résultats des simulations dans l’ensemble osmotique avec une contrainte
supplémentaire de 80 MPa. La ligne pleine correspond à une composition identique
dans les deux phases. La ligne pointillée est une approximation polynomiale des
résultats de simulation.
5.3 Mélanges binaires : résultats des simulations dans l’ensemble
osmotique 97
0.5
0.4
0.3
CCH4 (g/100g)
0.2
0.1
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
G
X CH4
4
CH2S (g/100g)
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
G
X H2S
effet, l’absorption du sulfure d’hydrogène dans le polyéthylène est plus favorisée par
rapport à l’absorption du méthane. Les figures 68 et 69 présentent les concentrations
5.3 Mélanges binaires : résultats des simulations dans l’ensemble
osmotique 98
de CH4 et d’H2 S en fonction des fractions molaires. Ces figures montrent que :
G
CCH4 ∝ XCH4
G (106)
CH2S ∝ XH2S
En s’appuyant sur l’équation (99), on peut également écrire :
G
CCH4 = (SCH4 P )XCH4
G (107)
CH2S = (SH2S P )XH2S
Ces équations indiquent que les solubilités du méthane (SCH4 ) et du sulfure d’hy-
drogène (SH2S ) sont indépendantes de la composition du mélange. Par conséquent,
elles sont égales aux solubilités du méthane et du sulfure d’hydrogène purs dans le
polyéthylène.
CH4-H2S (30-70) / PE
1.5 CH4-H2S (95-5) / PE
gCH4-CH2(r)
0.5
20
CH4-H2S (30-70) / PE
1.5 CH4-H2S (95-5) / PE
gH2S-CH2(r)
0.5
0
0 1 2 3 4 5 6 7 8 9 10
r(A)
Fig. 70 – Fonction de distribution radiale g(r) CH4 -CH2 (en haut) et H2 S-CH2
(en bas) pour deux compositions différentes du mélange. La phase gazeuse est un
mélange de méthane et de sulfure d’hydrogène à 3 MPa et 308 K. Les simulations
ont été effectuées dans l’ensemble osmotique avec une contrainte de 80 MPa.
Les fonctions de distribution radiale des gaz par rapport aux groupements CH2
et CH3 sont présentées sur les figures 70 et 71. Les fonctions g(r) montrent une autre
fois que les molécules de gaz ont plutôt tendance à se regrouper autour des groupe-
ments méthyle qu’autour des groupements méthylène. La répartition des différentes
espèces dans le polymère est pratiquement identique quelle que soit la composition
du mélange. Ainsi, la solubilité de méthane n’est affecté pas la présence de sulfure
d’hydrogène et vice versa.
5.4 Mélange ternaire 99
gCH4-CH3(r)
1.5
0.5
03
1.5
0.5
0
0 1 2 3 4 5 6 7 8 9 10
r(A)
Fig. 71 – Fonction de distribution radiale g(r) CH4 -CH3 (en haut) et H2 S-CH3
(en bas) pour deux compositions différentes du mélange. La phase gazeuse est un
mélange de méthane et de sulfure d’hydrogène à 3 MPa et 308 K. Les simulations
ont été effectuées dans l’ensemble osmotique avec une contrainte de 80 MPa.
Simulation NPT
0,45 f = XP
0,4
0,35
0,3
fH2S (MPa)
0,25
0,2
0,15
0,1
0,05
0
0 0,005 0,01 0,015 0,02 0,025 0,03 0,035 0,04 0,045
G
X H2S
Simulation NPT
f = XP
8,5
8
fCH4 (MPa)
7,5
6,5
6
0,6 0,65 0,7 0,75 0,8 0,85
G
X CH4
Simulation NPT
1,625
fCO2 (MPa)
1,62
1,615
1,61
0 0,005 0,01 0,015 0,02
G
X H2S
0.6
0.5
CH2S (g/100g)
0.4
0.3
0.2
0.1
0
0 0.01 0.02 0.03
G
X H2S
cette même valeur a été reprise pour calculer la solubilité du mélange CH4 -CO2 -H2 S à
cette température. La figure 75 montre les résultats de la simulation pour le sulfure
5.4 Mélange ternaire 102
2
1,9
(Sim σPE = 60 MPa)
1,8
1,7
1,6
1,5
1,4
CCH4(g/100g)
1,3
1,2
1,1
1
0,9
0,8
0,7
0,6
0,5
0,5 0,55 0,6 0,65 0,7 0,75 0,8 0,85 0,9 0,95 1
G
X CH4
ce qui conduit à une solubilité partielle de H2 S qui est constante. Ainsi, la solu-
bilité de H2 S,SH2S , ne dépend pas de la quantité d’H2 S dans le mélange. Cependant,
on ne peut pas relier cet effet à l’idéalité de la solubilité ; il faut étudier encore les
concentrations des autres espèces. La concentration du méthane en fonction de sa
fraction molaire est présentée sur la figure 76. Cette figure montre que la présence
d’une petite fraction de H2 S et de CO2 peut changer la solubilité du méthane. Dans
la gamme des compositions étudiées, la concentration de méthane n’est pas propor-
tionnelle à sa fraction molaire. La figure 77 présente la concentration de dioxyde
de carbone pour les trois compositions étudiées. Bien que les fractions molaires de
dioxyde de carbone dans le mélange gazeux soient identiques pour toutes les com-
positions, la concentration de CO2 est différente. La question est la suivante : cette
différence est-elle due à un comportement non-idéal du mélange ou bien elle est liée
5.5 Conclusion 103
σPE = 60 MPa
1,9
1,8
1,7
1,5
1,4
1,3
1,2
1,1
1
0 0,01 0,02 0,03 0,04
G
X H2S
aux erreurs statistiques des résultats de calcul dans cette fenêtre étroite de composi-
tion. Pour répondre à cette question, il faut étudier des compositions très différentes
du mélange qui dépasse la période de ce travail.
Dans l’intention de mieux visualiser la perturbation du système, plusieurs fonc-
tions de distribution radiales sont tracées sur les figures 78 et 79. Ces figures montrent
que les distributions des différentes espèces gazeuses autour des groupements ter-
minaux des chaı̂nes et autour des groupements CH2 ne sont pas influencées d’une
façon significative pour différentes compositions étudiées.
5.5 Conclusion
Dans ce chapitre, trois mélanges binaires et un mélange ternaire ont été étudiés.
La fugacité de chaque composant du mélange a été obtenue dans la phase gazeuse.
Ensuite, la solubilité de chacun des mélanges dans la phase amorphe du polyéthylène
semi-cristallin a été évaluée par simulation dans l’ensemble osmotique. À l’instar de
la solubilité des gaz purs, une contrainte effective sur la phase amorphe a été exercée
pour reproduire les données expérimentales. Ces travaux conduisent à la même valeur
de la contrainte effective. Ce résultat confirme l’unicité de la contrainte effective et
le fait qu’elle est seulement liée à la caractéristique de la phase polymère.
En ce qui concerne les mélanges binaires étudiés, aucun effet de non idéalité de
solubilité n’est observé : la solubilité en mélange est égale à la solubilité des corps
5.5 Conclusion 104
Ter (2-20-QS)
gCH4-CH2(r)
1,5
CH4 / PE
1
0,5
20
0,5
20
0,5
0
0 1 2 3 4 5 6 7 8 9 10
r(A)
Fig. 78 – Fonction de distribution radiale g(r) CH4 -CH2 (en haut), CO2 -CH2 (en
milieu) et H2 S-CH2 (en bas) pour le mélange ternaire de H2 S-CO2 -CH4 et pour
les gaz purs. Le mélange ternaire méthane, sulfure d’hydrogène dioxyde de carbone
dans la phase gazeuse est à 10 MPa. Le méthane pur, le dioxyde de carbone pur et
le sulfure d’hydrogène pur sont respectivement à 5 MPa, 5 MPa et 1,5 MPa. Les
simulations ont été effectuées dans l’ensemble osmotique avec une contrainte de 60
MPa à 333 K.
3
2,5 Ter (2-20-QS)
gCH4-CH3(r)
2 CH4 / PE
1,5
1
0,5
0
1,5 CO2 / PE
gCO2-CH3(r)
Ter (2-20-QS)
1
0,5
40
3 Ter (2-20-QS)
H2S / PE
gH2S-CH3(r)
0
0 1 2 3 4 5 6 7 8 9 10
r(A)
Fig. 79 – Fonction de distribution radiale g(r) CH4 -CH3 (en haut), CO2 -CH3 (en
milieu) et H2 S-CH3 (en bas) pour le mélange ternaire de H2 S-CO2 -CH4 et pour
les gaz purs. Le mélange ternaire méthane, sulfure d’hydrogène dioxyde de carbone
dans la phase gazeuse est à 10 MPa. Le méthane pur, le dioxyde de carbone pur et
le sulfure d’hydrogène pur sont respectivement à 5 MPa, 5 MPa et 1,5 MPa. Les
simulations ont été effectuées dans l’ensemble osmotique avec une contrainte de 60
MPa à 333 K.
5.5 Conclusion 106
107
τβ
P ea ' DCaoutchouc Sa = Pe (113)
α
L’hypothèse d’équivalence entre Da et DCaoutchouc n’est pas assez solide. Du point
de vu moléculaire, les structures des chaı̂nes de caoutchouc et de polyéthylène sont
différentes. De plus, cette hypothèse est fondée sur les comportements de ces po-
lymères à l’état fondu et n’est pas vérifiée dans l’état semi-cristallin.
Une troisième solution pour estimer le coefficient de diffusion dans la phase
amorphe a été proposée par Krishna et al. [106]. Ils ont considéré le polymère comme
une dispersion homogène des lamelles cristallines. Ils ont également supposé que la
diffusion dans la direction perpendiculaire aux lamelles est nulle. Puis, en établissant
une analogie entre le tenseur de diffusion et le tenseur diélectrique ils ont donné une
estimation de la limite inférieure de Da . Ils ont obtenu :
6.2 Diffusion 109
3
D1 = D (114)
2α
Où D1 est la limite inférieure du coefficient de diffusion dans la phase amorphe.
On a donc :
Da ≥ D1
Ainsi on peut également donner la limite inférieure de perméabilité :
3
P e1 = D1 Sa = Pe (115)
2α2
où
P ea ≥ P e1
Krishna et al. [106] ont corrigé les données expérimentales de Michaels [69] par
l’équation (114) et ont comparé ces données corrigées à des résultats de simulation.
Les simulations de Krishna et al. donnent un bon ordre de grandeur par rapport
aux valeurs ainsi corrigées.
Dans ce travail, nous allons comparer nos résultats de simulation avec des valeurs
estimées par les trois méthodes présentées dans cette section.
6.2 Diffusion
6.2.1 De l’autodiffusion à la diffusion fickienne
Concernant le processus de diffusion, plusieurs coefficients de diffusion peuvent
être définis. Dans l’annexe A trois coefficients sont présentés : le coefficient d’auto-
diffusion (D1∗ ), le coefficient de diffusion Maxwell-Stefan D12 MS
et celui de diffusion
fickienne D12 . Dans l’expression de la perméabilité et dans les phénomènes de trans-
port, c’est le coefficient de diffusion fickien qui intervient. En revanche, ce qui sera
calculé par la dynamique moléculaire est le coefficient d’autodiffusion. Dans cette
partie on reliera ces deux grandeurs a priori différentes.
Dans un système gaz-polymère, la mobilité des chaı̂nes de polymère est beaucoup
plus restreinte que celle des petites molécules de gaz18 . La diffusion du polyéthylène
est environ 100 fois plus petite que celle des molécules absorbées. En effet, les centres
de masse des chaı̂nes sont pratiquement fixes durant la période de diffusion des gaz. Si
on néglige la vitesse moyenne des chaı̂nes devant la vitesse des particules absorbées,
la diffusion de Maxwell-Stefan se réduit au coefficient d’autodiffusion. Cela peut être
démontré rigoureusement à partir de l’équation (137) de l’annexe A en négligeant
u2 devant u1 . Ainsi, pour la diffusion de gaz dans un polymère, on a :
18
Ce phénomène est aussi observé dans nos simulations. On reviendra sur ce sujet dans la pro-
chaine partie
6.2 Diffusion 110
4,4
Sim (σPE = 80 MPa)
4,2
Q=1
4
3,8
3,6
3,4
ln [fCO2 / Pst] 3,2
2,8
2,6
2,4
2,2
1,8
1,6
-0,4 -0,2 0 0,2 0,4 0,6 0,8 1 1,2 1,4 1,6 1,8 2
ln [CCO2 (g/100g)]
MS
D12 ' D1∗
Par conséquent, avec la notation introduite au chapitre 2 pour la diffusion dans
la phase amorphe du polymère semi-cristallin, on aura :
Da = Da∗ Q (116)
où Da∗ est le coefficient d’autodiffusion dans la phase amorphe et Q est le facteur
thermodynamique donné par :
∂ ln (f /Pst )
Q=
∂ ln (C)
où f est la fugacité du gaz, Pst est la pression standard et C est sa concentration
du gaz dissous en g/100g.
Le facteur thermodynamique peut être évalué par différentes méthodes. Afin
d’estimer Q dans nos calculs, il faut connaı̂tre la fugacité du gaz en fonction de
sa concentration dissous. Cette information est déjà fournie par les simulations de
Monte Carlo dans l’ensemble osmotique. La fugacité est un paramètre d’entrée de
la simulation tandis que la concentration est le résultat du calcul. La pente du loga-
rithme de la fugacité en fonction du logarithme de la concentration pour différentes
6.2 Diffusion 111
4,4
Sim (σPE = 80 MPa)
4,2 Q=1
3,8
3,6
3,4
ln [fCH4 / Pst] 3,2
2,8
2,6
2,4
2,2
1,8
1,6
-2 -1,8 -1,6 -1,4 -1,2 -1 -0,8 -0,6 -0,4 -0,2 0 0,2 0,4
ln [CCH4 (g/100g)]
4,5
4
ln [fH2S / Pst]
3,5
2,5
1,5
0,6 0,8 1 1,2 1,4 1,6 1,8 2 2,2 2,4 2,6 2,8 3
ln [CH2S (g/100g)]
molécule 1
1000 molécule 2
2
[(ri(t))-ri(0)] (A ) 800
600
2
400
200
0
0 500 1000 1500
t (ps)
1000
2
<x >
2
<y >
2
<z >
100
MSD (A )
2
10
1
10 100 1000
t (ps)
Fig. 85 – Déplacements carrés moyens du CO2 dans le polyéthylène suivant les trois
directions de l’espace à 293 K. La simulation a été effectuée sous une contrainte
effective de 80 MPa avec une pression de gaz de 5 MPa.
6.2 Diffusion 116
1000
CO2
PE
100
MSD (A )
2
10
0,1
10 100 1000 10000
t (ps)
3000
9 2 a)
2500 10 <D> = 0,1306 ± 0,0056 m /s
2000
MSD (A )
2
1500
1000
500
0
b)
1/6 ∂MSD/∂t (A /ps)
0,04
2
0,02
-0,02
efficient de diffusion sera plus petit. La figure 88 montre que la contrainte effective
donne des coefficients de diffusion proches des valeurs estimées d’après les données
de Michaels [69] (par les équations 112 et 114), d’après les données de Pino [171](par
l’équation 114) et d’après les données issues de l’expérience du temps retard de La-
guna [172] (par l’équation 114). Le déplacement minimale des particules mesurable
en RMN est compris entre 0,1 et 1 µm [173]. Or ce déplacement est de l’ordre des
tailles des hétérogénéités du matériau semi-cristallin. Le coefficient de diffusion ob-
tenu en RMN ne peut pas être corrigé facilement pour donner celui dans la phase
amorphe. Par conséquent, nous n’avons pas corrigé les données RMN de Laguna.
Par ailleurs, les mesures RMN [172] conduisent à un coefficient de diffusion très
grand. Pour expliquer cet écart, il faut noter que l’échantillon utilisé en RMN conte-
nait uniquement 26% de cristallinité contre une cristallinité d’environ 50% dans les
expériences du temps retard. En dépit de ces raisonnements, l’origine de cet écart
important n’a pas été expliquée par Laguna. Toutefois, ils avouent qu’il y a une large
distribution des valeurs mesurées du coefficient de diffusion et qu’un ordre de gran-
deur de différence n’est pas anormal. Sur la figure 88, l’extrapolation des données de
Durrill et al. [174] à haute température a également été représentée. La valeur ainsi
obtenue est proche de la valeur obtenue par simulation sous contrainte nulle. Cela
est tout à fait cohérent puisque la méthode d’extrapolation ne prend pas en compte
l’effet de cristallinité. Donc, à basse température, elle donne la diffusion dans un
état hypothétique 100% amorphe où il n’y a aucune contrainte supplémentaire.
6.2 Diffusion 118
0,6
9
0,4
0,2
0 5 10 15 20 25 30
Ca (g/100g)
10 *Da (m /s)
0,6
2
9
0,4
0,2
0
0 1 2 3 4 5
Ca (g/100g)
0,5
10 *Da (m /s)
2
0,4
9
0,3
0,2
0 2 4 6 8 10
Ca (g/100g)
σef f (MPa) P (MPa) 109 *Da (m2 /s) 109 *Da (m2 /s)
(σP E = 0) (σP E = σef f )
0,8 0,428 ± 0,014 0,131 ± 0,006
CO2 (293 K) 80 2,9 0,612 ± 0,009 0,150 ± 0,004
5 1,142 ± 0,015 0,222 ± 0,009
2,9 0,402 ± 0,038 0,126 ± 0,007
CH4 (298 K) 80 5 0,536 ± 0,014 0,139 ± 0,005
8,4 0,823 ± 0,025 0,167 ± 0,010
1,5 - 0,353 ± 0,009
H2 S (333 K) 60 3 - 0,402 ± 0,015
5 - 0,529 ± 0,007
6.3 Perméabilité
La détermination des coefficients de diffusion par dynamique moléculaire et les
solubilités par Monte Carlo, nous fournit les informations nécessaires pour évaluer la
perméabilité. La perméabilité s’obtient en multipliant le coefficient de diffusion par
la solubilité. Ces deux dernières grandeurs sont calculées sous contrainte effective afin
de prendre en compte l’effet des régions cristallines. Les perméabilités du CO2 et du
CH4 sous une contrainte effective de 80 MPa sont représentées sur la figure 91.a pour
trois pressions de gaz. Les estimations effectuées à partir des données de Michaels [69]
figurent également sur ce graphe. Elles sont assez proches des perméabilités obtenues
par simulation. Cet accord est prometteur car il assure la validité de la méthodologie
employée aussi bien en dynamique moléculaire qu’en Monte Carlo.
Les perméabilités de l’H2 S dans le polyéthylène amorphe à 333 K pour trois pres-
sions sont illustrées sur la figure 91.b. En raison du manque de données expérimental-
es, les perméabilités ainsi obtenues sont purement prédictives. Les valeurs de perméa-
bilité des gaz purs sont données dans le tableau 10.
Pour finir cette partie, nous commenterons un comportement général de la perm-
éabilité. Il s’agit de l’augmentation de la perméabilité avec la pression. Etant donnée
une solubilité constante dans cette gamme de pression, cette augmentation est plutôt
due à l’augmentation de la diffusion avec la pression (phénomène de plastification).
Cet aspect est probablement très important lors de l’utilisation de polymères sous
forte pression.
20
CH4 (Sim, σPE = 80 MPa)
a) CO2 (Sim, σPE = 80 MPa)
15 a
CH4 (d’après les données de Michaels )
a
CO2 (d’après les données de Michaels )
b
10 CH4 (d’après les données de Michaels )
b
700
H2S (Sim, σPE = 60 MPa)
3
b)
60
6
50
40
30
0 1 2 3 4 5 6 7 8 9 10
P (MPa)
cm [ST P ] 3
P (MPa) 106 *Pea ( cm.s.M Pa)
2,9 2,37 ± 1,06
CH4 (298 K) 5 3,11 ± 0,71
8,4 3,05 ± 0,66
0,8 7,13 ± 1,97
CO2 (293 K) 2,9 6,92 ± 1,07
5 11,95 ± 1,99
1,5 38,8 ± 3,3
H2 S 3 46,8 ± 4,2
5 61,1 ± 3,7
trapolées de Durill [174]. Ces mesures ont été acquises à pression atmosphérique
donc elles ne sont pas comparables avec celles de la simulation. En revanche, la ten-
dance des données de Durrill avec la température coı̈ncide bien avec celle obtenue
par dynamique moléculaire. Les limites inférieures de diffusion calculées à partir des
expériences de Michaels [69], Pino [171] et de Flaconnèche [34] sont aussi cohérentes
6.4 Effet de la température 123
-7
Sim
Ed = 21 kJ/mol
-7,5 D’après les données de Laguna
e
b
D’après les données de Laguna
a
-8 D’après les données de Michaels
b
D’après les données de Michaels
g
D’après les données de Durill
-8,5 b
D’après les données de Pino
log Da (m /s) b
D’après les données de Flaconneche
2
-9
-9,5
-10
-10,5
-11
2,6 2,8 3 3,2 3,4 3,6 3,8 4
-1
1000/T (K )
-7
Sim
-7,5 ED = 29 kJ/mol
a
D’après les données de Michaels
b
-8 D’après les données de Michaels
g
D’après les données de Lundberg
b
-8,5 D’après les données de Flaconneche
-9
log Da (m /s)
2
-9,5
-10
-10,5
-11
-11,5
-12
2,6 2,8 3 3,2 3,4 3,6 3,8 4
-1
1000/T (K )
-3
Sim
EPe = 17 kJ/mol
-3,5 D’après les données de Michaels
a
b
D’après les données de Michaels
b
D’après les données de Pino
-4 b
D’après les données de Flaconneche
log Pea (cm [STP]/cm/MPa/s)
-4,5
3
-5
-5,5
-6
-6,5
2,6 2,8 3 3,2 3,4 3,6 3,8 4
-1
1000/T (K )
-3
Sim
EPe = 31 kJ/mol
-3,5 a
D’après les données de Michaels
b
D’après les données de Michaels
-4 b
D’après les données de Flaconneche
-5
3
-5,5
-6
-6,5
-7
2,6 2,8 3 3,2 3,4 3,6 3,8 4
-1
1000/T (K )
95. Elle augmente avec la température en suivant une loi d’Arrhenius. La valeur
estimée à partir des expériences de Michaels [69] est très proche des valeurs simulées.
Les limites inférieures de perméabilité issues des expériences de Flaconnèche [34] et
de Michaels [69] sont aussi plus petites que les résultats de simulation. Ceci témoigne
d’un accord raisonnable entre l’expérience et la simulation.
Les énergies d’activation de la diffusion et de la perméabilité pour les cas de
CH4 et de CO2 sont représentées dans le tableau 11. Les valeurs issues des données
expérimentales de la littérature ont été également représentées dans ce tableau.
La comparaison montre que l’énergie d’activation donnée par la simulation est
légèrement plus petite que celle donnée par l’expérience. Cette différence vient du
fait que la simulation moléculaire de la phase amorphe ne prend pas en compte
l’effet d’immobilisation des chaı̂nes par les régions cristallines. Cette immobilisation
diminue le coefficient de diffusion par un facteur β. Le calcul plus précis des énergies
d’activation nécessite un modèle plus réaliste de la région amorphe.
6.5 Conclusion
Au cours de ce chapitre, nous avons utilisé la dynamique moléculaire pour calcu-
ler les coefficients de diffusion des gaz purs dans la phase amorphe du polyéthylène
à basse température. Les simulations de dynamique moléculaire ont été effectuées
sur les configurations préalablement obtenues par Monte Carlo sous une contrainte
effective. Cette contrainte effective, qui permet de reproduire la densité moyenne
6.5 Conclusion 126
7 Conclusion et perspectives
La perméabilité de gaz dans les polymères est un champ d’investissement aussi
bien expérimentalement qu’en modélisation. Elle fait l’objet d’enjeux industriels
importants liés à l’utilisation des polymères pour leurs propriétés barrières vis-à-vis
du transport des gaz ou comme membrane dans les procédés de séparation.
Les polymères, au dessous de la température de fusion, sont à l’état semi-cristallin.
Dans cet état la structure microscopique du polymère est hétérogène et elle dépend
de plusieurs facteurs [9–12, 16, 17, 21–23, 26, 27]. Le polymère à l’état semi-cristallin
est principalement constitué des deux phases : amorphe et cristalline. La morphologie
des polymères semi-cristallins a été étudiée par plusieurs techniques expérimentales
[19, 22, 28–32, 35–46, 48, 49]. Différents modèles théoriques ont été proposés pour
décrire la configuration des chaı̂nes dans l’interface entre les régions amorphe et
cristalline [12, 47, 50, 52, 53]. En outre, il a été montré que la morphologie du po-
lymère a une influence importante sur les paramètres de transport des gaz. De
nombreux auteurs se sont intéressés à l’impact de la cristallinité sur la solubi-
lité [33, 54, 56–68, 175–177]. D’autres recherches ont été consacrées à l’étude de
diffusion dans les polymères semi-cristallins [54, 59, 68–78, 84]. Plusieurs tentatives
pour découpler les paramètres de transport dans la phase amorphe et dans la phase
cristalline ont échoué [90, 91, 93]. Cela est dû au fait que les cristallites changent
les propriétés de la région amorphe. L’effet des régions cristallines sur la phase
amorphe peut être étudié par deux modèles différents. Le premier est le modèle
élastique [96–98] basé sur la présence des chaı̂nes pontantes entre les régions cris-
tallines. Ces chaı̂nes forment un réseau limitant le gonflement du polymère. Dans le
deuxième modèle, les régions cristallines changent la densité moyenne de la phase
amorphe à travers des régions interfaciales [31].
L’objectif de cette thèse a été de mettre au point un modèle prédictif de perm-
éation des gaz dans le polyéthylène semi-cristallin pour de larges gammes de pres-
sion, de température et de composition chimique de gaz. Nous avons appliqué des
simulations Monte Carlo pour calculer la solubilité des gaz simples dans la phase
amorphe du polyéthylène. Nous avons eu pour intérêt d’étudier l’influence de la mor-
phologie complexe du polymère semi-cristallin sur la solubilité. En nous basant sur
les évidences expérimentales [9, 22, 31, 46, 157], il est connu que les polymères semi-
cristallins contiennent une fraction importante de région interfaciale qui contribue
à la perméabilité des gaz.
Nos simulations en accord avec d’anciens travaux de la littérature confirment que
les solubilités calculées dans une phase purement amorphe sont plus grandes que les
solubilités expérimentales. Etant donné que pour les mêmes systèmes à l’état fondu,
simulation et expérience sont en bon accord [137], on conclut que ces différences
ne peuvent pas être attribuées au champ de force ou à la méthode de Monte Carlo
utilisés. Apparemment, le problème vient de la description de la phase perméable
dans l’état semi-cristallin. Ensuite, nous avons fait l’hypothèse que d’une part le
réseau formé par les chaı̂nes pontantes entre les cristallites, et d’autre part, le chan-
128
La comparaison avec les données expérimentales de perméabilité [34, 69, 171] est
satisfaisante.
Les coefficients de diffusion et perméabilité obtenues par simulation à différentes
températures peuvent être ajustés avec une loi d’Arrhenius. Il en résulte que les
énergies d’activation associées à la diffusion et à la perméabilité dans la phase
amorphe sont plus petites que celles obtenues [34,69] dans les vrais matériaux semi-
cristallins. Cela vient du fait que dans nos simulations l’immobilisation des chaı̂nes
de polymère par des régions cristallines n’est pas prise en compte.
La capacité prédictive de la simulation moléculaire à calculer la solubilité dans
l’état semi-cristallin est limitée par la connaissance de la densité effective. La méthod-
ologie appliquée durant cette thèse, basée sur l’évaluation de cette quantité par
comparaison des données expérimentales et des simulations de Monte Carlo dans
l’ensemble osmotique, est une voie possible. Cependant, une description plus précise
de l’interface cristal-amorphe à l’échelle moléculaire est souhaitable.
Pour la poursuite de ce travail, il est envisageable de réaliser une simulation ex-
plicite des régions cristallines en contact avec les régions amorphes. La simulation
de l’interface par la méthode de Monte Carlo implique d’améliorer l’échantillonnage
de l’espace des configurations. Introduire de la polydispersité dans les longueurs de
chaı̂ne de polymère pourrait être bénéfique et permettrait d’accepter les mouve-
ments de repontage plus fréquemment. D’autres mouvements tels que le End Brid-
ging [115,181] pourraient être introduits dans le code. Ces mouvements changent les
connectivités des chaı̂nes et donc ils relaxent aussi la topologie de la région inter-
faciale. L’utilisation de pas de dynamique moléculaire dans les simulations Monte
Carlo (Hybrid Monte Carlo [182]) est une technique qui permettrait d’associer les
atouts de la dynamique moléculaire à ceux du Monte Carlo.
La proportion des chaı̂nes pontantes, des boucles, des repliements serrés ainsi que
la concentration de groupements terminaux dans la phase amorphe dépendent de
plusieurs facteurs et ils peuvent être différents d’un polymère à l’autre. La solution
la plus appropriée est de générer différentes conformations de chaı̂nes d’une façon
aléatoire. La figure 96 illustre une configuration de polyéthylène contenant d’une
manière explicite la phase amorphe et la phase cristalline. La modélisation de cette
configuration par dynamique moléculaire permettra d’étudier la dynamique segmen-
tale des chaı̂nes dans la phase amorphe en fonction de la distance à l’interface avec
le cristal. La prochaine étape serait de mettre des molécules de gaz dans la phase
amorphe et d’estimer le facteur d’immobilisation des chaı̂nes β.
Une étude plus approfondie du polymère semi-cristallin peut être entreprise à
la lumière des résultats obtenus dans cette thèse. Le rôle de la phase cristalline
ainsi que son influence sur la phase amorphe restent encore à élucider aussi bien de
manière qualitative que quantitative.
Références 131
Références
[1] F. Bedoui, J. Diani, G. Régner, and W. Seiler. Acta Materialia, 54 :1513–1523,
2006.
[2] L. Mandelkern. Accounts of Chemical Research, 23 :380–386, 1990.
[3] Paul J. Flory. The Journal of Chemical Physics, 17(3) :223–240, 1949.
[4] J. Maxfield and L. Mandelkern. Macromolecules, 10 :1141, 1977.
[5] R. C. Allen. Ph.D. Dissertation, School of Materials Engineering Sciences,
Virginia Polytechnique Institute and State University, 1981.
[6] R. C. Allen and L. Mandelkern. Journal of Polymer Science. Polymer Physics
Edition, 20 :1465, 1982.
[7] L. Mandelkern, M. Glotin, and R. S. Benson. Macromolecules, 14 :22, 1981.
[8] I. G. Voigt-Martin, R. G. Alamo, and L. Mandelkern. Journal of Polymer
Science. Polymer Physics Edition, 24 :1283, 1986.
[9] C. Hedesiu, D. E. Demco, R. Kleppinger, A. A. Buda, B. Blumich, K. Remerie,
and V. M. Litvinov. Polymer, 48 :763–777, 2007.
[10] J. Dlugosz, C. Fraser, D. Grubb, A. Keller, J. Odell, and P. Goggin. Polymer,
17 :471, 1976.
[11] R. A. Shanks and G. Amarasinghe. J. Therm. Anal. Cal., 59 :471–482, 2000.
[12] E. Robelin-Souffaché and J. Rault. Macromolecules, 22 :3581–3594, 1989.
[13] J. D. Hoffman, L. J. Frolen, G. S. Ross, and J. I. Lauritzen. J. Res. Natl. Bur.
Stand. Sect., 59 :677, 1975.
[14] J. Point. Macromolecules, 12 :770, 1979.
[15] G. Allegra. J. Chem. Phys., 66 :5453, 1977.
[16] J. Rault. J. Macromol. Sci. Phys., 15 :567, 1978.
[17] S. Rastogi, Yefeng Yao, Dirk R. Lippits, Gunter, W. H. Hohne, R. Graf,
Hans W. Spiess, and P. J. Lemstra. Macromolecular Rapid Communication,
30 :826–839, 2009.
[18] B. Wunderlich and G. Czornyj. Macromolecules, 10 :906–913, 1977.
[19] G. M. Stack, L. Mandelkern, and I. G. Voigt-Martin. Macromolecules, 17 :321,
1984.
[20] J. Martinez-Salazar, P. S. Barham, and A. Keller. Journal of Polymer Science.
Polymer Physics Edition, 22 :1085, 1984.
[21] L. Mandelkern, K. W. McLaughlin, and R. G. Alamo. Macromolecules,
25 :1440–1444, 1992.
[22] R. Kitamaru, F. Horii, and K. Murayama. Macromolecules, 19 :636–643, 1986.
[23] E. Perez and D. L. VanderHart. Macromolecules, 20 :78–87, 1987.
Références 132
[49] X. Jing and S. Krimm. J. Polym. Sci. Polym. Lett. Ed., 21 :123, 1983.
[50] S. K. Kumar and D. Y. Yoon. Macromolecules, 22 :3458–3465, 1989.
[51] P. J. Flory. John Wiley, new york edition, 1969.
[52] P. J. Flory. Journal of American Chemistry Society, 84 :2857, 1962.
[53] S. Balijepalli and G. C. Rutledge. J. Chem. Phys., 109 :6523, 1998.
[54] A.S. Michaels and R.B. Parker Jr. J. of Polymer Science, 41 :53–71, 1959.
[55] P. C. Carman. Butterworth, london edition, 1956.
[56] R. B. Richards. Trans. Faraday. Soc., 41 :127, 1945.
[57] A. W. Myers, C. E. Rogers, V. Stannett, and M. Szwarc. Modern Plastics,
34 :157, 1957.
[58] C. E. Roger, V. Stannett, and M. Szwarc. J. Phys. Chem., 63 :1406, 1959.
[59] V. Vittoria. J. of Materials Sci., 30 :3954–3958, 1995.
[60] P. Paricaud, A. Galindo, and G. Jackson. Ind. Eng. Chem. Res., 43 :6871–
6889, 2004.
[61] R. M. Barrer, J. A. Barrie, and J. Slater. J. Polymer Sci., 27 :315, 1958.
[62] A.S. Michaels, W.R. Vieth, and J.A. Barrie. J. Appl. Phys., 34 :1–12, 1963.
[63] A. C. Puleo. The University of Texas, Austin, phd thesis edition, 1988.
[64] T. Suzuki, T. Tanaka, M. Nakajima, H. Yoshumizu, and Y. Tsujita. Polym J,
34 :891–896, 2002.
[65] C. Manfredi, M. A. D. Nobile, G. Mensitieri, G. Guerra, and M. J. Rapacciulo.
J. Polym. Sci. Polym. Phys., 35 :133–140, 1997.
[66] L. Guadagno, P. Baldi, V. Vittoria, and G. Guerra. Macromol. Chem. Phys.,
199 :2671–2675, 1998.
[67] R. Cowling and GS Park. Journal of Membrane Science, 5 :199–207, 1979.
[68] Z. Mogri and D. R. Paul. Polymer, 42 :2531–2542, 2001.
[69] A. S. Michaels and H. J. Bixler. Journal of Polymer Science, 50 :413–439,
1961.
[70] C. H. Klute. J. Polym. Sci., 41 :307–317, 1959.
[71] R. J. Pace and R. Datyner. J Polym Sci, Polym Phys Ed, 17 :437, 1979.
[72] R. J. Pace and R. Datyner. J Polym Sci, Polym Phys Ed, 17 :453, 1979.
[73] R. J. Pace and R. Datyner. J Polym Sci, Polym Phys Ed, 17 :465, 1979.
[74] C. Kiparissides, V. Dimos, T. Boultouka, A. Anastasiadis, and A. Chasiotis.
Journal of Applied Polymer Science, 87 :953, 2002.
[75] M. R. Wyllie and M. B. Spangler. Bull. Am. Assoc. Petrol. Geologists, 36 :359,
1952.
[76] A.S. Michaels, W.R. Vieth, and J.A. Barrie. J. Appl. Phys., 34 :13–20, 1963.
Références 134
[104] W. Yao, X. Hu, and Y. Yang. Journal of Applied Polymer Science, 103 :1737,
2006.
[105] E.F. Castro, E.E. Gonzo, and J.C. Gottifredi. Journal of Membrane Science,
113 :57, 1995.
[106] P.V. Krishna Pant and R.H. Boyd. Macromolecules, 26 :679–686, 1993.
[107] N. F. A. van der Vegt. Macromolecules, 33 :3153, 2000.
[108] P. Gestoso and N. C. Karayiannis. J. Phys. Chem., 112 :5646, 2008.
[109] H. Fujita. Hochpolym. Forsch., 3 :1, 1961.
[110] D. Frenkel and B. Smit. Academic Press, san diego edition, 1996.
[111] E. Bourasseau, P. Ungerer, and A. Boutin. J. Phys. Chem. B, 106 :5483–5491,
2002.
[112] JJ. de Pablo, M. Laso, and UW. Suter. J. Chem. Phys., 96 :2395–2403, 1992.
[113] RH. Boyd. Macromolecules, 22 :2477–2481, 1989.
[114] L.R. Dodd, T.D. Boone, and D.N. Theodorou. Molecular Physics, 78(4) :961–
996, 1993.
[115] V.G. Mavrantzas, T.D. Boone, E. Zervopoulou, and D.N. Theodorou. Macro-
molecules, 32 :5072–5096, 1999.
[116] N. C. Karayiannis, A.E. Giannousaki, V.G. Mavrantzas, and D.N. Theodorou.
J. Chem. Phys., 117(11) :5465, 2002.
[117] N. C. Karayiannis, V. G. Mavrantzas, and D. N. Theodorou. Phys. Rev. Lett.,
88(10) :105503–4, 2002.
[118] D.N. Theodorou. Molecular Physics, 102(2) :147–166, 2004.
[119] R. H. Swendsen and J. S. Wang. Phys. Rev. Lett., 57 :2607–2609, 1986.
[120] A.Z. Panagiotopoulos. Mol. Phys., 61 :813–826, 1987.
[121] S.K. Nath and J.J. de Pablo. J. Phys. Chem. B, 103 :3539–3544, 1999.
[122] M.P. Allen and D.J. Tildesley. Clarendon Press, oxford edition, 1987.
[123] S. Nose. J. Chem. Phys., 81 :511, 1984.
[124] W. G. Hoover. Phys. Rev., 31 :1695, 1985.
[125] H. C. Andersen. J. Comp. Phys., 52 :24, 1983.
[126] G. J. Martyna, M. L. Klein, and M. Tuckerman. J. Phys. Chem., 97 :2635,
1992.
[127] G. J. Martyna. Mol. Phys., 87 :1117, 1996.
[128] C.L. Kong. J. Chem. Phys., 59 :2464–2467, 1973.
[129] J. delhommelle and P. Millie. Mol. Phys., 99 :619, 2001.
[130] P. Ungerer, A. Boutin, and B. Tavitian. IFP publications, 2005.
[131] S. Toxvaerd. J. Chem. Phys., 93 :4290–4295, 1990.
Références 136
[182] B. Mehlig, D.W. Heermann, and B.M. Forrest. Physical Review B, 45 :679–685,
1991.
[183] Denis J. Evan and Gary P. Morris. Academic Press, london edition, 1990.
[184] R. Krishna and J. M. van Baten. Ind. Eng. Chem. Res., 44 :6939–6947, 2005.
[185] Dean R. Wheeler and John Newman. J. Phys. Chem. B, 108 :18353–18361,
2004.
[186] M.H. Klopffer and B. Flaconnèche. Oil and Gas Science and Technology,
56(3) :223–244, 2001.
[187] J. E. Mark. Springer, 2nd edition, us edition, 2007.
139
Fig. 97 – Inter diffusion des particules ”rouges” dans les particules ”bleues”. On
distingue deux cas différents. 1) Système à l’équilibre : Toutes les particules sont
réparties d’une manière homogène. Les couleurs sont imaginaires et elles sont uti-
lisées uniquement comme des étiquettes. Dans ce cas, le processus de diffusion est
caractérisé par le coefficient de «autodiffusion». 2) Système hors équilibre : Les
espèces rouges et bleues correspondent à deux types de molécule. Chaque espèce su-
bit un gradient de concentration. Dans ce cas, le processus de diffusion est caractérisé
par le coefficient de diffusion «fickien».
Système à l’équilibre
Supposons maintenant que les particules bleues et rouges dans la figure 97 cor-
respondent à deux espèces différentes. Cette situation correspond à un système hors
20 ∂C →
−
∂t + ∇. J = 0
141
∞ Ni Nj
−→ −−→
Z
1 X X
Λij = dt h( Vi,k (0)).( Vj,m (t))i (128)
3N 0 k=1 m=1
où N est le nombre total de particules dans le système. D’après Krishna et al. [184]
la relation du type Einstein de cette matrice est la suivante :
Ni Nj
1 1 1 X
h( (rl,i (t + ∆t) − rl,i (t))).( (−
−
→ −
→ r→ −→
X
Λij = − lim k,j (t + ∆t) − rk,j (t)))i (129)
6 ∆t→∞ N ∆t l=1 k=1
Les équations (128) et (129) montrent aussi que la matrice d’Onsager est symétrique.
C’est-à-dire :
Λij = Λji (130)
Cette relation s’appelle la relation de réciprocité d’Onsager (ORR).
Etant donné que la diffusion dans les systèmes hors équilibre conduit à un trans-
port de matière, à chaque espèce i peut être associé un flux correspondant. Le flux
est souvent défini dans le référentiel du centre de masse du système. Ainsi, le flux
de l’espèce i est défini comme :
→
−
Ji = Ci (−
→
ui − −
→
u CM ) (131)
Avec cette définition le flux total de toutes les espèces est annulé :
X− →
Ji = 0 (133)
i
Cela signifie que pour un système de n constituants, il existe n−1 flux indépendants.
Afin de relier les flux à la variation de concentration, les gradients de potentiel
chimique dans l’équation (126) doivent être exprimés en fonction des gradients de
concentration. Cette transformation est possible à travers la relation suivante :
Pour simplicité, on a remplacé la constante P0i par la pression standard Pst (=1
bar). Ceci ne change pas la valeur de dérivée puisque P0i et Pst ne sont que des
constantes numériques. En remplaçant les ∇µi ainsi obtenus dans l’équation (126),
les vitesses moyennes de chaque espèce ainsi que le flux associé peuvent être exprimés
en fonction des gradients de concentration :
n−1
−
→ X
M S ∂ ln (fi /Pst )
Ji = − Di,j ∇Cj (135)
j=1
∂ ln (Cj )
MS x2 x1
D12 = Λ11 + Λ22 − Λ12 − Λ21 (136)
x1 x2
La diffusion de Maxwell-Stefan peut être également obtenue par une relation de type
Green-Kubo. Pour un mélange binaire elle est égale à [185] :
N1 N2 ∞
Z
MS
D12 = dt < (−→
u2 (t) − −
→
u1 (t)).(−
→
u2 (0) − −
→
u1 (0)) > (137)
3N 0
où
∂ ln (f1 /Pst )
Q= (139)
∂ ln (C1 )
143
EP e = ED + ∆HS (145)
Les mouvements utilisés pour les molécules de gaz sont interdits aux molécules
de polymère.
– Mouvements pour les gaz : translation, rotation, insertion/destruction
– Valeurs maximales de déplacement (translation) : ajustée tous les 10000 pas
MC pour avoir 40% de taux d’acceptation du mouvement de translation.
– Valeurs maximales de rotation : ajustée tous les 10000 pas MC pour avoir 40%
de taux d’acceptation du mouvement de rotation.
– Nombre d’essais pour le biais du mouvement d’insertion : 5
Les mouvements utilisés pour les molécules de polymère sont interdits aux molécules
de gaz.
– Mouvements utilisés : Flip, Recroissance, Reptation, ConRot, DB, IDR
– Nombre de centres de force déplacés lors du mouvement de reptation : 2
D.2 Simulation de dynamique moléculaire 150
Changements de volume
– Valeur maximale pour le changement de volume : ajustée tous les 10000 pas
MC pour avoir 40% de taux d’acceptation pour le changement de volume.
Champ de force
– Rayon du coupure : 10 Å
– Rayon pour la mise à jour des listes d’interaction : 11 Å