@vdb: la Physique en Classe Prépa

L'ensemble des documents que je propose sur ce site sont placés sous la protection de la licence Creative Commons.

Atmosphère isotherme: simulation de Monte-Carlo

On présente ici une recherche de l'état d'équilibre d'un système en contact avec un thermostat. L'ensemble statistique dans lequel cette étude s'inscrit est l'ensemble canonique. La recherche de l'état d'équilibre est réalisé par une simulation de Monte-Carlo.

Le système étudié:

Il s'agit de molécules en contact d'un thermostat à la température . Le système admet un nombre fini () de niveaux d'énergie discrets. L'énergie du niveau ( et ) est donnée par: où est la masse d'une molécule, l'accélération de la pesanteur et un paramètre homogène à une longueur. Comme il est nécessaire de travailler avec des quantités adimensionnées pour faire la simulation numérique, on introduit le paramètre tel que .

On simule ainsi une atmosphère isotherme de taille finie découpée en cases. Le caractère fini du système entraîne d'inévitables effets de bords, que l'on peut minimiser en choisissant des conditions aux limites périodiques (la case est située au-dessus de la case ).

La simulation de Monte-Carlo

 

Pourquoi ?

Il est tout à fait possible de réaliser une étude théorique de ce système, elle est d'ailleurs présentée ci-dessous. Alors, pourquoi une simulation de Monte-Carlo pour rechercher l'état d'équilibre ?

