Équation de Poisson à deux dimensions

1. Introduction

Un exemple d'équation de Poisson est celle vérifiée par le potentiel électrostatique :

\frac{\partial^{2} V}{\partial x^{2}} + \frac{\partial^{2} V}{\partial y^{2}} = - \frac{ρ}{ε}

où $ρ (x, y)$ est la densité volumique de charge électrique.

Un autre exemple est l'équation de la chaleur en régime stationnaire vérifiée par la température :

\frac{\partial^{2} T}{\partial x^{2}} + \frac{\partial^{2} T}{\partial y^{2}} = - \frac{σ}{λ}

où $σ (x, y)$ est la densité volumique de puissance générée localement.

La forme générale de l'équation de Poisson est :

\frac{\partial^{2} u}{\partial x^{2}} + \frac{\partial^{2} u}{\partial y^{2}} = s (x, y)

où s(x,y) est une fonction connue, appelée la source.

Pour la discrétisation, on utilise aussi la forme intégrale de cette équation, pour un volume V délimitée par une surface S :

\iint_{S} \vec{g r a d} (u) \cdot \vec{n} d s = ∭_{V} s d v

La première intégrale est le flux sortant. Cette équation correspond au théorème de Gauss de l'électrostatique, à un bilan d'énergie pour un problème thermique.

On se limite ici à l'équation de Poisson bidimensionnelle en coordonnées cartésiennes. Dans ce cas, le gradient est dans le plan (x,y) et l'intégrale de surface se ramène à une intégrale sur une courbe fermée C; l'intégrale de volume se ramène à une intégrale sur la surface Σ délimitée par la courbe :

\int_{C} \vec{g r a d} (u) \cdot \vec{n} d l = \iint_{Σ} s d Σ

Voir aussi Discrétisation de l'équation de Poisson en coordonnées polaires et Discrétisation de l'équation de Poisson en géométrie axiale.

2. Maillage

Les nœuds du maillage sont définis par :

\begin{matrix} x_{i} = i Δ x \\ y_{j} = j Δ y \end{matrix}

L'indice i varie de 0 à N_x-1, l'indice j de 0 à N_y-1. On pose :

\begin{matrix} U_{i, j} = u (i Δ x, j Δ y) \\ S_{i, j} = s (i Δ x, j Δ y) \end{matrix}

Figure pleine page

3. Discrétisation de l'équation

La méthode des différences finies consiste à remplacer la dérivée seconde par une expression approchée à l'aide des valeurs au nœud et aux nœuds voisins. Pour la dérivée seconde par rapport à x, la différence finie la plus utilisée est :

\frac{U_{i + 1, j} - 2 U_{i, j} + U_{i - 1, j}}{Δ x^{2}}

En procédant de même pour l'autre dérivée, on obtient l'équation suivante pour chaque nœud du maillage :

Δ y^{2} (U_{i + 1, j} + U_{i - 1, j}) + Δ x^{2} (U_{i, j + 1} + U_{i, j - 1}) - 2 (Δ x^{2} + Δ y^{2}) U_{i, j} = Δ x^{2} Δ y^{2} S_{i, j}

Si U désigne la matrice colonne formée des N points du maillage, la relation précédente s'écrit sous forme matricielle :

A U = F

Il s'agit donc de résoudre un système linéaire à N inconnues.

On utilisera également la notation suivante :

D_{i, j} U_{i + 1, j} + G_{i, j} U_{i - 1, j} + H_{i, j} U_{i, j + 1} + B_{i, j} U_{i, j - 1} + C_{i, j} U_{i, j} = F_{i, j}

Les matrices D (droite), G (gauche), H (haut), B (bas), C (centre) et S (source) ont les mêmes dimensions que le maillage.

Pour un point non situé sur une frontière, on a donc :

\begin{matrix} D_{i, j} = Δ y^{2} \\ G_{i, j} = Δ y^{2} \\ H_{i, j} = Δ x^{2} \\ B_{i, j} = Δ x^{2} \\ C_{i, j} = - 2 (Δ x^{2} + Δ y^{2}) \\ F_{i, j} = S_{i, j} Δ x^{2} Δ y^{2} \end{matrix}

La discrétisation peut aussi être obtenue avec la méthode des volumes finis ([1]), qui consiste à discrétiser la forme intégrale exprimée pour un volume de contrôle. Pour un point situé dans le domaine (point A), le volume de contrôle est représenté sur la figure par le contour (a,b,c,d). Il s'agit d'un parallélépipède de côtés Δx,Δy,Δz. Le flux sortant du gradient sur la surface (ab) est discrétisé de la manière suivante :

\iint_{(a b)} \vec{g r a d} (u) \cdot \vec{n} d s ≃ \frac{U_{i, j + 1} - U_{i, j}}{Δ y} Δ x Δ z

