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

Nous allons maintenant étudier l'approximation par éléments finis $\mathcal{P}^{1}$, en particulier ses propriétés en terme de précision et de convergence. Nous donnerons en même temps le programme Maple correspondant.

Pour cela nous allons traiter le cas d'une équation générale du type (4.1) sur le domaine $\Omega=[0,L]$:

$$\left\{ \begin{array}{c} -\frac{d}{dx}(K(x)\frac{d}{dx}u(x))... ..., -K(L)\frac{du}{dx}(L)=\beta u(L)+\phi_{0}\end{array}\right.$$ (3.19)

Cette équation traduit un phénomène de diffusion avec un coefficient de diffusion $K(x)$ variable, couplé à un terme source fonction de la solution $u$ et d'un coefficient $\alpha(x)$ variable. Le second membre $f(x)$ traduit la partie du terme source indépendant de $u$. En $x=0$, on impose une condition de Dirichlet, et en $x=L$ une condition de Fourier (ou condition de Robin, ou condition mixte) qui impose que le flux de chaleur $\phi_{e}$en $x=L$ soit une fonction de la solution $u(L)$ (et non pas uniquement constant comme dans l'exemple précédent): $\phi_{e}=\beta u(L)+\phi_{0}$

3.2.1 Formulation faible

Pour obtenir la formulation faible de (4.19), on procède comme précédemment. On multiplie l'équation (4.19) par une fonction test $v(x)$, on intègre sur le domaine $\mathcal{D}$, puis on effectue une intégration par partie sur les termes de plus haut degré. Il vient:

$$\int_{0}^{L}K(x)\frac{du}{dx} \frac{dv}{dx} dx-[K(x)\frac... ...\int_{0}^{L}\alpha(x)u(x) v(x) dx=\int_{0}^{L}f(x)v(x) dx$$

Le calcul du terme de bord se fait en utilisant les conditions aux limites et en imposant à la fonction test $v(x)$ d'être une variation de la solution $u(x)$. Comme on impose la valeur de $u$ en $x=0$ ($u(0)=u_{e}$), sa variation doit être nulle en ce point. La fonction test $v$ doit donc s'annuler en $x=0$. et le terme de bord $K(0)\frac{du}{dx}(0)v(0)$ s'annulle. Le second terme de bord se calcule à l'aide de la seconde condition aux limites: $-K(L)\frac{du}{dx}(L)=\beta u(L)+\phi_{0}$. La formulation faible de (4.19) s'écrit alors:

$$\left\{ \begin{array}{c} \mbox{Trouvez }u(x) \mbox{ ave... ...n test }v(x)\mbox{ vérifiant }v(0)=0\end{array}\right.$$ (3.20)

Exercice: montrez que cette formulation faible est équivalente à la formulation variationnelle suivante:

$$\left\{ \begin{array}{c} \mbox{Trouvez }u(x) \mbox{ ave... ...L)^{2}-\int_{0}^{L}f(x)u(x) dx+\phi_{0}u(L)\end{array}\right.$$

3.2.2 Approximation par éléments finis

Pour construire l'approximation par éléments finis $\mathcal{P}^{1}$, on crée un maillage $\mathcal{M}^{h}$ du domaine de calcul constitué de $ne$ segments de coordonnées $[x_{k},x_{k+1}]$:

$$\mathcal{M}^{h}=\bigcup_{k=1}^{ne}[x_{k},x_{k+1}] \mbox{ avec } x_{1}=0, x_{ne+1}=L$$

On note $e_{k}=[x_{k},x_{k+1}]$, l'élément $k$ du maillage de longueur $h^{k}=x_{k+1}-x_{k}$.

Remarque: on a choisit de numéroter les noeuds du maillage de 1 à $ne+1,$ contrairement à l'exemple précédent, où on les avait numéroté de 0 à $N_{e}$. Ce choix sera justifier par la suite,lorsque l'on programmera la méthode. En effet dans de nombreux langage de programmation ( entre autre Maple et Matlab), les indices de tableaux commencent à 1 et non pas à 0. D'autre part, contrairement au “calcul à la main” précédent, nous ne chercherons pas à éliminer la valeur au premier noeud de façon à conserver une approche générale.

Figure: transformation: $\mathcal{T}_{4}:x=\frac{1-\xi}{2}x_{4}+\frac{1+\xi}{2}x_{5}$ et $\mathcal{T}_{7}:x=\frac{1-\xi}{2}x_{7}+\frac{1+\xi}{2}x_{8}$

Sur ce maillage, l'approximation $u^{h}(x)$ par éléments finis $\mathcal{P}^{1}$ de la solution est définie par ses valeurs $u^{h}(x_{i})=u_{i}$ aux $nn=ne+1$ points nodaux $\{ x_{i}\}_{i=1,nn}$, correspondant aux extrémités des segments du maillage. Cette approximation possède donc $nn=ne+1$ degrés de liberté $\{ u_{i}\}$ correspondant aux $nn$ noeuds d'interpolation $\{ x_{i}\}$.

Sur chaque élément $e_{k}$ cette approximation est définie à partir des 2 fonctions de forme $\{ N_{1},N_{2}\}$ de l'élément $\mathcal{P}^{1}$ (4.6) et du changement de variable $\xi=\xi(x)$ (4.5):

$$\begin{eqnarray*} u^{h}(x)=u_{k}N_{1}(\xi(x))+u_{k+1}N_{2}(\xi(x)) \mbox{ pou... ... e_{k}=[x_{k},x_{k+1}]\\ \xi(x)=2\frac{x-x_{k}}{x_{k+1}-x_{k}}-1\end{eqnarray*}$$

L'approximation globale est la somme de ces approximations élémentaires et s'écrit en fonction des fonctions de base $\Phi_{i}(x)$ :

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

La fonction de base $\Phi_{i}(x)$ est associée au noeud $S_{i}$ du maillage et est définie à partir des fonctions de forme $N_{q}$ par les relations:

$$\begin{eqnarray*} \Phi_{i}(x)=N_{1}(\xi) & \mbox{si} & x\in e_{i}\\ \Phi_{i}(x)... ...}(\xi) & \mbox{si} & x\in e_{i-1}\\ \Phi_{i}(x)=0 & \mbox{sinon}\end{eqnarray*}$$

Figure 4.13: Fonctions de base $\Phi _{i}$ et interpolation sur le maillage (4.12)

Le programme 4.2.2 implémente sous Maple l'interpolation éléments finis $\mathcal{P}^{1}$ sur un maillage de $N_{e}$ éléments.

programme Maple 4.2.2: Approximation éléments finis P1

# Approximation par éléments finis
> with(linalg): with(plots):
# Nombre d'éléments Ne 
> L:=3; Ne:=8; 
# Maillage et coordonnees des pts de maillage
> h:=L/Ne; Xp:=[seq(i*h,i=0..Ne)];
# d=degre d'interpolation, Nn=Nombre de noeuds, X=coordonnees des nds 
> d:=1; Nn:=Ne+1;X:=evalm(Xp):
# num=numérotation par élément et long=longueur de chaque element
> num:=k->[k,k+1]; long:=k->Xp[k+1]-Xp[k];
# Fonctions de forme de l'element P1
> N1:=unapply(interp([-1,1],[1,0],xi),xi);
> N2:=unapply(interp([-1,1],[0,1],xi),xi);
# Transformation x=x(xi) vers l'element de reference [-1,1]
> Xi:=unapply(N1(xi)*'Xp'[i]+N2(xi)*'Xp'[i+1],xi,i);
# Interpolation G d'une fonction g(x) sur le maillage
> Interpol:=proc(g)
> global Nn:
> local G,i;
> G:=vector(Nn,0):
> for i from 1 to Nn do G[i]:=g(X[i]); end;
> return evalm(G);
> end proc:
# Calcul de la norme de l'erreur d'interpolation: f(x)=fonction
# et F= interpolation (calcul sur l'elt de reference)
> Erreur:=proc(f,F)
> global Ne,Xp;
> local err,k,h;
> err:=0;
> for k from 1 to Ne do
>   h:=(Xp[k+1]-Xp[k])/2.0;
>   err:=err+h*
>      evalf(Int((f(Xi(xi,k))-(F[k]*N1(xi)+F[k+1]*N2(xi)))^2,
>          xi=-1..1));
> end:
> return sqrt(err);
> end:
> 
# Example: fonction cos(2x)
> f:=x->cos(2.0*x);
> F:=Interpol(f);
> P1:=plot(f(x),x=0..L):
> P2:=plot([seq([X[i],F[i]],i=1..Nn)],style=point,color=blue):
> display(P1,P2);
> err:=Erreur(f,F);

Le programme est écrit de façon modulaire et général, de telle sorte que l'on puisse très facilement implémenter une interpolation de degré $d$ quelconque. Pour cela on introduit comme variable le degré $d$ du polynôme d'interpolation (ligne 8), et on distingue les noeuds du maillage, que l'on note $Xp$ (ligne 6), des points d'interpolation, que l'on note $X$ (ligne 8). On introduit ensuite 2 fonctions (ligne 10): la fonction $num(k)$ qui pour un élément $k$ renvois les numéros des points d'interpolation de l'élément, et la fonction $long(k)$ qui calcul simplement la longueur $h_{k}$ de l'élément $k$. On calcule ensuite les fonctions de forme de l'élément $\{ N_{p}(\xi)\}$ comme polynômes de Lagrange en utilisant la fonction Maple interp (lignes 12,13). La transformation de l'élément $k$ vers l'élément de référence $x=x(\xi)$ est noté $Xi(\xi,k)$ (ligne 15):

$$x=Xp[k]*N_{1}(\xi)+Xp[k+1]*N_{2}(\xi)$$

On écrit ensuite une fonction $Interpol(f)$ qui calcul le vecteur $F$ des valeurs d'une fonction f(x) aux points d'interpolation: $F[i]=f(X[i])_{i=1,nn}$ (lignes 17 à 23).

Enfin on écrit une fonction $Erreur(f,F)$ qui calcule la norme de l'erreur entre une fonction $f(x)$ et son approximation $f^{h}$ sur le maillage élément finis. On choisit comme norme, l'intégrale du carré de la différence (norme $L^{2}$):

$$\Vert f-f^{h}\Vert^{2}=\int_{0}^{L}(f(x)-f^{h}(x))^{2} dx$$

Cette norme mesure l'erreur moyenne sur le domaine de calcul. Pour cela on calcule l'intégrale élément par élément, en éffectuant pour chaque élément la transformation vers l'élément de référence. Pour une approximation $\mathcal{P}^{1}$ sur un élément $e_{k}=[x_{k},x_{k+1}]$, on a:

$$f^{h}(\xi)=F[k]*\Phi_{1}(\xi)+F[k+1]*\Phi_{2}(\xi)$$

et l'erreur sur un élément s'écrit:

$$\Vert f-f^{h}\Vert_{e_{k}}^{2}=\frac{h_{k}}{2}\int_{-1}^{+1}(f(x(\xi)-(F[k]*\Phi_{1}(\xi)+F[k+1]*\Phi_{2}(\xi))^{2}dx$$

Le calcul de cette intégrale est laissé à Maple. Pour cela on a utilisé la fonction d'intégration $int$ sous sa forme directe d'intégration numérique i.e. $evalf(Int(..))$ (ligne 33). Cela évite que Maple cherche tout d'abord à effectuer une intégration analytique, puis calcule ensuite la valeur numérique. On gagne ainsi un factor 5 à 10 en temps de calcul, ce qui n'est pas négligeable pour des maillages importants.

La fin du programme est un exemple d'interpolation d'une fonction $f(x)=sin(2x)$ sur un maillage de 10 éléments de $[0,3].$ Les fonctions de base et l'interpolation sont données sur la figure (4.13). L'erreur d'interpolation vaut dans ce cas: $\Vert f-f^{h}\Vert=0.0606$

Exercice: calculer l'erreur d'interpolation pour des maillages de plus en plus fins.

3.2.3 Formulation faible discrète

La solution approchée $u^{h}$ par éléments finis $\mathcal{P}^{1}$ de (4.20) s'écrit sous la forme:

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

Elle doit vérifier les conditions aux limites fortes ( i.e. la condition de Dirichlet en $x=0$ : $u^{h}(0)=u_{e}$). En notant à nouveau que les fonctions de bases vérifient $\Phi_{i}(x_{j})=\delta_{ij}$, cette condition impose la valeur nodale de $u^{h}$ au noeud $x_{1}=0$ du maillage:

$$u^{h}(x)=\sum_{j=1}^{nn}u_{j}\Phi_{j}(x) \mbox{ avec } u_{1}=u_{e}$$ (3.21)

Les fonctions tests associées $v^{h}$ étant des variations de $u^{h}$, elles doivent donc s'annuller en $x=0$. Elles s'écrivent sous la forme générale suivante:

$$v^{h}(x)=\sum_{i=2}^{nn}v_{i}\Phi_{i}(x) \mbox{ car }v_{1}=0$$ (3.22)

Ces fonctions tests sont des combinaisons linéaires des $nn-1$ fonctions de base $\{\Phi_{i}(x)\}_{i=2,nn}$. En remplaçant dans la formulation faible (4.20) la solution exacte $u$ par la solution approchée $u^{h}$ donnée par (4.21) et la fonction test $v$ par une de ces $nn-1$ fonctions de base $\{\Phi_{i}(x)\}_{i=2,nn}$ , on obtient la formulation faible discrète:

$$\sum_{j=1}^{nn}u_{j}\left(\int_{0}^{L}K(x)\frac{d\Phi_{j}}{dx}\fr... ...i}(x) dx\right)+\underbrace{\beta u_{nn}\Phi_{i}(L)}_{terme C.L sur A.} =$$

$$\int_{0}^{L}f(x)\Phi_{i}(x) dx-\underbrace{\phi_{0}\Phi_{i}(L)}_{terme C.L. sur B}$$ (3.23)

Pour obtenir cette relation on a permuté la sommation $\sum_{i=1}^{nn}$ et l'intégration $\int_{0}^{L}$ et on a sortie les coefficients $u_{i}$ des intégrales. D'autre part on a remplacé $u^{h}(L)$ par sa valeur $u_{nn}$. L'équation (4.23) écrite pour $i=2,nn$ est un système linéaire de $nn-1$ inconnues $\{ u_{i}\}_{i=2,nn}$ (en notant que $u_{1}=u_{e}$ est fixé par la condition aux limites), qu'il suffit de résoudre pour obtenir la solution approchée $u^{h}$. De façon à construire un programme le plus général possible, on considérera que l'on a $nn$ inconnues $u_{i}$, qui sont données par les $nn-1$ équations (4.23) auxquelles on ajoute l'équation supplémentaire $u_{1}=u_{e}$. Dans cette approche on construit un système linéaire de $nn$ inconnues à $nn$ équations. Dans un premier temps cela permettra de construire la matrice $\mathbf{A}$ et le second membre $\mathbf{B}$ du système linéaire sans tenir compte des conditions aux limites et donc de façon générique. Puis dans un second temps, on appliquera les conditions aux limites sur le système linéaire: i.e. on remplacera la première équation par l'équation $u_{1}=u_{e}$ et on introduira le terme lié à la conditions aux limites en $x=L$

Les coefficients génériques de $\mathbf{A}$ et de $\mathbf{B}$ s'écrivent:

$$\begin{eqnarray*} \mathbf{A}_{ij} & = & \int_{0}^{L}K(x)\frac{d\Phi_{j}}{dx}\fra... ...(x) dx\\ \mathbf{B}_{i} & = & \int_{0}^{L}f(x)\Phi_{i}(x) dx\end{eqnarray*}$$

Pour calculer ces coefficients, on effectue un calcul élément par élément en déterminant les matrices élémentaires et les second membres élémentaires sur un élément $e_{l}$ . Les calculs des intégrales font intervenir des coefficients variables $K(x)$ et $\alpha(x)$ et nous utiliserons Maple pour calculer les intégrales faisant intervenir ces coefficients. Une autre approche consisterait à calculer une approximation éléments finis de ces coefficients sur chaque élément, ou à choisir une valeur moyenne par élément.

Exercice: comparer le calcul de la matrice élémentaire de raideur avec une approximation de $K(x)$ constante, $\mathcal{P}^{1}$ et exacte dans le cas où $K(x)$ est un polynôme de degré 1 et 2.

Pour le second membre $f(x)$ , nous utiliserons une approximation $\mathcal{P}^{1}$ sur le maillage éléments finis:

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

3.2.4 Matrice élémentaire

Sur un élément $e_{k}$ la matrice élémentaire est la somme d'une matrice de rigidité $\mathbf{K}^{k}$ et d'une matrice de masse $\mathbf{M}^{k}$, qui pour des éléments finis $\mathcal{P}^{1}$ sont des matrices 2*2 puisqu'il y a 2 fonctions de forme $\{ N_{1},N_{2}\}$ associées aux 2 points d'interpolations de numéros $\{ n_{1},n_{2}\}$ (figure4.6) :

$$N_{1}(\xi)=\frac{1-\xi}{2}, N_{2}(\xi)=\frac{1+\xi}{2} \„n_{1}=k, n_{2}=k+1$$

Figure: élément $\mathcal{P}^{1}$

De façon générale, pour des éléments finis $\mathcal{P}^{d}$ , les matrices élémentaires sont des matrices $(d+1)*(d+1)$ puisqu'il faut $d+1$ de points d'interpolation (de numéros $\{ n_{p}\}_{p=1,d+1}$ pour définir un polynôme de degré $d$ , et que l'on a donc $d+1$ fonctions de forme $\{ N_{p}(\xi)\}_{p=1,d+1}$.

Avec ces notations, de façon générique (i.e. valable pour un approximation $\mathcal{P}^{1}$, $\mathcal{P}^{2}$,$\ldots$ $\mathcal{P}^{d}$), les matrices élémentaires s'écrivent:

$$\mathbf{K}_{pq}^{k}=\frac{2}{h^{k}}\int_{-1}^{+1}K(\xi)\frac{dN_{p}}{d\xi}\frac{dN_{q}}{d\xi} d\xi ( p,q=1,d+1)$$ (3.24)

$$\mathbf{M}_{pq}^{k}=\frac{h^{k}}{2}\int_{-1}^{+1}\alpha(\xi)N_{p}(\xi)N_{q}(\xi) d\xi ( p,q=1,d+1)$$ (3.25)

Le programme 4.2.4 ci-dessous implémente le calcul de la matrice de rigidité, en programmant la relation (4.24) comme une procédure Maple. On laisse Maple effectuer les intégrations des fonctions de formes. On note enfin que pour renvoyer la valeur de la matrice élémentaire (et non son nom), on utilise la fonction evalm.

programme Maple 4.2.4: Calcul de la matrice de rigidité élémentaire

# Matrice élémentaire de rigidite
> MatRigidite:=proc(k,K)
> global d; # dimension de l'interpolation
> local Ke,p,q,h;
> h:=long(k);
> Ke:=matrix(d+1,d+1,0);
> for p from 1 to d+1 do
>   for q from 1 to d+1 do
>     Ke[p,q]:=2.0/h*
>  int(K(Xi(xi,k))*diff(N||p(xi),xi)*diff(N||q(xi),xi),xi=-1..1);
>   end;
> end;
> return evalm(Ke);
> end;
> 

Si la fonction $K(x)$ est contante, on retouve pour des éléments $\mathcal{P}^{1}$ la matrice élémentaire suivante:

$$\mathbf{K}^{k}=\frac{K}{h^{k}}\left[\begin{array}{cc} 1 & -1\\ -1 & 1\end{array}\right]$$

De même le programme 4.2.4 ci-dessous implémente le calcul de la matrice de masse, en programmant la relation (4.25) .

programme Maple 4.2.4: Calcul de la matrice de masse élémentaire

#  Matrice élémentaire de masse
> MatMasse:=proc(k,alpha)
> global d; # dimension de l'interpolation
> local Me,p,q,h
> h:=long(k);
> Me:=matrix(d+1,d+1,0);
> for p from 1 to d+1 do
>     for q from 1 to d+1 do
>      Me[p,q]:=h/2.0*
>      int(alpha(Xi(xi,k))*N||p(xi)*N||q(xi),xi=-1..1);
>     end;
> end;
> return evalm(Me);
> end;

Si la fonction $\alpha(x)$ est constante, on retouve pour des éléments $\mathcal{P}^{1}$ la matrice élémentaire suivante:

$$\mathbf{M}^{k}=\alpha h^{k}\left[\begin{array}{cc} 1/3 & 1/6\\ 1/6 & 1/3\end{array}\right]$$

3.2.5 Second membre élémentaire

Le calcul du second membre élémentaire procède de la même démarche. Pour un élément $\mathcal{P}^{1}$, de sommets $n_{1},n_{2}$, l'approximation de $f(x)$ sur un élément $e_{k}$ s'écrit en variable $\xi$:

$$f^{h}(\xi)=f_{n_{1}}N_{1}(\xi)+f_{n_{2}}N_{2}(\xi)$$

Pour un élément de degré $d$, l'approximation de $f(x)$ s'écrit:

$$f^{h}(\xi)=\sum_{q=1}^{d+1}f_{n_{q}}N_{q}(\xi)$$

d'où l'on déduit le second membre élémentaire générique :

$$\mathbf{B}_{p}^{l}=\sum_{q=1}^{d+1}f_{n_{q}}\frac{h_{k}}{2}\... ...}^{+1}N_{q}(\xi)N_{p}(\xi) d\xi \mbox{ pour } p=1,d+1$$ (3.26)

Le programme 4.2.5 ci-dessous implémente le calcul de ce vecteur élémentaire, en programmant la relation (4.26)

programme Maple 4.2.5: Vecteur second membre élémentaire

# Second membre elementaire
> SmbElem:=proc(k,F)
> global d; # dimension de l'interpolation
> local Be,q,p,h,noi,ni;
> h:=long(k); noi:=num(k);
> Be:=vector(d+1,0);
> for q from 1 to d+1 do
>   ni:=noi[q+1];
>   for p from 1 to d+1 do
>     Be[p]:=Be[p]+
>      F[ni]*h/2*int(N||q(xi)*N||p(xi),xi=-1..1);
>   end;
> end;
> return evalm(Be);
> end;

Dans le cas d'une approximation $\mathcal{P}^{1}$, on trouve l'expression:

$$\mathbf{B}^{k}=\frac{h_{k}}{6}\left[\begin{array}{c} 2f_{k}+f_{k+1}\\ f_{k}+2f_{k+1}\end{array}\right]$$

Exercice: démontrer cette dernière relation.

3.2.6 Assemblage

Le calcul de la matrice globale $\mathbf{A}$ et du second membre $\mathbf{B}$ s'effectue par une procédure générale d'assemblage, qui calcule les matrices élémentaires élément par élément et ensuite insère les coefficients de ces matrices à la bonne place dans la matrice globale.

Pour un élément $e_{k}$ de type $\mathcal{P}^{1}$ ayant comme numéro de sommets $\{ n_{1},n_{2}\}$, la matrice élémentaire $\mathbf{A}^{k}$ de dimension $(2,2)$ contribue aux coefficients de la matrice globale A suivants:

$$\begin{eqnarray*} A_{n_{1},n_{1}}\leftarrow A_{n_{1},n_{1}}+A_{1,1}^{k} \„... ... \„ & A_{n_{2},n_{2}}\leftarrow A_{n_{2},n_{2}}+A_{2,2}^{k}\end{eqnarray*}$$

Par exemple, sur le maillage de la figure (4.12), la matrice élémentaire $\mathbf{A}^{4}$sur l'élément $4$ contribue aux 4 coefficients suivants de $\mathbf{A}$:

$$\begin{eqnarray*} \mathbf{A}_{4,4}=\mathbf{A}_{4,4}+\mathbf{A}_{1,1}^{4}, & ... ...{4}, & \mathbf{A}_{5,5}=\mathbf{A}_{5,5}+\mathbf{A}_{2,2}^{4}\end{eqnarray*}$$

De façon générale, pour un élément $k$ de degré $d$ dont les numéros des noeuds d'interpolation sont notés $\{ n_{p}\}_{p=1,d+1}$, la matrice élémentaire $\mathbf{A}^{k}$ de dimension $(d+1,d+1)$ a une contribution dans les coefficients $A_{n_{p},n_{q}}$de la matrice globale $\mathbf{A}$ de dimension $(nn,nn)$. Plus précisément le coefficient $(p,q)$ de la matrice élémentaire $\mathbf{A}^{k}$ intervient dans le calcul du coefficient $(n_{p},n_{q})$ de la matrice globale $\mathbf{A}$, puisque $n_{p}$ est le numéro global du point d'interpolation $p$ de l'élément $k,$ et $n_{q}$ le numéro global du point d'interpolation $q$ .

$$\mathbf{A}_{n_{p},n_{q}}\leftarrow \mathbf{A}_{n_{p},n_{q}}+\mathbf{A}_{p,q}^{l}$$

L'algorithme d'assemblage général (1) est donnée ci dessous ainsi que le programme Maple 4.2.6 associé. On notera que l'on a fait varier les indices $p$ et $q$ à partir de 1 (et non de 0), pour tenir compte du fait que sous Maple (et Matlab) les indices des tableaux commencent à 1.

A l'issue de cet assemblage la matrice $\mathbf{A}$ et le second membre $\mathbf{B}$ calculée avec les paramètres de l'exemple précédent

$$L=3 m, D=2 cm, h=50 W/m^{2}/K, k=6000 W/m/K, T_{e}=60 C, T_{a}=20 C, \phi_{e}=32 W$$

s'écrivent pour un maillage de $ne=8$ éléments:

$$\mathbf{A}=\left[\begin{array}{ccccccccc} 507.0 & -501.0 & 0... ... 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & -501.0 & 507.0\end{array}\right]$$

$$\mathbf{B}=[118.0,236.0,236.0,236.0,236.0,236.0,236.0,236.0,118.0]$$

Cette matrice est symétrique et tri diagonale (car avec l'interpolation $\mathcal{P}^{1}$, une fonction de base associée à un noeud ou degré de liberté $x_{i}$ est non nulle sur l'intervalle $[x_{i-1},x_{i+1}]$).

Exercise: démontrer cette propriété

programme Maple 4.2.6: Assemblage de la matrice et du second membre

# Boucle d'assemblage
> A:=matrix(Nn,Nn,0): B:=vector(Nn,0):
> for k from 1 to Ne do
>   Ke:=MatRigidite(k,K);
>   Me:=MatMasse(k,alpha);
>   Be:=SmbElem(k,Vf);
>   noi:=num(k);
>   for p from 1 to d+1 do
>     ni:=noi[p];
>     for q from 1 to d+1 do
>       nj:=noi[q];
>       A[ni,nj]:=A[ni,nj]+Ke[p,q]+Me[p,q];
>     end;
>     B[ni]:=B[ni]+Be[p];
>   end;
> end:
> 

3.2.7 Prise en compte des conditions aux limites

Le calcul précédent est générique et ne tiens pas compte des conditions aux limites. Pour la condition aux limites de Dirichlet en $x=0$, on remplace la première équation par la condition:

$$T_{1}=T_{e}$$

ce qui reviens à annuler la première ligne de la matrice $\mathbf{A}$, puis à mettre un 1 sur le terme diagonale $\mathbf{A}_{11}$, et $T_{e}$ dans $\mathbf{B}_{11}$. Pour la condition aux limites en $x=L$ (condition mixte), il faut rajouter un terme dans la matrice $\mathbf{A}$ correspondant à :

$$\beta u_{Nn}\Phi_{i}(L)=\beta \delta_{i,Nn} u_{Nn}$$

qui est non nul uniquement pour $i=nn$. Il faut donc modifier uniquement le terme diagonale $A_{nn,nn}$:

$$A_{nn,nn}\leftarrow A_{nn,nn}+\beta$$

Pour le second membre $\mathbf{B}$, il faut rajouter le terme

$$-\phi_{0}\Phi_{i}(L)=-\phi_{0}\delta_{i,Nn}$$

qui est non nul uniquement pour $i=nn$. Il faut donc modifier uniquement le dernier terme de $\mathbf{B}:$

$$B_{nn}\leftarrow B_{nn}-\phi_{0}$$

Le petit programme Maple 4.2.7 correspondant s'écrit:

programme Maple 4.2.7: Application des conditions et résolution

# Conditions aux limites
> for j from 1 to Nn do A[1,j]:=0; end: A[1,1]:=1; B[1]:=ue;
> A[Nn,Nn]:=A[Nn,Nn]+beta; B[Nn]:=B[Nn]-phi;
# Résolution
> Uh:=evalf(linsolve(A,B));

Après application des conditions aux limites, la matrice $\mathbf{A}$ s'écrit:

$$\mathbf{A}=\left[\begin{array}{ccccccccc} 1.0 & 0.0 & 0.0 & ... ... 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & -501.0 & 507.0\end{array}\right]$$

et le second membre $\mathbf{B}$:

$$\mathbf{B}=[60.0,236.0,236.0,236.0,236.0,236.0,236.0,236.0,85.8]$$

On constate que la matrice $A$ n'est plus symétrique, à cause de la façon d'implémenter la condition aux limites de Dirichlet. Nous verrons dans le chapitre suivant comment imposer cette condition en conservant la symétrie.

3.2.8 Résultats

Avec les paramètres suivants donnés dans le programme Maple 4.2.8 ci-dessous, pour un maillage de $8$ éléments

programme Maple 4.2.8: Paramêtres du problème

# Parametres du problème
> L:=3; dd:=0.2;h:=50;Ta:=20;Te:=60; phie:=32;
> K:=evalf(6000.0*Pi*dd^2/4): alpha:=evalf(h*Pi*dd):
> f:=alpha*Ta: ue:=Te: beta:=0:phi:=phie:
# Vecteurs des parametres interpolés sur le maillage
> Vf:=Interpol(f):

la solution obtenue vaut:

$$\mathbf{X}=[60.0,55.3,51.3,48.2,45.7,43.8,42.4,41.6,41.3]$$

La comparaison de cette solution éléments finis avec la solution exacte est donnée sur la figure (4.15).

Figure 4.15: Solution éléments finis pour ne=8

On constate que l'erreur moyenne est très faible et vaut $0.083$ (soit $0.1\%$ en relatif), mais l'approximation de la dérivée en $x=L$ reste encore assez mauvaise: on trouve $158$ au lieu de $32$ (soit $390 \%$ d'erreur), ce qui est cependant meilleur qu'avec 3 éléments où on obtiens $374$. A titre de comparaison, avec le maillage de 3 éléments non régulièrement espacés du paragraphe 3.1, la solution est meilleure puisque l'on trouve un flux de $285$. On intuite ici l'intérêt en élément finis d'utiliser des maillages adaptés au problème.

3.2.9 Etude de la précision

Pour étudier la précision de la méthode des éléments finis $\mathcal{P}^{1}$, nous avons calculer l'erreur relative moyenne:

$$\frac{\vert\vert u-u^{h}\vert\vert}{\vert\vert u\vert\vert}=\sqrt{\frac{\int_{0}^{L}(u-u^{h})^{2}dx}{\int_{0}^{L}u^{2}dx}}$$

et l'erreur relative sur la condition aux limites en $x=L$:

$$\frac{\vert-K\left(\frac{du^{h}}{dx}\right)_{L}-\phi_{e}\vert}{\vert\phi_{e}\vert}$$

en fonction du nombre d'éléments $ne$ du maillage pour des maillages régulièrement espacés (i.e. la taille $h$ des éléments est proportionnelle à $1/ne$). Les résultats sont tracés en échelle logarithmique sur la figure (4.16). On constate que l'erreur relative moyenne varie en $\frac{1}{ne^{2}}$(i.e. en $\theta(h^{2})$) (comparaison avec une droite de pente $-2$), et que l'erreur relative sur la condition aux limites varie en $\frac{1}{ne}$(i.e. en $\theta(h)$ ). Ce résultat montre que l'approximation par éléments finis $u^{h}$ converge vers la solution exacte et que cette convergence est d'ordre 2. Ce résultat est cohérent avec l'erreur d'interpolation, qui comme nous l'avons vu est d'ordre 2 pour une approximation $\mathcal{P}^{1}$. On peut en fait démontrer que l'erreur par éléments finis est majoré par cette erreur d'interpolation.

Figure: Erreur relative en fonction de ne (elt $\mathcal{P}^{1}$)