5.3 Étude de la précision en éléments finis $\mathcal{P}^{1}$

Pour étudier la précision de l'approximation par éléments finis $\mathcal{P}^{1}$, nous choisissons un cas d'écoulement potentiel où l'on connaıt la solution analytique: l'écoulement autour d'un cylindre.

Soit un cylindre de rayon $a$, placé dans un écoulement uniforme de vitesse $U_{0}$ suivant x. L'écoulement potentiel autour du cylindre est donnée par la fonction de courant $\psi_{ex}$ en coordonnées polaires $(r,\theta)$:

$$\psi_{ex}(r,\theta)=U_{0}r \sin\theta-U_{0}\frac{a^{2}\sin\theta}{r}$$

qui vérifie l'équation de Laplace:

$$\Delta\Psi=0$$ (5.37)

et la condition aux limites de glissement sur le cylindre:

$$\psi(r=a)=0$$

Figure 6.17: écoulement potentiel autour d'un cylindre

Considérons le domaine $\Omega$ de la figure (6.17) de frontière ABCD (soit 1/4 du domaine complet). Sur ce domaine les conditions aux limites sont:

1. en entrée AD $(r=R)$ : $\psi_{DA}=U_{0}R \sin\theta-U_{0}\frac{a^{2}\sin\theta}{R}$,
2. sur l'axe AB $(y=0)$, on a une condition de symétrie $(v=0)$, et l'axe AB est la ligne de courant qui arrive sur l'obstacle. La valeur de $\psi$ vaut donc $\psi(A)$, soit $\psi_{AB}=0,$
3. sur l'obstacle BC $(r=a)$, la valeur de $\psi$ est constante et égale à celle sur AB, soit $\psi_{BC}=0$,
4. sur la sortie CD, on impose une condition de symétrie $(\frac{\partial\psi}{\partial n})_{CD}=(\frac{\partial\psi}{\partial x})_{x=0}=0.$

5.3.1 Maillage $\mathcal{P}^{1}$

Pour générer le maillage du domaine $\Omega$, nous utiliserons les programmes Matlab de génération de maillage décrits en annexe. Plus précisément, nous utiliserons la fonction quacou qui maille un domaine à l'aide d'une transformation “quadrilatères courbes” des frontières vers un carré unité que l'on maille (figure 6.18).

Figure 6.18: maillage par transformation quacou

On spécifie les points des cotés AB, BC , CD et DA en notant que le nombre de noeuds sur les faces opposées doivent être égaux. Le listing 6.3.1 ci dessous donne la fonction Matlab cylindremesh utilisée pour le maillage de $\Omega$, paramètré par le nombre de noeuds $n_{1}$et $n_{2}$ sur les cotés AB et BC.

programme matlab 6.3.1: Maillage éléments finis 2D $\mathcal{P}^1$: cylindremesh.m

function [G]=cylindremesh(n1,n2)
% maillage du cylindre
addpath ../Mailleur
% construction de la geometrie
R1=1; R2=3;
A=[-R2,0]; B=[-R1,0]; C=[0,R1]; D=[0,R2];
AB=[linspace(A(1),B(1),n1)',linspace(A(2),B(2),n1)'];
theta=linspace(pi,pi/2,n2)';
BC=[R1*cos(theta),R1*sin(theta)];
CD=[linspace(C(1),D(1),n1)',linspace(C(2),D(2),n1)'];
theta=linspace(pi/2,pi,n2)';
DA=[R2*cos(theta),R2*sin(theta)];
G=Quacou(AB,BC,CD,DA,[1,1,3,2]);
return;

Le maillage obtenu pour $n_{1}=20$ et $n_{2}=20$ est tracé sur la figure 6.19. Ce maillage comprend 400 noeuds et 722 éléments.

Figure 6.19: maillage du cylindre

5.3.2 Résolution

En utilisant les programmes précédents, on construit la matrice A et le second membre B du problème discrétisé par éléments finis. Avec l'approche utilisée, la matrice $A$ est une matrice $(nn,nn)$ qui contient $160000$ éléments dans le cas du maillage de la figure (6.19). Or la matrice A a une structure creuse, c'est à dire qu'il y a très peu d'éléments non nuls. Dans notre cas A possède uniquement $2412$ éléments non nuls, soit moins de $2\%$.

Exercice: montrer qu'en moyenne pour des éléments $\mathcal{P}^{1}$, le nombre d'éléments non nuls de A est $7nn$

Figure: Éléments non nuls de A, $\mathbf{A}^{-1}$ et $\mathbf{L}.\mathbf{U}$

