Sphères dures

1. Introduction

L'algorithme de Monte-Carlo Metropolis a été initialement introduit ([1]) pour simuler un système formé de sphères dures se déplaçant sur un plan de manière continue.

On considère ici une version simplifiée du modèle des sphères dures, dans laquelle les sphères se déplacent sur un pavage carré. Chaque cellule du pavage est occupée au maximum par une sphère.

2. Modèle des sphères dures

Le modèle des sphères dure consiste à représenter les atomes d'un fluide (monoatomique) par des sphères de rayon d. Ce modèle est utilisé en physique statistique pour obtenir une équation d'état ([2][3]).

Lorsque les sphères ne sont pas en contact, elles n'exercent aucune force l'une sur l'autre. Au contact, une force répulsive infinie apparaît, qui empêche l'interpénétration des sphères. L'énergie potentielle d'interaction entre deux sphères a donc la forme suivante :

\begin{matrix} r < d E (r) = \infty & (1) \\ r \geq d E (r) = 0 & (2) \end{matrix}

La distribution statistique des vitesses (distribution de Maxwell) est indépendante de la taille des sphères. Dans une simulation de Monte-Carlo, on peut donc se limiter à considérer les configurations du système dans l'espace et pas les vitesses des atomes : l'énergie cinétique n'est pas prise en compte. Nous avons donc repris la notation E comme utilisée dans Algorithme de Metropolis.

Puisque l'énergie potentielle est nulle, l'énergie interne du fluide de sphères dures se réduit à l'énergie cinétique des atomes (comme pour un gaz parfait) :

U = \frac{3}{2} N k_{B} T (3)

où N est le nombre de sphères.

L'équation d'état d'un fluide de sphères dures, qui relie la pression à la densité, est en revanche différente de celle du gaz parfait. On voit par exemple que si la densité est maximale (rangement compact des sphères), le système est dans un état figé, c'est-à-dire dans un état solide. C'est pour cette raison que le modèle des sphères dures est intéressant.

L'équation d'état s'obtient en utilisant la fonction de distribution radiale g(r). Cette fonction représente la probabilité pour que deux sphères (leur centre) soient distantes de r. Plus précisément, si l'on considère qu'une sphère se trouve centrée à l'origine, le nombre moyen de sphères centrées entre les rayons r et r+dr est

d N = 4 π \frac{N}{V} r^{2} g (r) d r (4)

