Developpement d’un code PIC/MCC 2d3v

Résolution des équations de champs

Le cadre de travail retenu étant celui d’un modèle électrostatique, la résolution des équations de Maxwell se limite `a la seule résolution de l’équation de Poisson ∇ · (ε∇φ(r, t)) = −ρ(r, t) sur Ω (4.1) associée aux conditions aux limites φ = f sur ΓD (4.2a) ∂φ ∂n = 0 sur ΓN (4.2b) o`u ΓD et ΓN sont les frontières correspondant respectivement aux conditions aux limites de type Dirichlet (potentiel électrique φ fixé) et Neumann (champ électrique normal fixé), n est la normale sortante au domaine de calcul Ω et f est une fonction de carré intégrable définie sur ΓD. On a par ailleurs ∂Ω = ΓD ∪ ΓN

Resolution des équations de champs

Méthodes de résolution aux éléments finis

La dérivation de la formulation faible du problème mixte pour l’équation de Poisson (4.1-4.2) réalisée en Annexe C conduit `a Z Ω (ε∇φ) · ∇vdΩ = Z Ω ρvdΩ ∀v ∈ H 1 0,ΓD (Ω). (4.3) La discrétisation de la formulation faible (4.3) est réalisée au moyen d’éléments de Lagrange quadrangles de degré un. Soit Nd le nombre de nœuds correspondant `a une condition aux limites de type Dirichlet et Nn le nombre de nœuds étant soit `a l’intérieur du domaine Ω, soit sur une frontière de type Neumann. En notant ψi : Ω → R la i-ème des Nt = Nd + Nn fonctions de base globales, le potentiel électrique φ est approché sur Ω par Φ = X Nn i=1 φiψi + X Nd j=1 fjψj (4.4) o`u fj est la valeur de la fonction f – potentiel électrique fixé – au nœud j et φi est la valeur du potentiel électrique φ au nœud i que l’on cherche `a déterminer. La substitution de l’approximation (4.4) dans (4.3) mène, après utilisation des fonctions de base globales ψi en tant que fonction test v, au système linéaire KΦ = Aρ + BF (4.5) o`u K, A et B sont respectivement des matrices de taille Nn × Nn, Nn × Nt et Nn × Nd, ρ est le vecteur ayant pour composantes les densités de charge `a chacun des Nt nœuds et F = (f1, · · · , fNd ) est le vecteur rendant compte du potentiel électrique imposé sur les frontières de type Dirichlet. Les matrices K, A et B sont obtenues par assemblage de matrices élémentaires connues élément par élément. En particulier, la matrice K a pour coefficient kij avec kij = Z Ω ε∇ψi · ∇ψjdΩ. (4.6) Puisque les coefficients kij , aij et bij des matrices K, A et B sont fixés une fois le maillage du domaine réalisé, le calcul en amont des matrices K−1A et K−1B permet de calculer directement le potentiel Φ `a chaque pas de temps sous la forme d’une somme de deux produits matriciels Φ = K −1Aρ + K −1BF, (4.7) dans laquelle seul ρ et éventuellement F sont fonction du temps. 

Validation

Afin de démontrer les capacités du solveur présenté au paragraphe précédent, celui ci est évalué sur différents cas test. La première configuration est un cas de charge simple pour lequel une solution analytique peut-ˆetre obtenue. Toutefois, afin de tester le code sur des configurations plus réalistes, le solveur est ensuite évalué en le confrontant `a un logiciel commercial éprouvé. 62 Developpement d’un code PIC/MCC 2d3v ´ Comparaison sur un cas analytique 0 0.5 1 0 0.5 1 0 500 1000 1500 z (m) r (m) φ (V) 0 200 400 600 800 1000 1200 1400 Figure 4.1: Potentiel électrique pour une distribution cylindrique uniforme de charge Dans un milieu homogène de permittivité ε0, le potentiel électrique peut ˆetre exprimé de manière analytique dans le cas d’un chargement cylindrique uniforme en volume. En effet, l’équation de Poisson s’écrit simplement dans ce cas 1 r ∂ ∂r  r ∂φ ∂r  = − ρ0 ε0 , (4.8) et le potentiel électrique φ a alors pour solution analytique φ˜(r) = ρ0 4ε0 (r 2 c − r 2 ) (4.9) o`u ρ0 est la densité uniforme de charge et rc le rayon du cylindre chargé. La figure 4.1 présente la solution numérique obtenue sur un maillage uniforme de dix mailles par dix mailles, avec pour conditions aux limites φ(rc, z) = 0 (4.10a) ∂φ ∂r (0, z) = 0 (4.10b) ∂φ ∂z (r, 0) = 0 et ∂φ ∂z (r, L) = 0. (4.10c) La densité uniforme de charge est ρ0 = 5 × 10−8 C.m−3 , et le rayon ainsi que la longueur du cylindre sont pris de longueur unitaire (rc = 1 m, L = 1 m). L’erreur relative commise par rapport `a la solution analytique φ˜ telle que rapportée dans le tableau 4.1 est au plus de 0.7 % sur l’axe, et décroit lorsque l’on s’éloigne de l’axe. On constate par ailleurs qu’un raffinement du maillage dans la direction radiale (20 mailles × 10 mailles) permet de ramener l’erreur relative `a des niveaux tout `a fait négligeables. r (m) 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 φ − φ˜ φ˜ ( o/oo) [10 × 10] 7.20 3.82 2.73 2.56 1.78 1.11 1.30 1.11 1.11 0.99 φ − φ˜ φ˜ ( o/oo) [20 × 10] 2.29 0.97 0.52 0.23 0.10 0.17 0.30 0.28 0.33 0.24 Table 4.1: Erreur relative en pour mille commise par le solveur relativement `a la solution analytique φ˜ pour différents raffinements de maillage : 10 mailles × 10 mailles et 20 mailles × 10 mailles. 

Confrontation avec un solveur

commercial Les performances du solveur sont évaluées sur le cas test présenté figure 4.2(a) pour différentes valeurs de permittivité relative εr et densité de charge ρ0. Pour ce faire, les résultats sont comparés `a des calculs réalisés au moyen du logiciel commercial Comsol R sur un maillage identique `a celui utilisé par le solveur (figure 4.2(b)), en utilisant également des éléments de Lagrange quadrangles d’ordre un. On s’intéresse plus particulièrement `a trois cas test : – ρ0 = 0 et εr = 1 : cas simple du calcul du potentiel électrique laplacien pour une électrode cylindrique portée `a un potentiel positif. – ρ0 = 5 × 10−8 C.m−3 et εr = 1 : mˆeme géométrie d’électrodes mais prise en compte du champ de charge d’espace induit par une distribution homogène de charge en volume. – ρ0 = 0 et εr = 4 : cas d’une électrode cylindrique placée `a l’intérieur d’un di- électrique, calcul du potentiel laplacien sur le domaine complet. La valeur de permittivité relative utilisée est représentative du plexiglas (εr = 3.3). Quel que soit le cas, les conditions aux limites sont les mˆemes, `a savoir un potentiel φ fixé sur la découpe du domaine (0 ≤ ξ ≤ 0.1 o`u ξ désigne r ou z), une condition type Neumann homogène sur le segment 0.1 < ξ < 0.3, et une condition de potentiel nul sur le reste des frontières du domaine.