Table des matières Python

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 :

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 :

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 :

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 nr(r), on effectue un échantillonnage régulier à p points de l'intervalle [-a/2,a/2]. On pose :

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 [rk,rk+1]=[k,k+1], on compte le nombre de particules gk 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 à rk=k est :

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

La fonction suivante calcule le tableau des valeurs gk, hk et nk. 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()
             
plotA.svgFigure 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()
             
plotA2.svgFigure 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()
             
plotB.svgFigure pleine page

La densité fluctue beaucoup pour r faible puis converge vers .

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()
             
plotC.svgFigure pleine page

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

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 1023).

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()
             
plotD.svgFigure 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()
             
plotE.svgFigure 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 nk, 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 π.

Creative Commons LicenseTextes et figures sont mis à disposition sous contrat Creative Commons.