Dynamique moléculaire : modèle des sphères molles à deux dimensions, équation d'état.

1. Introduction

Ce document fait suite à Dynamique moléculaire : modèle des sphères molles à deux dimensions.

On s'intéresse ici à l'obtention de l'équation d'état du système de sphères molles. On parle de sphères molles car le potentiel répulsif utilisé pour modéliser le contact des sphères présente une pente finie. Cependant, cette pente est très grande, et le comportement attendu est très voisin de celui d'un système de sphères dures.

L'équation d'état est obtenue à l'aide du viriel moyen et de l'énergie cinétique moyenne (moyennes temporelles) :

$$\frac{PA}{N{k}_{B}T}=1+\frac{6}{<Ec>}\epsilon <\sum _{i<j}({(d/r_{ij})}^{12}-{(d/r_{ij})}^6)>\text{(1)}$$

2. Courbes isothermes

import sys
sys.path.append("../spheresmolles2d") 
from spheresMolles2d import *
from matplotlib.pyplot import *
            

Il s'agit de faire plusieurs calculs de dynamique moléculaire, avec la même vitesse initiale (et donc la même température), en faisant varier la densité. Pour cela, on choisit de maintenir fixe le nombre de sphères et de faire varier la taille de la grille.

La fonction suivante effectue un calcul complet pour un nombre de sphères (par dimension d'espace) et une taille de domaine donnés. La fonction renvoit la densité, l'énergie cinétique moyenne par sphère (c.a.d. la température) et son écart-type, le rapport PA/Nk_BT et son écart-type.

def calcul(Nx,densite,vinit,h):
    sys = Systeme(Nx,densite)
    sys.initialiser(vinit)
    for k in range(100): # mise a l'equilibre
        sys.integration(sys.verlet,h,10)
    ec = numpy.zeros(0)
    rapport = numpy.zeros(0)
    iter = 100
    while iter > 0:
        iter -=1
        sys.integration(sys.verlet,h,10)
        sys.calculer_cinetique()
        ec = numpy.append(ec,sys.Ecinetique)
        rapport = numpy.append(rapport,sys.pression*sys.L**2/(sys.Ecinetique*Nx*Nx))
    return (numpy.mean(ec),numpy.std(ec),numpy.mean(rapport),numpy.std(rapport)) 
             

Voici le calcul pour 1600 sphères et pour différentes tailles de domaine :

Nx = 40
vinit = 1.0 
h = 0.005
densite = [0.1,0.2,0.3,0.4,0.5,0.6,0.7]
rapport = numpy.zeros(0)
dev_rapport = numpy.zeros(0)
ecinetique = numpy.zeros(0)
dev_ecinetique = numpy.zeros(0)
for dens in densite:
    (ec,dec,r,dr)=calcul(Nx,dens,vinit,h)
    print(ec,dec,r,dr)
    ecinetique = numpy.append(ecinetique,ec)
    dev_ecinetique =numpy.append(dev_ecinetique, dec)
    rapport = numpy.append(rapport,r)
    dev_rapport = numpy.append(dev_rapport,dr)
numpy.savetxt("spheres40-1.txt",[densite,rapport,dev_rapport,ecinetique,dev_ecinetique])
             

On trace le rapport PA/Nk_BT en fonction de la densité, définie comme le rapport de l'aire occupée par les disques sur l'aire A du domaine. Les barres d'incertitude correspondent aux écart-types des fluctuations.

[densite,rapport,dev_rapport,ecinetique,dev_ecinetique] = numpy.loadtxt("spheres40-1.txt")
figure()
errorbar(densite,rapport,yerr=dev_rapport,fmt=None)
plot(densite,rapport,'r')
xlabel("densite")
ylabel("PA/NkT")
title("T = %f"%numpy.mean(ecinetique))

             

figA.pdf