Au lieu de stocker la matrice A sous la forme d'un tableau à 2 indices, on utilise un stockage creux (sparse en anglais) pour lequel on ne conserve que les éléments non nuls. Pour cela on ne stocke que les valeurs $A_{i,j}$ et les indices $(i,j)$ des éléments non nuls. Considérons par exemple la matrice creuse $M$ suivante, qui possède 6 éléments non nuls sur 16:

$$\mathbf{M}=\left[\begin{array}{cccc} 2.0 & 0 & 0 & -1.0\\ 0... ...& 0\\ 0 & -3.0 & 1.0 & 0\\ 0 & 0 & 0 & -2.0\end{array}\right]$$

le stockage creux de cette matrice correspond au tableau de 6 lignes et 3 colonnes suivant:

$M_{ij}$

i

j 2.0

Dans le cas de la matrice A précédente ce stockage creux est très efficace puisque que l'on utilise un tableau de $3*2412$ éléments au lieu d'un tableau de $160 000$ éléments. Par contre si la matrice A est une matrice creuse, son inverse $\mathbf{A}^{-1}$ est en général une matrice pleine comme le montre la figure 6.20, où on a marqué par un point bleu les éléments non nuls de A et de $\mathbf{A}^{-1}$ . Pour la résolution d'un système linéaire linéaire avec une matrice creuse, il faut donc utiliser des méthodes numériques adaptées. Nous avons choisi d'utiliser une factorisation LU de la matrice A qui préserve au maximum le caractère creux de la matrice. Comme le montre la figure 6.20, cette factorisation génère $2*7 498$ éléments non nuls au lieu des $144 475$ de $\mathbf{A}^{-1}$ (soit $10\%$). Cette factorisation est relativement efficace puisque que la numérotation des noeuds est telle que les éléments non nuls sont proches de la diagonale (voir figure 6.20) . Dans le cas d'une numérotation quelconque, il faut appliquer une renumérotation des noeuds de façon à minimiser le remplissage de la matrice lors de la factorisation (par un algorithme de Cuthill-McKee ou un algorithme de degré minimum, tous les deux implémenter sous Matlab par les fonctions colrcm et colamd).

On peut encore améliorer le stockage en tenant compte du caractère symétrique de la matrice A et en ne stockant que la partie triangulaire supérieure de A. Dans ce cas, il faut symétriser la prise en compte des conditions aux limites de Dirichlet: pour un noeud i tel que $\psi_{i}=\psi^{e}$, il modifier la ligne i et la colonne i de A ainsi que le second membre B:

$$\begin{eqnarray*} \mathbf{B}\leftarrow\mathbf{B}-\psi^{e}[A_{i,k}]_{k=1,nn} & B_{1}=\psi^{e}\\ A_{i,k}=A_{k,i}=0 & (k=1,nn)\\ A_{i,i=1}\end{eqnarray*}$$

Exercice: modifier les programmes Matlab pour prendre en compte cette symétrie.

Pour utiliser un stockage “matrice creuse”, nous avons modifier le programme Matlab d'assemblage en déclarant la matrice A sparse (ligne 7 dans le programme 6.3.2 ci dessous).

programme matlab 6.3.2: Assemblage avec une matrice creuse: Assemblage.m

function [A,B]=Assemblage(G,K,alpha,F)
% assemblage de la matrice et du second membre
% K coef. diffusion, alpha coef. source, F terme source
% stocage de A sous forme creuse
% estimation du nbre d'element non nuls nzmax
nzmax=2*G.ddl*G.ne+G.nn;
A=sparse([],[],[],G.nn,G.nn,nzmax);
B=zeros(G.nn,1);
for k=1:G.ne
    n=G.Tbc(k,:); % numero des sommets
    Ke=MatriceRigidite(G,k);
    Me=MatriceMasse(G,k);
    Be=SmbElement(G,F,k);
    A(n,n)=A(n,n)+K*Ke+alpha*Me;
    B(n)=B(n)+Be;
end;
return;

Pour cela on fournit une estimation du nombre d'éléments non nuls de A à partir d'un nombre d'arêtes par élément (ligne 6).

Exercice: montrez que le coefficient $A_{ij}$ est non nul si $ij$ est une arête du maillage.