où V est le volume du système. La fonction de distribution radiale est non seulement un outil de calcul théorique, mais aussi une fonction accessible expérimentalement par diffraction de rayonnement (rayons X, neutrons). Elle apporte une information sur la structure locale du fluide (à l'échelle atomique), dont dépend son équation d'état. Pour un fluide de sphères dures, la théorie statistique ([3]) conduit à l'équation d'état suivante :

P V = N k_{B} T (1 + \frac{2 π}{3} \frac{N}{V} d^{3} g (d)) (5)

Il faut donc déterminer g(d), la valeur de la fonction de distribution radiale pour une distance égale au diamètre des sphères, c'est-à-dire lorsque les sphères sont en contact. Il s'agit d'une valeur moyenne qui peut être obtenue par une simulation de Monte-Carlo, en suivant l'algorithme de Metropolis. En introduisant n le nombre moyen de sphères par unité de volume en contact avec une sphère, l'expression (5) devient :

P V = N k_{B} T (1 + \frac{2 π}{3} d^{3} n) (6)

L'équation d'état peut être développée par rapport à la densité N/V, ce qui donne le développement du viriel :

P V = N k_{B} T (1 + B_{2} \frac{N}{V} + B_{3} (\frac{N}{V})^{2} \dots) (7)

Le calcul de n en fonction de la densité doit donc permettre de calculer les coefficients du viriel B₂,B₃,....

On se limite au modèle des sphères dures sur un plan (système bidimensionnel). Si A est l'aire du domaine carré, l'équation d'état s'écrit ([1]) :

P A = N k_{B} T (1 + \frac{π}{2} d^{2} n) (8)

où n est la densité surfacique moyenne des sphères en contact avec une sphère.

3. Algorithme

Les cellules du pavage sont carrées, avec une longueur de côté égal au diamètre d des sphères. La figure suivante représente une partie du pavage avec 4 sphères.

Figure pleine page

Les sphères se déplacent horizontalement, verticalement ou en diagonal, d'une case à la fois. Ce déplacement élémentaire respecte le principe d'ergodicité car il permet d'obtenir toutes les configurations possibles. L'algorithme de Metropolis appliqué à ce modèle consiste à générer une chaîne de configurations. Une configuration est générée à partir de la précédente en appliquant la procédure suivante :

Choix aléatoire d'une sphère et choix aléatoire d'un des 8 déplacements possibles.
Si le déplacement se fait vers une cellule vide, il est accepté.
Si le déplacement se fait vers une cellule déjà occupé, il est refusé. La nouvelle configuration est alors identique à la précédente.

Sur la figure, on voit les déplacements permis pour une sphère. Le déplacement vers la gauche est interdit à cause de la présence d'une autre sphère. Si ce déplacement interdit est sélectionné, le nouvel état est identique au précédent.

Ce modèle correspond à un potentiel répulsif infini entre les sphères. Lorsque les sphères ne sont pas en contact leur énergie d'interaction est nulle. Lorsqu'elle se chevauchent, l'énergie d'interaction est infinie. Le facteur de Boltzmann e^-βΔE peut donc prendre seulement deux valeurs : 1 ou 0. La première valeur correspond à l'acceptation du déplacement, le second à son refus. Il s'en suit que la température n'intervient pas dans ce modèle. L'énergie potentielle du système est constamment nulle.

Le seul paramètre est la densité. Si N est le nombre de sphères et N_c le nombre de cellules, on définit la densité par le rapport :

ρ = \frac{N}{N_{c}} (9)

Sa valeur maximale est 1. Dans ce cas, le système est figé puisque toutes les cases sont occupées.

On appliquera une condition limite périodique sur les bords. Comme condition initiale, on prendra un rangement compact des sphères, dans N cellules consécutives.

La chaîne de configurations obtenue permet de faire des calculs statistiques sur différentes grandeurs. On calculera la distance quadratique moyenne parcourue par les sphères et le nombre moyen de proches voisins.

4. Implémentation

On définit une classe dont le constructeur prend en argument le nombre de cellules dans une dimension d'espace et le nombre de sphères. Le constructeur crée les tableaux et initialise le système en plaçant les sphères dans l'ordre des cellules, de manière compacte.

spheres2d.py

import numpy
import random
import math

class Spheres2D:
    def __init__(self,nx,N):
        self.nx = nx
        self.Nc = nx*nx
        self.N = N
        self.rho = self.N*1.0/self.Nc
        self.pavage = numpy.zeros((nx,nx),dtype=numpy.uint8)
        self.spheres = numpy.zeros((N,2))
        self.dx = numpy.zeros(N)
        self.dy = numpy.zeros(N)
        self.dr2 = numpy.zeros(N)
        self.nvoisins = numpy.zeros(N)
        self.dir = [(1,0),(-1,0),(0,1),(0,-1),(-1,1),(1,1),(-1,-1),(1,-1)]
        i = j = 0
        for k in range(N):
            self.spheres[k][0] = i
            self.spheres[k][1] = j
            self.pavage[j][i] = 1
            i += 1
            if i==self.nx:
                i = 0
                j += 1

La fonction suivante implémente les conditions limites périodiques. Elle prend deux indices en argument et renvoit les indices effectifs pour le tableau du pavage.

    def indices(self,i,j):
        ii = i
        jj = j
        if ii == self.nx:
            ii = 0
        elif ii == -1:
            ii = self.nx-1
        if jj == self.nx:
            jj = 0
        elif jj == -1:
            jj = self.nx-1
        return (ii,jj)

La fonction suivante calcule et mémorise le nombre de proches voisins d'une sphère. Les proches voisins sont les sphères situées sur les 8 cases voisines.

    def voisins(self,s):
        nv = 0
        i = self.spheres[s][0]
        j = self.spheres[s][1]
        (ii,jj) = self.indices(i+1,j)
        nv += self.pavage[jj][ii]
        (ii,jj) = self.indices(i-1,j)
        nv += self.pavage[jj][ii]
        (ii,jj) = self.indices(i,j+1)
        nv += self.pavage[jj][ii]
        (ii,jj) = self.indices(i,j-1)
        nv += self.pavage[jj][ii]
        (ii,jj) = self.indices(i-1,j-1)
        nv += self.pavage[jj][ii]
        (ii,jj) = self.indices(i-1,j+1)
        nv += self.pavage[jj][ii]
        (ii,jj) = self.indices(i+1,j-1)
        nv += self.pavage[jj][ii]
        (ii,jj) = self.indices(i+1,j+1)
        nv += self.pavage[jj][ii]
        self.nvoisins[s] = nv

La fonction suivante effectue le pas élémentaire de l'algorithme de Metropolis. La distance au carrée parcourue depuis le début par la sphère déplacée est actualisée.

    def metropolis(self):
        s = random.randint(0,self.N-1)
        p = self.spheres[s]
        dir = self.dir[random.randint(0,7)]
        (ii,jj) = self.indices(p[0]+dir[0],p[1]+dir[1])
        if self.pavage[jj][ii]==0:
            self.pavage[jj][ii] = 1
            self.pavage[p[1]][p[0]] = 0
            self.spheres[s][0] = ii
            self.spheres[s][1] = jj
            self.dx[s] += dir[0]
            self.dy[s] += dir[1]
            self.dr2[s] = self.dx[s]*self.dx[s]+self.dy[s]*self.dy[s]

La fonction suivante effectue une boucle comportant n itérations. Une itération est constituée de N pas élémentaires. Ainsi une itération constitue en moyenne le déplacement (ou sa tentative) de toutes les sphères. Les deux listes renvoyées contiennent le déplacement quadratique moyen et le nombre de proches voisins moyen (une valeur pour chaque itération).

    def boucle(self,n):
        dr2 = numpy.zeros(n)
        nv = numpy.zeros(n)
        self.dx = numpy.zeros(self.N)
        self.dy = numpy.zeros(self.N)
        for k in range(n):
            dr2[k] = numpy.mean(self.dr2)
            for s in range(self.N):
                self.voisins(s)
            nv[k] = numpy.mean(self.nvoisins)
            for s in range(self.N):
                self.metropolis()
        return (dr2,nv)

5. Résultats

import numpy
import random
import math
from matplotlib.pyplot import *
import spheres2d

On choisit pour commencer 100 sphères sur un pavage de 400 cellules. On effectue 300 itérations.

nx = 20
N=100
spheres = spheres2d.Spheres2D(nx,N)
(dr2,nv)=spheres.boucle(500)
figure(figsize=(6,6))
plot(dr2)
xlabel('iteration')
ylabel('deltaR^2')
grid()

figA.pdf

Le déplacement quadratique moyen semble proportionnel au nombre d'itérations. Il s'agit d'une loi d'échelle caractéristique d'une marche au hasard. Voyons le nombre de proches voisins :

figure(figsize=(6,6))
plot(nv)
xlabel('iteration')
ylabel('n voisins')
grid()

figB.pdf

Lorsque l'équilibre est atteint, le nombre de voisins fluctue autour de 2.

Une simulation complète consiste à faire varier la densité. Pour chaque valeur de la densité, on procède à 300 itérations pour amener le système à l'équilibre, puis 300 itérations pour calculer les nombres de proches voisins.

rho = numpy.array([0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,0.95,0.98])
nvoisins = numpy.zeros(rho.size)
dnvoisins = numpy.zeros(rho.size)
Nc = nx*nx
for k in range(rho.size):
    spheres = spheres2d.Spheres2D(nx,int(rho[k]*Nc))
    (dr2,nv) = spheres.boucle(300)
    (dr2,nv) = spheres.boucle(300)
    nvoisins[k] = numpy.mean(nv)
    dnvoisins[k] = numpy.std(nv)
figure(figsize=(8,6))
errorbar(rho,nvoisins,yerr=dnvoisins,fmt=None)
xlabel('densite')
ylabel('n voisins')
grid()

figC.pdf

Ce modèle conduit à une variation linéaire du nombre de proches voisins en fonction de la densité. En terme de développement du viriel, cela signifie que l'on obtient le terme B₂ seulement, terme du premier ordre par rapport à la densité.

Sphères dures sur un pavage bidimensionnel

1. Introduction

2. Modèle des sphères dures

3. Algorithme

4. Implémentation

5. Résultats

Références