3.1 Première approche

On veut calculer la température $T(x)$ dans un barreau de longueur $L$, de section $S$ et de conductivité $k,$ dont l'extrémité gauche est isotherme (i.e $T(0)=T_{e}$), et l'extrémité droite reçoit un flux de chaleur $\phi_{e}$. En plus de la conduction dans le solide (le flux de chaleur par conduction dans une section s'écrit $\phi_{1}=-kS\frac{dT}{dx}$), le barreau échange de la chaleur par convection avec l'air ambiant à la température $T_{a}$ sur toute sa longueur. En notant $h$ le coefficient d'échange par convection par unité de surface, le flux de chaleur par convection s'écrit pour un élément de longueur $dx$ : $\phi_{2}=h dx p (T-T_{a})$ (où $p$ est le périmètre de la section du barreau)).

Figure 4.1: Température dans un barreau

En notant $K=kS$, et $\alpha=h p$, l'équation d'équilibre s'écrit: trouver $T(x)$ solution du système (4.1 )

$$\left\{ \begin{array}{c} -\frac{d}{dx}(K\frac{dT}{dx})+\alph... ...\\ T(0)=T_{e}, -K\frac{dT}{dx}(L)=\phi_{e}\end{array}\right.$$ (3.1)

3.1.1 Formulation faible

La formulation faible de ce problème s'obtient en multipliant l'équation par une fonction test $v(x)$, et en intégrant sur le domaine d'étude. Il vient:

$$\int_{0}^{L}(-\frac{d}{dx}(K\frac{dT}{dx}) v(x) dx+\int_{0}^{L}\alpha T(x) v(x) dx=\int_{0}^{L}\alpha T_{a}v(x) dx$$

Pour tenir compte des conditions aux limites et symétriser le problème, on effectue une intégration par partie sur les termes de plus haut degré:

$$\int_{0}^{L}K\frac{dT}{dx} \frac{dv}{dx} dx-[K\frac{dT}{d... ...0}^{L}\alpha T(x) v(x) dx=\int_{0}^{L}\alpha T_{a}v(x) dx$$

Si le problème est bien posé, le terme de bord $[K\frac{dT}{dx}v(x)]_{0}^{L}$ doit pouvoir se calculer en fonction des conditions aux limites. Pour cela on interprète la fonction test $v(x)$ comme une variation $\delta T(x)$ de la solution $T(x)$. Si on fixe la valeur de la solution en un point (condition de Dirichlet), sa variation est nulle et la fonction test $v(x)$ doit s'annuler en ce point. Pour notre problème, on doit donc imposer à $v(x)$ de s'annuller en $x=0$ puisque la valeur de $T(x)$ est fixée en ce point. Par contre en $x=L$, on impose aucune contrainte sur $v(x)$ et on utilise la condition de flux pour calculer le terme de bord. C'est une condition à la limite naturelle ou condition de Neuman. Avec ces conditions, le terme de bord se réduit à $\phi_{e}*v(L)$, et la formulation faible s'écrit:

$$\left\{ \begin{array}{l} \mbox{Trouvez }T(x) \mbox{ ave... ...n test }v(x)\mbox{ vérifiant }v(0)=0\end{array}\right.$$ (3.2)

On remarque que dans cette formulation, seule la condition au limite en $x=0$ est imposé de façon explicite. On parle alors de condition forte. La condition en $x=L$ n'est pas imposée explicitement, mais est prise en compte dans la formulation intégrale. C'est une condition naturelle.

Exercice: Montrer que la formulation faible (4.2) est équivalente au problème de minimisation (4.3)

$$\left\{ \begin{array}{c} \mbox{Trouvez }T(x) \mbox{ ave... ...nt_{0}^{L}\alpha T_{a}T(x) dx+\phi_{e}T(L)\end{array}\right.$$ (3.3)

3.1.2 Interpolation par éléments finis

Pour résoudre le problème (4.2) (dont il n'existe en général pas de solution analytique), on recherche une solution numérique approchée $T^{h}(x)$. En éléments finis, cette solution approchée est construite à partir de 2 données:

1. un maillage $\mathcal{M}^{h}$du domaine de calcul $\Omega$,
2. un choix d'interpolation $\mathcal{P}^{k}$ sur ce maillage.

Pour notre domaine de calcul unidimensionnel $\Omega=[0,L]$, le maillage correspond à un découpage de $\Omega$ en ne segments:

$$\Omega=\bigcup_{i=1}^{ne}[x_{i-1},x_{i}]$$

Sur chaque segment, on choisit une interpolation polynômiale. Le type de l'interpolation dépend du problème à traiter. Pour notre problème elle doit en particulier respecter les points suivants:

1. les intégrales dans la formulation faible (4.2) ou variationnelle (4.3) doivent exister, sans que le problème dégénère,
2. on cherche d'autre part une solution approchée globalement continue.

calcul de l'interpolation par éléments finis

La première condition impose que sur chaque élément le polynôme d'interpolation $p(x)$ soit au moins de degré 1 pour pouvoir calculer les dérivées premières de la solution qui interviennent dans la formulation faible (4.2). Sur chaque élément $e_{k}=[x_{k-1},x_{k}]$, la solution approchée est donc un polynôme $P(x)$ de degré $l\ge1$.

Figure 4.2: élément $e_{k}$

Pour un polynôme d'ordre 1 sur $[x_{k-1},x_{k}]$ (éléments finis $\mathcal{P}^{1}$), on utilise 2 points d'interpolation: les 2 extrémités du segment $\{ S_{1}=x_{k-1}, S_{2}=x_{k}\}$ (figure 4.2).

Exercice: déterminer les points d'interpolation pour un polynôme de degré 2, puis 3.

Figure: Maillage et interpolation $\mathcal{P}^{1}$

Considérons par exemple le maillage du domaine $\Omega=[0,1]$ de la figure (4.3) suivant:

$$\mathcal{M}^{h}=[0,\frac{1}{2}]\bigcup[\frac{1}{2},\frac{3}{4}]\bigcup[\frac{3}{4},1]$$

L'interpolation $f^{h}(x)$ par éléments finis $\mathcal{P}^{1}$ d'une fonction $f(x)$ s'écrit:

$$f^{h}(x)=\left\{ \begin{array}{l} 2(f_{1}-f_{0})x+f_{0} \mb... ...f_{2} \mbox{ si } \frac{3}{4}\le x\le1\end{array}\right.$$ (3.4)

où $\{ f_{0},f_{1},f_{2},f_{3}\}$ sont les valeurs nodales de $f$ aux points de maillage $\{0,\frac{1}{2},\frac{3}{4}1\}$. C'est une fonction linéaire par morceau et continue sur l'intervalle d'étude $\Omega$. Sa dérivée $\frac{df^{h}}{dx}$s'écrit:

$$\frac{df^{h}}{dx}(x)=\left\{ \begin{array}{l} 2(f_{1}-f_{0})... ...3}-4f_{2})x \mbox{ si } \frac{3}{4}<x<1\end{array}\right.$$

C'est une fonction constante par morceau et discontinue aux points de maillage.

Exercice: calculer pour ce même maillage l'interpolation par éléments finis $\mathcal{P}^{2}$

Figure 4.4: Interpolation éléments finis d'une fonction (à gauche) et de sa dérivée (à droite) sur un maillage de 3 segments avec des polynômes $P_{1}$

Sur la figure (4.4), on a tracé l'interpolation par éléments finis de la fonction $f(x)=sin(0.8\pi x)$ et de sa dérivée sur ce maillage de 3 segments en utilisant des polynômes de degré 1. On constate que l'approximation $\mathcal{P}^{1}$ est bien continue, mais que la dérivée est discontinue aux points de maillage.

3.1.2.1 erreur d'interpolation

On rappelle que l'erreur entre une fonction $f(x)$ et son polynôme d'interpolation $P(x)$ de degré $l$ sur l'intervalle $[a,b]$ passant par les $l+1$ points d'interpolation $\{ x_{j}\}_{j=1,l+1}$ s'écrit:

$$f(x)-P(x)=\frac{\prod_{j=1}^{l+1}(x-x_{j})}{(l+1)!}f^{(l+1)}(\chi) \mbox{ avec }\chi\in[a,b]$$

Exercice: en étudiant la fonction $W(t)=f(t)-P(t)-\frac{f(x)-P(x)}{\prod_{j=1}^{l+1}(x-x_{j})}\prod_{j=1}^{l+1}(t-x_{j})$, démontrer la formule précédente.

Pour une approximation linéaire sur un segment $e_{k}=[x_{k-1},x_{k}]$, l'erreur d'interpolation s'écrit

$$f(x)-P(x)=\frac{(x-x_{k})(x-x_{k-1})}{2}\frac{d^{2}f}{dx^{2}}(\chi)$$

On vérifie que l'erreur s'annulle aux points d'interpolation et est proportionnelle à la dérivée seconde de $f(x).$

A partir de cette relation locale, on peut déduire par calcul directe les majorations d'erreurs suivantes pour la norme du maximum et la norme moyenne de l'erreur sur l'intervalle $[x_{k-1},x_{k}]$ de longueur $h=x_{k}-x_{k-1}$.

$$\begin{eqnarray*} \max_{x\in[x_{k-1},x_{k}]}(\vert f(x)-p(x)\vert) & \le & \frac... ...}\sqrt{\int_{x_{k-1}}^{x_{k}}(\frac{d^{2}f}{dx^{2}}(\chi))^{2}dx}\end{eqnarray*}$$

Exercice: en intégrant l'erreur d'interpolation, démontrer ces formules

On vérifie que l'erreur moyenne d'approximation par éléments finis $\mathcal{P}^{1}$ est en $\theta(h^{2})$, i.e. est proportionnelle au carré de la longueur des éléments en un point de l'intervalle.

Exercice: montrez que l'erreur avec des éléments finis $\mathcal{P}^{2}$ est en $\theta(h^{3})$

3.1.3 Approximation par éléments finis:

L'approximation par éléments finis est définie de façon locale sur chaque élément (4.4). De façon à pouvoir la manipuler plus facilement, on va l'exprimer de façon globale. Pour cela on détermine une expression générique de l'approximation sur un élément $e_{k}=[x_{k-1},x_{k}]$.

3.1.3.1 Fonctions de forme

Pour ce faire on introduit une transformation géométrique qui permet de passer d'un élément quelconque $[x_{k-1},x_{k}]$ à un élément de référence $[-1,1]$, sur lequel on va définir l'approximation de manière générique.

Figure: Transformation $\mathcal{T}^{k}$ : $x=x_{k-1}\frac{1-\xi}{2}+x_{k}\frac{1+\xi}{2}$

Cette transformation affine $\mathcal{T}^{k}$ s'écrit:

$$\mathcal{T}^{k} \left\{ \begin{array}{c} [x_{k-1},x_{k}]\Lo... ...tarrow\xi=2\frac{x-x_{k-1}}{x_{k}-x_{k-1}}-1\end{array}\right.$$ (3.5)

Pour une approximation $\mathcal{P}^{1}$, on a 2 points d'interpolation et tout polynôme de degré 1 $P(\xi)$ s'écrit sur l'intervalle de référence comme combinaision linéaire des 2 polynômes de Lagrange $N_{1}$ et $N_{2}$ associés à ces 2 points d'interpolations:

$$P(\xi)=P(-1) N_{1}(\xi)+P(1) N_{2}(\xi)$$

avec

$$N_{1}(\xi)=\frac{1-\xi}{2} \mbox{ et } N_{2}(\xi)=\frac{1+\xi}{2}$$ (3.6)

Ce sont les deux “fonctions de forme” de l'élément $\mathcal{P}^{1}$ . Le tracé de ces fonctions de forme est donné sur la figure (4.6).

Exercice: déterminer les 3 fonctions de formes pour un élément $\mathcal{P}^{2}$

Figure: fonctions de forme $\mathcal{P}^{1}$ ($N_{1}$ en rouge, $N_{2}$ en vert)

En utilisant le changement de variable $\xi=\xi(x)$ (4.5), un polynôme de degré 1 en x sur un élément $[x_{k-1},x_{k}]$ s'écrit:

$$P(x)=P(x_{k-1}) N_{1}(\xi(x))+P(x_{k}) N_{2}(\xi(x))$$

Attention: lorsque l'on calcule la dérivée de l'approximation élément finis, la dérivée dans l'élément $[x_{k-1},x_{k}]$ n'est pas égale à la dérivée dans l'élément de référence. Il faut tenir compte du changement de variable:

$$\frac{dP}{dx}(x)=\frac{dP}{d\xi}(\xi)\frac{d\xi}{dx} \mbox{... ...,\frac{d\xi}{dx}=\frac{2}{h} \mbox{ et } h=x_{k}-x_{k-1}$$

Pour un polynôme de degré 1, on obtient:

$$\begin{eqnarray*} \frac{dP}{dx}(x) & = & P(x_{k-1})\frac{dN_{1}}{d\xi}\frac{2}{h... ..._{2}}{d\xi}\frac{2}{h}\\ & = & \frac{1}{h}(P(x_{k})-P(x_{k-1}))\end{eqnarray*}$$

Exercice: calculer la dérivée d'une approximation élément finis $\mathcal{P}^{2}$

3.1.3.2 Fonctions de base

Avec ces nouvelles notations, l'approximation élément finis $f^{h}$ de la fonction f(x) dans la relation (4.4) s'écrit

$$f^{h}(x)=\left\{ \begin{array}{l} f_{0}N_{1}(4x-1)+f_{1}N_{2... ...(8x-7) \mbox{ si } \frac{3}{4}\le x\le1\end{array}\right.$$

on constate que $f^{h}$ est une fonction linéaire des 4 valeurs nodales $\{ f_{0},f_{1},f_{2},f_{3}\}$. On peut donc écrire $f^{h}(x)$ comme une combinaison linéaire de ces valeurs:

$$f^{h}(x)=f_{0}\Phi_{0}(x)+f_{1}\Phi_{1}(x)+f_{2}\Phi_{2}(x)+f_{3}\Phi_{3}(x)$$

où les fonctions $\{\Phi_{0},\Phi_{1},\Phi_{2},\Phi_{3}\}$ sont définies par

$$\Phi_{0}(x)=N_{1}(4x-1) \mbox{ si } x\le\frac{1}{2} \Phi_{0}(x)=0 \mbox{ sinon}$$

$$\Phi_{1}(x)=N_{2}(4x-1) \mbox{ si } x\le\frac{1}{2} \Phi_{1}(x)=N_{1}(8x-5) \mbox{ si }\frac{1}{2}\le x\le\frac{3}{4} \Phi_{1}(x)=0 \mbox{ sinon}$$

$$\Phi_{2}(x)=N_{2}(8x-5) \mbox{ si }\frac{1}{2}\le x\le\frac{3}{4} \Phi_{2}(x)=N_{1}(8x-7) \mbox{ si }\frac{3}{2}\le x\le1 \Phi_{2}(x)=0 \mbox{ sinon}$$

$$\Phi_{3}(x)=N_{2}(8x-7) \mbox{ si }\frac{3}{2}\le x\le1 \Phi_{3}(x)=0 \mbox{ sinon}$$ (3.7)

Ces fonctions $\Phi _{i}$ sont appelées “fonctions de base” de l'approximation éléments finis $\mathcal{P}^{1}$. Elles sont tracées sur la figure (4.7) .

Exercice: pour ce même maillage, déterminer les fonctions de base $\mathcal{P}^{2}$

Figure: fonctions de base $\mathcal{P}^{1}$: $\varphi _{0}$ (rouge), $\varphi _{1}$ (vert), $\varphi _{3}$ (bleu), $\varphi _{4}$ (cyan)

Ces fonctions de base vérifient les propriétés importantes suivantes:

1. à chaque noeud $S_{i}$ du maillage, on associe une fonction de base $\Phi _{i}$
2. la fonction de base $\Phi _{i}$ vaut 1 au noeud $S_{i}$ du maillage et 0 sur tout les autres noeuds $\{ S_{j}\}_{j\ne i}$

$$\Phi_{i}(S_{j})=\delta_{ij}$$ (3.8)

3.1.4 Formulation faible discrète

Nous allons maintenant chercher une approximation par éléments finis de la solution du problème (4.2). Pour cela on choisit le maillage suivant du domaine de calcul $\Omega=[0,L]$ (c'est le même maillage que précédemment, mais sur le domaine de calcul de longueur $L$ du problème (4.2)):

$$\mathcal{M}^{h}=[0,\frac{L}{2}]\cup[\frac{L}{2},\frac{3L}{4}]\cup[\frac{3L}{4},L]$$

On choisit ensuite sur ce maillage une approximation par élément finis $\mathcal{P}^{1}$ de la solution approchée $T^{h}$ de (4.2). Pour ce maillage de $ne=3$ éléments, le nombre de points nodaux est égale à $nn=4$, et $T^{h}$ s'écrit:

$$T^{h}(x)=T(0)\Phi_{0}(x)+T(\frac{L}{2})\Phi_{1}(x)+T(\frac{3L}{4})\Phi_{2}(x)+T(L)\Phi_{3}(x)$$

On notera $\{ T_{i}\}_{i=0,3}$ les valeurs nodales $\{ T(0),T(\frac{L}{2}),T(\frac{3L}{4}),T(1)\}$, ce qui permet d'écrire $T^{h}$ sous la forme:

$$T^{h}(x)=\sum_{i=0}^{3}T_{i}\Phi_{i}(x)$$

La solution du problème (4.2) doit vérifier les conditions aux limites fortes ( Dirichlet), i.e.:

$$T^{h}(x=0)=T_{0}=T_{e}$$

La valeur nodale $T_{0}$ est donc fixée, et $T^{h}$ s'écrit:

$$T^{h}(x)=T_{e}\Phi_{0}(x)+\sum_{i=1}^{3}T_{i}\Phi_{i}(x)$$ (3.9)

Le problème discrétisé possède donc 3 degrés de liberté (ou inconnues): les 3 valeurs nodales $\{ T_{i}\}_{i=1,3}$. De façon générale après application des $nc$ conditions aux limites fortes, la solution élément finis possède $n=nn-nc$ degré de liberté. Pour un maillage de $ne$ éléments $\mathcal{P}^{1}$, on a $n=ne+1-nc$.

En remplaçant la solution exacte par la solution approchée (4.9) dans l'équation (4.2) il vient la formulation faible discrète:

$$\int_{0}^{L}K\left(T_{e}\frac{d\Phi_{0}}{dx}+\sum_{j=1}^{3}T_{j}\... ...t_{0}^{L}\alpha\left(T_{e}\Phi_{0}+\sum_{j=1}^{3}T_{j}\Phi_{j}\right)v^{h} dx$$

$$=\int_{0}^{L}\alpha T_{a}v^{h} dx-\phi_{e}v^{h}(L)$$ (3.10)

Ayant la forme (4.9) de la solution approchée $T^{h}=T^{h}(T_{1},T_{2},T_{3})$, on peut en déduire les fonctions tests $v^{h}(x)$ associées, puisque ce sont des variations $\delta T^{h}$ de $T^{h}$ qui s'écrivent:

$$v^{h}(x)=\delta T^{h}=\sum_{i=1}^{3}\delta T_{i}\Phi_{i}(x)$$ (3.11)

On vérifie que ces fonctions tests s'annullent en $x=0$ (condition de Dirichlet). Ces fonctions tests sont des combinaisons linéaires des 3 fonctions de base $\{\Phi_{i}(x)\}_{i=1,3}$ associées aux 3 degrés de liberté $\{ T_{i}\}_{i=1,3}$ .

En choisissant comme fonctions tests $v^{h}$ dans (4.10), respectivement les 3 fonctions de base $\{\Phi_{i}(x)\}_{i=1,3}$, on obtient les 3 équations qui permettent de calculer les 3 inconnues $\{ T_{i}\}_{i=1,3}$:

$$\begin{eqnarray*} \int_{0}^{L}K(T_{e}\frac{d\Phi_{0}}{dx}+\sum_{j=1}^{3}T_{j}\fr... ..._{3}dx\\ =\int_{0}^{L}\alpha T_{a}\Phi_{3}dx-\phi_{e}\Phi_{3}(L)\end{eqnarray*}$$

C'est un système linéaire de 3 équations à 3 inconnues, qui s'écrit après regroupement des termes sous la forme matricielle: $\mathbf{A.X=B}$, où la matrice $\mathbf{A,}$ le second membre $\mathbf{B}$ et le vecteur inconnu $\mathbf{X}$ sont donnés par:

$$\mathbf{A}=\left[\begin{array}{ccc} \int_{0}^{L}(K\frac{d\Ph... ...frac{d\Phi_{3}}{dx}+\alpha\Phi_{3}\Phi_{3})dx\end{array}\right]$$

$$\mathbf{X}=\left[\begin{array}{c} T_{1}\\ T_{2}\\ T_{3}\end{array}\right]$$ (3.12)

$$\mathbf{B}=\left[\begin{array}{c} \int_{0}^{L}\alpha T_{a}\P... ...frac{d\Phi_{3}}{dx}+\alpha\Phi_{0}\Phi_{3})dx\end{array}\right]$$

soit sous forme générique:

$$\mathbf{A}_{ij}=\int_{0}^{L}K\frac{d\Phi_{j}}{dx}\frac{d\Phi_{i}}{dx}dx+\int_{0}^{L}\alpha\Phi_{j}\Phi_{i} dx\mbox{ et } i,j=1,2,3$$

$$\mathbf{B}_{i}=\int_{0}^{L}\alpha T_{a}\Phi_{i}dx-\phi_{e}\Phi_{i... ...int_{0}^{L}(K\frac{d\Phi_{0}}{dx}\frac{d\Phi_{i}}{dx}+\alpha\Phi_{0}\Phi_{i})dx$$ (3.13)

On remarque que la matrice $\mathbf{A}$ est symétrique puisque la formulation faible est elle même symétrique en $T$ et $v$.

La solution de la formulation faible discrète (4.10) s'obtient donc par résolution d'un système linéaire. Il nous reste à calculer la matrice $\mathbf{A}$ et le second membre $\mathbf{B}$. Connaissant l'expression des fonctions de base (4.7), un calcul directe des intégrales permettrait d'obtenir la matrice $\mathbf{A}$ et le vecteur $\mathbf{B}$. Cette approche n'est possible que pour un petit nombre de degré de liberté, et on préfère utiliser une approche systématique pour le calcul de $\mathbf{A}$ et $\mathbf{B}$ , qui s'appelle l'assemblage.

3.1.5 Assemblage

Pour le calcul des coefficients $A_{ij}$ et $B_{i}$, nous allons tout d'abord donner quelques propriétés des fonctions de base $\Phi _{i}$ (voir figure 4.8):

1. la fonction de base $\Phi_{i}(x)$ associé au noeud $S_{i}$ est non nulle uniquement sur les éléments du maillage auxquels appartient le noeud $S_{i}$ , c'est à dire sur une toute petite partie du maillage. On dit que $\Phi_{i}(x)$ a un support compact. Pour une fonction de forme $\Phi_{i}(x)$ de type $\mathcal{P}^{1}$, ce support est le segment $[x_{i-1}, x_{i}]\cup[x_{i}, x_{i+1}]$.
2. en conséquence, le produit $\Phi_{i}(x)*\Phi_{j}(x)$ de 2 fonctions de formes est en générale toujours nulle, sauf si l'intersection des supports est non vide (ce qui est rare car les supports sont petits). Dans ce dernier cas le produit est non nul uniquement sur l'intersection des 2 supports. Pour les fonctions de forme $\mathcal{P}^{1}$, on vérifie que:
   
   $$N_{i}(x)*N_{j}(x)=0 \mbox{ si } \vert i-j\vert>1$$
   
   

   Figure 4.8: fonctions de base $\Phi _{i}$ et $\Phi _{j}$

   

3.1.5.1 Assemblage de la matrice

D'après la relation (4.13), la calcul de $A_{ij}$ fait intervenir 2 types d'intégrale. On décompose alors $\mathbf{A}$ en 2 matrices:

$$\mathbf{A}=\mathbf{K}+\mathbf{M} \mbox{ avec } K_{ij}... ...hi_{i}}{dx}dx, M_{ij}=\int_{0}^{L}\alpha\Phi_{j}\Phi_{i} dx$$

la matrice $\mathbf{K}$ est appelée matrice de rigidité (ou de raideur) et $\mathbf{M}$ matrice de masse.

Pour calculer ces intégrales on décompose l'intégrale sur le domaine $[0,L]$ en somme d'intégrales élémentaires sur chaque élément $[x_{k-1},x_{k}]$.

$$K_{ij}=\sum_{k=1}^{ne}\int_{x_{k-1}}^{x_{k}}K\frac{d\Phi_{j}... ...sum_{k=1}^{ne}\int_{x_{k-1}}^{x_{k}}\alpha\Phi_{j}\Phi_{i} dx$$

On est donc ramené à un calcul d'intégrales élémentaires sur chaque élément. En utilisant les propriétés des fonctions de base $\Phi _{i}$, chaque intégrale élémentaire est non nulle sur un élément $[x_{k-1},x_{k}]$ si et seulement si les noeuds $S_{i}$ et $S_{j}$ appartiennent tous les deux à cet élément. Pour des éléments $\mathcal{P}^{1}$, on a 4 cas possibles:

• soit $i=k-1$ et $j=k-1 \mbox{ ou } k$
• soit $i=k$ et $j=k-1 \mbox{ ou } k$

Pour un élément $\mathcal{P}^{1}$, on a donc à calculer uniquement 4 intégrales élémentaires sur un élément $k$, que l'on peut écrire sous forme matricielle:

$$K_{pq}^{k}=\int_{x_{k-1}}^{x_{k}}K\frac{d\Phi_{n_{q}}}{dx}\frac{d\Phi_{n_{p}}}{dx}dx \mbox{ pour } p=1,2 q=1,2$$ (3.14)

$$M_{pq}^{k}=\int_{x_{k-1}}^{x_{k}}\alpha\Phi_{n_{q}}\Phi_{n_{p}}dx \mbox{ pour } p=1,2 q=1,2$$ (3.15)

où $n_{1}=k-1$ et $n_{2}=k$ sont les numéros des 2 sommets de l'élément $k$.

Figure 4.9: maillage élements finis et code de couleur pour les éléments

Avec ces notations, et en utilisant les codes de couleur de la figure (4.9) pour les éléments du maillage, le premier élément de la matrice $\mathbf{A}$ (4.12) s'écrit

$$A_{11}=\color{red}K \color{black}_{22}^{1}+\color{red}M \col... ...}K \color{black}_{11}^{2}+\color{green}M \color{black}_{11}^{2}$$

puisque le noeud 1 a pour support les éléments $[x_{0},x_{1}]$ ($k=1$) et $[x_{1},x_{2}]$ ($k=2$ ). Par contre l'élément $A_{12}$s'écrit:

$$A_{12}=\color{green}K \color{black}_{12}^{2}+\color{green}M \color{black}_{12}^{2}$$

On a donc l'expression de $\mathbf{A}$ sous la forme:

$$_{\mathbf{A}=\left[\begin{array}{ccc} \color{red}K_{22}^{1} ... ...lack} & \color{blue}M_{22}^{3} \color{black}\end{array}\right]}$$

Il reste donc à calculer ces matrices élémentaires pour calculer $\mathbf{A}$. Pour ce faire on revient à la définition (4.7) des fonctions de base $\Phi_{i}(x)$ en fonction des fonctions de forme $N_{q}(\xi)$ en effectuant le changement de variable (4.5) sur l'élément de référence. Sur un élément $e_{k}=$ $[x_{k-1},x_{k}]$ (de longueur $h_{k}$), on a:

$$\Phi_{k-1}(x)=N_{1}(\xi)=\frac{1-\xi}{2} \mbox{ et }\Phi_{k}(x)=N_{2}(\xi)=\frac{1+\xi}{2}$$ (3.16)

En effectuant le changement de variables $x=x(\xi)$ (4.5) dans l'intégrale, il vient:

$$K_{pq}^{k}=\frac{h_{k}}{2}\int_{-1}^{+1}K(\frac{dN_{q}}{d\xi... ...{h_{k}}{2}\int_{-1}^{+1}\alpha N_{p}N_{q}d\xi p,q=1,2$$

Ces deux matrices élémentaires sont symétriques. De plus la somme des coefficients des lignes (et donc des colonnes) de $\mathbf{K}^{l}$ est nulle parce que la somme des fonctions de forme $N_{1}(\xi)+N_{2}(\xi)$ est égale à 1.

Exercice: démontrer que la somme des lignes (et des colonnes) de la matrice de rigidité est nulle

Pour des coefficients $K$ et $\alpha$ constants, ces 2 matrices élémentaires s'écrivent, après le calcul simple des intégrales:

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

Exercice: démontrer les expressions précédentes

En notant que les éléments du maillage ont pour longueur : $\color{red}h_{1}=\frac{L}{2} \color{black}, \color{green}h_{2}=\frac{L}{4} \color{black}, \color{blue}h_{3}=\frac{L}{4} \color{black}$, on obtient la matrice A suivante (en coloriant en rouge, vert et bleu les contributions des éléments 1,2 et 3):

$$\mathbf{A}=\frac{K}{L}\left[\begin{array}{ccc} \color{red}2 ... ...ck} & \color{blue}\frac{1}{12} \color{black}\end{array}\right]$$ (3.17)

On vérifie la symétrie de la matrice.

3.1.5.2 assemblage du second membre

Le calcul du second membre procède de la même démarche. Le calcul des intégrales se fait élément par élément en tenant compte des propriétés des fonctions de base. L'expression du second membre $\mathbf{B}$ de (4.13) contient 2 types de termes:

$$\mathbf{B}_{i}=\underbrace{\int_{0}^{L}\alpha T_{a}\Phi_{i}d... ...{d\Phi_{i}}{dx}+\alpha\Phi_{0}\Phi_{i})dx}_{cdt limite x=0}$$

1. un terme source à calculer sur tous les éléments,
2. deux termes liés aux conditions aux limites en $x=0$ et $x=L$, que l'on ne calcule que pour certains éléments.

Pour calculer le terme source, on le décompose en une somme d'intégrales élémentaires:

$$\int_{0}^{L}\alpha T_{a}\Phi_{i}dx=\sum_{k=1}^{ne}\int_{x_{k-1}}^{x_{k}}\alpha T_{a}\Phi_{i} dx$$

et on utilise les propriétés des fonctions de base $\Phi _{i}$. Pour des éléments $\mathcal{P}^{1}$, ces intégrales élémentaires sur un élément $e_{k}$ sont non nulles uniquement si $i=k-1$ ou $i=k$ . On a donc à calculer 2 intégrales élémentaires par élément, qui s'écrivent sous la forme d'un vecteur second membre élémentaire:

$$B_{p}^{k}=\int_{x_{k-1}}^{x_{k}}\alpha T_{a}\Phi_{n_{p}}dx \mbox{ pour } p=1,2$$

où $n_{1}=k-1$ et $n_{2}=k$ sont les 2 numéros des sommets de l'éélment $e_{k}$.

Avec cette notation, le terme source du second membre $\mathbf{B}$ s'écrit:

$$\mathbf{B}=\left[\begin{array}{c} \color{red}B_{2}^{1} \colo... ...}\\ \color{blue}B_{3}^{2} \color{black}\end{array}\right]+C.L.$$

Pour le calcul des seconds membres élémentaires, on utilise le changement de variable $x=x(\xi)$ (4.5) sur l'élément de référence. En utilisant la relation (4.16), il vient:

$$B_{p}^{k}=\frac{h_{k}}{2}\int_{-1}^{+1}\alpha T_{a}N_{p-1}d\xi$$

Si le coefficient $\alpha$ et la température de l'air $T_{a}$ sont constants, le second membre élémentaire s'écrit:

$$\mathbf{B}_{p}^{k}=\alpha T_{a}h_{k}\left[\begin{array}{c} \frac{1}{2}\\ \frac{1}{2}\end{array}\right]$$

Exercice: démontrer l'expression précédente.

On en déduit la contribution du terme source dans le second membre $\mathbf{B}$

$$\mathbf{B}=\alpha T_{a}L\left[\begin{array}{c} \color{red}\f... ...\ \color{blue}\frac{1}{8} \color{black}\end{array}\right]+C.L.$$

La contribution des termes liés aux conditions aux limites n'interviens que sur certaines composantes de $\mathbf{B}$. La contribution de la condition au limite en $x=0$ fait intervenir le produit de la fonction de base $\Phi_{0}$ par une fonction de base $\Phi _{i}$. Comme il a été indiqué dans les propriétés des fonctions de base, ce produit est non nul uniquement pour $i=1$ (et on a donc uniquement une contribution dans $\mathbf{B}_{1}$). L'intégrale se calcule sur l'élément $[x_{0},x_{1}]$:

$$T_{e}\int_{0}^{L}(K\frac{d\Phi_{0}}{dx}\frac{d\Phi_{1}}{dx}+... ...ac{d\Phi_{0}}{dx}\frac{d\Phi_{1}}{dx}+\alpha\Phi_{0}\Phi_{1})dx$$

Ce terme correspond justement aux intégrales élémentaires (4.14) et (4.15) pour l'élément $e_{0}=[x_{0},x_{1}]$ multiplié par $T_{e}$:

$$T_{e}\int_{0}^{L}(K\frac{\Phi_{0}}{dx}\frac{d\Phi_{1}}{dx}+\alpha\Phi_{0}\Phi_{1})dx=(K_{21}^{1}+M_{21}^{1})T_{e}$$

La contribution de la condition aux limites en $x=L$ fait intervenir $\Phi_{i}(L)$, qui est non nul uniquement pour $i=3$. On a donc uniquement une contribution dans $B_{3}$ qui s'écrit:

$$\phi_{e}\Phi_{3}(L)=\phi_{e}$$

Le second membre complet s'écrit alors:

$$\mathbf{B}=\alpha T_{a}L\left[\begin{array}{c} \color{red}\f... ...+\mathbf{M}_{21}^{1})T_{e}\\ 0\\ -\phi_{e}\end{array}\right]$$ (3.18)

3.1.6 Résolution

Pour la résolution numérique, on considère un barreau d'aluminium de longueur $L=3 m$ , de diamètre $D=2 cm$, dont le coefficient de conductivité thermique vaut $k=6000 W/m/K$. Il est maintenu à une température $T_{e}=60^{\circ}C$ en $x=0$ et on impose un flux de $\phi_{e}=32 W$ en $x=L$. Il est refroidit dans l'air à température ambiante $T_{a}=25^{\circ}C$ par convection forcée avec $h=50 W/m^{2}/K$ (pour de la convection naturelle $h\simeq10 W/m^{2}/K$).

Pour vérifier le calcul précédent, nous allons tout d'abord ne pas tenir compte de la convection (i.e. $\alpha=0$). Dans ce cas la solution exacte de l'équation (4.1) est triviale. La répartition de température est linéaire et vérifie:

$$T(x)=T_{e}+\frac{\phi_{e}}{K}x$$

Le système linéaire discret s'écrit:

$$\left[\begin{array}{ccc} 376.99 & -251.33 & 0\\ -251.33 & 5... ...\left[\begin{array}{c} 7539.82\\ 0\\ -32.00\end{array}\right]$$

ce qui donne après résolution la répartition de température suivante:

$$\left[\begin{array}{c} T_{1}\\ T_{2}\\ T_{3}\end{array}\ri... ...left[\begin{array}{c} 59.74\\ 59.61\\ 59.49\end{array}\right]$$

C'est exactement la solution analytique aux noeuds du maillage $(x=\frac{1}{2},\frac{3}{4},1)$. On constate ainsi que si la solution exacte est linéaire, la solution par éléments finis est égale à la solution exacte. Cela est naturelle, car l'approximation par élément finis $\mathcal{P}^{1}$ approche exactement une solution linéaire. On note aussi que la variation de température dans la barre est très faible dans ce cas (de l'ordre $0.5^{\circ}C$), car le flux de chaleur en sortie $\phi_{e}$ est faible .

Dans le cas $\alpha\ne0$ et avec des coefficients constants, on peut encore déterminer la solution analytique. Les calculs sont un peu plus complexes, et on donne le programme Maple 4.1.6 ci dessous qui permet de calculer cette solution.

programme Maple 4.1.6: Solution annalytique

# Programme Maple pour le calcul de la solution annalytique
> restart:
# Equation du problème
> -K*diff(T(x),x$2)+alpha*T(x)=alpha*Ta;eq:=%:
# Solution générale
> dsolve(eq,T(x)):Tex:=unapply(rhs(%),x);
# Determination des constantes avec les conditions aux limites
> Tex(0)=Te;
> K*D(Tex)(L)=phi;
> S:=simplify(solve({%,%%},{_C1,_C2})):
> Tex:=unapply(simplify(subs(S,Tex(x)),trig),x);
# Application
> L:=3; K:=evalf(6000*Pi*0.2^2/4); alpha:=evalf(50.0*Pi*0.2);
> Ta:=20; Te:=60; phi:=-32;
> plot(Tex(x),x=0..L,title="solution exacte");

Avec ce programme, on obtient l'expression suivante pour la solution exacte $T(x)$:

$$T(x)=20+3.066 e^{0.4082x}+36.93 e^{-0.4082x}$$

En remplaçant les valeurs numériques dans la matrice (4.17) et le second membre (4.18), on obtient le système suivant pour la solution approchée par éléments finis:

$$\left[\begin{array}{ccc} 400.56 & -247.40 & 0.0\\ -247.40 &... ...[\begin{array}{c} 7775.44\\ 471.24\\ 203.62\end{array}\right]$$

ce qui donne après résolution la répartition de température approchée $T^{h}$:

$$\left[\begin{array}{c} T_{1}\\ T_{2}\\ T_{3}\end{array}\ri... ...left[\begin{array}{c} 45.51\\ 42.26\\ 41.12\end{array}\right]$$

La solution analytique pour ces mêmes points de calcul vaut:

$$\left[\begin{array}{c} T_{1}\\ T_{2}\\ T_{3}\end{array}\ri... ...left[\begin{array}{c} 45.67\\ 42.42\\ 41.28\end{array}\right]$$

La solution exacte étant de type exponentielle, l'approximation par éléments finis ne peut pas donner la solution exacte. On constate cependant que dans ce cas la solution élément finis est très précise, puisque l'erreur nodale est inférieure à $0.1^{\circ}C$ (i.e. un écart relatif de $\simeq0.4\%$ ).

Figure 4.10: Solution exacte (en rouge) et solution approchée (en bleu)

Quelques remarques:

1. Même si la solution par élément finis approche la solution exacte aux noeuds du maillage, entre les noeuds l'approximation est linéaire. Donc l'erreur entre la solution exacte et la solution élément finis peut être grande comme le montre la figure (4.10).
2. Si on calcule la dérivée de la solution approchée, elle est constante par élément (voir figure (4.11)). On en déduit immédiatement que la condition au limite en $x=L$ n'est pas vérifiée par la solution approchée:

$$-K\frac{dT^{h}}{dx}=-K\frac{T_{4}-T_{3}}{L/4}=285\ne\phi_{e}=32$$

Figure 4.11: Dérivée de la solution exacte (en rouge) et approchée (en bleu)

On montrera dans la suite que les conditions aux limites naturelles sont vérifiées exactement par la solution approchée que lorsque le maillage devient très fin: i.e. à la limite quand la solution approchée tends vers la solution exacte.