On procède de même pour les surfaces (bc), (cd) et (da). L'intégrale de volume de la source est remplacé par :

S_{i, j} Δ x Δ y Δ z

4. Conditions limites

4.a. Condition de Dirichlet

Dans ce cas, le potentiel électrique, ou la température, est imposée sur une frontière. On écrira simplement pour tout point de la frontière :

\begin{matrix} D_{i, j} = 0 \\ G_{i, j} = 0 \\ H_{i, j} = 0 \\ B_{i, j} = 0 \\ C_{i, j} = 1 \\ F_{i, j} = U_{0} \end{matrix}

4.b. Condition de Neumann

Dans certains cas, le champ électrique peut être imposé sur une frontière. Par exemple sur une frontière perpendiculaire à l'axe Oy :

\frac{\partial V}{\partial y} = - E_{y}

Dans un problème de thermique, le flux thermique peut être imposé sur une frontière :

\frac{\partial T}{\partial y} = - \frac{φ_{y}}{λ}

On écrira donc la condition limite sous la forme générale :

\frac{\partial u}{\partial y} = w_{y}

Pour le point B situé sur le bord j=0, on écrit un développement de Taylor à l'ordre 2:

u (x_{i}, y_{j + 1}) ≃ u (x_{i}, y_{j}) + Δ y \frac{\partial u}{\partial y} (x_{i}, y_{j}) + \frac{Δ y^{2}}{2} \frac{\partial^{2} u}{\partial y^{2}} (x_{i}, y_{j})

La dérivée seconde par rapport à y est remplacée par une dérivée seconde par rapport à x à l'aide de l'équation de Poisson :

u (x_{i}, y_{j + 1}) ≃ u (x_{i}, y_{j}) + Δ y \frac{\partial u}{\partial y} (x_{i}, y_{j}) + \frac{Δ y^{2}}{2} (- \frac{\partial^{2} u}{\partial x^{2}} (x_{i}, y_{j}) + s (x_{i}, y_{j}))

Après avoir discrétisé la dérivée seconde comme plus haut, on obtient finalement pour le point de la frontière :

\begin{matrix} D_{i, j} = G_{i, j} = \frac{1}{2} Δ y^{2} \\ H_{i, j} = Δ x^{2} \\ B_{i, j} = 0 \\ C_{i, j} = - Δ x^{2} - Δ y^{2} \\ F_{i, j} = w_{y, i} Δ x^{2} Δ y + S_{i, j} \frac{1}{2} Δ x^{2} Δ y^{2} \end{matrix}

La même équation est obtenue plus simplement par la méthode des volumes finis, en appliquant la forme intégrale sur la surface (e,f,g,h) :

- w_{y, i} Δ x + \frac{U_{i + 1, j} - U_{i, j}}{Δ x} \frac{Δ y}{2} + \frac{U_{i - 1, j} - U_{i, j}}{Δ x} \frac{Δ y}{2} + \frac{U_{i, j + 1} - U_{i, j}}{Δ y} Δ x = \frac{1}{2} Δ x Δ y S_{i, j}

La condition limite de Neumann sur les points C et D s'obtient de manière analogue.

Pour le point F situé sur un coin, il y a une dérivée w_y imposée sur le côté parallèle à x, une dérivée w_y imposée sur le côté parallèle à y. La forme intégrale sur le volume de contrôle (x,y,z,w) conduit à :

- w_{x, j} \frac{1}{2} Δ y - w_{y, i} \frac{1}{2} Δ x + \frac{U_{i + 1, j} - U_{i, j}}{Δ x} \frac{Δ y}{2} + \frac{U_{i, j + 1} - U_{i, j}}{Δ y} \frac{Δ x}{2} = S_{i, j} \frac{1}{4} Δ x Δ y

On obtient donc pour le point F :

\begin{matrix} D_{i, j} = \frac{1}{2} Δ y^{2} \\ G_{i, j} = 0 \\ H_{i, j} = \frac{1}{2} Δ x^{2} \\ B_{i, j} = 0 \\ C_{i, j} = - \frac{1}{2} Δ x^{2} - \frac{1}{2} Δ y^{2} \\ F_{i, j} = w_{x, j} \frac{1}{2} Δ x Δ y^{2} + w_{y, j} \frac{1}{2} Δ x^{2} Δ y + S_{i, j} \frac{1}{4} Δ x^{2} Δ y^{2} \end{matrix}

La condition limite de Neumann sur le coin E s'obtient de la même manière. Les coins E et F sont similaires car le domaine de calcul se trouve à l'intérieur du coin. Il n'en est plus de même pour le coin G, sommet d'un objet situé dans le domaine. Le volume de contrôle pour ce point est (a,b,c,d,e,f). On obtient :

