@vdb: la Physique en Classe Prépa

@vdb / Ressources / Thèmes de physique / L'atmosphère isotherme

Contrat Creative Commons 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 [Maple Math] molécules en contact d'un thermostat à la température [Maple Math]. Le système admet un nombre fini ([Maple Math]) de niveaux d'énergie discrets. L'énergie du niveau [Maple Math] ( [Maple Math] et [Maple Math] ) est donnée par: [Maple Math][Maple Math] est la masse d'une molécule, [Maple Math] l'accélération de la pesanteur et [Maple Math] 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 [Maple Math] tel que [Maple Math].

On simule ainsi une atmosphère isotherme de taille finie découpée en [Maple Math] 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 [Maple Math] est située au-dessus de la case [Maple Math] ).

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 [Maple Math] 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 [Maple Math] une étape de cette chaîne et par [Maple Math] la configuration correspondante.

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

[Maple Math]

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

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:

[Maple Math]

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

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

[Maple Math] avec [Maple Math], d'où l'on déduit: [Maple Math] avec [Maple Math].

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

[Maple Math] si [Maple Math] et [Maple Math] si [Maple Math]. 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 [Maple Math] qui suit permet d'obtenir la température (en fait le paramètre [Maple Math]) 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 [Maple Math] est fixé à la valeur 0.2, ce qui correspond à une échelle de hauteur de pression [Maple Math], inférieure à la taille du système [Maple Math].

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

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

[Maple Math]
[Maple Math]
[Maple Math]

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

[Maple Math]

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 [Maple Math] et la case [Maple Math]. 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 à [Maple Math].

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 à [Maple Math]).

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

[Maple Math]

> U(0.2,10);

[Maple Math]

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 à [Maple Math] (pour L , on prend la valeur de [Maple Math] afin de minimiser l'influence des effets de taille finie):

> trouve(4.047236634,10);

[Maple Math]

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 [Maple Math] individuelle d'une molécule:

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

[Maple Math]

La fonction de partition du système des [Maple Math] molécules indiscernables est:

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

[Maple Math]

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

[Maple Math]

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

[Maple Math].

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

[Maple Math].

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

[Maple Math]

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

[Maple Math].

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

[Maple Math]

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));

[Maple Math]

L'énergie libre est obtenue directement par [Maple Math]. 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);

[Maple Math]

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;

[Maple Math]

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':

[Maple Math]

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

> 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);

[Maple Math]

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

[Maple Math]. 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:

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

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

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

[Maple Math]

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 [Maple Math], donc une fonction croissante de la température.

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

L'entropie est une fonction décroissante de [Maple Math], 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 [Maple Math], donc une fonction croissante de la température.

On représente ci-dessous la pression en fonction de [Maple Math] (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;

[Maple Math]
[Maple Math]

[Maple Math]
[Maple Math]

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

[Maple Math]
[Maple Math]

> 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 [Maple Math] petit (T faible), l'échelle de hauteur caractéristique du système ([Maple Math]) n'est plus négligeable devant la taille du système et les effets de bords sont importants. Ainsi, avec [Maple Math], on a [Maple Math]. Pour un système de taille [Maple Math], on trouve une distribution exponentielle sur [Maple Math].

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 [Maple Math] 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 à [Maple Math] et L

- une fois l'état d'équilibre initial atteint, on change L en 2L en maintenant [Maple Math] 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 [Maple Math] et de l'entropie du système, qui sont aussi des grandeurs fluctuantes. Pour ce calcul, le paramètre [Maple Math] 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. [Maple Math]. 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 [Maple Math] 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.