Le paradigme d'un système qu'il convient d'explorer via une simulation de Monte-Carlo est le modèle d'Ising. La richesse de ce système (existence possible d'une transition de phase à la température de Curie) ainsi que sa complexité théorique (intégrale d'échange et hamiltonien d'Heisenberg) font tout l'intérêt du modèle d'Ising mais aussi rendent impossible son approche par des étudiants des premières classes post-baccalauréat.

Le modèle présenté ici revêt un aspect théorique inaccessible à ces mêmes étudiants (notions de thermodynamique statistique) mais il est conceptuellement simple et permet d'illustrer ce qu'est la méthode de Monte-Carlo, de montrer ses performances (ce n'est pas pour rien qu'elle joue un rôle majeur dans la simulation numérique du comportement de systèmes physiques d'intérêt actuel comme les polymères).

 

Le principe de l'algorithme

Le but est de montrer que l'état d'équilibre du système étudié correspond à un minimum d'énergie libre du système et non pas à un minimum d'énergie interne. On veut aussi montrer que l'entropie est maximale à l'équilibre compte tenu de l'énergie (interne) du système.

Pour une valeur de la température donnée, il convient de faire évoluer le système à partir d'une configuration initiale arbitraire vers l'état d'équilibre. La configuration initiale retenue est celle où toutes les molécules sont condensées dans le premier niveau d'énergie: il s'agit de l'état d'équilibre à température nulle.

On génère donc des configurations successives du système, qui constituent une chaîne de Markov. Un processus de Markov est un processus pour lequel le futur, conditionné par le présent, est indépendant du passé. On désigne par une étape de cette chaîne et par la configuration correspondante.

Soit la probabilité que le système soit dans la configuration lors de l'étape . On désigne par la probabilité de transition de la configuration vers la configuration (on suppose cette probabilité de transition stationnaire). Alors, la probabilité de se retrouver dans la configuration à l'étape est:

Au bout d'un nombre d'itérations important ( grand), si le système arrive dans un état d'équilibre, alors on veut que .

On en déduit qu'une condition suffisante (mais pas nécessaire) pour que le système parvienne dans un état d'équilibre (donc stationnaire) est:

Cette condition porte le nom de bilan détaillé.

Comme on travaille dans l'ensemble canonique, on doit choisir la distribution de Boltzmann. Alors:

avec , d'où l'on déduit: avec .

Reste maintenant à choisir la probabilité de transition . Il existe plusieurs choix possibles, qui vérifient le principe du bilan détaillé. Nous retiendrons la fonction de Metropolis:

si et si . On peut vérifier que ce choix est compatible avec le principe du bilan détaillé et qu'il sera donc possible d'arriver à un état stationnaire du système grâce à cet algorithme. De plus, l'algorithme de Metropolis est en général ergodique: toute configuration peut être atteinte à partir d'une configuration quelconque en un nombre fini d'itérations. En conséquence, un des intérêts de l'algorithme de Metropolis est de ne pas laisser le système bloquer dans un état métastable puisque des transitions vers des configurations d'énergie plus élevée sont possibles !

 

Le problème des moyennes

En supposant que l'on ait fait suffisamment d'itérations pour atteindre l'état d'équilibre thermodynamique, il va falloir sélectionner parmi les configurations générées celles qui vont nous permettrent de calculer l'énergie interne, l'énergie libre et l'entropie du système. Clairement, les premières configurations générées n'ont aucun intérêt pour ce calcul car elles sont éloignées de l'état d'équilibre thermodynamique. On les écarte donc.

Précisons que lors de l'étude de systèmes plus complexes apparaissent d'autres difficultés. L'algorithme génère des configurations successives très ressemblantes. On peut alléger le calcul de la moyenne en ne retenant que certaines configurations régulièrement, pour avoir des mesures indépendantes ou du moins, le moins corrélées possible. Ensuite, pour s'affranchir des fluctuations, il est souvent profitable de calculer la moyenne des moyennes de plusieurs simulations correspondant à la même valeur de la température.

> restart:with(stats): with(statplots): with(plots):

> L:=15:N:=10000:

Simulation de Monte-Carlo: atmosphère isotherme

> montecarlo:=proc(alpha,kmax)
global L,N;
local n,i,k,nrj,compt_haut,compt_bas,j,m;
n[1]:=N:
for i from 2 to L do
n[i]:=0:
od:
nrj:=N:
m:=[seq(n[k],k=1..L)]:
for k from 1 to kmax do
compt_haut:=0:
compt_bas:=0:
if n[1]<>0 then
for j from 1 to n[1] do
if evalf(exp(-alpha))>rand()*1e-12 then compt_haut:=compt_haut+1:
fi:
od:
n[1]:=n[1]-compt_haut:
n[2]:=n[2]+compt_haut:
fi:
for i from 2 to L-1 do
compt_haut:=0:
compt_bas:=0:
if n[i]<>0 then
for j from 1 to n[i] do
if rand()*1e-12>0.5 then compt_bas:=compt_bas+1:
elif evalf(exp(-alpha))>rand()*1e-12 then compt_haut:=compt_haut+1:
fi:
od:
n[i]:=n[i]-compt_haut-compt_bas:
n[i+1]:=n[i+1]+compt_haut:
n[i-1]:=n[i-1]+compt_bas:
fi:
od:
if n[L]<>0 then
compt_bas:=0:
compt_haut:=0:
for j from 1 to n[L] do
if rand()*1e-12>0.5 then compt_bas:=compt_bas+1:
elif evalf(exp(-alpha))>rand()*1e-12 then compt_haut:=compt_haut+1:
fi:
od:
n[L]:=n[L]-compt_haut-compt_bas:
n[L-1]:=n[L-1]+compt_bas:
n[1]:=n[1]+compt_haut:
fi:
nrj:=nrj,add(j*n[j],j=1..L):
m:=m,[seq(n[j],j=1..L)];
od:
[[m],[nrj]];
end:

> moyenne:=proc(tab,rang)
global L,N;
local i,j,k,out,cpt;
for i from 1 to L do
out[i]:=0;
od:
for i from rang to nops(tab) do
for j from 1 to L do
out[j]:=out[j]+tab[i][j]:
od:
od:
[seq(evalf(out[k]/(nops(tab)-rang+1)),k=1..L)];
end:

La fonction qui suit permet d'obtenir la température (en fait le paramètre ) qui correspond à un état d'équilibre d'énergie donnée.

> trouve:=(x,L)->fsolve(U(y,L)=x,y):

Dans cette simulation, le paramètre est fixé à la valeur 0.2, ce qui correspond à une échelle de hauteur de pression , inférieure à la taille du système .

> res:=montecarlo(.2,300):

> popu:=moyenne(res[1],100);

> add(popu[i],i=1..L);

On vérifie que le nombre total de molécules reste constant au cours de la simulation: c'est l'avantage des conditions aux limites, qui connectent la case et la case . Sans cette précaution, le système se viderait par le haut au fur et à mesure que l'on établit de nouvelles configurations...à moins de pouvoir considérer un système de taille infinie, ce qui n'est bien sûr guère possible.

On se propose de vérifier que la répartition obtenue vérifie bien la distribution de Boltzmann (le contraire serait surprenant vu que l'échantillonage des configurations se fait au moyen de cette densité de probabilité !).

> plots[pointplot]([seq([i,popu[i]],i=1..L)]);

> plots[pointplot]([seq([i,log10(popu[i])],i=1..L)]);

Il semble que la distribution de Boltzmann soit vérifiée sur les 10 premières cases. On peut dire que le système simulé correspond à .

On représente maintenant l'évolution de l'énergie interne du système.

> E:=res[nops(res)]:

> plots[pointplot]([seq([i,E[i]],i=1..nops(E))]);

Ce graphe montre l'évolution de l'énergie du système à partir de la configuration initiale. On observe en fait la relaxation du système, préparé à température nulle, et mis en contact à l'instant initial avec un thermostat de température non nulle. L'énergie du système finit par atteindre une valeur stationnaire, abstraction faite des fluctuations.

Il y a là une différence majeure avec le cas d'un système isolé (ensemble microcanonique): l'énergie interne du système thermostaté à l'équilibre est libre de fluctuer alors que sa température est maintenue constante.

On calcule la valeur moyenne de cette énergie rapportée à l'effectif du système pour comparer la valeur obtenue par simulation à la valeur théorique (système de volume égal à ).

> nrj:=(evalf(add(E[i],i=100..nops(E))/(nops(E)-99)))/N;

> U(0.2,10);

La simulation donne une valeur qui diffère de 2% par rapport à la valeur théorique. Ce bon accord permet de valider la procédure écrite pour rechercher l'état d'équilibre.

Vérifions que cette valeur de l'énergie, obtenue par simulation, correspond bien à (pour L , on prend la valeur de afin de minimiser l'influence des effets de taille finie):

> trouve(4.047236634,10);

Etude théorique: fonctions d'état et équation d'état du système

L'expression des différentes fonctions d'état est obtenu à partir du calcul de la fonction de partition. Comme les molécules sont considérées comme indiscernables, on peut commencer par le calcul de la fonction de partition individuelle d'une molécule:

> L:='L':z:=Sum('exp(-alpha*k)','k=1..L');

La fonction de partition du système des molécules indiscernables est:

> N:='N':Z:=z^N/N!;

On peut maintenant calculer les différentes fonctions d'état thermodynamiques:

Cette relation s'écrit aussi, sous la forme:

.

On définit ensuite une énergie adimensionnée et rapportée à l'effectif :

.

> EE:=(alpha,L)->-Sum(-k*exp(-alpha*k),k = 1 .. L)/Sum(exp(-alpha*k),k = 1 .. L);

Lorsqu'on fait un calcul à la main, en utilisant l'approximation de Stirling pour évaluer , on arrive à la fonction suivante (notée ):

.

> U:=(alpha,L)->1/(1-exp(-alpha))-L*exp(-L*alpha)/(1-exp(-L*alpha));

On peut vérifier que le calcul effectué à la main (moyennant une approximation bien légitime) ou le calcul exact effectué par Maple donnent les mêmes résultats numériques:

> evalf(EE(1,15)),evalf(U(1,15));

L'énergie libre est obtenue directement par . On utilise une énergie libre adimensionnée pour la simulation numérique.

Calcul exact de l'énergie libre par Maple

> N:='N':

> FF:=(alpha,L)->-ln(Sum(exp(-alpha*k),k = 1 .. L)^N/N!)/(N*alpha);

Calcul de l'énergie libre adimensionnée (à la main, moyennant l'approximation de Stirling):

> F:=(alpha,L)->1+1/alpha*log((1-exp(-alpha))/(1-exp(-L*alpha)))-1/alpha+ln(N)/alpha;

On vérifie l'accord numérique des deux calculs (surtout pour valider le calcul à la main !):

> N:=10000:evalf(FF(1,15)),evalf(F(1,15));N:='N':

Pour le calcul de l'entropie adimensionnée, on l'obtient simplement en écrivant . On peut aussi la calculer à partir de la fonction de partition: , puis l'adimensionner en divisant par .

> S:=(alpha,L)->alpha*(exp(-alpha)/(1-exp(-alpha))-L*exp(-L*alpha)/(1-exp(-L*alpha)))-log((1-exp(-alpha))/(1-exp(-L*alpha)))+1-ln(N);

On peut aussi établir l'expression de la pression canonique du système, sachant que le volume du système est sa taille L :

. Cette quantité est homogène à une force. On peut donc l'adimensionner en la divisant par le poids total du système Nmg . Le calcul de p conduit au résultat suivant:

. Il s'agit de l'équation d'état du système reliant les variables d'état p , L et (lié à la température).

Il est intéressant de calculer la limite haute température, c'est-à-dire la limite <<1 (ce qui équivaut à mgLa<< ):
, soit en réintroduisant la température, , ce qui s'écrit encore sous la forme . On retrouve la loi des gaz parfaits.

> pression:=(alpha,L)->1/(1-exp(-alpha*L));

Et maintenant, les différents graphes:

> GU:=plot(U(alpha,15),alpha=0..2):display(GU);

L'énergie du système est une fonction décroissante de , donc une fonction croissante de la température.

> plot(S(alpha,15),alpha=0..2);

L'entropie est une fonction décroissante de , et par conséquent une fonction croissante de la température: plus la température est importante et plus le nombre de configurations accessibles au système augmente.

> plot(FF(alpha,15),alpha=0..2);

L'énergie libre est une fonction décroissante de , donc une fonction croissante de la température.

On représente ci-dessous la pression en fonction de (en trait gras) ainsi que la loi du gaz parfait (en trait fin).

> Gp:=plot(pression(alpha,15),alpha=0..1,thickness=2):
GGP:=plot(1/(15*alpha),alpha=0..1,view=[0..1,0..3.5]):
display(Gp,GGP);

A basse température, la pression réduite du système tend vers 1, ce qui signifie que la pression du système tend vers Nmg qui représente le poids total du système.

Séquences de simulations de température croissante

Afin de valider les simulations, on veut retrouver par la simulation la valeur de l'énergie interne du système et la comparer aux résultats théoriques. Cela permet aussi de visualiser l'influence de la taille finie du système sur la simulation.

> sim:=proc(alpha1,alpha2,pas,kmax,borne)
global L,N;
local alpha,nrj,E,res,cpt,k;
cpt:=1:
for alpha from alpha1 to alpha2 by pas do
res:=montecarlo(alpha,kmax):
E:=res[nops(res)]:
nrj[cpt]:=(evalf(add(E[i],i=borne..nops(E))/(nops(E)-borne+1)))/N;
cpt:=cpt+1;
od:
[seq(nrj[k],k=1..cpt-1)];
end:

> EE:=sim(0.2,1.8,0.2,300,60):

> EE;L;N;

> EE:=[4.058504545, 2.914687603, 2.206392149, 1.819161570, 1.584344215, 1.434422314, 1.328484711, 1.252577273, 1.197214050];

> GE:=plots[pointplot]([seq([0.2+(i-1)*0.2,EE[i]],i=1..nops(EE))],symbol=cross):display(GU,GE);

L'accord entre les valeurs obtenues par simulation et les résultats théoriques sont compatibles: ceci valide la procédure utilisée.

Pour petit (T faible), l'échelle de hauteur caractéristique du système () n'est plus négligeable devant la taille du système et les effets de bords sont importants. Ainsi, avec , on a . Pour un système de taille , on trouve une distribution exponentielle sur .

Evolution monotherme du système

On considère le système étudié dans son état d'équilibre à la température imposée par le thermostat. On imagine une transformation qui consiste à doubler le volume du système (on multiplie par deux) sans fournir de travail au système (on retire par exemple une paroi escamotable): il s'agit d'une détente face au vide...mais la présence du thermostat fait qu'elle est monotherme. Grâce à la procédure montecarlo , on recherche le nouvel état d'équilibre du système. Comment évoluent les différentes fonctions thermodynamiques au cours de cette transformation ?

On écrit une procédure qui réalise les opérations suivantes:

- recherche de l'état d'équilibre initial correspondant à et L

- une fois l'état d'équilibre initial atteint, on change L en 2L en maintenant constant.

- on recherche le nouvel état d'équilibre

> relaxation:=proc(alpha,kmax,A)
global L,N;
local r1,popu,n,i,k,nrj,compt_haut,compt_bas,j,m;
r1:=montecarlo(alpha,kmax);
for i from 1 to L do
n[i]:=r1[1][kmax][i]:
od:
for i from L+1 to 2*L do
n[i]:=0:
od:
nrj:=r1[2][kmax]:
L:=A*L:
m:=[seq(n[k],k=1..L)]:
for k from 1 to kmax do
compt_haut:=0:
compt_bas:=0:
if n[1]<>0 then
for j from 1 to n[1] do
if evalf(exp(-alpha))>rand()*1e-12 then compt_haut:=compt_haut+1:
fi:
od:
n[1]:=n[1]-compt_haut:
n[2]:=n[2]+compt_haut:
fi:
for i from 2 to L-1 do
compt_haut:=0:
compt_bas:=0:
if n[i]<>0 then
for j from 1 to n[i] do
if rand()*1e-12>0.5 then compt_bas:=compt_bas+1:
elif evalf(exp(-alpha))>rand()*1e-12 then compt_haut:=compt_haut+1:
fi:
od:
n[i]:=n[i]-compt_haut-compt_bas:
n[i+1]:=n[i+1]+compt_haut:
n[i-1]:=n[i-1]+compt_bas:
fi:
od:
if n[L]<>0 then
compt_bas:=0:
compt_haut:=0:
for j from 1 to n[L] do
if rand()*1e-12>0.5 then compt_bas:=compt_bas+1:
elif evalf(exp(-alpha))>rand()*1e-12 then compt_haut:=compt_haut+1:
fi:
od:
n[L]:=n[L]-compt_haut-compt_bas:
n[L-1]:=n[L-1]+compt_bas:
n[1]:=n[1]+compt_haut:
fi:
nrj:=nrj,add(j*n[j],j=1..L):
m:=m,[seq(n[j],j=1..L)];
od:
m:=op(r1[1]),m;
nrj:=op(r1[2]),nrj;
[[m],[nrj]];
end:
trouveL:=proc(dist)
local i;
i:=nops(dist);
while dist[i]=0 do
i:=i-1;
od:
i;
end:

> L:=15:res:=relaxation(0.4,200,2):

> E:=res[nops(res)]:

> plots[pointplot]([seq([i,E[i]],i=1..nops(E))]);

Ce graphe permet de montrer l'évolution de l'énergie interne du système au cours de cette transformation. On constate que la variation d'énergie interne est positive.

Pour suivre l'évolution de l'énergie libre et de l'entropie, on commence par rechercher le volume occupé par le système pour chaque configuration générée par l'algorithme. Ce volume est une variable fluctuante, qui diffère du volume total L offert au système, que l'on définit comme étant le "numéro" du niveau de plus haute énergie occupé par au moins une molécule.

> LL:=[seq(trouveL(res[1][i]),i=1..nops(res[1]))]:

Pour chacune de ces configurations intermédiaires, caractérisée par son volume, on calcule, à l'aide de l'expression théorique, la valeur du potentiel thermodynamique et de l'entropie du système, qui sont aussi des grandeurs fluctuantes. Pour ce calcul, le paramètre reste fixé à sa valeur 0,4 car la température n'est pas une variable fluctuante: elle est imposée par le milieu extérieur au système.

Commençons par visualiser l'évolution du potentiel thermodynamique énergie libre:

> f:=[seq(F(0.4,LL[i]),i=1..nops(LL))]:

> plots[pointplot]([seq([i,f[i]],i=2..nops(f))],view=[0..nops(LL),18.7..18.8]);

On constate qu'au cours de l'évolution monotherme (spontanée) du système, le potentiel thermodynamique énergie libre du système décroît. . Il est minimum à l'équilibre thermodynamique compte tenu des contraintes qui s'imposent au système . Au cours de la transformation, la variation de l'énergie libre du système est négative.

Ceci illustre que le second principe de la thermodynamique est un principe d'évolution: en modifiant une contrainte initiale imposée au système, on force celui-ci à relaxer vers un nouvel état d'équilibre, où le potentiel thermodynamique trouve une nouvelle valeur minimale, compatible avec les contraintes que le milieu extérieur impose au système.

Comment évolue l'entropie du système ?

> s:=[seq(S(0.4,LL[i]),i=1..nops(LL))]:

> plots[pointplot]([seq([i,s[i]],i=2..nops(s))],view=[0..nops(LL),-6.4..-6.2]);

Au cours de la transformation, l'entropie du système augmente. Qualitativement, on peut comprendre qu'en augmentant l'espace disponible, on augmente aussi le nombre de configurations microscopiques possibles.
Récapitulons, l'énergie interne et l'entropie du système augmentent au cours de cette évolution. La décroissance du potentiel thermodynamique énergie libre résulte de la compétition entre les deux termes énergétiques et entropiques. L'augmentation de l'énergie interne a tendance à faire croître l'énergie libre...mais c'est sans compter sur l'entropie: la transformation se traduit par une augmentation considérable du nombre de configurations microscopiques possibles, l'entropie augmente. C'est l'augmentation importante du terme entropique qui fait finalement pencher la balance dans le sens d'une décroissance de l'énergie libre.

Références

A propos des simulations de Monte-Carlo

D.P. Landau et R. Alben, "Monte-Carlo calculations as an aid in teaching statistical mechanics", Am. J. Phys., 41 , 394-400 (1973).

K. Binder et D.W. Heermann, "Monte-Carlo Simulation in Statistical Physics", Springer-Verlag.

D. Chandler, "Introduction to Modern Statistical Mechanics", Oxford University Press.

Ouvrages de physique statistique

B. Diu, C. Guthmann, D. Lederer et B. Roulet, "Eléments de physique statistique", Hermann (1989).

M. Le Bellac, "Des phénomènes critiques aux champs de jauge", Inter-Editions, éditions du CNRS.

S. Vauclair, "Eléments de physique statistique", Inter-Editions, éditions du CNRS (vraiment, il s'agit d'un livre tout à fait remarquable: concis, précis, simple, clair avec d'évidentes qualités pédagogiques !).

M. Moreau, Cours de physique stochastique du DEA de Physique des Liquides (Université Paris VI), 1997.