Méthode particle-in-cell

Pour de nombreux types de problèmes, la méthode particle-in-cell (PIC) classique est relativement intuitive et simple à mettre en œuvre. Cela explique probablement son succès, notamment pour la simulation des plasmas, pour laquelle la méthode comprend généralement les procédures suivantes :

• Intégration des équations du mouvement.
• Interpolation des termes sources de charge et de courant sur le maillage du champ.
• Calcul des champs aux points du maillage.
• Interpolation des champs du maillage aux emplacements des particules.

Les modèles qui incluent les interactions des particules uniquement par le biais des champs moyens sont appelés « PM » (particule-maillage). Ceux qui incluent les interactions binaires directes sont appelés « PP » (particule-particule). Les modèles avec les deux types d'interactions sont appelés « PP-PM » ou « P³M ».

La méthode PIC est sensible aux erreurs dues à ce qu'on appelle le « bruit particulaire discret »^[5]. Cette erreur est de nature statistique et reste aujourd'hui assez mal comprise.

Les algorithmes PIC géométriques modernes reposent sur un cadre théorique très différent. Ces algorithmes utilisent des outils de variétés discrètes, de formes différentielles interpolantes^[6] et d'intégrateurs symplectiques canoniques ou non canoniques pour garantir l'invariance de jauge et la conservation de la charge, de l'énergie-impulsion et, plus important encore, de la structure symplectique de dimension infinie du système particule-champ^[7],[8]. Ces caractéristiques recherchées sont attribuées au fait que les algorithmes PIC géométriques reposent sur le cadre théorique des champs plus fondamental et sont directement liés au principe variationnel de la physique.

Au sein de la communauté de recherche sur les plasmas on étudie des systèmes composés de différentes particules (électrons, ions, neutres, molécules, particules de poussière, etc.). L'ensemble des équations associées aux codes PIC comprend donc la force de Lorentz qui décrit le mouvement et est résolue par le module de déplacement ou « pousseur » du code, ainsi que les équations de Maxwell qui déterminent les champs électrique et magnétique et sont calculées par le solveur de champs.

Les systèmes réels étudiés sont souvent extrêmement grands en termes de nombre de particules qu'ils contiennent. Afin de rendre les simulations possibles on utilise des « super-particules ». Une super-particule (ou « macroparticule ») est une particule de calcul qui représente de nombreuses particules réelles. Il peut s'agir d'un grand nombre d'électrons ou d'ions dans le cas d'une simulation de plasma ou de tourbillons élémentaires dans une simulation de fluide. Il est possible de faire varier le nombre de particules car l'accélération due à la force de Lorentz dépend uniquement du rapport charge/masse. Une super-particule suivra donc la même trajectoire qu'une particule réelle (aux erreurs numériques près).

Le nombre de particules réelles correspondant à une super-particule doit être choisi de manière à pouvoir recueillir des statistiques suffisantes sur le mouvement de la particule. S'il existe une différence significative entre la densité des différentes espèces du système (entre ions et neutres par exemple), des rapports particules réelles/super-particules distincts peuvent être utilisés pour chacune d'elles.

Même avec des super-particules, le nombre de particules simulées est généralement très élevé (> 10⁵), et le module de déplacement des particules est souvent l'étape la plus gourmande en temps de calcul. Par conséquent, il doit doit être rapide. Des efforts considérables ont été consacrés à l'optimisation des différents schémas.

