4.2 Première approche avec Matlab

On reprend l'exemple précédent, mais avec un cylindre de section carré. Les dimensions du problème sont indiquées sur la figure (4.1).

Figure 4.1: écoulement potentiel dans un canal

L'équation d'équilibre régissant la fonction de courant $\psi(x,y)$ n'est autre que l'équation de Laplace

$$\Delta\psi=\frac{\partial^{2}\psi}{\partial x^{2}}+\frac{\partial^{2}\psi}{\partial y^{2}}=0$$

où $\psi(x,y)$ est la fonction de courant telle que:

$$u=\frac{\partial\psi}{\partial y}, v=-\frac{\partial\psi}{\partial x}$$

$u$ et $v$ étant les composantes de la vitesse. Les conditions aux limites en vitesse pour ce problème sont:

1. en entrée, une vitesse constante suivant x : $u=U_{0}, v=0$
2. en sortie, la même vitesse qu'en entrée : $u=U_{0}, v=0$
3. une condition de glissement sur les parois: $\overrightarrow{u}.\overrightarrow{n}=0$

En utilisant la définition de $\psi$, on en déduit les conditions aux limites pour la fonction de courant (en choisissant une origine pour $\psi$ en $y=0$):

1. en entrée en $x=-5H$ : $\psi(-5H,y)=\int_{0}^{y}u dy=U_{0}y$
2. en sortie en $x=5H$ : $\psi(5H,y)=\int_{0}^{y}u dy=U_{0}y$
3. sur les parois : $\psi=cste$

Par raison de symétrie du domaine, des conditions aux limites et de l'équation du problème, on peut limiter le domaine d'étude au domaine $\Omega$ de frontière ABCDEF (soit $1/4$ du domaine initial). Les conditions aux limites sur ce domaine sont:

1. en entrée AF $(x=-5H)$ : $\psi_{AF}=U_{0}y$,
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(-5H,0)$, soit $\psi_{AB}=0,$
3. sur l'obstacle BCD, la valeur de $\psi$ est constante et égale à celle sur AB, soit $\psi_{BCD}=0$,
4. sur la paroi supérieure EF $(y=2H)$, la valeur de $\psi$ est constante et égale à $\psi(-5H,2H)$, soit $\psi_{EF}=2U_{0}H$
5. sur la sortie DE, on impose une condition de symétrie $(\frac{\partial\psi}{\partial n})_{DE}=(\frac{\partial\psi}{\partial x})_{x=0}=0.$

4.2.1 Formulation faible

La formulation faible du problème s'obtiens en multipliant l'équation par une fonction test $v(x,y)$, et en intégrant sur le domaine $\Omega$ . En utilisant la formule de Green, on intègre par partie le terme en Laplacien pour tenir compte des conditions aux limites et symétriser le problème. En utilisant les conditions aux limites, et en interprétant la fonction test $v(x,y)$ comme une variation $\delta\psi$ de la solution $\psi(x,y)$, la formulation faible s'écrit:

$$\left\{ \begin{array}{cc} \mbox{Trouver } \psi(x,y)\mbox{... ...,y) \mbox{ t.q.} v_{ABCD}=0, v_{EF}=0\end{array}\right.$$ (4.5)

Dans la suite, pour effectuer les calculs, nous choisirons comme valeurs numériques $H=1$ et $U_{0}=\frac{1}{2}$, ce qui impose un débit unité en entrée.

La formulation faible (4.5) est équivalente à la formulation variationnelle suivante:

$$\left\{ \begin{array}{cc} \mbox{Trouver } \psi(x,y) & \mb... ...}=U_{0}y, \psi_{ABCD}=0,\:\psi_{EF}=2U_{0}H\end{array}\right.$$ (4.6)

4.2.2 Interpolation par éléments finis $\mathcal{P}^{1}$

Pour résoudre numériquement le problème (4.5), on cherche une solution approchée $\psi^{h}(x,y)$. En éléments finis cette solution est définie par la donnée:

1. d'un maillage $\mathcal{M}^{h}$ du domaine de calcul $\Omega$,
2. d'une interpolation sur chaque élément du maillage.

Pour des éléments finis $\mathcal{P}^{1}$, le maillage est un maillage triangulaire et l'interpolation est une interpolation polynômial de degré 1 sur chaque élément, et continue globalement. Pour le problème considéré, nous choisissons le maillage de la figure (4.2) qui contient $ne=12$ éléments et $nn=11$ sommets.

Figure 4.2: maillage du canal de 11 noeuds et 12 éléments

La description d'un maillage comprend deux informations principales:

1. une information géométrique: les coordonnées de chacun des noeuds du maillage,
2. une information topologique: le numéro des 3 sommets de chacun des éléments du maillage, appelé table de connection.

Ces informations sont données par 2 tableaux: un tableau $\mathbf{X}$ de nombres réels de dimension $2nn$ pour les coordonnées, et un tableau $\mathbf{Tbc}$ de nombres entiers de dimension $3ne$ pour la table de connection. Pour le maillage de la figure (4.2), ces valeurs des deux tableaux sont données ci dessous (table 4.1 ).

Tableau 4.1: Coordonnées et Table de connection pour le maillage (4.2)

X	1	2 1

Tbc	1	2	3 1

Le domaine de calcul $\Omega$ est donc discrétisé en $ne$ éléments triangulaires $e_{k}$ :

$$\Omega=\mathcal{M}^{h}=\bigcup_{k=1}^{ne}e_{k} \mbox{ ... ...e } \mbox{ de sommets }\{ Tbc(k,1),Tbc(k,2),Tbc(k,3)\}$$

4.2.2.1 Interpolation sur un élément finis $\mathcal{P}^{1}$

Sur chaque élément $e_{k}$ du maillage, l'approximation $f^{h}$ par élément finis $\mathcal{P}^{1}$ d'une fonction $f(x,y)$ est un polynôme de degré 1 en $x$ et $y$, qui s'écrit:

$$f_{\vert e_{k}}^{h}=ax+by+c$$

Pour déterminer ce polynôme, il faut 3 points d'interpolation $\{ P_{j}\}_{j=1,3}$ et la valeur de la fonction $\{ F_{j}=f(P_{j})\}_{j=1,3}$ en ces 3 points. Les 3 équations $\{ P(x_{j},y_{j})=F_{j}\}_{j=1,3}$ permettent alors de déterminer les 3 coefficients de $f_{\vert e_{k}}^{h}$. Pour assurer la condition de continuité globale de la solution, on choisit comme points d'interpolation les 3 sommets $\{ S_{j}^{k}\}_{j=1,3}$ du triangle $k$ . L'approximation est donc définie par ses valeurs aux $nn$ noeuds du maillage.

Considérons par exemple la fonction $f(x,y)=(1+\frac{x}{10})*(\frac{y}{2})^{4}$ sur le maillage (4.2). Elle est définie par le tableau de valeurs nodales $F$ suivants:

$$F=[0.0437, 0.2848, 0.306, 0.3164, 1.0, 0.9, 0.5, 0.0, 0.0, 0.0563, 0.0625]$$

et elle est tracée en 3D et en iso-valeurs sur la figure (4.3).

Figure 4.3: Interpolation $\mathcal{P}^{1}$ en 2D

Sur l'élément 5 qui a pour sommets les noeuds $\{ n_{1}=2,n_{2}=6,n_{3}=1\}$, cette fonction est le polynôme de degré 1 qui prend les valeurs $F_{n_{1}}=0.2848, F_{n_{2}}=0.9, F_{n_{3}}=0.047$ sur ces 3 sommets. Dans l'espace c'est un plan qui est représenté en perspective sur la figure (4.4). Son expression analytique $f_{\vert e_{5}}^{h}=ax+by+c$ est obtenue en résolvant le système de 3 équations à 3 inconnues $\{ a,b,c\}$:

$$\left\{ \begin{array}{ccc} ax_{n_{1}}+by_{n_{1}}+c & = & F_{... ...}\\ ax_{n_{3}}+by_{n_{3}}+c & = & F_{n_{3}}\end{array}\right.$$ (4.7)

où $\{ x_{n_{1}},y_{n_{1}}\}$ sont les coordonnées du premier sommet $S_{1}$ (de numéro $n_{1}=2$), $\{ x_{n_{2}},y_{n_{2}}\}$ les coordonnées du second sommet $S_{2}$ (de numéro $n_{2}=6$), $\{ x_{n_{3}},y_{n_{3}}\}$ les coordonnées du troisième sommet $S_{3}$ (de numéro $n_{3}=1$).

La résolution du système (4.7) permet d'obtenir les expressions suivantes pour $\{ a,b,c\}$

$$\left\{ \begin{array}{c} a=\frac{(y_{n_{2}}-y_{n_{3}})F_{n_{... ...-(x_{n_{3}}-x_{n_{1}})(y_{n_{2}}-y_{n_{1}})}\end{array}\right.$$ (4.8)

Les relations (4.8) montrent que l'interpolation $\mathcal{P}^{1}$ de $f(x,y)$ sur l'élément $e_{k}$ ( $f_{\vert e_{k}}^{h}$) est une combinaison linéaire des valeurs nodales $\{ F_{n_{1}},F_{n_{2}},F_{n_{3}}\}$, qui s'écrit:

$$f_{\vert e_{k}}^{h}(x,y)=F_{n_{1}} p_{1}(x,y)+F_{n_{2}} p_{2}(x,y)+F_{n_{3}} p_{3}(x,y)$$ (4.9)

Les fonctions $\{ p_{1}(x,y), p_{2}(x,y), p_{3}(x,y)\}$ sont les 3 polynômes de degré 1, associés aux 3 sommets $\{ S_{1},S_{2},S_{3}\}$ de l'élément, suivants:

$$\left\{ \begin{array}{c} p_{1}(x,y)=\frac{(y_{n_{2}}-y_{n_{3... ...-(x_{n_{3}}-x_{n_{1}})(y_{n_{2}}-y_{n_{1}})}\end{array}\right.$$ (4.10)

Ces polynômes de degré 1 (dont on vérifie les symétries par permutation d'indices) vérifient :

$$p_{i}(S_{j})=\delta_{ij} \mbox{ pour } i,j=1,3$$

Figure 4.4: polynômes $p_{1},p_{2},p_{3}$ et interpolation $f^{h}$ sur l'élément 5 (2,6,1)

Ainsi le polynôme $p_{1}(x,y)$ vaut 1 sur le sommet $S_{1}$ et 0 sur les 2 autres sommets $S_{2}$ et $S_{3}$; de même pour $p_{2}(x,y)$ et $p_{3}(x,y)$ avec une permutation des indices. La représentation en perspective de ces 3 fonctions est donnée sur la figure (4.4) (pour le tracé on a remplacer les coordonnées des sommets par leurs valeurs obtenues dans le tableau 4.1: i.e.

$$\begin{eqnarray*} x_{n_{1}}=X(2,1)=-1.0, y_{n_{1}}=X(2,2)=1.5,\\ x_{n... ...2)=2.0,\\ x_{n_{3}}=X(1,1)=-3.0, y_{n_{3}}=X(1,2)=1.0,\end{eqnarray*}$$

Ces polynômes vérifient en outre la relation:

$$p_{1}(x,y)+p_{2}(x,y)+p_{3}(x,y)=1$$

Dans le cas d'une interpolation $\mathcal{P}^{1}$, ces polynômes ont une interprétation géométrique. Ce sont les coordonnées barycentriques. Ces coordonnées sont définies de la façon suivante: pour chaque point $M$ de coordonnées $(x,y)$, le vecteur $\overrightarrow{OM}$ s'écrit en fonction des sommets du triangle, comme combinaison des vecteurs $\overrightarrow{OS_{1}}$, $\overrightarrow{OS_{2}}$, $\overrightarrow{OS_{3}}$. Les coefficients sont les coordonnées barycentriques par rapport au triangle considéré, i.e.

$$\overrightarrow{OM}=\lambda_{1}\overrightarrow{OS_{1}}+\lamb... ... \mbox{ avec } \lambda_{1}+\lambda_{2}+\lambda_{3}=1$$ (4.11)

Les valeurs de $\lambda_{1},\lambda_{2},\lambda_{3}$ vérifient les relations (voir la figure 4.5 pour les notations):

$$\lambda_{1}=\frac{Aire(MS_{2}S_{3})}{Aire(S_{1}S_{2}S_{3})}=\frac{\frac{1}{2}\overline{MH_{1}}\overline{.S_{2}S_{3}}}{Aire(S_{1}S_{2}S_{3})}$$

$$\lambda_{2}=\frac{Aire(MS_{3}S_{1})}{Aire(S_{1}S_{2}S_{3})}=\frac{\frac{1}{2}\overline{MH_{2}}\overline{.S_{3}S_{1}}}{Aire(S_{1}S_{2}S_{3})}$$ (4.12)

$$\lambda_{3}=\frac{Aire(MS_{1}S_{2})}{Aire(S_{1}S_{2}S_{3})}=\frac{\frac{1}{2}\overline{MH_{3}}\overline{.S_{1}S_{2}}}{Aire(S_{1}S_{2}S_{3})}$$

qui sont justement les expressions (4.10) des polynômes $\{ p_{1},p_{2},p_{3}\}$.

On a aussi une propriété remarquable des coordonnées barycentriques: un point $M$ est à l'intérieur du triangle si et seulement si ses 3 coordonnées barycentriques sont positives, il est sur un des cotés du triangle si la coordonnée barycentrique par rapport au coté opposé est nulle, et il est à l'extérieur du triangle si une au moins des coordonnées barycentriques est négative.

Figure 4.5: coordonnées barycentriques

4.2.3 Approximation par éléments finis $\mathcal{P}^{1}$

L'approximation par éléments finis est donc définie de façon locale sur chaque élément, en calculant des formules d'interpolation du type (4.9) et (4.10). De façon a obtenir une expression générique pour l'interpolation, on va, comme dans le chapitre précédent introduire une transformation d'un élément $e_{k}$ vers un élément de référence. Cet élément de référence est le triangle rectangle unité dans le plan de référence $(\xi,\eta)$ (figure 4.6).

4.2.3.1 interpolation sur l'élément de référence et fonctions de forme

Figure 4.6: transformation $\mathcal{T}_{k}: (x,y)\Leftrightarrow(\xi,\eta)$ vers l'élément de référence $\mathcal{P}^{1}$

Cette transformation $\mathcal{T}^{k}$ est une transformation affine:

$$\mathcal{T}_{k} \left\{ \begin{array}{ccc} e^{k} & \Longlef... ...egin{array}{c} \xi\\ \eta\end{array}\right]\end{array}\right.$$ (4.13)

dont l'expression de $(\xi,\eta)$ en fonction de $(x,y)$ s'écrit :

$$\left[\begin{array}{c} \xi\\ \eta\end{array}\right]=\left[\... ...\right]+\left[\begin{array}{c} b_{1}\\ b_{2}\end{array}\right]$$

Les 6 coefficients de la transformation sont déterminés par les conditions de transformation des 3 sommets $\{ S_{i}\}_{i=1,3}$ de $e_{k}$ vers les 3 sommets $\{\hat{S}_{i}\}_{i=1,3}$ de l'élément de référence $\hat{e}$:

$$\left\{ \hat{S}_{i}=T_{k}(S_{i})\right\} _{i=1,3}$$

On remarque que $\xi(x,y)$ est un polynôme de degré 1 en $(x,y)$, qui vaut 0 en $S_{1}$ (car $\xi(\hat{S}_{1})=0$) , 1 en $S_{2}$ (car $\xi(\hat{S}_{2})=1$), et 0 en $S_{3}$ (car $\xi(\hat{S}_{3})=0$). C'est donc le polynôme d'interpolation $p_{2}(x,y)$ (ou la coordonnée barycentrique $\lambda_{2}(x,y)$) associé au sommet $S_{2}$, dont l'expression est donnée en (4.10). De même $\eta(x,y)$ est un polynôme de degré 1 en $(x,y)$, qui vaut 0 en $S_{1}$ (car $\eta(\hat{S}_{1})=0$), 0 en $S_{2}$ (car $\eta(\hat{S}_{2})=0$), et 1 en $S_{3}$ (car $\eta(\hat{S}_{3})=1$). C'est donc le polynôme d'interpolation de Lagrange$p_{3}(x,y)$ (ou la coordonnée barycentrique $\lambda_{3}(x,y)$) associé au sommet $S_{3}$. On a donc:

$$\xi=\lambda_{2}(x,y), \eta=\lambda_{3}(x,y)$$

et

$$1-\xi-\eta=\lambda_{1}(x,y)$$

En notant $(x_{n_{q}},y_{n_{q}})$ les coordonnées du sommet $S_{q}$ de l'élément $e_{k}$, on obtiens l'expression générique de la transformation de $e_{k}$ vers l'élément de référence $\hat{e}$:

$$\xi=\frac{(y_{n_{3}}-y_{n_{1}}) x+(x_{n_{1}}-x_{n_{3}}) y+(x_... ...2}}-x_{n_{1}})(y_{n_{3}}-y_{n_{1}})-(x_{n_{3}}-x_{n_{1}})(y_{n_{2}}-y_{n_{1}})}$$ (4.14)

$$\eta=\frac{(y_{n_{1}}-y_{n_{2}}) x+(x_{n_{2}}-x_{n_{1}}) y+(x... ...2}}-x_{n_{1}})(y_{n_{3}}-y_{n_{1}})-(x_{n_{3}}-x_{n_{1}})(y_{n_{2}}-y_{n_{1}})}$$

Sur cet élément de référence, les 3 polynômes d'interpolation $\{ p_{1},p_{2},p_{3}\}$ ont une expression très simple: ce sont les fonctions forme $\{ N_{1},N_{2},N_{3}\}$ de l'approximation $\mathcal{P}^{1}$ sur l'élément de référence associées aux 3 sommets $\{\hat{S}_{1},\hat{S}_{2},\hat{S}_{3}\}$:

$$N_{1}(\xi,\eta)=1-\xi-\eta, N_{2}(\xi,\eta)=\xi, N_{3}(\xi,\eta)=\eta$$ (4.15)

Pour les calculs d'intégrale, on aura besoin du changement de variables $x=x(\xi \eta)$ et $y=y(\xi,\eta)$, qui correspond donc à la transformation réciproque. Cette transformation étant affine, $x(\xi \eta)$ est une combinaison linéaire des valeurs aux noeuds, i.e. des abscisses $\{ x_{n_{q}}\}$ des sommets de l'élément $e_{k}$, et de même pour $y(\xi,\eta)$. On a donc

$$x=x_{n_{1}}N_{1}(\xi,\eta)+x_{n_{2}}N_{2}(\xi,\eta)+x_{n_{3}}N_{3}(\xi,\eta)$$ (4.16)

$$y=y_{n_{1}}N_{1}(\xi,\eta)+y_{n_{2}}N_{2}(\xi,\eta)+y_{n_{3}}N_{3}(\xi,\eta)$$

La matrice jacobienne $\mathbf{J}_{k}$= $\frac{D(x,y)}{D(\xi,\eta)}$ de cette transformation se calcule simplement:

$$\mathbf{J}_{k}=\left[\begin{array}{cc} \frac{\partial x}{\pa... ...y_{n_{2}}-y_{n_{1}}) & (y_{n_{3}}-x_{n_{1}})\end{array}\right]$$ (4.17)

De même la matrice jacobienne de la transformation inverse $\frac{D(\xi,\eta)}{D(x,y)}$ est égale à l'inverse de $\mathbf{J}^{k}$:

$$\mathbf{J}_{k}^{-1}=\left[\begin{array}{cc} \frac{\partial\x... ...y_{n_{1}}-y_{n_{2}}) & (x_{n_{2}}-x_{n_{1}})\end{array}\right]$$ (4.18)

en notant $aire_{k}=\frac{1}{2}((x_{n_{2}}-x_{n_{1}})(y_{n_{3}}-y_{n_{1}})-(x_{n_{3}}-x_{n_{1}})(y_{n_{2}}-y_{n_{1}}))$ l'aire du triangle $e^{k}$

En notant $\{ n_{1},n_{2},n_{3}\}$ les numéros des 3 sommets de l'élément $e_{k}$ qui sont donnés par la table de connection ( $n_{1}=Tbc(k,1)$, $n_{2}=Tbc(k,2)$, $n_{3}=Tbc(k,3)$), l'approximation $f^{h}$ s'écrit sur l'élément $e_{k}$:

$$f_{\vert e_{k}}^{h}(\xi,\eta)=F_{n_{1}}N_{1}(\xi,\eta)+F_{n_{2}}N_{2}(\xi,\eta)+F_{n_{3}}N_{3}(\xi,\eta)$$ (4.19)

Cette expression a la même forme que la relation (4.9), mais l'expression (4.15) des fonctions de forme est beaucoup plus simple que l'expression (4.10) des polynômes d'interpolation, ce qui va nous permettre en particulier un calcul plus simple des intégrales dans la formulation faible.

Attention: la dérivée d'une fonction dans l'élément de référence n'est évidemment pas égale à la dérivée dans le plan physique $(x,y)$. Il faut tenir compte du changement de variable:

$$\frac{\partial f}{\partial x}=\frac{\partial f}{\partial\xi}... ...+\frac{\partial f}{\partial\eta}\frac{\partial\eta}{\partial x}$$

qui s'écrit en fonction de la transposée de l'inverse de la matrice jacobienne $\mathbf{J}^{k}$. Le vecteur gradient s'écrit:

$$\overrightarrow{\nabla_{x,y}}f=(\mathbf{J}_{k}^{-1})^{t}\overrightarrow{\nabla_{\xi,\eta}}f$$

$$\left[\begin{array}{c} \frac{\partial f}{\partial x}\\ \fra... ...rtial\xi}\\ \frac{\partial f}{\partial\eta}\end{array}\right]$$ (4.20)

4.2.3.2 fonctions de base

Nous avons montré que l'approximation $f^{h}(x,y)$ par éléments finis $\mathcal{P}^{1}$ d'une fonction $f(x,y)$ sur le maillage (4.2) était déterminée par les valeurs nodales $\{ F_{i}\}$ de $f$ aux $nn=11$ noeuds $\{ M_{i}\}$ du maillage. Sur chaque élément, $f^{h}$ est un polynôme de degré 1 donnée par l'expression (4.19), qui est une fonction linéaire des 3 valeurs aux sommets de l'élément.

On peut donc écrire l'approximation $f^{h}(x,y)$ comme une combinaison linéaire des valeurs nodales $\{ F_{i}=f^{h}(M_{i})\}_{i=1,nn}$ :

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

Les fonctions $\Phi_{i}(x,y)$ sont les fonctions de base de l'approximation. Elles sont telles que sur chaque élément $e_{k}$ on retrouve l'approximation (4.19), ce sont donc des polynômes de degré 1 en $(x,y)$. De plus elles vérifient la propriété suivante pour chaque noeud $M_{j}$ du maillage de coordonnées $(x_{j},y_{j})$:

$$\Phi_{i}(M_{j})=\Phi_{i}(x_{j},y_{j})=\delta_{ij}$$

Autrement dit la fonction de base $\Phi_{i}(x,y)$ est une fonction continue qui vaut 1 au noeud $i$ du maillage, 0 sur tout les autres noeuds, et qui sur chaque élément est un polynôme de degré 1. On en déduit que sur un élément $e_{k}$ dont le noeud $i$ n'est pas un sommet, la fonction $\Phi_{i}(x,y)$ est identiquement nulle (car elle est nulle sur les trois sommets). Le support de la fonction $\Phi_{i}(x,y)$ (i.e. le lieu des points où la fonction est non nulle) se limite donc aux éléments du maillage ayant le noeud $i$ pour sommet:

$$\Phi_{i}(x,y)=\left\{ \begin{array}{ccc} 0 & \mbox{sur\, l'\... ..., l'\'{e}l\'{e}ment}\, k\, tq. & M_{i}=S_{q}\end{array}\right.$$ (4.21)

Ainsi la fonction de base $\Phi _{2}(x,y)$ est non nulle uniquement sur les éléments $\{5,6,7,8\}$ et vaut

$$\Phi_{2}(x,y)=\left\{ \begin{array}{cc} N_{1}(\xi,\eta) & \m... ...eta) & \mbox{sur} e_{8}\\ 0 & \mbox{sinon}\end{array}\right.$$

Figure 4.7: fonction de base $\Phi _{2}(x,y)$

En effet le sommet 2 est le premier sommet sur l'élément $e_{5}$ dans la table de connection (4.1), (i.e. Tbc(5,1)=2); le troisième sur l'élément $e_{6}$ (i.e. Tbc(6,3)=2), le second sur l'élément $e_{7}$ (i.e. Tbc(7,2)=1), et enfin le premier sur l'élément $e_{8}$ (i.e. Tbc(8,1)=2). Cette fonction de base $\Phi _{2}(x,y)$ est tracée en perspective sur la figure (4.7), et sa forme est une forme pyramidale.

4.2.4 Formulation faible discrète

La solution approchée $\psi^{h}$ du problème (4.5) est donc définie à partir de ces $nn=11$ valeurs nodales $\{\psi_{i}\}_{i=1,nn}$ aux sommets du maillage de la figure (4.2).

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

La solution du problème (4.5) doit de plus vérifiée les conditions aux limites fortes, i.e.:

$$\psi_{ABCD}^{h}=0,\:\psi_{EF}^{h}=2U_{0}H, \Psi_{AF}^{h}=U_{0}y$$

donc la valeur de $\psi^{h}$ est fixé sur les noeuds du maillage se trouvant sur ces frontières:

$$\psi_{8}=\psi_{9}=\psi_{10}=\psi_{11}=0, \psi_{5}=\psi_{6}=\psi_{7}=\psi^{e} \mbox{ avec } \psi^{e}=2U_{0}H$$

Grâce à un choix judicieux de la numérotation des noeuds du maillage, la solution approchée $\psi^{h}$ s'écrit:

$$\psi^{h}(x,y)=\sum_{j=1}^{4}\psi_{j}\Phi_{j}(x,y) + \psi^{e}*(\Phi_{5}(x,y)+\Phi_{6}(x,y)+\Phi_{7}(x,y))$$ (4.22)

Après application des conditions aux limites de Dirichlet, le problème discrétisé possède 4 degrés de liberté. En remplaçant la solution exacte par la solution approchée (4.22) dans (4.5) , il vient la formulation faible discrète:

$$\int_{\Omega}\left\{ \frac{\partial}{\partial x}\left(\sum_{j=1}^... ...i^{e}\sum_{j=5}^{7}\Phi_{j}(x,y)\right)\frac{\partial v^{h}}{\partial x}\right. +$$

$$\left.\frac{\partial}{\partial y}\left(\sum_{j=1}^{4}\psi_{j}\Phi... ...sum_{j=5}^{7}\Phi_{j}(x,y)\right)\frac{\partial v^{h}}{\partial y}\right\} dxdy = 0$$ (4.23)

Les fonctions tests $v^{h}$sont déduites de la forme de la solution approchée (4.22), puisque ce sont des variations $\delta\psi^{h}$ de $\psi^{h}$:

$$v^{h}(x,y)=\delta\psi^{h}=\sum_{i=1}^{4}\delta\psi_{i} \Phi_{i}(x,y)$$

On vérifie que ces fonctions s'annulent sur les frontières de Dirichlet. Ce sont des combinaisons linéaires des 4 fonctions de base $\{\Phi_{i}\}_{i=1,4}$ associées aux 4 degrés de liberté $\{\psi_{i}\}_{i=1,4}$.

En choisissant comme fonctions tests $v^{h}$ dans (4.23), respectivement ces 4 fonctions de base $\{\Phi_{i}\}_{i=1,4}$, on obtiens les 4 équations qui vont permettre de calculer les 4 inconnues $\{\psi_{i}\}_{i=1,4}$:

$$\begin{eqnarray*} \int_{\Omega}\left[\frac{\partial}{\partial x}\left(\sum_{j=1}... ...ght)\frac{\partial\Phi_{i}}{\partial y}\right]dxdy & =0 & (i=1,4)\end{eqnarray*}$$

Après regroupement des termes, en permutant l'intégration et la sommation, les équations s'écrivent:

$$\sum_{j=1}^{4}\psi_{j}\int_{\Omega}\left(\frac{\partial\Phi_... ...}}{\partial y}\frac{\partial\Phi_{i}}{\partial y}\right)d\omega$$

C'est un système linéaire de 4 équations à 4 inconnues, qui s'écrit sous la forme matricielle: $\mathbf{AX}=\mathbf{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}{cccc} \int_{\Omega}\overright... ...hi_{4}.\overrightarrow{\nabla}\Phi_{4}d\omega\end{array}\right]$$

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

$$\mathbf{B}=-\psi^{e}\left[\begin{array}{c} \sum_{j=5}^{7}\in... ...hi_{j}.\overrightarrow{\nabla}\Phi_{4}d\omega\end{array}\right]$$

soit sous forme générique:

$$\mathbf{A}_{ij}=\int_{\Omega}\left(\frac{\partial\Phi_{j}}{\parti... ...rac{\partial\Phi_{j}}{\partial y}\frac{\partial\Phi_{i}}{\partial y}\right)dxdy i=1,4 j=1,4$$

$$\mathbf{B}_{i}=-\sum_{j=5}^{7}\psi^{e}\int_{\Omega}\left(\frac{\p... ...rac{\partial\Phi_{j}}{\partial y}\frac{\partial\Phi_{i}}{\partial y}\right)dxdy i=1,4$$ (4.25)

On remarque que la matrice $\mathbf{A}$ est symétrique due à la symétrie de la formulation faible.

La recherche d'une solution approchée par éléments finis dans la formulation faible se ramène donc à la résolution du système linéaire (4.24). Il nous reste donc à calculer la matrice $\mathbf{A}$ et le second membre $\mathbf{B}$. A nouveau, comme dans le chapitre précédent, nous n'allons pas calculer directement les intégrales dans les relations (4.25), mais suivre une approche systématique pour le calcul de $\mathbf{A}$ et $\mathbf{B}$

4.2.4.1 assemblage de la matrice

Le calcul des coefficients de $\mathbf{A}$ se fait élément par élément, en notant que l'intégrale sur le domaine $\Omega$ est la somme d'intégrales élémentaires sur chacun des triangles $e_{k}$ du maillage:

$$\mathbf{A}_{ij}=\sum_{k=1}^{ne}\int_{e_{k}}\left(\frac{\part... ...rtial\Phi_{i}}{\partial y}\right)dxdy=\sum_{l=1}^{ne}A_{ij}^{k}$$

En utilisant la propriété des fonctions de base $N_{i}$ qui sont non nulles uniquement sur le support du noeud i, on constate que les intégrales élémentaires $A_{ij}^{k}$ sont presque toujours nulles sauf si l'élément $e_{k}$ est dans le support du noeud i et du noeud j, c'est à dire si i et j sont des sommets de l'élément $e_{l}$.

On a donc en réalité $3*3=9$ intégrales élémentaires à calculer par élément $e_{k}$, ce sont, en notant $\{ n_{1},n_{2},n_{3}\}$ les numéros des 3 sommets de l'élément k:

$$\mathbf{A}_{pq}^{k}=\int_{e_{k}}\left(\frac{\partial\Phi_{n_... ...{q}}}{\partial y}\right)dxdy \mbox{ pour }p=1,3 q=1,3$$ (4.26)

Avec ces notations, le premier coefficient de la matrice $\mathbf{A}$ s'écrit

$$\mathbf{A}_{11}=\mathbf{A}_{11}^{1}+\mathbf{A}_{33}^{2}+\mat... ...{3}+\mathbf{A}_{11}^{4}+\mathbf{A}_{33}^{5}+\mathbf{A}_{11}^{6}$$

puisque le noeud 1 a pour support les éléments $\{ e_{1},e_{2},e_{3},e_{4},e_{5},e_{6}\}$ et correspond au premier sommet sur l'élément $e_{1}$, au troisième sur $e_{2}$, .... De même le second coefficient de $\mathbf{A}$ s'écrit:

$$\mathbf{A}_{12}=\mathbf{A}_{31}^{5}+\mathbf{A}_{13}^{6}$$

puisque les seuls éléments ayant comme sommet les noeuds 1 et 2 sont les éléments $\{ e_{5},e_{6}\}$. Sur l'élément $e_{5}$ le noeud 1 correspond au troisième sommet et le noeud 2 au premier, alors que sur l'élément $e_{6}$ le noeud 1 correspond au premier sommet et le noeud 2 au troisième.

L'assemblage complet de la matrice $\mathbf{A}$ donne donc:

$$\mathbf{A}=\left[\begin{array}{cccc} A_{11}^{1}+A_{33}^{2}+A... ...}^{10}+A_{21}^{12} & A_{22}^{10}+A_{22}^{12}\end{array}\right]$$ (4.27)

Pour calculer les intégrales élémentaires (4.26), on effectue la transformation vers l'élément de référence.

Pour calculer l'intégrale d'une fonction $f(x,y)$ sur un élément $e_{k}$, on effectue le changement de variable $x=x(\xi \eta)$ et $y=y(\xi,\eta)$ pour se ramener à un calcul d'intégrale sur l'élément de référence. Le calcul de cette intégrale sur l'élément de référence s'effectue par partie, en intégrant d'abord en $\eta$ puis en $\xi$ (figure 4.8). On a donc:

$$\int_{e_{k}}f(x,y) dxdy=\int_{0}^{1}\int_{0}^{1-\xi}f(\xi,\eta) \det(\mathbf{J}_{k}) d\eta d\xi$$

puisque l'on a la relation suivante entre les éléments de surface:

$$dxdy=\det(\mathbf{J}_{k}) d\xi d\eta$$

où $\mathbf{J}_{k}$ est le jacobien (4.17) de la transformation $T_{k}^{-1}$ de l'élément de référence $\hat{e}$ vers l'élément $e_{l}$ vers . En utilisant la relation (4.20) pour le calcul des dérivées, la matrice élémentaire s'écrit après ce changement de variable:

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

Figure 4.8: Intégration sur l'élément de référence

En notant que la matrice jacobienne $\mathbf{J}_{k}$ est constante ainsi que les gradients des fonctions de forme, et que la surface de l'élément de référence vaut:

$$\int_{0}^{1}\int_{0}^{1-\xi}d\eta d\xi=\mbox{aire du triangle}=\frac{1}{2},$$

l'intégrale se simplifie:

$$\mathbf{A}_{pq}^{k}=\frac{\det(\mathbf{J}_{k})}{2}\left[\beg... ...l\xi}\\ \frac{\partial N_{q}}{\partial\eta}\end{array}\right]$$ (4.29)

Un calcul directe du déterminant du Jacobien fournit la valeur:

$$\det(\mathbf{J}^{k})=2 aire_{k}$$

que l'on peut vérifier en notant que ce déterminant est le rapport entre l'aire du triangle $e_{k}$ : $aire_{k}$, et l'aire du triangle de référence $\hat{e}$ : $1/2$ et. De même par un calcul directe le produit matriciel $(\mathbf{J}_{k}^{-1}).(\mathbf{J}_{k}^{-1})^{t}$ s'écrit:

(

$$(\mathbf{J}_{k}^{-1}).(\mathbf{J}_{k}^{-1})^{t}=\frac{1}{(2a... ...\ (y_{3}^{k}-y_{1}^{k})(y_{1}^{k}-y_{2}^{k})\end{array}\right]$$

Compte tenu de l'expression (4.15) des fonctions de forme, le calcul de leur gradient est trivial:

$$\overrightarrow{\nabla}N_{1}=\left[\begin{array}{c} -1\\ -1... ...ow{\nabla}N_{3}=\left[\begin{array}{c} 0\\ 1\end{array}\right]$$

Pour calculer la matrice élémentaire, on peut effectuer le calcul directe des 9 coefficients à partir de la relation (4.29), mais on peut aussi remarquer que la matrice élémentaire est symétrique, et que la somme des lignes (et des colonnes) est nulle (car la somme des gradients des fonctions de forme est nulle). Il suffit donc de calculer 3 coefficients: $A_{22}^{k}, A_{33}^{k}, A_{23}^{k}$ , les autres étant déduits comme indiqué ci dessous:

$$\mathbf{A}^{k}=\left[\begin{array}{ccc} \mathbf{A}_{22}^{k}+... ... & \mathbf{A}_{23}^{k} & \mathbf{A}_{33}^{k}\end{array}\right]$$ (4.30)

Le calcul de ces 3 coefficients donne

$$\mathbf{A}_{22}^{k}=\frac{(x_{3}^{k}-x_{1}^{k})^{2}+(y_{3}^{k}-y_... ...{33}^{k}=\frac{(x_{1}^{k}-x_{2}^{k})^{2}+(y_{1}^{k}-y_{2}^{k})^{2}}{4aire_{k}},$$ (4.31)

$$\mathbf{A}_{23}^{k}=\frac{(x_{3}^{k}-x_{1}^{k})(x_{1}^{k}-x_{2}^{k})+(y_{3}^{k}-y_{1}^{k})(y_{1}^{k}-y_{2}^{k})}{4aire_{k}}$$

d'où l'expression de la matrice élémentaire $\mathbf{A}^{k}$:

$$\frac{1}{4aire_{k}}\left[\begin{array}{ccc} (x_{2}^{k}-x_{3}... ...& (y_{3}^{k}-y_{1}^{k})(y_{1}^{k}-y_{2}^{k})\end{array}\right]$$ (4.32)

Pour calculer la matrice de notre système, il suffit donc de calculer les 11 matrices élémentaires correspondants aux $Ne=11$ éléments du maillage en utilisant les relations (4.30) et (4.31), et de reporter les coefficients dans la matrice globale (4.27).

On obtiens ainsi les valeurs des coefficients de la matrice du système:

$$\mathbf{A}=\left[\begin{array}{cccc} 0.5938 & -0.25 & 0.0 & ... ...0 & 40.0 & -16.0\\ 0.0 & 0.0 & -16.0 & 32.0\end{array}\right]$$ (4.33)

qui est bien entendu symétrique.

4.2.4.2 assemblage du second membre

Le calcul du second membre (4.25) procède de la même démarche. On remarque aussi que, dans notre cas, le second membre ne contient que des termes provenant des conditions aux limites. L'intégrale à calculer est exactement la même que pour la matrice $\mathbf{A}$, et on écrit donc le second membre sous la forme:

$$\mathbf{B}_{i}=-\sum_{j=5}^{7}\psi^{e}\left(\sum_{k=1}^{ne}\... ...partial y}\frac{\partial\Phi_{i}}{\partial y}\right)dxdy\right)$$

C'est la somme de coefficients de matrices élémentaires $\mathbf{A}_{pq}^{k}$ (4.26): pour calculer le second membre $\mathbf{B}_{i}$, on doit prendre en compte les matrices élémentaires des éléments du maillage ayant le noeud $i$ comme sommet et un des noeuds $j$ sur le bord $EF$ (i.e. $j=5,6,7$). Ainsi pour le second membre $\mathbf{B}_{1}$ on a:

$$\mathbf{B}_{1}=-\mathbf{A}_{13}^{4}\psi_{7}-\mathbf{A}_{12}^... ...psi_{e}-\mathbf{A}_{12}^{4}\psi_{e}-\mathbf{A}_{32}^{5}\psi_{e}$$

puisque que l'élément $e_{4}$ a pour premier sommet le noeud 1 et pour troisième sommet le noeud 7 (qui sur $EF),$ mais aussi pour second sommet le noeud 6 (qui sur $EF),$ et que l'élément $e_{5}$ a pour troisième sommet le noeud 1 et pour second sommet le noeud 6.

Le second membre complet s'écrit donc:

$$\mathbf{B}=\left[\begin{array}{c} -\mathbf{A}_{13}^{4}\psi_{... ...{12}\psi_{e}\\ -\mathbf{A}_{23}^{12}\psi_{e}\end{array}\right]$$

Le calcul précédent nous a fournit les matrices élémentaires, et on obtiens comme valeurs de $\mathbf{B}$ :

$$\mathbf{B}=\left[\begin{array}{c} 0.0156\\ 10.25\\ 4.0\\ 8.0\end{array}\right]$$ (4.34)

compte tenue de la valeur de $\psi_{e}=1.0$

4.2.4.3 résolution

La résolution du système linéaire avec la matrice (4.33) et le second membre (4.34), nous fournit la valeur de la solution approchée pour les 4 degrés de liberté du système:

$$X=\left[\begin{array}{c} 0.2377\\ 0.5021\\ 0.5010\\ 0.5005\end{array}\right]$$

Compte tenu des conditions aux limites, on obtiens la solution approchée sur tous les noeuds du maillage:

$$\psi_{e}=\left[\begin{array}{ccccccccccc} 0.2377 & 0.5021 & ... ...005 & 1.0 & 1.0 & 1.0 & 0.0 & 0.0 & 0.0 & 0.0\end{array}\right]$$

Cette solution est représentée sur la figure (4.9) sous la forme d'iso-valeurs en couleur. Une ligne $\psi=cste$ corresponds à une couleur fixe, dont la valeur est fournit par la palette de couleurs à droite de la figure. Ces lignes iso-valeurs $\psi=cste$ sont justement les lignes de courant de l'écoulement . Compte tenu de la petitesse du maillage, ces lignes de courant approchées ne sont pas très régulières, mais on retrouve le comportement global de l'écoulement qui est défléchi par l'obstacle.

Figure 4.9: Iso-valeurs de la solution approchée

On peut comparer cette solution à la solution calculée sur un maillage beaucoup plus fin de $3200$ éléments et $1691$ noeuds. Les iso-valeurs de cette solution, qui est très proche de la solution exacte sont tracées sur la figure (4.10).

Figure 4.10: Iso valeurs de la solution sur un maillage très fin

Figure 4.11: Comparaison des profils calculés (x) avec les profils de référence

En comparant les deux figures, on constate que l'allure de l'écoulement est bien la même, par contre la déviation des lignes de courants sur le maillage grossier (4.9) est beaucoup trop important en amont de l'obstacle, comparé à la solution de référence (4.10). Par contre sur l'obstacle, la solution est quasiment identique et correspond à une répartition linéaire. Cela est confirmé par la figure suivante (4.11), où on a tracé les profils de la solution approchée (en traits pointillés) pour deux abscisses $x=-1$ (droit de l'obstacle) et $x=-3$ (amont de l'obstacle), comparés à la solution de référence (en traits pleins).