La diffusion des particules

Qu'est-ce que la diffusion des particules

La loi de Fick

La diffusion des particules est un phénomène courant en physique, en chimie et en biologie. Il est encore plus courant que cela, vous pouvez le constater tous les matins. Vous versez du café, du lait ou du thé dans votre bol puis vous y jetez un ou plusieurs morceaux de sucre. Le sucre se dissout et les molécules de sucre se dispersent dans le breuvage. Oui, mais comment?

Il y a deux possibilités:

le liquide est au repos, c'est à dire qu'il n'existe aucun mouvement de convection. La convection, souvent d'origine thermique mais pas seulement, est un mouvement global du liquide. En l'absence de convection, les molécules de sucre se dispersent par diffusion. Elles bougent, au hasard de leur énergie cinétique et des collisions avec les molécules constitutives du liquide (de l'eau essentiellement) et les autres molécules de sucre (éventuellement, cela dépend de la concentration). Ce mode de dispersion est très lent. Pour un bol de 10 cm de diamètre, un simple calcul montre que la diffusion complète du sucre dans le bol prendrait environ 33 semaines (ce temps est égal à L²/D, L étant le diamètre du bol et D le coefficient de diffusion, valant 0,52*10^-9, le tout en USI évidemment). 33 semaines !! Vous allez être en retard pour le cours de physique!
la diffusion étant trop lente, on ajoute un mouvement de convection. Et dans notre cas, le plus souvent, il est créé par une petite cuillère ! A laquelle s'ajoute très souvent la convection thermique dans un café chaud... C'est un réflexe de remuer son café avec une petite cuillère, maintenant vous savez pourquoi!

Le phénomène de diffusion des particules fut étudié par un physiologiste allemand Adolph Fick, qui proposa en 1855, après beaucoup de mesures, une loi empirique, qui relie linéairement le courant volumique de particules (le nombre de particules qui traversent une unité de volume) et le gradient de concentration (nombre de particules par unité de volume). C'est la loi de Fick, que l'on écrit \( \vec{J}(\vec{r},t) = -D\vec{grad}(n(\vec{r},t))\). Je suppose que les lecteurs connaissent la signification du gradient. A noter le signe - affectantle gradient, qui traduit le fait que les molécules diffusent vers les zones de faible concentration.