Ces schémas se divisent en deux catégories : les solveurs implicites et les solveurs explicites. Alors que les solveurs implicites (par exemple, le schéma d'Euler implicite) calculent la vitesse des particules à partir des champs déjà mis à jour, les solveurs explicites utilisent uniquement la force de l'étape de temps précédente. Ils sont donc plus simples et plus rapides, mais nécessitent un pas de temps plus petit. On utilise ainsi la méthode saute-mouton, une méthode explicite du second ordre, ou bien l'algorithme de Boris qui annule le champ magnétique dans l'équation de Newton-Lorentz^[9],[10].

Pour les applications plasma, la méthode saute-mouton prend la forme suivante :

${\displaystyle {\frac {\mathbf {x} _{k+1}-\mathbf {x} _{k}}{\Delta t}}=\mathbf {v} _{k+1/2}}$

${\displaystyle {\frac {\mathbf {v} _{k+1/2}-\mathbf {v} _{k-1/2}}{\Delta t}}={\frac {q}{m}}\left(\mathbf {E} _{k}+{\frac {\mathbf {v} _{k+1/2}+\mathbf {v} _{k-1/2}}{2}}\times \mathbf {B} _{k}\right)}$

où l'indice ${\displaystyle k}$ se réfère aux quantités de l'étape de temps précédente, ${\displaystyle k+1}$ aux quantités mises à jour de l'étape de temps suivante (c.-à-d. ${\displaystyle t_{k+1}=t_{k}+\Delta t}$) et les vitesses sont calculées entre les étapes de temps habituelles ${\displaystyle t_{k}}$.

Les équations du schéma de Boris, substituées aux équations ci-dessus, sont :

${\displaystyle \mathbf {x} _{k+1}=\mathbf {x} _{k}+{\Delta t}\,\mathbf {v} _{k+1/2}}$

${\displaystyle \mathbf {v} _{k+1/2}=\mathbf {u} '+q'\mathbf {E} _{k}}$

avec :

${\displaystyle \mathbf {u} '=\mathbf {u} +(\mathbf {u} +(\mathbf {u} \times \mathbf {h} ))\times \mathbf {s} }$

${\displaystyle \mathbf {u} =\mathbf {v} _{k-1/2}+q'\mathbf {E} _{k}}$

${\displaystyle \mathbf {h} =q'\mathbf {B} _{k}}$

${\displaystyle q'=\Delta t\,q/(2m)}$

${\displaystyle \mathbf {s} =2\mathbf {h} /(1+h^{2})}$

Grâce à son excellente précision aux temps longs, l'algorithme de Boris est devenu la norme de facto pour la propagation d'une particule chargée. On a constaté que l'excellente précision de l'algorithme non relativiste de Boris, est due à la conservation du volume de l'espace des phases, bien qu'il ne soit pas symplectique. La borne supérieure globale sur l'erreur d'énergie, généralement associée aux algorithmes symplectiques, reste valable pour cet algorithme, ce qui en fait un algorithme efficace pour la dynamique multi-échelle des plasmas. Il a également été démontré qu'il est possible d'améliorer la poussée relativiste de Boris pour qu'elle conserve le volume et admette une solution à vitesse constante dans des champs E et B croisés^[11].

Les méthodes les plus couramment utilisées pour résoudre les équations de Maxwell (ou plus généralement, les équations aux dérivées partielles) appartiennent à l'une des trois catégories suivantes :

• Méthodes différences finies (MDF)
• Méthodes éléments finis (MEF)
• Méthodes spectrales

Avec la MDF, le domaine continu est remplacé par une grille discrète de points sur laquelle sont calculés les champs électrique et magnétique. Les dérivées sont ensuite approchées par les différences entre les valeurs des points voisins de la grille, transformant ainsi les équations aux dérivées partielles en équations algébriques.

Avec la MEF, le domaine continu est divisé en un maillage discret d'éléments. Les équations aux dérivées partielles sont traitées comme un problème aux valeurs propres et une solution d'essai est initialement calculée à l'aide de fonctions de base localisées dans chaque élément. La solution finale est ensuite obtenue par optimisation jusqu'à l'obtention de la précision requise.

Les méthodes spectrales, telles que la transformée de Fourier rapide (FFT), transforment également les EDP en un problème de valeurs propres, mais cette fois-ci les fonctions de base sont d'ordre élevé et définies globalement sur tout le domaine. Le domaine lui-même n'est pas discrétisé ; il reste continu. Là encore, une solution d'essai est trouvée en insérant les fonctions de base dans l'équation aux valeurs propres, puis optimisée afin de déterminer les meilleures valeurs des paramètres initiaux.

L'appellation « particles-in-cell » provient de la manière dont les macro-quantités du plasma (densité numérique, densité de courant, etc.) sont attribuées aux particules de la simulation (pondération des particules). Les particules peuvent se situer n'importe où sur le domaine continu, mais les macro-quantités, tout comme les champs, ne sont calculées qu'aux points du maillage. Pour obtenir les macro-quantités, on suppose que les particules ont une forme donnée, déterminée par la fonction de forme ${\displaystyle S(\mathbf {x} -\mathbf {X} )}$ où ${\displaystyle \mathbf {x} }$ est la coordonnée de la particule et ${\displaystyle \mathbf {X} }$ celle du point d'observation.

Le choix le plus simple et le plus courant pour la fonction de forme est sans doute le schéma dit « cloud-in-cell » (CIC), qui est un schéma de pondération linéaire du premier ordre. Quel que soit le schéma choisi, la fonction de forme doit satisfaire les conditions suivantes^[12] : les champs obtenus par le solveur de champs sont déterminés uniquement aux points de la grille et ne peuvent être utilisés directement dans le module de déplacement des particules pour calculer la force agissant sur celles-ci, ils doivent être interpolés par pondération :

${\displaystyle \mathbf {E} (\mathbf {x} )=\sum _{i}\mathbf {E} _{i}S(\mathbf {x} _{i}-\mathbf {x} )}$

où l'indice ${\displaystyle i}$ désigne le point de la grille. Afin de garantir la cohérence des forces agissant sur les particules, le calcul des macro-quantités à partir des positions des particules sur la grille et l'interpolation des champs des points de la grille vers les positions des particules doivent également être cohérents, car les deux méthodes apparaissent dans les équations de Maxwell. De plus, le schéma d'interpolation des champs doit conserver la quantité de mouvement. Ceci peut être réalisé en choisissant le même schéma de pondération pour les particules et les champs et en assurant simultanément la symétrie spatiale appropriée (c'est-à-dire l'absence de force auto-induite et le respect des lois de Newton) du solveur de champ.

Simuler l'interaction pour chaque paire d'un grand système serait trop coûteux en calcul ; c'est pourquoi plusieurs mméthodes de type Monte-Carlo ont été développées. Une méthode largement utilisée est le « modèle de collision binaire »^[13] dans lequel les particules sont regroupées selon leur cellule, puis appariées aléatoirement, et enfin les paires entrent en collision.

Dans un plasma réel, de nombreuses autres interactions peuvent intervenir, allant des collisions élastiques, telles que les collisions entre particules chargées et neutres, aux collisions inélastiques, telles que la collision d'ionisation électron-neutre, en passant par les réactions chimiques ; chacune d'elles nécessitant un traitement distinct. La plupart des modèles de collision gérant les collisions chargées-neutres utilisent soit le schéma « Monte-Carlo direct », dans lequel toutes les particules portent une information sur leur probabilité de collision, soit le schéma « collision nulle »^[14],[15] qui n'analyse pas toutes les particules mais utilise plutôt la probabilité de collision maximale pour chaque espèce chargée.

Comme pour toute méthode de simulation le pas de temps et la taille de la grille doivent être soigneusement choisis afin que les phénomènes d'échelle temporelle et spatiale étudiés soient correctement résolus. De plus, le pas de temps et la taille de la grille influent sur la vitesse et la précision du code.

Pour une simulation de plasma électrostatique utilisant un schéma d'intégration temporelle explicite (par exemple, la méthode saute-mouton, la plus courante), deux conditions importantes concernant la taille de la grille Δx et le pas de temps Δt doivent être respectées pour garantir la stabilité de la solution :

Δx < 3,4 λ_D

Δt ≤ 2 ω_pe^-1

Ces conditions peuvent être obtenus en considérant les oscillations harmoniques d'un plasma unidimensionnel non magnétisé. Cette dernière condition est strictement requise, mais des considérations pratiques liées à la conservation de l'énergie suggèrent d'utiliser une contrainte beaucoup plus stricte où le facteur 2 est remplacé par un nombre d'un ordre de grandeur inférieur. L'utilisation de ${\displaystyle \Delta t\leq 0.1\omega _{pe}^{-1}}$ est typique^[12],[16]. Sans surprise, l'échelle de temps naturelle dans le plasma est donnée par l'inverse de la fréquence plasma ${\displaystyle \omega _{pe}^{-1}}$ et l'échelle de longueur par la longueur de Debye ${\displaystyle \lambda _{D}}$.

Pour une simulation explicite de plasma électromagnétique, le pas de temps doit également satisfaire la condition de Courant-Friedrichs-Lewy (condition CFL) :

Δt < Δx / c

où Δx ≈ λ_D et c est la vitesse de la lumière.

En physique des plasmas la simulation PIC a été utilisée avec succès pour étudier les interactions laser-plasma, l'accélération des électrons et le chauffage ionique dans l'ionosphère aurorale, la magnétohydrodynamique, la reconnexion magnétique ainsi que le gradient de température ionique et d'autres micro-instabilités dans les tokamaks, les décharges sous vide et les plasmas poussiéreux.

Les modèles hybrides peuvent utiliser la méthode PIC pour le traitement cinétique de certaines espèces, tandis que d'autres espèces (de distribution maxwellienne) sont simulées par un modèle fluide.

Les simulations PIC ont également été appliquées en dehors de la physique des plasmas à des problèmes de mécanique des solides et de mécanique des fluides^[17],[18] comme la méthode particle-in-cell multiphasique (en).