Modèles de Structures Aléatoires de Type Réaction-Diffusion - Thèse de Morphologie Mathématique - Luc Decker, Ecole des Mines de Paris (1999)

suivant: 2.3 Automates cellulaires   remonter: 2. Modèles de diffusion   précédent: 2.1 Systèmes à base   [Table des matières]

2.2 Marcheurs aléatoires

2.2.1 Principes

Dans la théorie du Mouvement Brownien, on considère qu'un fluide est constitué de particules élémentaires, qui disposent de leur propre liberté de mouvement. A l'échelle la plus réduite, il s'agit de molécules. Mais ce concept reste valide à plus grande échelle, et sert de base à une classe de modèles de fluides. Les particules engendrent un comportement collectif (ou moyen), qui est caractérisé au niveau macroscopique par des statistiques sur les grandeurs mesurables pour chaque particule isolément. On considère ainsi la résultante des contributions des particules élémentaires, qui forment une population représentative du fluide, dans toute sa variabilité. On peut y voir une application du célèbre théorème de la Limite Centrale : la somme de variables aléatoires indépendantes suit une loi normale.

Selon une approche lagrangienne, on va s'intéresser à la trajectoire d'une particule individuelle en fonction du temps. Dans son mouvement aléatoire désordonné, la particule participe au processus de diffusion du fluide. Dans le cas présent, on ne cherchera pas à reproduire la cause de ce mouvement brownien, c'est à dire les collisions ou interactions avec d'autres particules. Du point de vue des simulations, introduisons la notion de marcheur aléatoire : une particule dont la direction de déplacement change aléatoirement à chaque pas de temps. Cette approche ne relève pas du même principe que les systèmes de particules. D'une part, parce les trajectoires des particules sont simulées une par une. D'autre part parce que l'objectif final est d'obtenir des mesures quantitatives sur le fluide simulé. Les conditions aux limites influencent le comportement du fluide, comme par exemple dans le cas d'un milieu poreux. Inversement, les mesures portant sur la diffusion du fluide - la dispersion des particules - constituent un outil pour caractériser le milieu traversé.

2.2.2 Mise en oeuvre

Considérons un milieu poreux représenté sous la forme d'une image binaire tri-dimensionnelle. Dans ce domaine discret, les voxels sont répartis entre la matrice (les parties solides du milieu) et les pores (le vide). Soit un fluide baignant ce milieu poreux, avec l'hypothèse que sa diffusion est nulle à l'intérieur de la matrice. Pour simuler ce fluide, on considère une population de $N$ particules réparties uniformément dans les pores assez éloignés des limites du domaine, à l'instant initial $t_{0} = 0$. La trajectoire de chaque particule dans le domaine est alors simulée indépendamment, par déroulement du procédé itératif suivant:
A chaque pas de temps élémentaire,