\begin{matrix} D_{i, j} = \frac{1}{2} Δ y^{2} \\ G_{i, j} = Δ y^{2} \\ H_{i, j} = \frac{1}{2} Δ x^{2} \\ B_{i, j} = Δ x^{2} \\ C_{i, j} = - \frac{3}{2} Δ x^{2} - \frac{3}{2} Δ y^{2} \\ F_{i, j} = - w_{x, j} \frac{1}{2} Δ x Δ y^{2} - \frac{1}{2} w_{y, i} Δ x^{2} Δ y + S_{i, j} \frac{3}{4} Δ x^{2} Δ y^{2} \end{matrix}

4.c. Relation de passage entre deux milieux différents

Sur une frontière entre deux milieux différents, il faut tenir compte des différences de conductivité dans un problème de thermique, ou des différences de constantes diélectriques ou magnétiques dans un problème d'électromagnétisme. Cela revient à écrire la forme intégrale de la manière suivante :

\int_{C} a \vec{g r a d} u \cdot \vec{n} d l = \iint_{Σ} s d Σ

où le coefficient a représente la propriété relative du milieu (conductivité, permittivité, etc), égale à a₁ ou a₂. On introduit ici deux constantes pour obtenir un traitement symétrique des deux milieux, mais une des deux constantes sera égale à 1.

En un point sur une frontière parallèle à y, la discrétisation de la forme intégrale s'écrit alors :

a_{1} \frac{U_{i - 1, j} - U_{i, j}}{Δ x} Δ y + a_{2} \frac{U_{i + 1, j} - U_{i, j}}{Δ x} Δ y + \frac{a_{1} + a_{2}}{2} \frac{U_{i, j - 1} - U_{i, j}}{Δ y} Δ y x + \frac{a_{1} + a_{2}}{2} \frac{U_{i, j + 1} - U_{i, j}}{Δ y} Δ x = Δ x Δ y S_{i, j}

On obtient finalement :

\begin{matrix} D_{i, j} = a_{2} Δ y^{2} \\ G_{i, j} = a_{1} Δ y^{2} \\ H_{i, j} = \frac{a_{1} + a_{2}}{2} Δ x^{2} \\ B_{i, j} = \frac{a_{1} + a_{2}}{2} Δ x^{2} \\ C_{i, j} = - (a_{1} + a_{2}) (Δ x^{2} + Δ y^{2}) \\ F_{i, j} = S_{i, j} Δ x^{2} Δ y^{2} \end{matrix}

Dans le cas $a_{1} = a_{2} = 1$ , on retrouve la discrétisation de l'équation de Poisson.

5. Méthode itérative de Gauss-Seidel

Dans une méthode itérative (ou méthode de relaxation), le système AU=F est résolu par itérations successives, jusqu'à convergence vers la solution. Au départ, les valeurs données aux $U_{i, j}$ sont quelconques (généralement nulles). À chaque itération, les N points du maillage sont parcourus dans un ordre déterminé, par exemple par rangées successives. Dans la méthode de Gauss-Seidel, les valeurs $U_{i, j}$ sont stockées dans un seul tableau. Le noeud courant est modifié par la relation :

U_{i, j}^{*} = \frac{1}{C_{i, j}} (F_{i, j} - D_{i, j} U_{i + 1, j} - G_{i, j} U_{i - 1, j} - H_{i, j} U_{i, j + 1} - B_{i, j} U_{i, j - 1})

Les noeuds du membre de droite sont soit des noeuds modifiés dans la précédente itération, soit des noeuds modifiés dans la même itération. Par exemple, dans le cas d'un balayage par lignes successives de gauche à droite et de bas en haut, les termes G et B ont été actualisés dans la même itération.

L'inconvénient de cette méthode est sa lenteur de convergence. Elle peut être accélérée par la méthode de sur-relaxation, qui consiste à adopter une moyenne pondérée de la valeur calculée par le schéma des différences finies et de l'ancienne :

U_{i, j}^{n e w} = ω U_{i, j}^{*} + (1 - ω) U_{i, j}^{o l d}

Le paramètre de sur-relaxation doit vérifier :

1 < ω < 2

Pour un problème donné, il faut chercher le paramètre optimal.

Pour plus de détails voir Méthode de Gauss Seidel .

Pour les très grands maillages, la convergence de la méthode de Gauss-Seidel peut être accélérée de manière spectaculaire par les méthodes multigrilles, qui consistent à utiliser successivement des maillages de différentes tailles lors du processus d'itération. Pour plus de détails, voir Méthode multigrilles.

Références

[1] R.H. Pletcher, J.C. Tannehill, D.A. Anderson, Computational Fluid Mechanics and Heat Transfer, (CRC Press, 2013)