5.2 Éléments finis triangulaire $\mathcal{P}^{1}$

Après l'étude du paragraphe précédent, nous allons maintenant étudier l'approximation par éléments finis $\mathcal{P}^{1}$ d'un problème générique sur un domaine bidimensionnel $\Omega$ quelconque (figure 6.12) .

Figure 6.12: domaine $\Omega$ de calcul

La frontière $\Gamma=\partial\Omega$ se décompose en 4 sous frontières disjointes (et éventuellement nulle):

$$\Gamma=\Gamma_{1}\cup\Gamma_{2}\cup\Gamma_{3}\cup\Gamma_{4}$$

Chacune de ses frontières correspond à un type de conditions aux limites appliquées à la solution $u(x,y):$

1. $\Gamma_{1}$ est une frontière de type Dirichlet homogène: i.e. $u_{\Gamma_{1}}=0,$
2. $\Gamma_{2}$ est une frontière de type Dirichlet non-homogène: i.e. $u_{\Gamma_{2}}=u_{e},$
3. $\Gamma_{3}$ est une frontière de type Neuman homogène: i.e. $(\frac{\partial u}{\partial n})_{\Gamma_{3}}=0,$
4. $\Gamma_{4}$ est une frontière de type Fourier: i.e. $-K (\frac{\partial u}{\partial n})_{\Gamma_{4}}=\beta u_{\Gamma_{4}}+\phi_{0}$.

Le problème générique que nous allons étudier s'écrit alors:

$$\left\{ \begin{array}{cc} -\frac{\partial}{\partial x}\left(... ...)_{\Gamma_{4}}=\beta u_{\Gamma_{4}}+\phi_{0}\end{array}\right.$$ (5.32)

L'équation d'équilibre est une équation de convectionconduction avec un terme source fonction de la solution. $K$ représente le coefficient de diffusion $(K>0)$, $\alpha$ le coefficient du terme source $(\alpha\ge0)$ linéaire en $u(x,y)$, et $f(x,y)$ le terme source indépendant de $u(x,y)$.

5.2.1 Formulation faible

La formulation faible de (6.32) s'obtiens comme toujours en multipliant par une fonction test $v(x,y),$ puis en intégrant sur le domaine $\Omega$. On intègre ensuite par partie le terme en dérivée seconde en utilisant la formule de Green. En utilisant les conditions aux limites et en interprétant la fonction test $v(x,y)$ comme une variation $\delta u$ de la solution $u(x,y)$, la formulation faible s'écrit:

$$\left\{ \begin{array}{cc} \mbox{Trouver } u(x,y) \mbox{... .... } v_{\Gamma_{1}}=0, v_{\Gamma_{2}}=0\end{array}\right.$$ (5.33)

Pour calculer la solution approchée $u^{h}$, nous avons besoin d'un maillage de $\Omega$, et d'une approximation de la solution sur ce maillage.

5.2.2 Maillage éléments finis et interpolation $\mathcal{P}^{1}$

Le domaine de calcul $\Omega$ est découpé en triangles comme sur la figure (6.13). Ce maillage est construit à l'aide d'outils informatiques ou de logiciels spécialisés. On trouvera en annexe les outils de maillage écrit avec Matlab, qui nous ont servi à générer les différents maillages utilisés dans ce livre.

Figure: Maillage éléments finis $\mathcal{P}^{1}$

De façon général, un outil de maillage doit fournir les informations suivantes:

1. le nombre de noeuds nn, le nombre d'éléments ne,
2. les coordonnées $(x_{i},y_{i})$ de chaque noeud $M_{i}$ du maillage,
3. les numéros des sommets de chaque élément $e_{k}$ (ou table de connection): $tbc_{k,1}, tbc_{k,2}\„tbc_{k,3}$
4. pour chaque noeud $M_{i}$, une information (i.e. un entier) $frt_{i}$ qui précise si le noeud est interne ($frt_{i}=0$), ou sur une frontière $\Gamma_{l}$ ($frt_{i}=l$)
5. pour chaque élément $e_{k}$, une information de région $reg_{k}$, indiquant le numéro du domaine auquel appartiens l'élément. Par défaut il n'y a qu'un seul domaine et $reg_{k}=1$ pour tous les éléments.

Avec les outils utilisés, ses informations sont stockés sous forme de tableau dans un fichier de géométrie, ayant par convention une extension .msh. Ce fichier est écrit en ASCII, et est donc lisible sur n'importe quel système et avec n'importe quel langage. Il est compatible avec le format utilisé par FREEFEM^5.2, un logiciel éléments finis du domaine public. Le fichier de maillage est structuré par ligne de la façon suivante:

ligne

champ 1

champ 2

champ 3

champ 4 1

Le fichier de maillage de la figure 6.13 contiens par exemple les valeurs suivantes :

30 40 maillage_triangle_P1
0.000000 0.000000 2
0.550000 -0.125000 3
1.100000 -0.250000 3
...................
0.850000 0.900000 4
1.175000 0.950000 4
1.500000 1.000000 1
1 2 7 1
6 1 7 1
.........
24 25 30 1
29 24 30 1

Sous Matlab, nous utiliserons une structure de données pour la géométrie de type structure, qui contiendra les champs suivants (en notant G le nom variable géométrie):

• [G.nn] nbre de noeuds de la géométrie G : entier
• [G.ne] nbre d'éléments de la géométrie G : entier
• [G.dim] la dimension d'espace (=2 en 2D) : entier
• [G.ddl] nbre de degré de liberté par éléments (=3 pour des éléments $\mathcal{P}^{1}$ en 2D) : entier
• [G.X] tableau des coordonnées $(x,y)$ des noeuds: tableau G.nn*G.dim réels
• [G.Tbc] table de connection: tableau G.ne*G.ddl entiers
• [G.Frt] numéro de frontière par noeuds: tableau G.nn entiers
• [G.Reg] numéro de région par éléments: tableau G.ne entiers

La fonction Matlab Lecture.m (programme 6.2.2) lit une géométrie éléments finis sur fichier, et initialise la variable structure géométrie G.

programme matlab 6.2.2: Lecture de la géometrie: Lecture.m

function [G]=Lecture(fichier)
% lecture du maillage 
% ouverture du fichier
fid=fopen(fichier,'r');
% lecture de la dimension 
[L,count]=fscanf(fid,'%d %d %s',3);
G.nn=L(1); G.ne=L(2);
% initialisation structure geometrie
G.dim=2; G.ddl=3;
G.X=zeros(G.nn,G.dim); G.Frt=zeros(G.nn,1);
G.Tbc=zeros(G.ne,G.ddl); G.Reg=zeros(G.ne,1);
% lecture des coordonnees
for i=1:G.nn
    [L,count]=fscanf(fid,'%f %f %d',3);
    G.X(i,1:2)=[L(1) L(2)];
    G.Frt(i)=L(3);
end;
% table deconnection
for l=1:G.ne
    [L,count]=fscanf(fid,'%d %d %d %d',4);
    G.Tbc(l,1:3)=[L(1) L(2) L(3)];
    G.Reg(l)=L(4);
end;
fclose(fid);
return;

Ensuite sous Matlab, pour lire la géométrie stockée dans le fichier maillage.msh, on écrit simplement:

>>  G1=Lecture('maillage.msh')

G1 =

     nn: 30
     ne: 40
    dim: 2
    ddl: 3
      X: [30x2 double]
    Frt: [30x1 double]
    Tbc: [40x3 double]
    Reg: [40x1 double]

5.2.2.1 Approximation éléments finis

On a montré précédemment qu'une approximation $u^{h}(x,y)$ sur un maillage éléments finis est une combinaison linéaire des valeurs nodales $\{ u_{i}\}_{i=1,nn}$ sur les noeuds du maillage. les coefficients de cette combinaison linéaire sont les fonctions de base $\Phi_{i}(x,y)$ :

$$u^{h}(x,y)=\sum_{i=1}^{nn}u_{i} \Phi_{i}(x,y) \mbox{ avec } u_{i}=u^{h}(x_{i},y_{i})$$

Ces fonctions $\{\Phi_{i}\}$ forment une base locale, i.e. une fonction $\Phi_{i}(x,y)$ est non nulle uniquement sur une petite partie du maillage: le support du noeud $M_{i}$, i.e. l'ensemble des éléments $e_{l}$ ayant le noeud $M_{i}$ comme sommet ( $M_{i}\in e_{l}$). Elles vérifient en outre les propriétés suivantes:

1. la fonction de base $\Phi_{i}(x,y)$ vaut 1 au noeud i et 0 sur tous les autres noeuds:
   
   
   $$\Phi_{i}(x_{j},y_{j})=\delta_{ij}$$
   
   

2. la fonction de base $\Phi_{i}(x,y)$ est non nulle uniquement sur son support $Sup_{i}=\{\cup e_{k} /M_{i}\in e_{k}\}$
   
   
   $$N_{i}(x,y)\neq0 { si } (x,y)\in Sup_{i}$$
   
   

3. la fonction de base $\Phi_{i}(x,y)$ est orthogonale à presque toutes les autres fonctions de base $\Phi _{j}$, en particulier celles qui sont associées à des noeuds j non voisin de i (t.q. $Sup_{i}\cap Sup_{j}=\emptyset$) :
   
   
   $$\Phi_{i}(x,y)*\Phi_{j}(x,y) \left\{ \begin{array}{cc} =0 & ... ...mbox{ si } (x,y)\in(Sup_{i}\cap Sup_{j})\end{array}\right.$$
   
   

Un exemple de fonctions de base (associées aux noeuds $7,12,24$ du maillage de la figure (6.13)) est tracé sur la figure (6.14), ainsi que leurs supports. On constate sur cette figure que les fonctions de base $\Phi_{7}$ et $\Phi_{24}$ sont orthogonales (i.e. le produit $\Phi_{7}*\Phi_{24}$ est toujours nul), ainsi que $\Phi_{12}$ et $\Phi_{24}$ , et que le produit $\Phi_{7}*\Phi_{12}$ est non nul uniquement sur les 2 éléments ($e_{9},e_{12}$).

Figure 6.14: fonctions de base $\Phi _{7},\Phi _{12},\Phi _{24}$ et leurs supports

Pour calculer une fonction de base sur un élément $e_{k}$ du maillage, on utilise une transformation (6.10,6.11,6.13) de cet élément vers un élément de référence $\hat{e}$ . Sur cet élément de référence, la fonction de base associée à l'un des sommets $\{ S_{q}\}_{q=1,3}$ de l'élément $e_{k}$ coıncide alors avec l'une des fonctions de forme $N_{q}(\xi,\eta)$ (6.12) associées aux sommets $\{\hat{S}_{q}\}$ de l'élément de référence $\hat{e}$ (6.18).

5.2.3 Formulation faible discrète

La solution approchée $u^{h}$ de l'équation (6.32) s'écrit sous la forme:

$$u^{h}(x,y)=\sum_{i=1}^{nn}u_{i} \Phi_{i}(x,y)$$

Elle doit de plus vérifier les conditions aux limites de Dirichlet, i.e. $u_{\Gamma_{1}}^{h}=0, u_{\Gamma_{2}}^{h}=u_{e}$. Ces conditions de Dirichlet imposent la valeur nodale de $u^{h}$ sur tous les noeuds se trouvant sur la frontière $\Gamma_{1}$ ou $\Gamma_{2}$, i.e:

$$u_{l}=0 \mbox{ si } M_{l}\in\Gamma_{1}, u_{l}=u_{e} \mbox{ si } M_{l}\in\Gamma_{2}$$

Soit $nf_{1}$ le nombre de noeuds sur la frontière $\Gamma_{1}$, et $nf_{2}$ le nombre de noeuds sur la frontière $\Gamma_{2}$, le nombre de degré de liberté $nl$ de $u^{h}$ est égale à $nl=nn-nf_{1}-nf_{2}$. Cependant la numérotation du maillage étant quelconque, on ne peut comme dans l'exemple précédent avoir une numérotation simple des degrés de liberté. En notant $\overline{\Omega}_{d}=\Omega-\Gamma_{1}-\Gamma_{2}$, on écrit la solution sous la forme:

$$u^{h}(x,y)=\sum_{j\in\overline{\Omega}_{d}}u_{j} \Phi_{j}(x,y)+\sum_{l\in\Gamma_{2}}u_{e} \Phi_{l}(x,y)$$ (5.34)

Les fonctions tests associées $v^{h}(x,y)$ sont des variations de $u^{h}$, qui doivent donc s'annuller sur les frontières de Dirichlet $\Gamma_{1}$ et $\Gamma_{2}$. Elles s'écrivent sous la forme générale suivante:

$$v^{h}(x,y)=\sum_{i\in\overline{\Omega}_{d}}v_{i} \Phi_{i}(x,y)$$

Ces fonctions tests sont une combinaison linéaire des $n_{l}$ fonctions de base $\Phi _{i}$ ( $i\in\overline{\Omega}_{d}$). En remplaçant la solution exacte $u$ par l'expression (6.34) de $u^{h}$ et la fonction test $v$ par une des fonctions de base précédentes $\Phi _{i}$, on obtiens la formulation faible discrète:

$$\sum_{j\in\overline{\Omega}_{d}}u_{j}\left(\underbrace{\int_{\Ome... ...ma_{4}}\beta\Phi_{j}.\Phi_{i} d\gamma}_{\mbox{terme C.L. sur} A}\right) =$$

$$\underbrace{\int_{\Omega}f\Phi_{i} d\omega}_{\mbox{terme géné... ...ce{\int_{\Gamma_{4}}\phi_{0}\Phi_{i} d\gamma}_{\mbox{terme C.L. sur} B} -$$

$$\sum_{l\in\Gamma_{2}}u_{l}\underbrace{\left(\int_{\Omega}K(\frac{... ...ma_{4}}\beta\Phi_{l}.\Phi_{i} d\gamma\right)}_{\mbox{terme C.L. sur }B}$$ (5.35)

La relation (6.35) écrite pour les $nl$ fonctions de base $N_{i}$ appartenant à $\overline{\Omega}_{d}$, est un système linéaire de $nl$ inconnues $u_{j}$ ( $j\in\overline{\Omega}_{d}$), qu'il suffit de résoudre pour obtenir la solution approchée $u^{h}$. Le calcul directe de ce système linéaire nécessiterait une renumérotation des noeuds, de façon à avoir la même numérotation pour les inconnues et degrés de liberté (i.e. de 1 à $nl$ pour les noeuds $j\in\overline{\Omega}_{d}$), et ce chaque fois que l'on change les conditions aux limites du problème. De façon à écrire un programme qui soit le plus général possible, on construira tout d'abord un système linéaire générique pour toutes les valeurs nodales $\{ u_{i}\}_{i=1,nn}$, puis on appliquera les conditions aux limites du problème sur ce système linéaire de dimension $nn$.

Les coefficients génériques de la matrice $\mathbf{A}$ et du second membre $\mathbf{B}$ s'écrivent:

$$\begin{eqnarray*} \mathbf{A}_{ij} & = & \int_{\Omega}K(\frac{\partial\Phi_{j}}{\... ... d\omega\\ \mathbf{B}_{i} & = & \int_{\Omega}f\Phi_{i} d\omega\end{eqnarray*}$$

Pour calculer ces coefficients, on effectue un calcul élément par élément en déterminant les matrices et les second membres élémentaires sur chaque élément $e_{k}$:

$$\begin{eqnarray*} \mathbf{A}_{ij} & = & \sum_{k=1}^{ne}\left(\underbrace{\int_{e... ...\Phi_{i} d\omega}_{\mbox{second membre élémentaire}}\right)\end{eqnarray*}$$

5.2.4 Calcul des matrices élémentaires $\mathcal{P}^{1}$

En tenant compte des propriétés des fonctions de base, sur un élément $e_{k}$ on doit calculer uniquement $3*3$ intégrales élémentaires faisant intervenir les 3 fonctions de base associées aux 3 sommets de l'élément. En notant $\{ n_{1},n_{2},n_{3}\}$ les numéros de ces 3 sommets, on a donc à calculer sur un élément $e_{k}$ la matrice élémentaire $3*3$ suivante:

$$\mathbf{A}_{pq}^{k}=\int_{e_{k}}K(\frac{\partial\Phi_{n_{q}}... ...{e_{k}}\alpha\Phi_{n_{q}}.\Phi_{n_{p}} d\omega (p,q=1,3)$$

En considérant que les coefficients $K$ et $\alpha$ sont constants sur l'élément $e_{k}$, cette matrice est la somme d'une matrice élémentaire de rigidité $\mathbf{K}^{k}$ et d'une matrice de masse $\mathbf{M}^{k}$:

$$\mathbf{A}_{pq}^{k}=K \mathbf{K}_{pq}^{k}+\alpha \mathbf{M}_{pq}^{k}$$

Le calcul de ces 2 matrices s'effectue à l'aide du changement de variables vers l'élément de référence $\hat{e}$.

5.2.4.1 matrice élémentaire de rigidité

Nous rappelons l'expression générique (6.25) de la matrice de rigidité calculée sur l'élément de référence:

$$\mathbf{K}_{pq}^{k}=\int_{0}^{1}\left[\int_{0}^{1-\xi}\left(... ...end{array}\right]\right)\det(\mathbf{J}_{k}) d\eta\right]d\xi$$

Le calcul effectué au paragraphe 6.1.4 montre que cette matrice dépend de 3 coefficients (6.28), que l'on peut écrire sous forme vectorielle:

$$\mathbf{K}_{22}^{k}=\frac{\left\Vert \overrightarrow{S_{1}S_... ...arrow{S_{1}S_{3}}\otimes\overrightarrow{S_{2}S_{1}}\right\Vert$$

que l'on reporte dans l'expression de $\mathbf{K}^{k}$

$$\mathbf{K}^{k}=\left[\begin{array}{ccc} K_{22}^{k}+K_{33}^{k... ...{23}^{k}-K_{33}^{k} & K_{23}^{k} & K_{33}^{k}\end{array}\right]$$

En utilisant les notations matricielles de Matlab, ces formulent se programment très simplement. La fonction MatriceRigidite (programme 6.2.4) implémente directement sous Matlab les relations précédentes. On a aussi utiliser le fait que le produit vectoriel $\overrightarrow{S_{1}S_{3}}\otimes\overrightarrow{S_{2}S_{1}}$ a une seule composante non nulle, qui est suivant l'axe $z$, et qui est positive si les sommets sont donnés dans l'ordre trigonométrique dans la table de connection (ce qui est le cas).

programme matlab 6.2.4: Matrice élémentaire de rigidité: MatriceRigidite.m

function [Ke]=MatriceRigidite(G,k)
% calcul de la matrice de rigidite de l'element k
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));
K22=S13*S13'/(4*Aire);
K33=S21*S21'/(4*Aire);
K23=(S13*S21')/(4*Aire);
Ke=[K22+K33+2*K23,-K23-K22,-K23-K33;
    -K23-K22     ,  K22   , K23    ;
    -K23-K33     ,  K23   , K33    ];
return;

5.2.4.2 matrice élémentaire de masse

De manière identique, le changement de variable permet d'obtenir l'expression générique de la matrice de masse:

$$\mathbf{M}_{pq}^{k}=\int_{0}^{1}\int_{0}^{1-\xi}N_{q}(\xi,\eta)N_{p}(\xi,\eta)\det(\mathbf{J}_{k}) d\eta d\xi$$

Sur l'élément de référence, le calcul de cette matrice est simple, puisque le déterminant du jacobien est constant, $\det(\mathbf{J}^{k})=2 aire_{k}$ et les fonctions de formes $L_{q}$ et $L_{p}$ sont des polynômes très simples en ($\xi,\eta$). Compte tenu de l'expression de ces polynômes (6.12) et des propriétés de symétrie, il suffit de calculer 2 intégrales:

$$\int_{0}^{1}\int_{0}^{1-\xi}\eta\xi d\eta d\xi=\frac{1}{24}, \int_{0}^{1}\int_{0}^{1-\xi}\xi^{2} d\eta d\xi=\frac{1}{12}$$

En effet les propriétés de symétrie imposent:

$$\int_{0}^{1}\int_{0}^{1-\xi}(N_{1})^{2} d\eta d\xi=\int_{0... ...} d\eta d\xi=\int_{0}^{1}\int_{0}^{1-\xi}\xi^{2} d\eta d\xi$$

$$\int_{0}^{1}\int_{0}^{1-\xi}N_{1}N_{2} d\eta d\xi=\int_{0}... ...} d\eta d\xi=\int_{0}^{1}\int_{0}^{1-\xi}\eta\xi d\eta d\xi$$

d'où l'expression de la matrice de masse élémentaire:

$$\mathbf{M}^{k}=\frac{aire_{k}}{12}\left[\begin{array}{ccc} 2 & 1 & 1\\ 1 & 2 & 1\\ 1 & 1 & 2\end{array}\right]$$

dont la fonction Matlab MatriceMasse est donné (programme ) ci dessous.

programme matlab 6.2.4: Matrice élémentaire de masse: MatriceMasse.m

function [Me]=MatriceMasse(G,k)
% calcul de la matrice de rigidite de l'element k
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));
Me=Aire/12*[2,1,1; 1,2,1; 1,1,2];
return;

5.2.5 Calcul du second membre élémentaire $\mathcal{P}^{1}$

Pour chaque élément $e_{k}$, on a à calculer 3 intégrales élémentaires associées à chacun des 3 sommets de l'élément. Avec les mêmes notations que précédemment, ces intégrales s'écrivent:

$$\mathbf{B}_{p}^{k}=\int_{e_{k}}f.\Phi_{n_{p}} d\omega (p=1,3)$$

Pour calculer ces intégrales, on remplace la fonction $f(x,y)$ par son interpolation $f^{h}(x,y)$ sur le maillage éléments finis:

$$f^{h}(x,y)=\sum_{i=1}^{nn}f_{i} \Phi_{i}(x,y)$$

ce qui permet d'écrire les intégrales sous la forme:

$$\mathbf{B}_{p}^{k}=\sum_{q=1}^{3}f_{n_{q}}\left(\int_{e_{k}}\Phi_{n_{q}}.\Phi_{n_{p}} d\omega\right) (p=1,3)$$

puisque sur l'élément $e_{k}$, l'interpolation $f^{h}$ de $f$ s'écrit: $f^{h}(x,y)=\sum_{q=1}^{3}f_{n_{q}}.N_{n_{q}}(x,y)$.

C'est donc le produit du vecteur $\mathbf{F}^{k}=[f_{n_{1}},f_{n_{2}},f_{n_{3}}]$ des valeurs nodales par la matrice élémentaire de masse $\mathbf{M}^{k}$. Le second membre élémentaire s'écrit donc:

$$\mathbf{B}^{k}=\mathbf{M}^{k}\mathbf{F}^{k}=\frac{aire_{k}}{... ...array}{c} f_{n_{1}}\\ f_{n_{2}}\\ f_{n_{3}}\end{array}\right]$$

et la fonction Matlab SmbElement est donné (programme 6.2.5) ci dessous.

programme matlab 6.2.5: Second membre élémentaire: SmbElement.m

function [Be]=SmbElement(G,F,k)
% calcul du second membre elmentaire (F valeurs nodales de f)
n=G.Tbc(k,:); % numero des sommets de l'element k
Fk=F(n);
Mk=MatriceMasse(G,k);
Be=Mk*Fk;
return;

5.2.6 Assemblage

L'assemblage de la matrice $\mathbf{A}$ et du second membre $\mathbf{B}$ consiste à passer en revu les éléments, à calculer les matrices élémentaires de masse et de rigidité ainsi que le second membre élémentaire, puis à mettre ses valeurs aux bons endroits dans la matrice et le second membre globale.

L'algorithme d'assemblage (2) ci dessous donne le principe de l'assemblage.

Le programme Matlab ci dessous implémente cet algorithme, en utilisant les notations matricielles de Matlab qui permettent d'éviter l'écriture de boucles. Ainsi les boucles internes 5 à 12 sont remplacées par 2 lignes de codes matriciels Matlab (lignes 11 et 12 dans le programme 6.2.6)

La fonction Matlab Assemblage est donné (programme 6.2.6) ci dessous.

programme matlab 6.2.6: Assemblage de A et B: 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;

5.2.7 Conditions aux limites

L'application des conditions aux limites sur le système linéaire obtenu après l'assemblage dépend du type de conditions aux limites.

5.2.7.1 Conditions de Dirichlet

L'imposition des conditions aux limites de Dirichlet (frontières $\Gamma_{1}$ et $\Gamma_{2}$) consiste à fixer la valeur de la solution aux noeuds $M_{i}$ se trouvant sur la frontière de Dirichlet ( $\Gamma_{1}\cup\Gamma_{2}$). Pour cela on remplace simplement dans le système linéaire l'équation $i$ par l'équation:

$$\begin{eqnarray*} u_{i}=0 & \mbox{si } M_{i}\in\Gamma_{1}\\ u_{i}=u_{e} & \mbox{si } M_{i}\in\Gamma_{2}\end{eqnarray*}$$

Dans la matrice $\mathbf{A}$, cela revient à annuler la ligne $i$ et à mettre un 1 sur la diagonale, et dans le second membre $\mathbf{B}$, à remplacer la composante i par 0 ou $u_{e}$ selon le cas:

$$\mathbf{A}_{ij}=0, \mathbf{A}_{ii}=1, \mathbf{B}_{i}=0 \mbox{ ou } u_{e}$$

Les termes supplémentaires dues à cette conditions aux limites sont automatiquement pris en compte dans les autres équations.

5.2.7.2 Conditions de Neuman

La condition aux limites de Neuman homogène (frontière $\Gamma_{3}$) est déjà prise en compte dans la formulation et ne nécessite pas de modification supplémentaire. Par contre la condition mixte sur la frontière $\Gamma_{4}$ nécessite le calcul d'intégrales de bords:

$$\mathbf{A}_{ij}^{cl}=\int_{\Gamma_{4}}\beta\Phi_{j}.\Phi_{i}... ...,\mathbf{B}_{i}^{cl}=-\int_{\Gamma_{4}}\Phi_{0}N_{i} d\gamma$$ (5.36)

Compte tenu de la propriété des fonctions de base, ces contributions interviennent que pour des fonctions de base associées à des noeuds $M_{i}$ sur la frontière $\Gamma_{4}$. Le calcul de ces intégrales se décompose en calcul élémentaire sur les arêtes de la frontière $\Gamma_{4}$. Pour cela nous allons tout d'abord déterminer les arêtes frontières de la géométrie.

Soit AF le tableau des arêtes frontières, i.e. l'ensemble des arêtes des éléments $e_{k}$ du maillage se trouvant sur la frontière $\Gamma=\partial\Omega$ de la géométrie:

$$AF=\bigcup_{k=1}^{ne}(e_{k}\cap\Gamma)$$

Chaque arête frontière $AF_{l}$ est définie par le numéro du premier et du second sommet de l'arête (parcourue dans le sens trigonométrique: i.e. avec une normale extérieure à droite). Le tableau AF est donc un tableau de $2*naf$ entiers (en notant $naf$ le nombre total d'arêtes frontières). Par exemple pour la maillage de la figure (6.13), le nombre d'arêtes frontières vaut $naf=18$, et le tableau AF vaut:

arêtes AF

numéro 1er sommet

numéro 2nd sommet 1

Remarque: l'ordre des arêtes est totalement arbitraire.

Pour déterminer ces arêtes, nous utiliserons l'algorithme (3) suivant:

Il faut noté qu'une arête frontière n'est pas forcément une arête ayant ses 2 sommets sur la frontière. Ainsi l'élément 33 du maillage de la figure (6.13) a une arête $[27,21]$ dont les deux sommets sont sur la frontière, mais qui n'est pas frontière. Donc après avoir déterminé la liste des arêtes ayant les 2 sommets sur la frontière (lignes 1 à 3), il faut éliminer les arêtes internes, i.e. apparaissant 2 fois dans la liste. On peut en effet vérifier la propriété suivante:

$ij \mbox{ arete interne} \Leftrightarrow\exists k,l \mbox{ t.q. } ij\in e_{k} \mbox{ et } ji\in e_{l}$

Le programme Matlab 6.2.7 correspondant est donné ci dessous. La fonction Matlab AretesFrt utilise la fonction find de Matlab qui permet de déterminer la liste des indices d'un tableau, dont les valeurs vérifient une condition. Pour tester si un noeud est sur une frontière, on utilise le tableau Frt .

programme matlab 6.2.7: détermination des arêtes fronières: AretesFrt.m

function [AF]=AretesFrt(G)
% determine la liste des aretes frontieres d'un maillage G
% liste des elts frontieres (i.e au moins 1 nds frontiere)
if (G.ddl==3)
 Count=G.Frt(G.Tbc(:,1))+G.Frt(G.Tbc(:,2))+G.Frt(G.Tbc(:,3));
 Efrt=find(Count>0);
 % liste des aretes des elements frontieres
 AF0=[G.Tbc(Efrt,1:2); G.Tbc(Efrt,2:3); G.Tbc(Efrt,3), G.Tbc(Efrt,1)] ;
elseif (G.ddl==4)
 Count=G.Frt(G.Tbc(:,1))+G.Frt(G.Tbc(:,2))+G.Frt(G.Tbc(:,3))+G.Frt(G.Tbc(:,4));
 Efrt=find(Count>0);
 % liste des aretes des elements frontieres
 AF0=[G.Tbc(Efrt,1:2); G.Tbc(Efrt,2:3); G.Tbc(Efrt,3:4); G.Tbc(Efrt,4), G.Tbc(Efrt,1)] ;

else
 disp('Erreur type de geometrie'); return;   
end;
% recherche des aretes frontieres parmi cette liste i.e. les 2 nds sur la frontiere
LAF0=find(G.Frt(AF0(:,1))>0 & G.Frt(AF0(:,2))>0) ;
% liste des aretes possiblement frontieres
AF1=AF0(LAF0,:);
% recherche des aretes doubles (non frontiere)
na=size(AF1,1); 
for k=1:na-1
    ij=AF1(k,:);
    % test si arete ji est dans la liste
    I=(find(AF1(k+1:na,1)==ij(2) &  AF1(k+1:na,2)==ij(1)));
    if (~isempty(I))    AF1(k,:)=[0 0]; AF1(k+I,:)=[ 0 0];    end;
end;
% elimination
LAF1=find(AF1(:,1)~=0);
AF=AF1(LAF1,:);
% fin
return;

Ayant la liste AF des arêtes frontières, une intégrale sur $\Gamma_{4}$ se décompose en intégrale élémentaire:

$$\int_{\Gamma_{4}}f d\gamma=\sum_{l=1, AF_{l}\in\Gamma_{4}}^{naf}\int_{AF_{l}}f d\gamma$$

Pour calculer les intégrales de bord (6.36), il suffit donc de calculer des intégrales élémentaires sur les arêtes frontières de $\Gamma_{4}$, soit pour une arête $AF_{l}$ de sommets $n_{1}=AF(l,1)$ et $n_{2}=AF(l,2)$ les 6 intégrales suivantes:

$$\mathbf{Af}_{pq}^{l}=\int_{M_{n_{1}}}^{M_{n_{2}}}\beta\Phi_{... ...{n_{2}}}\phi_{0}\Phi_{n_{p}} d\gamma \mbox{pour }p,q=1,2$$

puisque sur cette arête, seules 2 fonctions de base sont non nulle: les 2 fonctions de base $\Phi_{n_{1}}$ et $\Phi_{n_{2}}$ associées aux 2 sommets. Pour calculer ces intégrales on effectue un changement de variable de l'arête $AF_{l}$ vers le segment de référence $[-1,1]$ (figure 6.15).

Figure 6.15: calcul des intégrales frontières

Sur cet élément, l'expression des 2 fonctions de base est très simple: ce sont les 2 polynômes de Lagrange $\mathcal{P}^{1}$en $\xi$

$$\Phi_{n_{1}}(x,y)\mid_{EF_{l}}=N_{1}(\xi)=\frac{1-\xi}{2} \... ...et }\Phi_{n_{2}}(x,y)\mid_{EF_{l}}=N_{2}(\xi)=\frac{1+\xi}{2}$$

puisque $\Phi_{n_{1}}$ est une fonction affine qui vaut 1 au noeud $n_{1}$ et 0 au noeud $n_{2}$, et idem pour $\Phi_{n_{2}}$. La variable $\xi$ a une interprétation géométrique, puisque c'est l'abscisse curviligne sur le segment $[n_{1},n_{2}]$.

Avec ces notations, les intégrales élémentaires de bord s'écrivent:

$$\mathbf{Af}_{pq}^{l}=\int_{-1}^{+1}\beta N_{q}(\xi).N_{p}(\x... ...phi_{0}N_{p}(\xi) \frac{h^{l}}{2}d\xi \mbox{pour }p,q=1,2$$

en notant $h^{l}$ la longueur de l'arête $AF_{l}$ . Un calcul élémentaire fournit la valeur de ces intégrales dans le cas où $\beta$ et $\Phi_{0}$ sont constants par arêtes:

$$\mathbf{Af}^{l}=\beta h^{l}\left[\begin{array}{cc} \frac{1}{... ...t[\begin{array}{c} \frac{1}{2}\\ \frac{1}{2}\end{array}\right]$$

ayant ses 2 sommets sur la frontière. Ayant ces intégrales élémentaires, il suffit ensuite d'appliquer une procédure d'assemblage pour insérer ces contributions dans la matrice $\mathbf{A}$ et le second membre $\mathbf{B}$.

Le programme Matlab ci dessous implémente cet assemblage. La fonction Matlab Climite applique les conditions aux limites en modifiant la matrice et le second membre générique (i.e. calculés sans tenir compte des conditions aux limites). On utilise aussi la convention suivante: les arêtes de la frontière $\Gamma_{4}$ sont les arêtes ayant au moins un noeud dont le numéro de frontière est égale à 4. Ainsi pour le maillage de la figure (6.13), la frontière $\Gamma_{4}$ va du noeud 30 au noeud 26 (et non de 29 à 27). On calcule ainsi toutes les intégrales de bords, mais on modifie aussi l'équation pour les 2 noeuds extrémités 30 et 26. Si ces noeuds sont sur une frontière de Dirichlet (ce qui est le cas), il faudra imposer la condition de Dirichlet après le calcul de ces intégrales de bords. C'est ce qui est fait dans la fonction Matlab Climite , où on impose d'abord les conditions faibles, puis ensuite les conditions fortes.

programme matlab 6.2.7: Application des conditions aux limites: Climite.m

function [A1,B1]=Climite(G,A,B,beta,phi0,Ue)
% application des C.L. sur la matrice A et B
A1=A; B1=B;
% liste des aretes frontieres
AF=AretesFrt(G);
% application des cdts mixtes sur gamam4
AF4=AF(find((G.Frt(AF(:,1))==4)|(G.Frt(AF(:,2))==4)),:);
na=size(AF4,1);
for k=1:na
    n1=AF4(k,1); n2=AF4(k,2); % numero des 2 sommets
    dx=norm(G.X(n2,:)-G.X(n1,:)); % longueur de l'arete
    A1(n1,n1)=A1(n1,n1)+beta*dx/3; A1(n2,n2)=A1(n2,n2)+beta*dx/3;
    A1(n1,n2)=A1(n1,n2)+beta*dx/6; A1(n2,n1)=A1(n1,n2)+beta*dx/6;
    B1(n1)=B1(n1)-phi0*dx/2;
    B1(n2)=B1(n2)-phi0*dx/2;
end;
% applications des cdts fortes
ND1=find(G.Frt(:)==1);
for i=ND1'
 A1(i,:)=0; A1(i,i)=1.0; B1(i)=0;
end;
ND2=find(G.Frt(:)==2);
% pour les cdts de Dirichlet on peut passer une fonction 
if isa(Ue,'double') 
  for i=ND2'
    A1(i,:)=0; A1(i,i)=1.0; B1(i)=Ue;
  end;
else
% Ue est une fonction inline de (x,y)
  for i=ND2'
    A1(i,:)=0; A1(i,i)=1.0; 
    B1(i)=Ue(G.X(i,1),G.X(i,2));
  end;
end;
return;

5.2.8 Résolution

Après imposition des conditions aux limites, la solution approchée $u^{h}$ s'obtiens par résolution du système linéaire. Le script Matlab 6.2.8 qui enchaıne la suite des opérations pour calculer cette solution approchée est donné ci-dessous.

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

% resolution du probleme model
G=Lecture('maillage.msh');
% parametres
alpha=0.; K=1.0; beta=0.0; phi0=-1.0; Ue=2;
% second membre
F=zeros(G.nn,1);
% assemblage
[A,B]=Assemblage(G,K,alpha,F);
% conditions aux limites
[A,B]=Climite(G,A,B,beta,phi0,Ue);
% resolution
% renumerotation pour optimisation
m = symamd(A);
[LA,UA]=lu(A(m,m));
U=UA\(LA\B);

Figure 6.16: solution $u^{h}$ pour $\phi _{0}=0$ (à gauche) et $\phi _{0}=-1$ (à droite)

Pour les valeurs des paramètres du programme ci dessus, on a tracé la solution obtenue pour 2 valeurs de $\phi_{0}$: $\phi _{0}=0$ et $\phi _{0}=-1$ (figure 6.16). On constate, sur ce problème simple de diffusion sans terme source, l'influence de la condition aux limites sur $\Gamma_{4}$ , qui pour le second cas dévie les iso-valeurs de la solution vers la droite par rapport au cas homogène. Ceci correspond bien à l'imposition d'une valeur de dérivée normale positive. On note aussi que les conditions de Neuman ne sont pas vérifiées exactement par la solution approchée, car les iso-valeurs ne sont pas exactement perpendiculaires aux frontières de Neuman. C'est seulement à la limite, que ces conditions seront vérifiées exactement.