Densité volumique

1. Introduction

Ce document introduit la notion de densité volumique. Une simulation de Monte-Carlo est effectuée pour montrer l'influence de la taille du volume utilisé pour calculer la densité. Voir aussi la simulation javascript qui reprend ces calculs.

2. Définition du modèle

On considère N particules dans un volume cubique de côté a. Ces particules peuvent être des charges dans un milieu chargé, des molécules particulières dans une solution ou dans un gaz, etc. Le nombre N est très grand. La densité volumique moyenne des particules sur le volume total est :

$$n_{\infty }=\frac{N}{a^3}\text{(1)}$$

Les particules sont réparties aléatoirement dans le cube, avec une densité de probabilité uniforme.

Soit O le centre du cube. On note h(r) le nombre de particules contenues dans la sphère de centre O et de rayon r. La fonction h est croissante.

On définit alors la fonction suivante :

$$n\left(r\right)=\frac{h\left(r\right)}{\frac{4}{3}\pi r^3}\text{(2)}$$

qui représente le nombre de particules par unité de volume dans la sphère de rayon r. Cette fonction doit vraisemblablement tendre vers la densité moyenne :

$$\underset{r\to \infty }{lim}n\left(r\right)=n_{\infty }\text{(3)}$$

Nous allons faire une simulation pour mettre en évidence la vitesse de convergence vers la densité moyenne.

3. Simulation

Les coordonnées x,y,z des N particules sont obtenues avec un générateur de nombres pseudo-aléatoires dans l'intervalle [-a/2,a/2].

Le rayon de la plus grande sphère incluse dans le cube est a/2. Pour obtenir la fonction n_r(r), on effectue un échantillonnage régulier à p points de l'intervalle [-a/2,a/2]. On pose :

$$r_k=k\Delta r\text{(4)}$$

avec k un entier variant de 0 à p-1. Pour simplifier les calculs, on pose Δr=1. Le nombre de points p est donc la valeur entière de a/2.

Pour chaque intervalle [r_k,r_k+1]=[k,k+1], on compte le nombre de particules g_k dont la distance r au centre se trouve dans l'intervalle. Il s'agit en fait du nombre de particules comprises entre les sphères de rayon k et k+1.

Le nombre de particules dont la distance est inférieure à r_k=k est :

$$h_k=\sum _{j=0}^k{g}_j\text{(5)}$$

On obtient enfin la valeur de n(r) pour la distance r_k par :

$$n_k=\frac{h_k}{\frac{4}{3}\pi r_k^3}\text{(6)}$$

La fonction suivante calcule le tableau des valeurs g_k, h_k et n_k. La taille du cube et le nombre de particules sont données en argument.

import math
import numpy
import random

def densite(a,N):
    ninf = N/math.pow(a,3)
    p = int(math.floor(a/2))
    g = numpy.zeros(p)
    h = numpy.zeros(p)
    n = numpy.zeros(p)
    random.seed()
    for i in range(N):
        x = random.uniform(-a/2,a/2)
        y = random.uniform(-a/2,a/2)
        z = random.uniform(-a/2,a/2)
        r = math.sqrt(x*x+y*y+z*z)
        k = int(math.floor(r))
        if k<p:
            g[k] += 1
    for k in range(1,p):
        h[k] = h[k-1]+g[k]
    for k in range(1,p):
        n[k] = h[k]/(4.0/3*math.pi*math.pow(k,3))
    return[ninf,g,h,n]
             

Un exemple :

a=1000.0
N = int(1e5)
[ninf,g,h,n] = densite(a,N)
             

La densité moyenne est :

print(ninf)
--> 0.0001

Tracé du nombre de particules comprises entre r et r+Δr:

from matplotlib.pyplot import *             

figure(figsize=(8,6))
plot(g)
xlabel("r")
ylabel("g")
grid()
             

Figure pleine page

Tracé du nombre de particules dont le rayon est inférieur à r:

figure(figsize=(8,6))
plot(h)
xlabel("r")
ylabel("h")
grid()
             

Figure pleine page

Tracé de la densité dans la sphère de rayon r :

figure(figsize=(8,6))
plot(n)
xlabel("r")
ylabel("n")
axis([0,500,0,ninf*5])
grid()
             

Figure pleine page

La densité fluctue beaucoup pour r faible puis converge vers $n_{\infty }$.

Il est intéressant de faire plusieurs simulations et de tracer la densité sur la même figure :

figure(figsize=(8,6))
for i in range(5):
    [ninf,g,h,n] = densite(a,N)
    plot(n)
xlabel("r")
ylabel("n")
axis([0,500,0,ninf*5])
grid()
             

Figure pleine page

On constate que la valeur limite est pratiquement atteinte vers r=300. Cette figure permet de définir deux échelles de distance :

• L'échelle microscopique (r<300) où la densité fluctue avec la taille de la sphère considérée et d'un point à l'autre du milieu.
• L'échelle mésoscopique, petite par rapport à l'échelle macroscopique, mais sur laquelle la densité de fluctue pas. À cette échelle, les grandeurs thermodynamiques peuvent être définies.

Dans le cas présent, la sphère minimale qui définit l'échelle mésoscopique a un rayon r=300. Elle contient environ 10000 particules, ce qui est grand mais infime comparé aux nombres de particules dans un volume macroscopique (de l'ordre de 10²³).

Voyons les résultats pour une densité plus grande :

a=1000.0
N = int(1e6)
[ninf,g,h,n] = densite(a,N)
figure(figsize=(8,6))
for i in range(5):
    [ninf,g,h,n] = densite(a,N)
    plot(n)
xlabel("r")
ylabel("n")
axis([0,500,0,ninf*5])
grid()
             

Figure pleine page

La convergence est plus rapide, ce qui signifie que l'échelle mésoscopique est plus petite. Voyons à l'inverse le cas d'un milieu beaucoup moins dense :

a=1000.0
N = int(1e3)
[ninf,g,h,n] = densite(a,N)
figure(figsize=(8,6))
for i in range(5):
    [ninf,g,h,n] = densite(a,N)
    plot(n)
xlabel("r")
ylabel("n")
axis([0,500,0,ninf*5])
grid()
             

Figure pleine page

La convergence est beaucoup plus lente. On remarque aussi que pour chaque simulation la densité est nulle en dessous d'un rayon de l'ordre de 50 à 100.

Remarque : si la division par π est omise dans le calcul de n_k, on obtient une suite qui converge vers π multipliée par la densité moyenne. Il s'agit d'une méthode de Monte-Carlo pour le calcul de π.