• Le temps discret $t$ est incrémenté.
• Soit $(x(t), y(t), z(t))$ la position d'une particule. Un tirage aléatoire donne la nouvelle position $(x',y',z')$ vers laquelle la particule est susceptible de se diriger. Les déplacements selon chacune des trois directions sont séparables; ils sont unitaires ou nuls, et définis à l'aide des probabilités suivantes:

  $$\hspace{-1cm} \left\{ \begin{array}{l} P[x' = x(t) + 1] = p_{x}\ ... ...z}}\ \ P[z' = z(t) - 1] = 1 - p_{z} - p_{o_{z}}\\ \end{array} \right. \ \ \ \$$ (II.-7)

  
  

• Si la position $(x',y',z')$ appartient aux pores du milieu, la particule peut s'y déplacer. D'où:

  $$x(t+1)=x',\ \ y(t+1)=y',\ \ z(t+1)=z'$$ (II.-8)

  
  

• Dans le cas contraire, la particule se heurte à un obstacle et reste donc immobile durant un intervalle de temps, soit:

  $$x(t+1)=x(t),\ \ y(t+1)=y(t),\ \ z(t+1)=z(t)$$ (II.-9)

  
  

La procédure s'achève par l'un des deux cas suivants:

• Un temps $t = T_{max}$ est atteint; celui-ci est fixé arbitrairement. Plus ce temps est court, plus les distances parcourues par les particules sont faibles. Il en résulte que $T_{max}$ définit l'échelle de la simulation, et par extension, l'échelle spatiale de validité des grandeurs mesurées à partir des trajectoires.
• La particule quitte le domaine au temps $t_{f} < T_{max}$, en atteignant l'un des bords. On mesure ainsi son temps de séjour ou temps de traversée du milieu poreux.

Les marcheurs aléatoires ainsi définis se déplacent sur une trame cubique. Par combinaison des mouvements autorisés dans chaque direction, il existe $26$ déplacements unitaires possibles pour une particule à condition que les probabilités d'immobilisation $p_{o_{x}}$, $p_{o_{y}}$ et $p_{o_{z}}$ soient non nulles. Cette valeur correspond à la connexité de la trame: chaque point admet 26 voisins. Une particule peut aussi rester immobile avec la probabilité $p_{o_{x}} p_{o_{y}} p_{o_{z}}$. Par ailleurs, les probabilités $p_{x}$, $p_{y}$ et $p_{z}$ permettent d'introduire une dérive dans les déplacements des particules. Pour que ces déplacements soient symétriques dans chaque direction, il faut que les probabilités de déplacements positifs $(+1)$ et négatifs $(-1)$ soient égales, d'où

$p_{x} = \frac{1 - p_{o_{x}}}{2}$, $p_{y} = \frac{1 - p_{o_{y}}}{2}$, $p_{z} = \frac{1 - p_{o_{z}}}{2}$

Dans le cas contraire, le déséquilibre s'apparente à l'action d'une force d'entraînement des particules: il est ainsi possible de simuler un écoulement à travers le milieu poreux. La Figure II-1 montre quelques exemples de trajectoires; un histogramme des temps de séjours de particules à l'intérieur d'un domaine parfaitement homogène est donné sur la Figure II-2.

La vitesse moyenne $v_{i}$ d'une particule $i$ au cours de sa trajectoire s'achevant au temps $t_{f_{i}}$ est estimée par

$$v_{x_{i}} = \frac{1}{t\ (t_{f_{i}} - t)} \sum_{h = 0}^{t_{f_{i}} - t} (x_{i}(h+t) - x_{i}(h))$$ (II.-10)

dans la direction $x$. Des expressions semblables sont établies pour les directions $y$ et $z$. On en déduit les projections de la vitesse moyenne de la population de particules:

$$\displaystyle{v_{x} = \frac{1}{N} \sum_{i = 0}^{N} v_{x_{i}}} \ \... ... v_{y_{i}}} \ \ \ \displaystyle{v_{z} = \frac{1}{N} \sum_{i = 0}^{N} v_{z_{i}}}$$ (II.-11)

Cette vitesse moyenne tend à être nulle en l'absence de dérive. Pour des marcheurs aléatoires dans un milieu homogène, on a:

$$v_{x} = 2\ p_{x} + p_{0_{x}} - 1 \ \ \ \ v_{y} = 2\ p_{y} + p_{0_{y}} - 1 \ \ \ \ v_{z} = 2\ p_{z} + p_{0_{z}} - 1$$ (II.-12)

Deux exemples de distributions expérimentales des vitesses moyennes sont présentés sur la Figure II-3.

Parallélisation. La simulation des marcheurs aléatoires a été implantée sur le calculateur parallèle IBM-SP2 de l'Ecole des Mines de Paris. Une gestion dynamique des tâches attribuées aux différents processeurs a été mise en place. Lorsque les trajectoires de centaines de milliers de particules doivent être générées, on définit la taille d'une tranche de calcul élémentaire (1000 particules par exemple) que chaque processeur traite indépendamment en chargeant la totalité du milieu poreux. En cours de simulation, les processeurs communiquent entre eux pour s'informer du nombre total de particules restant à simuler, et le cas échéant, entreprennent de nouvelles tranches de calcul. Grace à ce système, les calculs s'achèvent presque au même moment sur tous les processeurs, chacun en ayant accompli une part proportionnelle à sa disponibilité, en cas de concurrence avec d'autres applications. L' efficacité d'un tel programme parallèle est proche de 100% puisque les communications entre processeurs sont extrêmement réduites : l'utilisation de $N$ processeurs permet de diviser les temps de simulation par $N$.

2.2.3 Estimation des coefficients de diffusion

La modélisation d'un fluide par marcheurs aléatoires permet d'obtenir une estimation de ses coefficients de diffusion dans un milieu poreux. Il s'agit de la principale application de cette classe de modèle. Ces coefficients de diffusion peuvent participer à la caractérisation de la morphologie de structures aléatoires quelconques, lorsqu'elles sont considérées comme un milieu poreux binaire. Nous effectuerons de telles mesures sur les structures générées par réaction-diffusion au cours de la partie III.

Soient $D_{x}$, $D_{y}$ et $D_{z}$ les coefficients de diffusion du fluide à travers le milieu poreux, pour les trois directions. Ces coefficients macroscopiques sont définis pour un milieu homogène équivalent. Leur estimation est obtenue par l'une ou l'autre des méthodes suivantes, qui exploitent les statistiques des trajectoires des marcheurs aléatoires: [Matheron79,Jeulin92a,Akoulenko95]

• Méthode du variogramme des positions des particules
  Dans ce cas, on enregistre les positions $(x(t), y(t), z(t))$ de chaque particule du temps $t = 0$ au temps $t = T_{max}$. Il est indispensable que toutes les particules restent dans le domaine jusqu'au temps $T_{max}$, à fixer en conséquence. D'après la loi de Fick et les résultats de la théorie du mouvement Brownien, on a :

  $$D_{x} \ t = \frac{1}{2} E[X_{h+t} - X_{h} - v_{x} t]^{2}$$ (II.-13)

  
  
  où $v_{x}$ est la vitesse moyenne des particules (sa projection selon la direction $x$). Cette expression est à comparer avec le variogramme d'ordre deux des positions des particules au cours du temps, défini par :

  $$\gamma(t) = \frac{1}{2} E[X_{h+t} - X_{h}]^{2}$$ (II.-14)

  
  
  Si $v_{x} = 0$, en l'absence de dérive, on estime directement $D_{x}$ à partir de la pente de la moyenne des variogrammes expérimentaux $\gamma(t)$ sur $N$ particules. Dans le cas contraire, on a :

  $$\frac{\gamma(t)}{t} = \frac{v_{x}^{2}}{2} t + D_{x}$$ (II.-15)

  
  
  et il en résulte que $D_{x}$ est l'ordonnée à l'origine de la droite approximant la courbe expérimentale $\frac{\gamma(t)}{t}$. En général, pour une population de particules suffisamment importante, les courbes théoriques et expérimentales sont parfaitement confondues. Les temps de séjour minimaux sont connus approximativement en appliquant la méthode des temps de séjour décrite ci-dessous, ce qui permet de définir la durée $T_{max}$ des simulations. Lorsque la taille du domaine est réduite dans une direction, ces temps de séjours sont très courts et cette méthode conduit à une estimation du coefficient de diffusion à une échelle spatiale relativement faible.
  

  Dans le cas particulier où les particules se déplacent dans un espace homogène, on obtient les valeurs de référence des coefficients de diffusion en fonction des probabilités de la marche aléatoire:

  $$D_{x} = p_{x} \ \ \ \ D_{y} = p_{y} \ \ \ \ D_{z} = p_{z}$$ (II.-16)

  
  

  Diffusion non-fickienne. Pour certaines catégories de milieux poreux, la diffusion n'obéit pas à la loi de Fick. C'est même une situation couramment observée dans l'étude de problèmes concrets, pour lequels la loi de Fick est une simple approximation. La morphologie complexe d'un milieu poreux est à l'origine de non-linéarités, par exemple lorsque ce dernier renferme des cavités sans issues, qui jouent le rôle de "pièges" à une échelle donnée et peuvent être contournées à une échelle plus grande. Si l'estimation de la diffusion par marcheurs aléatoires reste valable, la relation qui relie le variogramme des trajectoires et le coefficient de diffusion prend la forme suivante pour certains types de diffusion :

  $$D_{x} \ t^{\alpha} = \frac{1}{2} E[X_{h+t} - X_{h} - v_{x} t]^{2}\ \ \ \ \ \alpha < 2$$ (II.-17)

  
  
  L'exposant $\alpha$ caractérise le mode de diffusion; pour $\alpha = 1$, on retrouve la loi de Fick.
  

• Méthode des temps de séjour
  Le temps de séjour $T_{a}$ d'une particule correspond au temps que celle-ci passe dans le domaine tri-dimensionnel, avant d'en atteindre un bord fixé situé à une distance $a$ de sa position de départ moyenne, soit $t_{f} = T_{a}$ pour $x(t_{f}) = a$. Pour mesurer ces temps de séjour, on imposera donc une dérive dans le mouvement des particules afin d'être certain qu'elles quittent effectivement le domaine par le bord choisi et l'on fixe $T_{max} = \infty$. Dans le cas où le milieu suit une loi de Fick de coefficient de diffusion $D$, la fonction de répartition (ou distribution) des temps de séjour des particules est donnée par:

  $$F_{T_{a}}(t) = P[X(t) \geq a]\ P[X(t) \geq a\ \vert\ t \geq T_{a}]$$ (II.-18)

  
  
  Lorsque la force effective est suffisamment élevée, on a $P[X(t) \geq a\ \vert\ t \geq T_{a}] \approx 1$ et on peut alors montrer que:

  $$P[X(t) \geq a] = \frac{1}{\sqrt{4 \pi D_{x} t}} \int_{a}^{\infty} e^{- \frac{(x - v_{x} t)^2}{4 D_{x} t}} dx$$ (II.-19)

  
  

  soit

  $$F_{T_{a}}(t) = \Phi(\infty) - \Phi(\frac{a - v_{x} t}{\sqrt{2 D_{x} t}}) = \frac{1}{2} - \Phi(\frac{a - v_{x} t}{\sqrt{2 D_{x} t}})$$ (II.-20)

  
  
  où

  $$\Phi(x) = \frac{1}{\sqrt{2 \pi}} \int_{0}^{x} e^{-\frac{k^{2}}{2}} dk$$ (II.-21)

  
  
  est la fonction de Laplace; on a aussi $\Phi(x) = \frac{1}{2}\ erf(\frac{x}{\sqrt{2}}))$ où la fonction $erf$ donne l'intégrale de la distribution de Gauss (fonction erreur). Les coefficients de diffusion seront alors déterminés par ajustement de la fonction de distribution théorique avec l'histogramme cumulé des temps de séjours résultant de la simulation des trajectoires de plusieurs milliers de particules.
  

  Remarquons que cette méthode ne s'applique pas dans le cas où une partie des temps de séjour des particules sont infinis, lorsque certaines particules restent captives dans les porosités qui ne communiquent pas avec l'extérieur. Notons aussi que cette méthode permet d'estimer un coefficient de séjour pour des temps longs et pour de grandes échelles de porosité, ce qui correspond à la méthode physique de mesure d'un coefficient de diffusion par l'intermédiaire d'un traceur traversant l'échantillon.

  Enfin, la mise en place de parois étanches sur les bords du domaine, dans une ou deux directions, permet de mesurer les temps de séjour dans les directions restantes sans être confronté au problème des sorties des particules par un autre bord que celui qui est fixé pour la mesure, même quand la dérive de la marche aléatoire est réduite. L'intéret de fermer ainsi une partie des bords du domaine apparait surtout lorsque la taille du domaine dans une direction est très inférieure à ses dimensions dans les autres directions. On constate que dans le cas de parois latérales étanches, les distributions théoriques et expérimentales des temps de séjour restent en correspondance.