Sur chaque élément on a 3 cotés donc 6 arêtes , soit au total une estimation de $6ne$ arêtes et donc $6ne+nn$ éléments non nuls. Cette estimation est grossière puisque l'on compte deux fois les arêtes internes. Une estimation exacte nécessiterait de déterminer les arêtes de la géométrie. Dans notre cas le calcul exacte donne $2242$ arêtes pour une estimation de $4332,$ soit un nombre d'éléments non nuls estimé de $4732$ pour $2642$.

Exercice: écrire l'algorithme et le programme Matlab permettant de déterminer les arêtes de la géométrie.

C'est la seule modification apportée au programme d'assemblage, puisque les calculs des matrices et second membres élémentaires sont inchangés.

Avec ces modifications, la solution approchée $\psi ^{h}$ s'obtiens par résolution du système linéaire par factorisation LU. Le script Matlab qui enchaıne la suite des opérations pour calculer cette solution approchée est donné ci-dessous.

programme matlab 6.3.2: Résolution par éléments finis 2D $\mathcal{P}^1$

% maillage du cylindre
G=cylindremesh(20,20);
% second membre
F=zeros(G.nn,1);
% assemblage
[A,B]=Assemblage(G,1.0,0.0,F);
% conditions aux limites
Psie=inline('y-y./(x.^2+y.^2)','x','y');
[A,B]=Climite(G,A,B,0.0,0.0,Psie);
% resolution
[LA,UA]=lu(A);
Psi=UA\(LA\B);
% calcul de l'erreur
err=Erreur(G,Psie,Psi)

c4prog11

Le solution obtenue est tracé sur la figure (6.21).

Figure 6.21: solution $\psi ^{h}$ par éléments finis $\mathcal{P}^{1}$ avec le maillage de la figure 6.19

5.3.3 Intégration numérique

Pour déterminer la précision de cette solution numérique $\psi ^{h}$, il faut calculer la norme de l'erreur entre la solution exacte $\psi_{ex}$ et la solution approchée $\psi ^{h}$:

$$\Vert\psi_{ex}-\psi^{h}\Vert=\sqrt{\int_{\Omega}(\psi_{ex}-\psi^{h})^{2}d\omega}$$

Le calcul de l'intégrale se fait par assemblage, avec un calcul élément par élément:

$$\int_{\Omega}(\psi_{ex}-\psi^{h})^{2}d\omega=\sum_{k=1}^{ne}\int_{e_{k}}(\psi_{ex}-\psi^{h})^{2}d\omega$$

et chaque intégrale élémentaire est calculé par transformation sur l'élément de référence:

$$\int_{e_{k}}(\psi_{ex}-\psi^{h})^{2}d\omega=\det(\mathbf{J}_... ...m_{p=1}^{3}\psi_{n_{p}}N_{p}(\xi,\eta)\right)^{2} d\eta d\xi$$ (5.38)

La fonction $\psi_{ex}$ n'étant pas polynômiale, on ne sait pas calculer analytiquement cette intégrale. On choisit une intégration numérique par points de Gauss. On approxime l' intégrale d'une fonction $f(\xi,\eta)$ sur le triangle de référence par une formule de quadrature de Gauss: i.e par une somme pondérée (de poids $w_{i})$ de valeurs de $f(\xi,\eta)$ en $ng$ points de Gauss de coordonnées $(\xi_{i},\eta_{i})$:

$$\int_{0}^{1}\int_{0}^{1-\xi}f(\xi,\eta) d\eta d\xi\simeq \sum_{i=1}^{ng}w_{i} f(\xi_{i},\eta_{i})$$

Les poids $w_{i}$ et les coordonnées $(\xi_{i},\eta_{i})$ des points de Gauss sont tels que la formule de quadrature soit exacte pour des polynômes de degré le plus élevé possible. Cette formule possède $3ng$ degrés de liberté, on peux donc imposer $3ng$ contraintes. Ainsi pour intégrer exactement des polynômes de degré 1 en $(\xi,\eta)$ (ce qui impose 3 conditions) , il suffit d'un seul point de Gauss, l'intersection des médianes:

$$\int_{0}^{1}\int_{0}^{1-\xi}f(\xi,\eta) d\eta d\xi\simeq\frac{1}{2}f(\frac{1}{3},\frac{1}{3})$$

Exercice: démontrer cette formule de quadrature.

Pour un triangle, les points de Gauss doivent vérifier des propriétés de symétrie: i.e. leurs coordonnées barycentriques sont obtenues par permutations circulaires et les poids sont équivalents. Les points de Gauss sont donc à une même distance sur les médianes du triangle.

Ainsi la formule de quadrature utilisant 3 points de Gauss sur les médianes s'écrit:

$$\int_{0}^{1}\int_{0}^{1-\xi}f(\xi,\eta) d\eta d\xi\simeq w... ...,\xi_{1})+w_{1}f(1-2\xi_{1},\xi_{1})+w_{1}f(\xi_{1},1-2\xi_{1})$$

En écrivant que cette formule doit être exacte pour les polynômes de degré $\le2$, i.e. pour $f=1, \xi, \eta, \xi\eta, \xi^{2}, \eta^{2}$, on en déduit les relations (non linéaires) que doivent vérifier $w_{1}$ et $\xi_{1}$.Parmi ces relations, deux sont indépendantes, ce qui permet de déterminer $w_{1}$ et $\xi_{1}$.

Le programme Maple 6.3.3 ci-dessous détermine ces relations et résout le problème non-linéaire. Dans ce cas on trouve 2 solutions:

$$w_{1}=\frac{1}{6}, \xi_{1}=\frac{1}{2} \mbox{ ou }w_{1}=\frac{1}{6}, \xi_{1}=\frac{1}{6}$$

qui permette d'intégrer exactement des polynômes de degré $\le2$. La seconde solution est la plus précise, puisque l'erreur d'intégration pour des polynômes de degré 3 est la plus faible.

programme Maple 6.3.3: détermination des points de Gauss sur un triangle

> restart:
# Formule de quadrature sur un triangle avec 3 pts
> F:='F':Ie:=int(int(F(chi,eta),eta=0..1-chi),chi=0..1);
> F:='F':Ih:=w1*F(xi1,eta1)+w2*F(xi2,eta2)+w3*F(xi3,eta3);
# Cdts de symetrie
> w2:=w1; w3:=w1; eta1:=xi1; xi2:=1-2*xi1; eta2:=xi1; xi3:=xi1;
> eta3:=1-2*xi1;
> Ih;
# 6 conditions
> F:=(xi,eta)->1;eval(Ih-Ie=0);rel1:=%:
> w1:=solve(rel1,w1);
> F:=(xi,eta)->xi;eval(Ih-Ie=0);rel2:=%:
> F:=(xi,eta)->eta;eval(Ih-Ie=0);rel3:=%:
> F:=(xi,eta)->xi*eta;expand(eval(Ih-Ie=0));rel4:=%:
> F:=(xi,eta)->xi^2;expand(eval(Ih-Ie=0));rel5:=%:
> F:=(xi,eta)->eta^2;expand(eval(Ih-Ie=0));rel6:=%:
# Une seule relation indépendante: 2 solutions
> solve(rel6,xi1);
> xi1:=1/2; F:='F': Ih;
> F:=(xi,eta)->xi^3;eval(Ih-Ie=0);
> xi1:=1/6; F:='F': Ih;
> F:=(xi,eta)->xi^3;eval(Ih-Ie=0);

Avec cette approche, on peut déterminer les points de Gauss pour des intégrations exactes de polynômes. Dans le tableau (6.2) sont donnés les formules de quadrature de Gauss précises de l'ordre 1 à 5. Les paramètres exactes de ces formules de Gauss, obtenus avec Maple en suivant la même démarche que précédemment, sont donnés dans le tableau (6.3), ainsi que leurs valeurs numériques approchées^5.3..

Tableau 6.2: points de Gauss pour l'intégration numérique sur un triangle

ordre de précision

formule

figure

poids

$\xi_{i}$

$\eta_{i}$ 1

On constate que pour l'ordre 3 (i.e. avec 4 points de Gauss), le poids $p_{2}$ associé au noeud central est négatif. La formule de quadrature ne préserve donc pas la positivité de l'intégrale lorsque l'intégrante $f(\xi,\eta)$ est positif: i.e.

$$f(\xi,\eta)\ge0\Longrightarrow\int_{0}^{1}\int_{0}^{1-\xi}f(\xi \eta) d\eta d\xi>0$$

Cette formule n'est donc pas utilisable en pratique et seul les ordres 2,4 ou 5 sont utilisés.

Tableau 6.3: paramètres pour l'intégration numérique sur un triangle

ordre

paramêtres

valeur approchée 1

Pour terminer sur l'intégration sur le triangle de référence, nous donnons les formules d'intégration exacte pour les polynômes. Un calcul directe montre que pour tous les entiers positifs n et m on a :

$$\int_{0}^{1}\int_{0}^{1-\xi}\xi^{m}\eta^{n} d\eta d\xi=\frac{m!n!}{(n+m+2)!}$$

d'où l'on déduit la formule d'intégration pour les fonctions de forme:

$$\int_{0}^{1}\int_{0}^{1-\xi}N_{1}^{l}N_{2}^{m}N_{3}^{n} d\eta d\xi=\frac{l!m!n!}{(l+n+m+2)!}$$

Exercice: démontrer cette formule (à l'aide de Maple).

Pour calculer l'intégrale élémentaire dans l'expression de l'erreur (6.38), nous avons utilisé une intégration numérique d'ordre 2 avec 3 points de Gauss. La fonction Matlab Erreur 6.3.3 correspondante qui calcul l'erreur globale (6.38),est donnée ci dessous.

programme matlab 6.3.3: calcul de l'erreur d'intégration : Erreur.m

function [err]=Erreur(G,f,F)
% calcul de l'erreur en norme L2
% pts de Gauss
Xsi=[1.0/6.0,2.0/3.0,1.0/6.0];
Eta=[1.0/6.0,1.0/6.0,2.0/3.0];
P=[1-Xsi-Eta;Xsi;Eta];
err=0.0;
for k=1:G.ne
    % calcul de l'erreur / elt
    n=G.Tbc(k,:); % numero des sommets de l'element k
    S21=G.X(n(1),:)-G.X(n(2),:);
    S13=G.X(n(3),:)-G.X(n(1),:);
    Aire=0.5*(S13(1)*S21(2)-S13(2)*S21(1));
    % calcul des coordonnees des pts de Gauss
    X=P*G.X(n,:);
    % erreur sur l'elt k
    errk=(f(X(:,1),X(:,2))-P*F(n));
    err=err+errk'*errk*Aire/3;
end;
err=sqrt(err);
return;

5.3.4 Détermination de la précision

La précision d'une solution par éléments finis dépend du degré de l'élément et de la taille $h$ des éléments. Pour un triangle $e_{k}$, on définit la taille $h_{k}$ comme la longueur du plus grand coté (figure6.22 ). On définit aussi le diamètre $d_{i}$ du cercle inscrit et le diamètre $d_{e}$ du cercle circonscrit, ainsi que le rapport d'aspect du triangle $R=\frac{d_{i}}{d_{e}}$. Pour un triangle équilatéral, on a $d_{e}=h$ et $d_{i}=\frac{\sqrt{2}}{3}h$, donc $R=\frac{\sqrt{2}}{3}$. Par contre pour un triangle de plus en plus aplati, le rapport d'aspect tend vers 0 (car $d_{i}\rightarrow0$) , ce qui indique la dégénérescence du triangle.

Figure 6.22: longueurs caractéristiques d'un triangle

Avec ces définitions, on montre que l'erreur d'interpolation $\mathcal{P}^{1}$ sur un triangle vérifie une relation analogue au cas 1D:

$$\left\Vert f-f^{h}\right\Vert =\sqrt{\int_{e_{k}}(f-f^{h})^{... ...2}})^{2}+(\frac{\partial^{2}f}{\partial y^{2}})^{2}\right)dxdy}$$

Cependant la constante $C$ dépend du rapport d'aspect $R$ et l'estimation d'erreur dégénère lorsque $R\rightarrow0$, i.e. lorsque le triangle devient de plus en plus aplatis.

L'erreur d'approximation étant borné par l'erreur d'interpolation, on déduit une estimation de l'erreur d'approximation $\left\Vert \psi_{ex}-\psi^{h}\right\Vert$ sur un maillage non dégénéré de taille caractéristique $h$:

$$\left\Vert \psi_{ex}-\psi^{h}\right\Vert =\sqrt{\int_{\Omega... ...)^{2}+(\frac{\partial^{2}f}{\partial y^{2}})^{2}\right)d\omega}$$

Figure: Norme de l'erreur d'approximation $\mathcal{P}^{1}$ en fonction de $n\propto h^{-1}$

Nous avons donc construit une série de maillage de $\Omega,$ avec une taille $h\rightarrow0$ et un rapport d'aspect $R$ constant en fixant un nombre de noeuds identiques sur les cotés lors de la génération du maillage: i.e. $n_{1}=n_{2}=n$. Dans ce cas la taille des éléments $h$ est proportionnelle à $\frac{1}{n}$. On a donc tracé sur la figure 6.23, la norme de l'erreur en fonction de $n$ et on vérifie que l'on obtiens une décroissance de l'erreur en $n^{-2}$, soit en $\theta(h^{2})$, ce qui montre que l'approximation par éléments finis $\mathcal{P}^{1}$ est d'ordre 2.