Modèle des sphères dures à deux dimensions : équation d'état

1. Introduction

Ce document fait suite à Modèle des sphères dures à deux dimensions, où est présenté l'algorithme de Metropolis appliqué au modèle des sphères dures.

On s'intéresse ici à l'obtention de la pression en fonction de la densité, en suivant la méthode utilisée dans [1].

Il s'agit de calculer la fonction de distribution radiale g(r) pour des valeurs proches de r=d (diamètre des sphères), puis d'obtenir g(d) par extrapolation.

2. Extrapolation de la fonction de distribution radiale

import sys
sys.path.append("../spheres2d") 
from spheres2dMetropolis import *
from matplotlib.pyplot import *
            

On calcule la fonction de distribution jusqu'à r=2d :

N = 7
densite = 0.4
systeme = Systeme(N,densite)
systeme.initialiser()
systeme.boucle_metropolis(100000)
figure(figsize=(8,6))
(r,g) = systeme.calculer_distribution_radiale(3,20,10000,100)
plot(r,g,'o')
xlabel('r')
ylabel('g')
title('densite = %f'%densite)
axis([0,3,0,g.max()])
              

Figure pleine page

Pour obtenir la valeur extrapolée de g(d), on utilise une interpolation polynomiale avec l'algorithme de Neville. La fonction suivante fait ce calcul, en prenant les P premiers points :

def extrapoler(r,g,P):
    def recurrence(i,m,x):
        if m==0:
            return g[i]
        else:
            return ((x-r[i+m])*recurrence(i,m-1,x)+(r[i]-x)*recurrence(i+1,m-1,x))/(r[i]-r[i+m])
    return recurrence(0,P-1,1.0)
                

print(extrapoler(r,g,5))
--> 2.2831652360918731

3. Équation d'état

L'équation d'état est donnée par la relation suivante :

$$\frac{PA}{N{k}_{B}T}=1+\frac{\pi }{2}\frac{d^{2}N}{L^2}g\left(d\right)\text{(1)}$$

On trace donc ce rapport en fonction de la densité (aire occupée par les sphères divisée par l'aire totale).

import math
densite = [0.1,0.2,0.3,0.4,0.5]
pression = []
for k in range(len(densite)):
    systeme = Systeme(N,densite[k])
    systeme.initialiser()
    systeme.boucle_metropolis(100000)
    (r,g) = systeme.calculer_distribution_radiale(2,20,10000,100)
    pression.append(1.0+math.pi/2*systeme.N/(systeme.L**2)*extrapoler(r,g,5))
figure(figsize=(8,6))
plot(densite,pression,'o')
xlabel('densite') 
ylabel('PA/NkT') 
axis([0,0.5,0.0,max(pression)])
              

Figure pleine page