2.2.4 Application: évaluation de la diffusion au travers d'un échantillon réel

Dans le cadre d'une étude commandée par le Commissariat à l'Energie Atomique, nous avons eu à caractériser les propriétés de transport d'un matériau poreux [Tovena98,Decker98a]. Les données qui nous été fournies consistent en une séquence de 50 photographies représentant les sections d'un échantillon cylindrique du matériau. Ces prises de vue ont été réalisées au cours de sa découpe en tranches successives. Il s'agit d'un milieu hétérogène particulièrement intéressant, essentiellement constitué de feuilles repliées et compressées. L'étude porte en particulier sur la diffusion dans la direction verticale $(y)$. En effet, l'observation des sections a permis de constater qu'il n'existe aucun passage vertical à deux dimensions.

La détection des porosités par traitement d'image constitue la première étape de l'étude, que nous ne détaillerons pas ici. L'exploitation des clichés a posé de nombreuses difficultés, en raison de leur qualitité médiocre. L'échantillon a ensuite été reconstitué dans une image binaire tri-dimensionnelle, par mise en correspondance des sections consécutives. La fraction volumique totale des pores est de l'ordre de 22%. Dans ce milieu poreux virtuel, un test de percolation a tout d'abord été réalisé afin de vérifier qu'il existait bien un chemin entre faces permettant de le traverser. Dans ce but, un algorithme de propagation géodésique a été appliqué [Jeulin92b,Demarty96]. Les porosités ouvertes sur l'extérieur de l'échantillon sont également détectées par cet algorithme; leur volume représente 88% de l'ensemble des porosités : le milieu poreux est fortement interconnecté (Figure II-5).