Le coefficient D est le coefficient de diffusion. Sa valeur est très petite (de l'ordre de 10^-6 à 10^-9 voire moins), toujours positive et de dimension L²T^-1.

La loi de Fick est une approximation linéaire, qui est vérifiée dans des conditions particulières mais relativement courantes. En particulier, la concentration doit être faible pour pouvoir négliger les interactions entre les molécules de sucre (dans notre cas..).

L'équation de diffusion

Pour établir l'équation de diffusion, on considère un élement de volume et l'on établit l'équation de continuité sur cet élément. Pour simplifier, en restant dans le cas de la diffusion du sucre, on considère qu'aucune molécule de sucre n'est créée dans ce volume. Ayant établi cette équation de continuité, on l'injecte dans la loi de Fick et le tour est joué. On obtient une équation qui décrit la variation de concentration de nos molécules de sucre en fonction du temps et de la distance à l'origine du morceau de sucre.

Pour un volume cylindrique de hauteur dx et de base dS, l'équation de continuité s'écrit, avec nos hypothèses \( \frac{\partial n}{\partial t} = -\frac{\partial J(x)}{\partial x}\). En reportant cette équation dans la loi de Fick, j'obtiens donc \( \frac{\partial n}{\partial t} = D \Delta n(r,t) \) avec \( \Delta n(r,t) \) le laplacien de la fonction concentration n(r,t).

Plusieurs remarques concernant cette équation:

le terme de gauche est du premier ordre et le terme de droite du deuxième ordre (un laplacien...) : cette équation n'est donc pas invariante au renversement du temps (t ou -t). La diffusion est un phénomène irréversible par construction.
En régime stationnaire, c'est à dire quand \( \frac{\partial n}{\partial t} = 0 \), on retrouve l'équation \( \Delta n(r,t) = 0 \) également connue sous le nom d'équation de Laplace, que nous avons étudié ici.

La recherche de la solution pour cette équation dans le cas d'un régime quelconque passe par une analyse physique du phénomène. Comme nous l'avons vu plus haut, la diffusion d'une particule est essentiellement caractérisée par les collisions de la particule avec les molécules environnantes. On peut considérer qu'entre deux collisions, la vitesse de la particule est constante. Cela ne vous rappelle rien ? Allez voir ici... La probabilité qu'une particule se trouve, à un instant t, à une distance d de l'origine est donnée par une loi gaussienne. La variance de cette distribution gaussienne est proportionnelle à la racine carrée du temps.

Je vous propose de vérifier ces différents points à l'aide d'une simulation de la diffusion dans un plan.

Comment simuler la diffusion de particules

Pour simuler la diffusion d'une population de particules, on utilisera un modèle qui reproduit au mieux leur mouvement fait de trajectoires rectilignes uniformes (la vitesse des particlues est supposée constante) et de collisions.Il s'agit bien sur du modèle de la marche au hasard, déjà vu lors de notre simulation du mouvement brownien. J'ai choisi de construire une simulation 2D, c'est à dire que notre population de particules évoluera dans un plan.

A l'origine des temps, toutes les particules sont positionnées à l'origine du repère. Pour chaque génération, espacée d'un temps τ unitaire, la simulation calculera les cordonnées (x,y)_j de chaque particule, j étant le numéro de la génération. Après le calcul des coordonnées de la génération j, la simulation affichera l'ensemble des particules dans le repère. Il sera ainsi possible de suivre visuellement la diffusion génération après génération.

Reste à déterminer comment calculer les coordonnées (x,y)_j à partir des coordonnées (x,y)_j-1 de telle sorte qu'elles soient distribuées selon une loi uniforme. Dans mon programme, j'ai construit deux lois d'évolution qui répondent au critère de distribution attendu:

une marche 2D à pas constant. Selon cette loi, on calcule x_j = x_j-1 + L0*cosθ(j) et y_j = y_j-1 + L0*sinθ(j) . L0 est une constante, le pas de la marche au hasard. θ(j) est un angle tiré au hasard sur le domaine [0,2*π].
une marche 2D à variables aléatoires indépendantes. On calcule x_j = x_j-1 + a_x(j) et y_j = y_j-1 + a_y(j). Les variables aléatoires a_x et a _y sont distribuées selon uune loi uniforme dans le domaine [0,1].

Le programme de simulation permettra de choisir entre ces deux lois au lancement de la simulation.

Selon ces modèles, et en appliquant le théorème de la limite centrale, on démontre que le rayon quadratique moyen des particules pour une génération donnée est proportionnel à la racine carrée du temps, chaque génération étant séparée par le temps τ caractéristique unitaire de la simulation. Nous essayerons de vérifier cette propriété par la simulation.

Le programme

Le programme de simulation est écrit en C, développé comme d'habitude dans l'environnement Dev-Cpp. Il utilise la librairie graphique DisLin, comme toujours dans mes programmes en C, mais aussi, et c'est plus exceptionnel, la librairie GSL (GNU Scientific Library) pour la régression linéaire. Bien sur, vous pouvez ne pas utiliser cette librairie, mais plutôt la routine disponible sur la page Méthode des moindres carrés en TS.

Quelques commentaires sur le source du programme:

- les structures de données du programme:il en y en a deux. La structure SPoint permet de stocker les données propres à chaque particule, ses coordonnées et sa distance par rapport à l'origine. La structure SGeneration décrit les données propres à une génération (population de particules à un instant donné) : le numéro de la génération, le tableau de stockage des particules pour la génération et la distance moyenne par rapport à l'origine de l'ensemble des particules de la génération.

//*************************************************************

// Définition des structures de données

//*************************************************************

struct SPoint {

double x;

double y;

double d;

};

struct SGeneration {

int numero;

SPoint population[MAXPAR];

double dbarre;

};

- les routines de génération des nombres aléatoires. Ce sont celles déjà utilisées ici :

double Random(void)

{

const double a = 16807.0; // 7^5

const double c = 2147483647.0; // 2^31 - 1

double x,y;

x = a*seed;

seed = (long)(fmod(x,c));

y = seed/c; // normalisation par c

return (y);

}

//******************************************************************************

// Génération de la constante d'initialisation du générateur Random

// à partir de la date courante selon Anderson (1990)

//******************************************************************************

long GenerationSeed(void)

{

struct tm *temps;

time_t maintenant;

int t0,t1,t2,t3,t4,t5;

long lseed;

// récupération de l'heure actuelle

time(&maintenant); // enregistrement de l'heure actuelle

temps = localtime(&maintenant); // conversion en structure tm

t0 = temps->tm_sec;

t1 = temps->tm_min;

t2 = temps->tm_hour;

t3 = temps->tm_mday;

t4 = temps->tm_mon;

t5 = temps->tm_year;

// calcul du seed d'après l'algo d'Anderson (seed varie entre 0 et 2^31-1

lseed = t5 + 70*(t4 + 12*(t3 + 31*(t2 + 23*(t1 + 59*t0))));

if ((lseed%2) == 0) lseed = lseed-1; // s'assurer que seed est impair

return lseed;

}

- les lois d'évolution que nous avons décrites ci-dessus :

// Loi 2D à pas constant

void Evolution1(const double x0, const double y0, double *x1, double *y1)

{

const double L0 = 0.1;

*x1 = x0 + L0*cos(DEUXPI*Random());

*y1 = y0 + L0*sin(DEUXPI*Random());

}

// Marche au hasard selon une loi de probabilité uniforme

void Evolution2(const double x0, const double y0, double *x1, double *y1)

{

*x1 = x0 + Random();

*y1 = y0 + Random();

}

- la boucle de simulation: elle comprend deux boucles. La première pour chaque chaque génération, qui calcule la nouvelle génération et l'affiche. Et une boucle imbriquée, qui pour chaque génération calcule la position de chaque particule selon la loi d'évolution sélectionnée au lancement du programme.

// Boucle de simulation

for(j=1;j < NbGenerations;j++)

{

// calcul de l'évolution de la génération j

for(i=0;i < NbParticules;i++)

{

// calcul des coordonnées de la nouvelle génération

x0 = Simulation[j-1].population[i].x;

y0 = Simulation[j-1].population[i].y;

switch (NLoi)

{

case 1: Evolution1(x0,y0,&x1,&y1);

case 2: Evolution2(x0,y0,&x1,&y1);

default : Evolution1(x0,y0,&x1,&y1);

}

//normalisation des coordonnées

if (x1 <= XMIN) x1 = XMIN;

if (x1 >= XMAX) x1 = XMAX;

if (y1 <= YMIN) y1 = YMIN;

if (y1 >= YMAX) y1 = YMAX;

// calcul de la distance entre la particule et le centre

Simulation[j].population[i].d = sqrt(x1*x1 + y1*y1);

//enregistrement des coordonnées de la nouvelle génération

Simulation[j].population[i].x = x1;

Simulation[j].population[i].y = y1;

}

// affichage de la génération i

color("red"); for(i=0;i<NbParticules;i++)

rlsymb(CODESYMB,Simulation[j].population[i].x,Simulation[j].population[i].y);

//effacement de la génération avec une tempo

Sleep(TEMPO);

endgrf();

color("black");

graf(XMIN,XMAX,XMIN,GRADX,YMIN,YMAX,YMIN,GRADY);

color("white");

for(i=0;i<NbParticules;i++)

rlsymb(CODESYMB,Simulation[j].population[i].x,Simulation[j].population[i].y);

}

- enfin, après des tracés divers, le calcul de régression linéaire entre le carré de la distance moyenne des particules, soit le rayon quadratique moyen, et le temps. Ce calcul est réalisé par la routine gsl_fit_linear de la libraire GSL. Les coefficients de la régression sont affichés dans la fenêtre console du programme:

// Calcul de la régression entre t et dbarre*dbarre au moyen de la GSL GSL

X = (double *)calloc(NbGenerations, sizeof(double));

Y = (double *)calloc(NbGenerations, sizeof(double));

for(j=1;j<NbGenerations;j++)

{

X[j] = j;

Y[j] = Simulation[j].dbarre*Simulation[j].dbarre;

}

gsl_fit_linear(X,1,Y,1,NbGenerations,&c0,&c1,&cov00,&cov01,&cov11,&sumsq);

printf("r^2 = %g*t + %g\n",c1,c0);

free(X);

free(Y);

Le source est téléchargeable ici, ainsi que le fichier projet, ce qui vous évitera de chercher les librairies à charger...

Les résultats

Evolution de la population de particules

Le programme affiche deux courbes:

la courbe supérieure trace la répartition spatiale de la population de particules au fil du temps. Chaque point rouge figure une particule
la courbe inférieure est la courbe de variation du rayon quadratique moyen en fonction du temps. Dans la simulation, discrète par essence, le temps est discrétisé par le numéro de génération.

Les trois figures ci-dessous montrent l'évolution d'une population de 1 000 particules qui diffusent selon une loi uniforme (choix 2) pendant 200 générations. La première figure montre l'état à la génération 50, la seconde à la génération 100 et la troisième à la génération 200.

Sur les courbes inférieures, on note la forme de plus en plus caractéristique d'une courbe représentative de la fonction y = sqrt(x).

La diffusion, qui répond à une loi gaussienne, est isotrope dans le plan. Cela n'est pas évident sur les copies d'écran, qui sont comme vous le constatez, allongées dans le sens horizontal. C'est un artefact de la librairie graphique. Je n'ai pas encore trouvé comme tracer une figure isotrope avec DisLin...

Vérification de la loi de diffusion

J'ai reproduit ci-dessous la console du programme, qui m'indique pour une expérience donnée (1 000 particules sur 200 générations) l'équation obtenue par la régression linéaire. Vous pouvez également utiliser Excel pour calculer cette régression en utilisant les données stockées dans le fichier Diffusion.csv, produit lors de chaque expérience. Attention toutefois : le programme stocke les nombres avec un point décimal, alors qu'Excel, en version française, attend des virgules décimales. Vous devrez donc convertir au préalable les points en virgule.

Vous constatez qu'il existe une relation linéaire entre le carrée du rayon quadratique moyen et le temps, c'est à dire entre le rayon quadratique moyen et la racine carrée du temps. Le coefficient de corrélation est très bon, proche de 0,999 pour 200 générations. Evidement, si vous augmentez le nombre de particules, dans la limite de la mémoire de votre PC, le coefficient de corrélation sera encore meilleur. Pour ce faire, modifiez la constante MAXPAR dans le source du programme.