Finalement, des simulations de marcheurs aléatoires sont effectuées dans les porosités ouvertes afin de déterminer les coefficients de diffusion d'un milieu homogène équivalent, dans les directions $x$ et $y$. D'une part, la méthode du variogramme donne une estimation à petite échelle, avec des distances moyennes parcourues par les marcheurs aléatoires de l'ordre de $2 mm$ pour $T_{max} = 1200$ itérations; pour des temps plus longs, une partie des marcheurs sortent en effet du domaine. Dans la direction horizontale $x$, la diffusion respecte la loi de Fick et on obtient un rapport $D_{x}/D_{ref} = 0.32$, $D_{ref}$ étant le coefficient de diffusion de référence pour un milieu totalement poreux. La Figure II-6 présente le variogramme des trajectoires dans la direction verticale; comme ce dernier n'est pas linéaire, on formule l'hypothèse que la loi de Fick n'est pas vérifiée. En conséquence, une loi de diffusion de type $D\ t^{\alpha}$ est appliquée pour obtenir un ajustement avec une loi théorique; la méthode usuelle des "moindres carrés" est mise en oeuvre pour estimer simultanément $D$ et $\alpha$. D'autre part, la méthode des temps de séjour donne des coefficients de diffusion plus élevés à grande échelle, lorsque l'interconnection tri-trimensionnelle du milieu poreux intervient davantage. L'analyse des résultats permet aussi de vérifier l'anisotropie du milieu poreux étudié: le coefficient de diffusion est cinq fois plus faible dans la direction $y$ orthogonale aux porosités lamellaires.

suivant: 2.3 Automates cellulaires   remonter: 2. Modèles de diffusion   précédent: 2.1 Systèmes à base   [Table des matières]

Decker, Luc. "Modèles de structures aléatoires de type réaction-diffusion". PhD diss. (191 p.), Paris, ENSMP-CMM, 1999.