# Solution de l'équation de Poisson dans un domaine bidimensionnel par la méthode des éléments finis

Mohamed Boun-jad* et Toufik Zebbiche

Laboratoire des Sciences de l'Aéronautique, Institut d'Aéronautique
Université Saâd Dahleb, B.P. 270, Blida 09000, Algeria

(reçu le 15 Avril 2012 – accepté le 29 Septembre 2013)

**Résumé** - Ce travail présente une modeste contribution sur le domaine des éléments finis. On a essayé de résoudre l'équation de Poisson avec les conditions aux limites de type Dirichlet dans un domaine quelconque bidimensionnel simplement connexe. Après une brève théorie sur la résolution de cette équation, on a trouvé les solutions analytiques exactes que pour les sections circulaires, rectangulaires et elliptiques. D'où notre intérêt est orienté vers la recherche des solutions numériques approchées. La méthode utilisée est celle des éléments finis. Deux programmes de calcul numérique ont été réalisés dans ce contexte. Le premier est consacré pour la génération du maillage triangulaire dans les domaines concernés. On a développé ici une technique permettant de générer un maillage de type 'H' dans n'importe quel domaine simplement connexe. Le deuxième programme réalisé est celui pour la résolution numérique de l'équation de Poisson dans un domaine simplement connexe avec les conditions aux limites de Dirichlet. Le programme utilise un seul type d'élément fini est le triangle à trois nœuds. Le problème résolu a un grand intérêt physique et pratique dans pas mal de disciplines, telles que, le transfert de chaleur dans des ailettes de moteur, l'écoulement à travers les conduites quelconques et d'autres problèmes tels que, l'étude de la torsion des poutres de sections non circulaires arbitraires comme la torsion des pales d'hélicoptères et les ailettes dans les refroidisseurs ainsi que les aubes dans les compresseurs. Les résultats présentés sont choisis pour des sections simples pour but de comparaison et pour d'autres formes complexes.

**Abstract** – This work presents a modest contribution to the field of finite elements. We tried to solve the Poisson equation with boundary conditions of Dirichlet in any two-dimensional simply connected domain. After a brief theory on solving this equation, we found the exact analytical solutions for circular, rectangular and elliptical sections. Hence our interest is directed towards the search of approximate numerical solutions. The method used is that of finite elements. Two numerical calculation programs have been made in this context. The first is devoted to the generation of triangular mesh in the areas concerned. We have developed here a technique to generate a mesh type 'H' in any simply connected domain. The second program is realized that for the numerical solution of the Poisson equation in a simply connected domain with Dirichlet boundary conditions. The program uses a single type of finite element is the triangle with three nodes. The problem solved is a great physical and practical interest in a lot of disciplines, such as the heat transfer fins engine, the flow through the pipes and any other problems such as the study of the torsion beams of arbitrary non-circular sections such as the twist of the blades of helicopters and coolers in the vanes and blades in the compressor. The results shown are chosen for simple sections for comparison purposes, and also for other complex shapes.

<footnote>* z_toufik270169@yahoo.fr</footnote>
**Mot clés:** Méthode des éléments finis - Elément triangulaire - Equation de Poisson - Conditions aux limites type Dirichlet - Domaine simplement connexe - Génération de maillage type H - Méthode de résolution directe de Khaletski - Matrice bande - Stockage des matrices sous forme de vecteur - Erreur de calcul - Valeur moyenne - Valeur maximale - Dérivée de la solution.

# 1. INTRODUCTION

L'évolution actuelle de la technologie amène le chercheur à réaliser des projets de plus en plus complexes, coûteux, et soumis à des contraintes de sécurité de plus en plus sévères. Nous pensons aux projets spatiaux, aéronautiques et nucléaires dans lesquels la sécurité est vitale. Pour dominer ces projets et autres, le chercheur a besoin de modèles qui lui permettent de simuler le comportement de systèmes physiques complexes.

Les sciences de l'ingénieur (mécanique des solides et des fluides, thermique, ...) permettent de décrire le comportement de systèmes physiques, grâce à des équations aux dérivées partielles [1-7]. La méthode des éléments finis est l'une des méthodes les plus utilisées aujourd'hui pour résoudre effectivement ces équations. Elle nécessite l'utilisation intensive de l'ordinateur. C'est une méthode très générale qui s'applique à la majorité des problèmes rencontrés dans la pratique; problèmes stationnaires ou non stationnaires, linéaires ou non linéaires, définis dans un domaine géométrique quelconque à une, deux ou trois dimensions [1], [3], [8], [14].

La méthode des éléments finis consiste à utiliser une approximation simple des variables inconnues pour transformer les équations aux dérivées partielles en équations algébriques [5], [8-12]. Elle fait appel aux trois domaines suivants:

* Sciences de l'ingénieur pour construire les équations aux dérivées partielles [8-9].

* Méthodes numériques pour construire et résoudre les équations algébriques [8], [13-21].

* Programmation et informatique pour exécuter efficacement les calculs sur l'ordinateur.

La méthode des éléments finis est maintenant très répandue dans les industries, en particulier en construction aéronautique, aérospatiale, navale et nucléaire et en général en mécanique des fluides. Elle nécessite l'écriture des programmes informatiques professionnels répondant à la résolution de ces problèmes. Pour que la méthode des éléments finis soit efficace dans les applications industrielles, il faut utiliser des programmes d'assistance à la préparation des données et à l'interprétation des résultats [12], [16], [22]. La compréhension de la méthode exige en effet des connaissances dans des domaines variés [1-3], [8], [10], [23].

* Compréhension du problème physique étudié et connaissance intuitive de la nature de la solution cherchée.

* Approximation des inconnues par sous-domaines et construction de fonctions d'interpolations.

* Construction des équations du système étudié sous forme variationnelle, soit à partir de méthodes énergétiques, soit à partir d'équations aux dérivées partielles.

* Technique d'organisation matricielle des données.

* Méthodes numériques d'intégration, de résolution de systèmes d'équations algébriques et différentielles, linéaires et non linéaires.
La présente étude consiste alors à la résolution numérique de l'équation de Poisson dans un domaine simplement connexe quelconque et complexe avec les conditions aux limites de type Dirichlet. En effet, la complexité de ces domaines induit que la solution analytique (exacte) du problème n'existe pas, et évidemment notre intérêt s'oriente vers la recherche des solutions numériques approchées. La méthode des éléments finis sera introduite et appliquée avec succès, puisqu'elle s'adapte sur n'importe quelle frontière du domaine [8], [9], [11], [23-27].

La discrétisation des équations, gérant le phénomène physique seront donc fondées sur l'utilisation des éléments finis bi dimensionnel (2D). L'élément utilisé est celui le triangle à trois nœuds. C'est l'élément le plus simple. Il donne une bonne présentation de la frontière de la géométrie complexe [2], [8], [11], [13], [16].

Cependant, deux difficultés principales se rencontrent, si la méthode des éléments finis est utilisée.

* La première difficulté majeure, c'est la manière optimale de génération de maillage, ainsi que l'adaptation des équations au modèle mathématique de la résolution.

* L'autre difficulté est le choix de la méthode numérique adéquate à la résolution du système d'équations algébriques, vu le nombre important de degré de liberté alloués. Pour cela, il faut subdiviser le domaine en petits sous-domaines de géométrie simple connue, telle que la géométrie triangulaire par l'intermédiaire des fonctions d'interpolations. Un nombre important d'éléments finis sera établi et appliqué pour le développement d'un programme numérique, afin d'avoir une convergence vers la solution exacte.

Vu la dimension importante de la matrice du système, le stockage dans la mémoire de l'ordinateur sera réalisé en utilisant la technique des matrices bandes afin d'éviter le stockage des cases nulles. Pour avoir un nombre important de zéros, il faut faire une bonne numérotation des nœuds de maillage. De même, la résolution du système d'équations sera faite de manière directe en utilisant l'algorithme de Khalatski, qui sera modifié spécialement pour des matrices bandes et symétriques [2], [8], [28-32].

La solution du problème traité représente la distribution du champ de solution dans le domaine. Puisque l'équation de Poisson peut représenter pas mal des phénomènes physiques, il est parfois très intéressant de calculer la distribution de la dérivée de la solution à travers le domaine, ainsi que l'intégrale de la solution.

## 2. ANALYSE DE L'EQUATION DE POISSON PAR LA METHODE DES ELEMENTS FINIS

L'analyse des éléments finis est assez compliquée, par le fait que les problèmes à deux dimensions sont décrits par des équations aux dérivées partielles. Le contour $Γ$ du domaine $Ω$ à deux dimensions est en général une courbe. Cependant, les éléments finis ont des formes géométriques simples, que l'on peut utiliser comme une approximation donnée au domaine à deux dimensions. D'autre part, pour les problèmes à 2D, on ne cherche pas uniquement l'approximation de la solution pour les équations aux dérivées partielles données, mais on approxime aussi le domaine donné par un maillage convenable en éléments finis [9].

Par conséquent, on obtient une bonne approximation de l'erreur (due à l'approximation de la solution) comme l'erreur de la discrétisation (due à l'approximation du domaine) dans l'analyse des problèmes des éléments finis à deux
dimensions. Le maillage des éléments finis concerne les éléments simples 2D, tels que, les triangles, rectangles et quadrilatères qu'ils soient connectés entre eux par des points nodaux sur le contour de l'élément. La possibilité de représenter la géométrie du domaine par la collection des éléments finis rend la méthode valable et pratique et un outil pour la résolution des problèmes aux valeurs aux limites touchent une variété des domaines de l'ingénieur [2], [8-12].

Considérons le problème de la recherche de la solution u de l'équation différentielle de second ordre de Poisson avec un deuxième membre suivante [6], [8-9]:

$$ \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} = f \quad ; \quad u = 0 \text{ sur } \Gamma \qquad (1) $$

La valeur de f est considérée égale à l'unité. Le plan apporté à ce domaine est celui (x0y). La forme variationnelle de l'équation (1) est donnée par:

$$ 0 = \int_{\Omega_e} \left[ \frac{\partial v}{\partial x} \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} \frac{\partial u}{\partial y} - v f \right] d x \times dy \qquad (2) $$

L'équation (2) est la base de la forme variationnelle du modèle en éléments finis de l'équation (1).

## 2.1 Formulation en éléments finis

La forme variationnelle de l'équation (1) indique que l'approximation de u choisie doit être au moins bilinéaire en u et en y. Supposons que u est approximé par l'expression suivante:

$$ u = \sum_{j=1}^{j=n} u_j \times \Psi_j \qquad (3) $$

Où, $u_j$ sont les valeurs de u aux points $(x_j, y_j)$ et $\Psi_j$ sont les fonctions d'interpolations linéaires avec la propriété:

$$ \Psi_i(x_j, y_j) = \delta_{ij} \qquad (4) $$

La forme spécifiée de $\Psi_j$ sera dérivée pour les éléments triangulaires linéaires. Substituons l'équation (3) dans la forme variationnelle (2), on obtient:

$$ \sum_{j=1}^{j=n} K_{ij}^{(e)} \times u_j^{(e)} = F_i^{(e)} \qquad (5) $$

Telle que:

$$ K_{ij}^{(e)} = \int_{\Omega_e} \left[ \frac{\partial \Psi_i}{\partial x} \frac{\partial \Psi_j}{\partial x} + \frac{\partial \Psi_i}{\partial y} \frac{\partial \Psi_j}{\partial y} \right] d x \times dy \qquad (6) $$

$$ F_i^{(e)} = \int_{\Omega_e} f \times \Psi_i \times d x \times dy \quad i = 1, 2, \dots, n \qquad (7) $$

L'équation (8) représente le modèle mathématique en éléments finis de l'équation (1). Pour les problèmes à deux dimensions, la relation entre le nombre de nœuds (qui est égal au nombre de termes dans le polynôme de l'approximation de u) et le degré du
polynôme n'est pas unique. Par exemple, le polynôme, contient trois termes linéairement indépendants, et il est linéaire pour x et y :

$$u(x,y) = c_1 + c_2 \times x + c_3 \times y \quad (8)$$

L'approximation (8) exige un élément avec trois nœuds (avec un seul degré de liberté pour chaque nœud). L'élément bidimensionnel avec trois nœuds est un triangle avec nœuds aux sommets du triangle. Ici nous dérivons les fonctions d'interpolations linéaires pour l'élément triangulaire.

## 2.2 Fonctions d'interpolations d'un élément triangulaire à trois nœuds

Considérons l'approximation linéaire (8). L'ensemble $(1, x, y)$ est linéairement indépendant et complet. Nous devons réécrire l'approximation (8), pour qu'elle satisfasse les conditions suivantes:

$$u(x_i, y_i) = u_i \quad i=1,2,3 \quad (9)$$

où $(x_i, y_i)$ $(i = 1,2,3)$ sont les coordonnées du trois sommets du triangle. Autrement, on détermine les trois constantes $c_i$ dans l'équation (8) en fonctions de $u_i$ à partir de l'équation (9). Sous la forme matricielle, on obtient:

$$\begin{cases} u_1 \\ u_2 \\ u_3 \end{cases} = \begin{bmatrix} 1 & x_1 & y_1 \\ 1 & x_2 & y_2 \\ 1 & x_3 & y_3 \end{bmatrix} \times \begin{cases} c_1 \\ c_2 \\ c_3 \end{cases} \quad (10)$$

Notons que les nœuds sont numérotés dans le sens inverse des aiguilles d'une montre. Résolvant le système d'équations (10) pour $c_i$ $(i = 1,2,3)$, ensuite substituant les résultats dans l'équation (8), on obtient:

$$u(x,y) = \sum_{i=1}^{i=3} u_i \times \Psi_i^{(e)} \quad (11)$$

où $\Psi_i^{(e)}$ sont les fonctions d'interpolations linéaires pour l'élément triangulaire.

$$\Psi_i^{(e)} = \frac{1}{2A_e} \times (\alpha_i + \beta_i x + \gamma_i y) \quad i=1,2,3 \quad (12)$$

avec $A_e$, l'aire du triangle. Elle est donnée par:

$$2A_e = \begin{vmatrix} 1 & x_1 & x_1 \\ 1 & x_2 & y_2 \\ 1 & x_3 & y_3 \end{vmatrix} \quad (13)$$

Et $\alpha_i, \beta_i, \gamma_i$ sont des constantes données par:

$$\left. \begin{array}{l} \alpha_i = x_j y_k - x_k y_j \\ \beta_i = y_j - y_k \\ \gamma_i = x_k - x_j \end{array} \right\} \quad i \neq j \neq k \quad (14)$$

Et les indices $i, j, k$ permutent dans l'ordre naturel.
## 2.3 Calcul des matrices de rigidités [2], [8-10], [12], [14], [33]

Le calcul des matrices de l'élément $[K^{(e)}]$ et le vecteur force $\{F^{(e)}\}$ dans l'équation (7) par les méthodes conventionnelles (par intégration exacte) est, en général, pas facile [34]. Nous récrivons $[K^{(e)}]$ dans l'équation (7) comme la somme de deux matrices de base $[S^{11}]$ et $[S^{22}]$ données par:

$$[K^{(e)}] = [S^{11}] + [S^{22}] \qquad (15)$$

où

$$S_{ij}^{11} = \int_{\Omega_e} \frac{\partial \Psi_i}{\partial x} \frac{\partial \Psi_j}{\partial x} dx dy \qquad (16)$$

$$S_{ij}^{22} = \int_{\Omega_e} \frac{\partial \Psi_i}{\partial y} \frac{\partial \Psi_j}{\partial y} dx dy \qquad (17)$$

On peut écrire aussi:

$$f_i^{(e)} = \int_{\Omega_e} f \times \psi_i \, dx dy \qquad (18)$$

## 2.4 Matrice de rigidité d'un élément triangulaire à trois nœuds [8, 9], [34]

Pour un triangle, les formules des intégrales exactes suivantes sont nécessaires pour l'évaluation des éléments des matrices de l'élément [8, 9], [26], [34].

$$I_{mn} = \int_{\Delta} x^m \times y^n dx dy \qquad (19)$$

Alors, nous avons:

$$I_{00} = A \qquad (20)$$

$$I_{10} = A \times x_G \quad ; \quad x_G = \frac{1}{3} \sum_{i=1}^{i=3} x_i \qquad (21)$$

$$I_{01} = A \times y_G \quad ; \quad y_G = \frac{1}{3} \sum_{i=1}^{i=3} y_i \qquad (22)$$

$$I_{11} = \frac{A}{12} \left[ \sum_{i=1}^{i=3} x_i y_i + 9 x_G y_G \right] \qquad (23)$$

$$I_{20} = \frac{A}{12} \left[ \sum_{i=1}^{i=3} x_i^2 + 9 x_G^2 \right] \qquad (24)$$

$$I_{02} = \frac{A}{12} \left[ \sum_{i=1}^{i=3} y_i^2 + 9 y_G^2 \right] \qquad (25)$$

En utilisant les fonctions d'interpolations linéaires (12) dans les équations (6) et (17), on obtient:

$$S_{ij}^{11} = \frac{1}{4A} \beta_i \beta_j \quad ; \quad S_{ij}^{22} = \frac{1}{4A} \gamma_i \gamma_j \qquad (26)$$
Connaissant que: $\alpha_i + \beta_i x_G + \gamma_i y_G = 2 A/3$ [qui peut être obtenu qu'à partir des équations (14), (20), (21) et (22)], on a alors:

$$f_i^{(e)} = \frac{f}{2}(\alpha_i + \beta_i x_G + \gamma_i y_G) = \frac{f \times A}{3} \quad (27)$$

Une fois les coordonnées des nœuds de l'élément connues, on peut calculer facilement $\alpha_i$, $\beta_i$ et $\gamma_i$ à partir de l'équation (14) et en les substituant ensuite dans l'équation (26) pour obtenir les matrices de l'élément.

## 2.5 Assemblage des matrices d'éléments

L'assemblage est l'opération qui consiste à construire la matrice globale $[K]$ et le vecteur global $\{F\}$ du domaine complet à partir des matrices élémentaires $[K^{(e)}]$ et des vecteurs élémentaires $\{F^{(e)}\}$ des éléments finis. Pour représenter l'assemblage des matrices d'éléments finis, nous illustrons la procédure en considérant que le maillage est construit par deux éléments triangulaires.

Soient $[K_{ij}^{(1)}]$, $[K_{ij}^{(2)}]$ $(i, j = 1, 2, 3)$, les coefficients des matrices de rigidité correspondant à ces deux éléments $\{F_i^{(1)}\}$ et $\{F_i^{(2)}\}$ $(i = 1, 2, 3)$ vecteurs forces des deux éléments respectivement.

Le maillage d'éléments finis est présenté dans la figure 1. A partir de ce maillage, nous notons la correspondance entre les valeurs globales et les valeurs nodales des éléments.

$$u_1 = u_1^{(1)}; u_2 = u_2^{(1)}; u_3 = u_3^{(2)}; u_4 = u_4^{(1)} = u_3^{(2)} \quad (28)$$

<figure><img src="image_1.png" /><figcaption>Fig. 1: Assemblage des deux éléments triangulaire à 3 nœuds</figcaption></figure>

Pour l'élément (1)

$$[K^{(1)}] = \begin{bmatrix} 1 & 4 & 2 \\ K_{11}^{(1)} & K_{12}^{(1)} & K_{13}^{(1)} \\ K_{21}^{(1)} & K_{22}^{(1)} & K_{23}^{(1)} \\ K_{31}^{(1)} & K_{32}^{(1)} & K_{33}^{(1)} \end{bmatrix}$$
$$ \{ F^{(1)} \} = \begin{array}{c|c} 1 & F_1^{(1)} \\ 4 & F_2^{(1)} \\ 2 & F_3^{(1)} \end{array} $$

Pour l'élément (2)

$$ [K^{(2)}] = \begin{array}{c|c|c} 1 & K_{11}^{(2)} & K_{12}^{(2)} \\ 3 & K_{21}^{(2)} & K_{22}^{(2)} \\ 4 & K_{31}^{(2)} & K_{32}^{(2)} \end{array} \begin{array}{c|c} 1 & 3 & 4 \\ K_{11}^{(2)} & K_{12}^{(2)} & K_{13}^{(2)} \\ K_{21}^{(2)} & K_{22}^{(2)} & K_{23}^{(2)} \\ K_{31}^{(2)} & K_{32}^{(2)} & K_{33}^{(2)} \end{array} $$

$$ \{ F^{(2)} \} = \begin{array}{c|c} 1 & F_1^{(2)} \\ 3 & F_2^{(2)} \\ 4 & F_3^{(2)} \end{array} $$

La matrice $[K]$ et le vecteur force $\{F\}$ du domaine complet sont obtenus respectivement comme assemblage des matrices de rigidité et vecteurs forces des deux éléments constituants le domaine complet.

$$ [K] = \begin{bmatrix} K_{11}^{(1)} + K_{11}^{(2)} & K_{13}^{(1)} & K_{12}^{(2)} & K_{12}^{(1)} + K_{13}^{(2)} \\ K_{41}^{(1)} & K_{33}^{(1)} & 0 & K_{32}^{(1)} \\ K_{21}^{(2)} & 0 & K_{22}^{(2)} & K_{23}^{(2)} \\ K_{12}^{(1)} + K_{31}^{(2)} & K_{13}^{(1)} & K_{32}^{(2)} & K_{22}^{(1)} + K_{33}^{(2)} \end{bmatrix} \quad (29) $$

$$ \{F\} = \begin{cases} F_1^{(1)} + F_1^{(2)} \\ F_3^{(1)} \\ F_2^{(2)} \\ F_2^{(1)} + F_3^{(2)} \end{cases} \quad (30) $$

La matrice (29) représente la matrice de rigidité globale du domaine complet. Elle est de rang $4 \times 4$, vu que le nombre de nœuds dans le domaine est égal à 4. On note ici que chaque nœud contient un seul degré de liberté (une inconnue). L'équation (30) représente le vecteur force du domaine complet. Il est de rang 4 pour cet exemple.

### 3. GENERATION DE MAILLAGE

La mise en œuvre effective de la méthode des éléments finis décrite dans le paragraphe précédent nécessite l'utilisation des méthodes numériques puissantes pour résoudre le système algébrique, qui en résulte.
Dans ce paragraphe, on traite le problème de la génération de maillage pour les domaines simplement connexes bidimensionnels. Il existe plusieurs méthodes et techniques de génération de maillage que se soit triangulaires ou quadrilatères [1], [23-27].

Dans notre étude, on va développer une méthode algébrique qui donne un maillage structuré, dont les cellules sont de forme quadrilatères qu'on peut facilement les transformer en des cellules triangulaires. Le maillage obtenu est de la forme 'H'.

## 3.1 Procédure de génération de maillage

Nous sommes intéressés aux domaines de type simplement connexes (2D). La forme générale de ce domaine peut être présentée par la figure 2.

Pour faire la génération de maillage à l'intérieur de ce domaine, et parmi les plusieurs méthodes et techniques employées et suivant la méthode de résolution du système algébrique, ainsi que la manière de considération la forme de la matrice de rigidité, il est très intéressant d'appliquer l'idée de génération de maillage dans les régions rectangulaires pour ce type de problème.

<figure><img src="image_2.png" /><figcaption>Fig. 2: Frontière d'un domaine simplement connexe</figcaption></figure>

On obtient par conséquent quatre cotés qui sont respectivement AD, AB, BC et CD. Si on fait des transformations continues pour chaque coté, on peut rendre ces lignes sous forme de droites comme le montre la figure adjacente. Alors le coté AB du domaine joue le rôle de la partie basse du rectangle. Le coté CD est celui de la partie haute du rectangle.

Les cotés AD et BC sont ceux de la partie gauche et droite du rectangle. Pour le rectangle lui même, nous considérons la subdivision par $N_x$ et $N_y$ nœuds sur l'axe horizontal des x et vertical des y respectivement. Par analogie, les cotés AB et CD du domaine seront subdivisés en $N_x$ nœuds et les cotés AD et BC seront subdivisés en $N_y$ nœuds.

### 3.1.1 Génération des nœuds sur le contour du domaine

Avant de commencer la génération des nœuds internes, il faut commencer par ceux de la frontière de chaque coté. Plusieurs procédures peuvent être considérées et cela suivant la disposition de la courbure de chaque coté. Dans le cas général, il n'est pas facile de procéder à une méthode universelle qui sera applicable pour n'importe quelle courbe, mais on essaye d'adapter pour chaque forme complexe, une procédure convenable, sans oublier de dire que le choix des points A, B, C et D influe sur la forme des quatre cotés du domaine.
Parmi les critères de choix, des positions des points A, B, C et D, est que les fonctions sur les cotés AB, BC, CD et AD doivent être univoque en x ou y, c'est-à-dire, pour chaque valeur de x, il y a une image sur les cotés AB et CD. En ce qui concerne les fonctions qui passent par les points AD et BC et qu'elles doivent être aussi univoques, mais en y.

Dans ce qui suit, on suppose que les fonctions des cotés choisies du domaine sont univoques. Considérons les cotés AB et CD, les positions des points A, B, C et D sont données. On regarde ensuite est-ce qu'on divise le segment $[x_A, x_B]$ en $N_x$ nœuds, ensuite déterminer les $y_i$ ($i = 1, 2, \dots, N_x$) correspondantes ou, de diviser le segment $[y_A, y_B]$ en $N_x$ nœuds, ensuite déterminer les $x_i$ ($i = 1, 2, \dots, N_x$). Pour prendre décision, on essaye de calculer les valeurs suivantes:

$$x_{dis} = |x_B - x_A| \quad ; \quad y_{dis} = |y_B - y_A| \qquad (31)$$

On calcule la valeur maximale entre $x_{dis}$ et $y_{dis}$ par:

$$d_{max} = \text{Max}(x_{dis} \cdot y_{dis}) \qquad (32)$$

Si $d_{max} = x_{dis}$, on propose les valeurs des $x_i$ ($i = 1, 2, \dots, N_x$), ensuite déterminer les $y_i$ en utilisant la fonction qui passe par les points A et B du coté AB. Dans certains cas, pour prendre la décision, il est mieux de voir encore si la fonction est donnée par $y = f_{AB}(x)$ ou $x = g_{AB}(y)$ pour ne pas calculer la fonction inverse. On suppose qu'on a fait le choix des $x_i$, ensuite calculer les $y_i$, alors, on divise le segment $[x_A, x_B]$ en $N_x$ nœuds dont les points A et B sont inclus. Pour simplifier, on choisit un pas constant de subdivision et si le coté est courbe au voisinage du point A ou B dans n'importe quelle région, il est recommandé d'utiliser une procédure de condensation des nœuds, dont on propose sa procédure dans le prochain paragraphe. On obtient, par conséquent ($N_x - 1$) sous-intervalles. La longueur de chaque intervalle est donnée par:

$$\Delta x_{AB} = \frac{|x_B - x_A|}{(N_x - 1)} \qquad (33)$$

Les abscisses des points de ce côté sont données par:

$$x_{AB}(i) = x_A + (i - 1)\Delta x_{AB} \quad i = 1, 2, \dots, N_x \qquad (34)$$

Les ordonnées des points sont données par:

$$y_{AB}(i) = f_{AB}(i) \qquad (35)$$

Où $f_{AB}(x)$, est la fonction du côté AB supposée donner en coordonnées cartésiennes.

Concernant le côté et par analogie avec le côté AB, on aura le suivant en changeant uniquement A par D, B par C et $f_{AB}(x)$ et $f_{DC}(x)$. Alors,

$$\Delta x_{CD} = \frac{|x_C - x_D|}{(N_x - 1)} \qquad (36)$$
Les abscisses et les ordonnées des points de ce côté sont données par:

$$x_{\text{DC}}(i) = x_{\text{D}} + (i - 1)\Delta x_{\text{CD}} \quad i = 1, 2, \dots, N_x \quad (37)$$

$$y_{\text{DC}}(i) = f_{\text{DC}} \times (x_{\text{DC}}(i)) \quad i = 1, 2, \dots, N_x \quad (38)$$

Où $f_{\text{DC}}(x)$, est la fonction du côté DC.

La génération des nœuds sur les contours AD et BC se fait de la même démarche que celle des cotés AB et DC, ici suivant la disposition de ces cotés. Si $x_A \neq x_D$ et $x_D \neq x_C$ en plus les deux fonctions sont univoques, on peut suivre exactement la démarche précédente pour les cotés AB et DC. On aura:

$$\Delta x_{\text{AD}} = \frac{|x_A - x_D|}{(N_y - 1)} \quad (39)$$

$$x_{\text{AD}}(i) = x_{\text{D}} + (j - 1)\Delta x_{\text{AD}} \quad j = 1, 2, \dots, N_y \quad (40)$$

$$y_{\text{AD}}(i) = f_{\text{AD}} \times (x_{\text{AD}}(j)) \quad j = 1, 2, \dots, N_y \quad (41)$$

Mais si $x_A = x_D$, l'algorithme (40) et (41) ne fonctionne plus et il faut considérer de la manière suivante. On divise le segment $[y_A, y_D]$ en $N_y$ nœuds dont les points A et B sont inclus, on aura:

$$\Delta y_{\text{AD}} = \frac{|y_D - y_A|}{(N_y - 1)} \quad (42)$$

Les coordonnées des points de ce coté sont données par:

$$y_{\text{AD}}(j) = y_A + (j - 1)\Delta y_{\text{AD}} \quad j = 1, 2, \dots, N_y \quad (43)$$

$$x_{\text{AD}}(j) = f_{\text{AD}}^{-1} \times (y_{\text{AD}}(j)) \quad j = 1, 2, \dots, N_y \quad (44)$$

Où $f_{\text{AD}}^{-1}(y)$, est la fonction inverse du côté AD. Pour déterminer la racine de la fonction inverse et si $f_{\text{AD}}(y)$ est implicite, on utilise l'algorithme de dichotomie par exemple.

Pour le côté BC, on aura:

$$\Delta y_{\text{BC}} = \frac{|x_B - x_C|}{(N_y - 1)} \quad (45)$$

$$x_{\text{BC}}(j) = x_C + (j - 1)\Delta y_{\text{BC}} \quad j = 1, 2, \dots, N_y \quad (46)$$

$$y_{\text{BC}}(j) = f_{\text{BC}} \times (x_{\text{BC}}(j)) \quad j = 1, 2, \dots, N_y \quad (47)$$

Mais si $x_B = x_C$, l'algorithme (46) et (47) ne fonctionne plus, et il faut une autre procédure. On divise le segment $[y_B, y_C]$ en Finalement, après qu'on génère les nœuds sur la frontière on obtient le schéma de la figure 3.

On note ici qu'entre chaque deux nœuds successifs sur un coté, on assimile la courbure de la frontière à un segment de droite. L'exemple présenté sur cette figure est pour $N_x = 10$ et $N_y = 6$, dont les nœuds B et C sont inclus. On obtient:
$$ \Delta y_{BC} = \frac{y_C - y_B}{(N_y - 1)} \qquad (48) $$

<figure><img src="image_3.png" /><figcaption>Fig. 3: Génération des nœuds sur la frontière du domaine</figcaption></figure>

$$ \begin{gathered} y_{BC}(j) = y_B + (j-1)\Delta y_{BC} \\ x_{BC}(j) = f_{BC}^{-1}(y_{BC}(j)) \end{gathered} \qquad (49) $$

Où, $f_{BC}^{-1}(y)$, représente la fonction inverse du côté BC.

### 3.1.2 Fonction de condensation 1D

La fonction de condensation unidimensionnelle (stretching function) est généralement utilisée pour distribuer les points, le long d'une frontière particulière ou une région d'un domaine nécessitant une résolution précise. Par exemple, pour étudier un écoulement visqueux sur un profil d'aile (les aérodynamiciens s'intéressent à ce problème) nécessite l'introduction de cette fonction ou équivalente pour renforcer la précision de la solution au niveau de la surface du profil d'aile et par conséquent, les points de la frontière se condensent volontairement vers une région voulue. Le gradient de la solution est fort au niveau de cette région [1, 23].

Mathématiquement, il est intéressant d'introduire une variable indépendante dans la fonction pour la normalisation. Pour la fonction de condensation unidimensionnelle appliquée sur le côté EA de la figure 4, la variable indépendante normalisée est donnée par:

$$ \eta^* = \frac{\eta - \eta_A}{\eta_E - \eta_A} \qquad (50) $$

Avec, $0 \le \eta^* \le 1$ et $\eta_A \le \eta \le \eta_E$

Où, $\eta$ peut représenter x ou y.

On peut même donner la distribution sur l'intervalle $[0, 1]$ par $\eta^*$ avec des sous intervalles égaux. La fonction de condensation développée par Robert en 1971, ensuite modifier par Eisemann en 1979 est donnée par [1], [23]:

$$ s = P \times \eta^* + (1 - P) \times \left( 1 - \frac{\tanh[Q(1 - \eta^*)]}{\tanh[Q]} \right) \qquad (51) $$

Où P et Q sont deux paramètres de contrôle de points de maillage. Si la valeur de P est proche de l'unité, le deuxième membre de la formule disparait rapidement. Une
fois ‘s’ obtenu, il est demandé de spécifier la distribution de x et y. Par exemple, en définissant,

$$ \frac{x - x_A}{x_A - x_E} = f(s) \quad ; \quad \frac{y - y_A}{y_A - y_E} = g(s) \qquad (52) $$

Ces formules génèrent $x(s)$ et $y(s)$ directement. Un simple choix est que $f(s) = g(s) = s$, tels que les relations (51) donnent:

$$ x = x_A + s \times (x_A - x_E) \quad , \quad y = y_A + s \times (y_A - y_E) \quad (53) $$

Des distributions typiques des points sur le segment EA en utilisent la relation (52) pour différentes valeurs de P et Q, sont présentées sur la figure 4 suivante:

<figure><img src="image_4.png" /><figcaption>Fig. 4: Distribution des nœuds selon la relation (51)</figcaption></figure>

Pour des valeurs de $P > 1.0$, il est possible de condenser les nœuds vers le point A. De même, pour condenser vers les nœuds vers le point E, il suffit de prendre des valeurs de P proche de zéro. On peut remarquer que si on choisit la fonction $f(s) = 1 - s$, on trouve la configuration inverse, c'est-à-dire, pour des valeurs de $P > 1.0$ les points se condensent vers l'extrémité E au lieu de l'extrémité A.

### 3.1.3 Génération des nœuds internes

Après qu'on ait déterminé les positions des nœuds de la frontière du domaine, on procède maintenant à la présentation de la méthode développée pour la génération des nœuds internes. Le maillage obtenu est de type ‘H’, c'est-à-dire que la lettre ‘H’ figure dans le maillage.

On obtient des petits cellules ou domaines de formes quadrilatères. Chaque nœud est caractérisé par deux valeurs $(x, y)$ qui permet de déterminer la position exacte de ce point. La procédure qu'on présente est divisée en deux parties. Elle consiste en premier lieu, de déterminer les abscisses $x_i$ de tous les points en ignorant les ordonnées, et cela par une interpolation entre les cotés gauche et droite. On obtient une infinité de solution de position, puisque les $y_i$ sont absentes. Une fois, les abscisses sont obtenus, on procède immédiatement à la détermination les ordonnées $y_i$ de ces nœuds par une interpolation entre les cotés bas et haut du domaine, (Fig. 6).

Comme le nombre de nœuds sur les cotés gauche et droite sont égaux à $N_y$, alors, à chaque nœud d'un coté sera relié avec un nœud de l'autre coté par une courbe qui contient $N_x$ nœuds, dont les extrémités sont inclus, (Fig. 5).

Pour déterminer les abscisses des nœuds sur la courbe FG, on s'intéresse uniquement aux abscisses des nœuds des extrémités gauche et droite qui sont connues.
La projection de la courbe reliant les points F et G sur l'axe horizontal donne un segment de droite.

<figure><img src="image_5.png" /><figcaption>Fig. 5: Illustration de la procédure de détermination des abscisses</figcaption></figure>

On divise le segment FG en Nₓ points dont ceux des extrémités sont inclues. On obtient par conséquent, (Nₓ - 1) intervalles égaux. L'abscisse de chaque point sur cet intervalle et par suite sur la courbe FG est obtenue par:

$$x_i = x_F + (i - 1) \times \Delta x \quad i=1,2,3,\dots,N_x \quad (54)$$

$$\text{Avec, } \Delta x = \frac{x_G - x_F}{(N_x - 1)} \quad (55)$$

En faisant varier le point F sur le coté gauche en allant du point A jusqu'au point D et le point G du coté droit en allant du point B jusqu'au point C, on obtient les abscisses des nœuds des autres courbes et cela en faisant l'extension de la relation (54) par:

$$\left. \begin{aligned} \Delta x(j) &= \frac{x_{\text{côté2}}(j) - x_{\text{côté4}}(j)}{(N_x - 1)} & \quad j=1,2,3,\dots,N_y \\ x(i,j) &= x_{\text{côté4}}(j) + (i-1) \times \Delta x(j) & \quad i=1,2,3,\dots,N_x \end{aligned} \right\} \quad (56)$$

La deuxième étape consiste à déterminer les ordonnées des nœuds internes en utilisant par conséquent l'interpolation linéaire entre le coté 3 du haut et le coté 1 du bas par:

<figure><img src="image_6.png" /><figcaption>Fig. 6: Illustration de la procédure de détermination des ordonnées</figcaption></figure>

Comme chaque courbe contient Nᵧ nœuds, alors, chaque segment contient Nᵧ nœuds. En divisant ce dernier en (Nᵧ - 1) intervalles égaux, l'ordonnée de chaque nœud est donnée par:

$$y_i = y_R + (j - 1) \times \Delta y \quad j=1,2,3,\dots,N_y \quad (57)$$
$$ \text{Avec } \Delta y = \frac{y_H - y_R}{(N_y - 1)} \qquad (58) $$

En faisant varier maintenant le point R sur le coté 1 du bas et le point H sur le coté 3 du haut, on obtient une autre courbe dont les ordonnées de ces nœuds sont données par l'extension de la relation (57).

$$ \left. \begin{aligned} \Delta y(i) &= \frac{y_{\text{côté3}}(i) - y_{\text{côté1}}(i)}{(N_x - 1)} & \quad & i=1,2,3,\dots,N_x \\ y(i,j) &= y_{\text{côté1}}(i) + (j-1) \times \Delta y(i) & \quad & j=1,2,3,\dots,N_y \end{aligned} \right\} \qquad (59) $$

Finalement, les relations (56) et (59) déterminent rapidement les positions des nœuds internes du maillage. Une fois, on s'intéresse à un raffinement du maillage dans une région qui contient un fort gradient. On peut utiliser la procédure de condensation des nœuds présentée dans le paragraphe précédent 3.1.2.

## 3.2 Influence de la numérotation des nœuds sur la forme de la matrice de rigidité

La numérotation des nœuds joue un rôle très important pour la mémorisation des matrices, surtout, si le rang de la matrice devient de plus en plus élevé [17-19].

Pour une matrice carrée quelconque d'ordre $N \times N$, il faut stocker $N^2$ cases, ce qui devient impossible pour les valeurs de N élevée. Le besoin en mémoire de l'ordinateur pour notre problème est important, ce qui nous oblige à penser à une technique pour résoudre un système d'équations d'ordre très élevé afin de faire la résolution du problème physique avec une bonne précision voulue.

Si on regarde le domaine du problème physique, on peut rendre la matrice de rigidité correspondante pleines de zéros qui convergent vers une région connue et les éléments non nuls se concentrent vers la diagonale principale, appelées matrices bandes [17-19], [21], [28-32]. En plus la formulation mathématique en éléments finis du problème physique donne une matrice symétrique.

Ces deux propriétés (bande + symétrie) donnent un avantage important et répondent aussi bien à la question de mémorisation des éléments de la matrice de rigidité. On peut diviser les éléments de la matrice [K] en deux parties: - Eléments nuls, - Eléments non nuls.

Le résultat de l'application des opérations mathématiques entre deux nombres réels ou entiers, dont l'un est nul, est connu au départ sans refaire un calcul. Donc, on n'a pas besoin de donner à l'ordinateur de refaire ce calcul, puisque le zéro est connu et il n'est pas nécessaire de réserver toute une case pour stocker un zéro. Comme le nombre des zéros peut être important, suivant la numérotation, il reste d'indiquer uniquement leurs positions exactes dans la matrice de rigidité [K].

Alors, il est très important de rendre la matrice de rigidité pleine de zéro afin de minimiser la mémorisation. Ce souhait ne sera réalisable qu'après avoir fait une très bonne numérotation des nœuds dont le but de serrer au maximum les zéros vers le coin supérieur de droite de la matrice. Ce type de matrice est appelée matrice bande. Elle est caractérisée par la valeur de la demi-bande. Le problème se pose donc au niveau de la valeur de cette bande. Le nombre des zéros de la matrice augmente si la valeur de la demi-bande diminue [17-19], [21], [29-32].

La partie hachurée contient les éléments non nuls (partie bande). Il est à noter qu'on peut trouver dans la partie bande des éléments nuls.
Mais si on veut travailler avec les matrices bandes, on ne peut pas éliminer ces éléments de la mémorisation puisque le travail avec ce type de matrice, ainsi que la décomposition en deux matrices triangulaires et aussi bandes.

<figure><img src="image_7.png" /><figcaption>Fig. 7: Illustration de la matrice bande</figcaption></figure>

On note encore que si on calcule le nombre des zéros qui se trouvent dans la bande, on peut trouver pour certaines numérotations un nombre qui dépasse la moitié. Pour illustrer ce qu'on a dit, on considère quelques exemples en faisant varier la taille de la matrice [K] pour présenter uniquement les positions des cases non nulles (sans calcul de la valeur réelle) et de calculer le nombre des zéros.

La demi-largeur de la demi-bande est égale à la différence maximale plus un (problème des intervalles) entre les numéros de deux nœuds d'un même élément, multiplié par le nombre de degré de liberté des nœuds (pour notre cas, chaque nœud contient un degré de liberté). Si on a un domaine plan rectangulaire, formé par NE éléments triangulaires, la valeur de la demi bande B de la matrice, avant l'application des conditions aux limites, est donnée par:

Soit $B_L$ la valeur de la demi-bande donnée par l'élément numéro $L$ suivant. Alors,

<figure><img src="image_8.png" /><figcaption>Fig. 8: Elément fini triangulaire à trois nœuds</figcaption></figure>

$$B_L = \text{Max}(i, j, k) - \text{Min}(i, j, k) + 1 \qquad (60)$$

Par conséquent, la valeur de B est donnée par:

$$B = \text{Max}(B_1, B_2, B_3, \dots, B_{\text{NE}}) \qquad (61)$$

Considérons un domaine de forme rectangulaire qui sera discrétisé en trente (30) éléments triangulaires à trois nœuds, comme indique la figure suivante.

La numérotation des éléments ne pose aucun problème sur la valeur de la demi bande de la matrice de rigidité, c'est-à-dire, elle est indépendante de la bande. Chaque élément du domaine contient trois nœuds, alors la matrice de rigidité élémentaire (3×3).
<figure><img src="image_9.png" /><figcaption>Fig. 9: Numérotation des éléments</figcaption></figure>

On considère pour cet exemple, trois cas de numérotation des nœuds et pour chaque cas. On va construire la matrice de rigidité du domaine complet et cela indépendamment des conditions aux limites. Il est à noter que pour cet exemple, le nombre de nœuds est égal à 24, alors la matrice de rigidité $[K]$ est de rang $(24×24)$. On va considérer trois cas de numérotations.

Cas 1-

On va numéroter les nœuds suivant l'axe horizontal de gauche vers la droite en montant suivant l'axe vertical des y comme le montre la figure suivante.

<figure><img src="image_10.png" /><figcaption>Fig. 10: Exemple sur la numérotation des nœuds</figcaption></figure>

La matrice de rigidité de ce domaine par cette numérotation sur la base d'un élément triangulaire est donnée par la figure 11 suivante:

<figure><img src="image_11.png" /><figcaption>Fig. 11: Forme de la matrice de rigidité pour une numérotation choisie suivant le cas 1</figcaption></figure>

On peut remarquer que:

1. La largeur de la demi-bande B est égale à 8.
2. Le nombre des zéros à l'extérieur de la bande sur la partie supérieure est égal à 136 cases.

3. Il y a deux nœuds internes sont reliés avec sept (07) éléments finis.

4. Chaque nœud interne de maillage est commun au maximum entre six (06) éléments finis. Les nœuds de frontière, sauf les deux coins sont liés entre trois (03) éléments finis. Les deux points des coins sont liés en un seul élément fini.

### Cas 2-

On va numéroter maintenant les nœuds suivant l'axe vertical des y de bas en haut en déplaçant suivant l'axe horizontal x de gauche vers la droite, (Fig. 12).

<figure><img src="image_12.png" /><figcaption>Fig. 12: Exemple sur la numérotation des nœuds</figcaption></figure>

La matrice de rigidité de ce domaine par cette numérotation sur la base d'un élément triangulaire est donnée par la figure 13 suivante:

<figure><img src="image_13.png" /><figcaption>Fig. 13: Forme de la matrice de rigidité pour une numérotation choisie suivant le cas 2</figcaption></figure>

On peut remarquer dans ce cas que:

1. La largeur de la demi-bande B est égale à 6.

2. Le nombre des zéros à l'extérieur de la bande sur la partie supérieure est égal à 171 cases.

3. Il y a deux nœuds internes, qui sont reliés avec sept (07) éléments finis.

4. Chaque nœud interne de maillage est commun au maximum entre six (06) éléments finis. Les nœuds de frontière, sauf les deux coins sont liés entre trois (03) éléments finis. Les deux points des coins sont liés en un seul élément fini.

### Cas 3-

On va numéroter maintenant les nœuds d'une façon aléatoire, (Fig. 14). Notant ici, que nous avons plusieurs choix.

La matrice de rigidité de ce domaine par cette numérotation sur la base d'un élément triangulaire est donnée par la figure 15 suivante.
<figure><img src="image_14.png" /><figcaption>Fig. 14: Exemple sur la numérotation des nœuds</figcaption></figure>

On procède maintenant au calcul de la valeur de la bande de cette matrice de la figure 14. Comme la numérotation est aléatoire, on essaye de calculer la bande donnée pour chaque élément par la formule (60), ensuite la bande de la matrice complète par la formule (61). On aura:

$$B_1 = \text{Max}(1,2,13) - \text{Min}(1,2,13) + 1 = 13$$

$$B_2 = \text{Max}(1,18,2) - \text{Min}(1,18,2) + 1 = 13$$

$$B_3 = \text{Max}(18,20,2) - \text{Min}(18,20,2) + 1 = 19$$

$$B_4 = \text{Max}(18,19,20) - \text{Min}(18,19,20) + 1 = 3$$

$$B_5 = \text{Max}(19,11,20) - \text{Min}(19,11,20) + 1 = 10$$

Pour les éléments des numéros 6 jusqu'à 24, on trouve respectivement les valeurs suivantes 17, 9, 12, 10, 20, 9, 12, 6, 19, 17, 10, 6, 5, 5, 17, 11, 11, 16, 19, 19, 13, 18, 12, 12, 13.

<figure><img src="image_15.png" /><figcaption>Fig. 15: Forme de la matrice de rigidité pour une numérotation choisie suivant le cas3</figcaption></figure>

Alors, la valeur de la demi-bande de la matrice de rigidité [ K ] est calculée comme la valeur maximale de ces valeurs qui peut être donnée par la formule (61) donnant la valeur B=20. Pour ce cas, le nombre des zéros gagné est égal à 10 cases.

En conclusion, on déduit entre les trois cas de numérotation que le deuxième cas correspond à une numérotation optimale pour laquelle les cases non nulles se condensent le plus possible vers la diagonale principale et par conséquent le nombre des zéros sera très grand par rapport aux deux autres et que la bande est la plus petite.

Dans le cas général d'un domaine quelconque, il n'est pas facile d'obtenir une numérotation optimale. On est amené alors à faire une estimation pour obtenir la plus petite numérotation possible.

Finalement, si on veut généraliser à un domaine simplement connexe, quelconque, bidimensionnel avec conditions aux limites de type Dirichlet, qui peut être considéré ou
transformé vers un rectangle, avec une discrétisation de $N_x$ nœuds suivant l'axe horizontal des x et de $N_y$ nœuds suivant l'axe vertical des y, on aura les résultats suivants:

1. Le nombre total de nœuds NN du domaine est égal à $NN = N_x \times N_y$.

2. Le nombre total d'éléments finis triangulaires à trois nœuds dans le domaine est égal à $NE = 2(N_x - 1) \times (N_y - 1)$.

3. On commence la numérotation suivant la direction, qui possède le nombre de nœuds minimal:

Si $N_x \le N_y$, on suit la numérotation comme dans le cas 1.

Si $N_x > N_y$, on suit la numérotation suivant le cas 2.

4. La largeur de la demi bande avant l'introduction des conditions aux limites est égale à $B= \text{Min}(N_x, N_y)+2$. On propose l'algorithme suivant qui permet de calculer la largeur de la demi bande de la matrice de rigidité pour n'importe quel type de numérotation des nœuds.

B=0
POUR I=1 JUSQU'À NET
FAIRE
IP=MIN(NUMNT(I,1),NUMNT(I,2),NUMNT(I,3))
IG=MAX(NUMNT(I,1),NUMNT(I,2),NUMNT(I,3))
BI=IG-IP+1
SI BI>B ALORS B=BI
FIN FAIRE

5. Le nombre de nœuds, qui se trouvent sur la frontière du domaine, et où la solution du problème est connue, c'est-à-dire à éliminer de la résolution du système (conditions aux limites) est égal à $NDDLE = 2N_x + 2N_y - 4$.

6. Le nombre des inconnues restant du problème qui entrent dans la résolution du système d'équations, c'est-à-dire le rang de la matrice de rigidité est égal à $NDDLE = NN-NDLE = (N_x-2) \times (N_y-2)$.

7. La valeur de la demi bande finale de la matrice de rigidité pour la résolution du système d'équations, après l'introduction des conditions aux limites, est égal à:

$$B = \text{Min}(N_x, N_y) \qquad (62)$$

8. Le nombre des zéros qui se trouvent à l'extérieur de la bande (partie supérieure) de la matrice de rigidité [K] dans la phase finale après l'introduction des conditions aux limites est donné par:

$$NZERO = \frac{1+2+\cdots+N-B}{N-B} = \frac{(N-B) \times (N-B+1)}{2} \qquad (63)$$

### 3.3 Filtrage des nœuds

Il est parfois très efficace d'utiliser une option appelée filtrage des nœuds. Dans certain cas, la taille des cellules de maillage est non homogène. On trouve des cellules de grandes et de très petites tailles, surtout pour les points au voisinage de la frontière. Pour rendre la même grandeur aux tailles des cellules, on repositionne les points internes de maillage en utilisant la procédure présentée par la relation (64) suivante [1, 8]:
$$x_i^{\text{nouveau}} = x_i^{\text{ancien}} - \frac{\omega}{m} \sum_{k=1}^{k=m} (x_k^{\text{ancien}} - x_i^{\text{ancien}}) \quad (64)$$

$$y_i^{\text{nouveau}} = y_i^{\text{ancien}} - \frac{\omega}{m} \sum_{k=1}^{k=m} (y_k^{\text{ancien}} - y_i^{\text{ancien}})$$

Avec, $(x_k, y_k)$, $k = 1, 2, 3, \dots$, m. Les positions des ‘m’ nœuds reliés ave le nœud n° i.

L'algorithme (64), est à refaire entre 25 jusqu'à 100 itérations. On peut même aller à 1000 itérations. Pour notre type de maillage, chaque nœud interne est relié avec six (06) nœuds, sauf deux, les adjacents du coin droite bas et du haut gauche. Les nœuds de la frontière ne sont pas concernés par le filtrage.

### 3.4 Connexion des éléments finis du maillage

Après avoir fait la discrétisation du domaine en petits éléments, on obtient un ensemble de points et d'éléments finis. Chaque élément contient sa propre matrice de rigidité $[K^{(e)}]$ et le vecteur force $\{f^{(e)}\}$. Le problème consiste à faire l'assemblage de ces matrices élémentaires pour obtenir celle du domaine complet. Pour arriver aux résultats, il faut connaître les numéros des nœuds de chaque élément, alors, il faut déclarer une matrice d'ordre ($NE \times 3$), si les éléments obtenus sont triangulaires à trois nœuds. On désigne cette matrice par NUMNT. Par exemple, pour les trois cas de numérotations précédents, le vecteur NUMNT ($NE, 3$) sera rempli de la façon suivante.

Les trois valeurs de chaque ligne du vecteur NUMNT ($NE, 3$) représentent les numéros des lignes et des colonnes de la matrice de rigidité globale du domaine complet avant l'introduction des conditions aux limites.

<table><thead><tr><th></th><th colspan="3">Cas 1</th><th colspan="3">Cas 2</th><th colspan="3">Cas 3</th></tr><tr><th></th><th>1</th><th>2</th><th>3</th><th>1</th><th>2</th><th>3</th><th>1</th><th>2</th><th>3</th></tr></thead><tbody><tr><th>1</th><td>1</td><td>8</td><td>7</td><td>1</td><td>6</td><td>2</td><td>1</td><td>2</td><td>13</td></tr><tr><th>2</th><td>1</td><td>2</td><td>8</td><td>2</td><td>5</td><td>6</td><td>2</td><td>18</td><td>2</td></tr><tr><th>3</th><td>2</td><td>3</td><td>8</td><td>3</td><td>5</td><td>6</td><td>3</td><td>18</td><td>2</td></tr><tr><th>4</th><td>2</td><td>9</td><td>9</td><td>4</td><td>5</td><td>9</td><td>4</td><td>19</td><td>20</td></tr><tr><th>5</th><td>3</td><td>10</td><td>9</td><td>5</td><td>9</td><td>14</td><td>5</td><td>19</td><td>11</td></tr><tr><th>6</th><td>3</td><td>4</td><td>10</td><td>6</td><td>9</td><td>13</td><td>6</td><td>3</td><td>8</td></tr><tr><th>7</th><td>4</td><td>11</td><td>10</td><td>7</td><td>13</td><td>18</td><td>7</td><td>3</td><td>11</td></tr><tr><th>8</th><td>4</td><td>5</td><td>9</td><td>8</td><td>13</td><td>17</td><td>8</td><td>3</td><td>14</td></tr><tr><th>9</th><td>5</td><td>6</td><td>11</td><td>9</td><td>17</td><td>21</td><td>9</td><td>14</td><td>5</td></tr><tr><th>10</th><td>6</td><td>12</td><td>11</td><td>10</td><td>21</td><td>22</td><td>10</td><td>5</td><td>24</td></tr><tr><th>11</th><td>7</td><td>14</td><td>13</td><td>11</td><td>2</td><td>7</td><td>11</td><td>13</td><td>7</td></tr><tr><th>12</th><td>7</td><td>8</td><td>14</td><td>12</td><td>2</td><td>6</td><td>12</td><td>13</td><td>2</td></tr><tr><th>13</th><td>8</td><td>15</td><td>14</td><td>13</td><td>6</td><td>11</td><td>13</td><td>2</td><td>4</td></tr><tr><th>14</th><td>8</td><td>9</td><td>15</td><td>14</td><td>6</td><td>10</td><td>14</td><td>2</td><td>20</td></tr><tr><th>15</th><td>9</td><td>16</td><td>15</td><td>15</td><td>10</td><td>15</td><td>15</td><td>2</td><td>20</td></tr><tr><th>16</th><td>9</td><td>10</td><td>16</td><td>16</td><td>10</td><td>14</td><td>16</td><td>16</td><td>4</td></tr><tr><th>17</th><td>10</td><td>17</td><td>16</td><td>17</td><td>14</td><td>19</td><td>17</td><td>11</td><td>12</td></tr><tr><th>18</th><td>10</td><td>11</td><td>17</td><td>18</td><td>14</td><td>18</td><td>18</td><td>11</td><td>8</td></tr><tr><th>19</th><td>11</td><td>18</td><td>17</td><td>19</td><td>18</td><td>23</td><td>19</td><td>19</td><td>9</td></tr><tr><th>20</th><td>11</td><td>12</td><td>18</td><td>20</td><td>18</td><td>22</td><td>23</td><td>8</td><td>24</td></tr><tr><th>21</th><td>13</td><td>14</td><td>19</td><td>21</td><td>3</td><td>7</td><td>4</td><td>21</td><td>7</td></tr><tr><th>22</th><td>14</td><td>20</td><td>19</td><td>22</td><td>7</td><td>8</td><td>4</td><td>15</td><td>17</td></tr><tr><th>23</th><td>14</td><td>21</td><td>20</td><td>23</td><td>7</td><td>12</td><td>8</td><td>23</td><td>20</td></tr><tr><th>24</th><td>14</td><td>15</td><td>21</td><td>24</td><td>7</td><td>11</td><td>12</td><td>24</td><td>7</td></tr><tr><th>25</th><td>15</td><td>22</td><td>21</td><td>25</td><td>11</td><td>16</td><td>12</td><td>25</td><td>4</td></tr><tr><th>26</th><td>15</td><td>16</td><td>22</td><td>26</td><td>11</td><td>15</td><td>16</td><td>26</td><td>4</td></tr><tr><th>27</th><td>16</td><td>23</td><td>22</td><td>27</td><td>15</td><td>20</td><td>16</td><td>27</td><td>16</td></tr><tr><th>28</th><td>16</td><td>17</td><td>23</td><td>28</td><td>15</td><td>19</td><td>20</td><td>28</td><td>16</td></tr><tr><th>29</th><td>17</td><td>24</td><td>23</td><td>29</td><td>19</td><td>24</td><td>20</td><td>29</td><td>12</td></tr><tr><th>30</th><td>17</td><td>18</td><td>24</td><td>30</td><td>19</td><td>23</td><td>24</td><td>30</td><td>12</td></tr></tbody></table>

Par exemple, pour la ligne 24 du premier cas, on a les valeurs 14, 15 et 21, c'est-à-dire, l'élément triangulaire numéro 24 a les trois nœuds de numéros 14, 15 et 21 dans le sens antihoraire. Autrement dit, le nœud local 1 de ce triangle est égal à 14. On peut prendre le nœud local 1 égale à 15 ou 21, mais il faut que les trois numéros permutent dans le sens antihoraire.
Si la numérotation des nœuds est aléatoire, comme le cas 3, le remplissage de cette matrice sera faite d'une façon manuelle, vu l'inexistence d'une formule de récurrence.

### 3.5 Préparation du système d'équations [9]

Dans ce paragraphe, on va présenter les étapes nécessaires et la technique développée pour faire la préparation du système d'équations algébriques, plus précisément la matrice de rigidité [K] globale du domaine complet et le vecteur force {$f$} correspondant. Notre problème ici est de préparer le système d'équations avant et après l'introduction des conditions aux limites, en plus, sous forme de vecteur. Le rang de la matrice de rigidité du système avant l'introduction des conditions aux limites est égal à NN, et après l'introduction des conditions aux limites, le rang de la matrice sera réduit à la valeur NDDLR inférieure à NN. On remarque que le rang sera diminué par la valeur NDDLE. On note ici que la matrice de rigidité est singulière, c'est-à-dire, son déterminant est nul et cela avant l'introduction des conditions aux limites et que l'application de ces dernières enlève la singularité de la matrice et elle la rend régulière. Autrement dit, il existe une solution mathématique au système d'équations à résoudre.

Le problème qui réside ici est que les numéros des lignes et des colonnes à enlever de la matrice de rigidité avant l'introduction aux limites durant la phase de l'application de ces dernières sont les numéros des nœuds de la frontière du domaine. Pour arriver à l'enlèvement de ces lignes et colonnes, il faut déclarer un vecteur unicolonne de type entier de rang NN, nommé par exemple 'ICL' afin d'indiquer chaque nœud du domaine par sa présence ou l'absence dans la matrice de rigidité finale par un code. Le code utilisé ici est les deux chiffres 0 et 1.

- Si la case ICL (I) = 1, alors, le nœud de N°'1' sera présent dans la matrice finale.

- Si la case ICL (I) = 0, alors, le nœud de N°'1' sera absent dans la matrice finale.

Par exemple, pour les trois cas de numérotations des nœuds présentés auparavant, le vecteur ICL sera rempli de la façon suivante:

**Cas 1**

<table><thead><tr><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th><th>10</th><th>11</th><th>12</th></tr></thead><tbody><tr><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>1</td><td>1</td><td>1</td><td>1</td><td>0</td></tr><tr><td>13</td><td>14</td><td>15</td><td>16</td><td>17</td><td>18</td><td>19</td><td>20</td><td>21</td><td>22</td><td>23</td><td>24</td></tr><tr><td>0</td><td>1</td><td>1</td><td>1</td><td>1</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr></tbody></table>

**Cas 2**

<table><thead><tr><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th><th>10</th><th>11</th><th>12</th></tr></thead><tbody><tr><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>1</td><td>1</td><td>0</td><td>0</td><td>1</td><td>1</td><td>0</td></tr><tr><td>13</td><td>14</td><td>15</td><td>16</td><td>17</td><td>18</td><td>19</td><td>20</td><td>21</td><td>22</td><td>23</td><td>24</td></tr><tr><td>0</td><td>1</td><td>1</td><td>1</td><td>0</td><td>0</td><td>1</td><td>1</td><td>0</td><td>0</td><td>0</td><td>0</td></tr></tbody></table>

**Cas 3**

<table><thead><tr><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th><th>10</th><th>11</th><th>12</th></tr></thead><tbody><tr><td>0</td><td>1</td><td>0</td><td>1</td><td>0</td><td>0</td><td>1</td><td>1</td><td>0</td><td>0</td><td>1</td><td>1</td></tr><tr><td>13</td><td>14</td><td>15</td><td>16</td><td>17</td><td>18</td><td>19</td><td>20</td><td>21</td><td>22</td><td>23</td><td>24</td></tr><tr><td>0</td><td>0</td><td>0</td><td>1</td><td>0</td><td>0</td><td>0</td><td>1</td><td>0</td><td>0</td><td>0</td><td>0</td></tr></tbody></table>
On peut tirer de ce remplissage que le nombre des cases contenant le chiffre '1' est le même pour les deux cas et est égal à 8 pour ces exemples. Autrement dit, ce nombre est le rang de la matrice de rigidité finale, après l'application des conditions aux limites. Par conséquent, le nombre des zéros contenant dans le vecteur ICL est celui NDDLE.

Si la numérotation est aléatoire, le remplissage du vecteur 'ICL' devient d'une façon manuelle vue l'absence d'une formule de récurrence, mais pour notre numérotation adoptée dans le calcul, on peut trouver une formule de cette dernière et le remplissage devient automatique. On propose la procédure suivante suivant le cas 2. On aura:

DO I = 1, NN
ICL (I) = 0
ENDDO
DO J = 2, NY - 1
DO I = 2, NX - 1
K = (J-1) * NX + I
ICL(K) = 1
ENDDO
ENDDO

Après qu'on ait fait le remplissage du vecteur 'ICL', pour la connaissance des lignes et des colonnes enlevées de la matrice de rigidité pendant l'introduction des conditions aux limites, il nous reste à connaître le cas contraire, c'est-à-dire, le numéro du nœud dans le maillage correspondant à la ligne ou à la colonne de la matrice de rigidité finale (valable pour la résolution du système d'équations).

La forme de la matrice de rigidité des trois cas de numérotation après l'introduction des conditions aux limites possède la forme suivante, et cela après l'enlèvement des lignes et des colonnes possédant les numéros des nœuds de la frontière du domaine. Rappelons ici que le rang de la matrice devient (8 × 8). On peut regarder exactement les numéros des lignes les colonnes restants. On aura:

Cas 1

<figure><img src="image_16.png" /><figcaption>Fig. 16: Forme de la matrice [K] du cas 1 après application des conditions aux limites</figcaption></figure>

On peut remarquer à titre d'information que le nombre de case non nul à l'intérieur de la matrice [K] est le même, quelque soit le type de numérotation. La différence réside au niveau de leurs positions. A l'extérieur des matrices dans les figures 16, 17 et 18, on voit des petits chiffres. Normalement ces numéros commencent de 1 jusqu'à NDDLR, puisque la solution mathématique d'un système d'équations à la base d'une matrice d'ordre NDDLR donne les NDDLR solutions numérotées de 1 jusqu'à NDDLR. Mais, pour notre problème, ces solutions représentent les résultats aux nœuds internes du domaine qui possèdent des numéros pas forcement de 1 jusqu'à NDDLR.
Par exemple, la première solution du système sur la base des matrices de ces figures donne mathématiquement la valeur de $u(1)$ pour les trois cas, mais cette valeur représente la solution au nœud 8 pour le cas 1 et la solution au nœud 6 pour le cas 2.

Pour arriver à faire ce transfert, on déclare un vecteur de type entier, nommé par exemple par ‘IPOS’ de dimension NDDLAR, contient pour chaque case par ordre croissant, le numéro du nœud interne. Pour expliciter ce paragraphe, considérons le remplissage du vecteur ‘IPOS’ pour les trois cas.

### Cas 2

<figure><img src="image_17.png" /><figcaption>Fig. 17: Forme de la matrice [ K ] du cas 2 après application aux conditions aux limites</figcaption></figure>

### Cas 3

<figure><img src="image_18.png" /><figcaption>Fig. 18: Forme de la matrice [ K ] du cas 3 après application des conditions aux limites</figcaption></figure>

### Cas 1

$$ \text{IPOS} = \begin{bmatrix} 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 \\ 8 & 9 & 10 & 11 & 14 & 15 & 16 & 17 \end{bmatrix} $$

### Cas 2

$$ \text{IPOS} = \begin{bmatrix} 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 \\ 6 & 7 & 10 & 11 & 14 & 15 & 16 & 17 \end{bmatrix} $$

### Cas 3

$$ \text{IPOS} = \begin{bmatrix} 1 & 2 & 3 & 4 & 5 & 6 & 7 & 8 \\ 2 & 4 & 7 & 8 & 11 & 12 & 16 & 20 \end{bmatrix} $$

Le remplissage de ce vecteur peut être fait à partir du vecteur ‘ICL’, quelque soit le type de numérotation. On propose l’algorithme suivant pour faire ce remplissage.

J = 1
DO K = 1, NN
IF (ICL (K).EQ.1) THEN
IPOS (J) = K
J = J + 1
ENDIF
ENDDO

Maintenant, si on fait une petite comparaison sur le changement de la forme de la
matrice de rigidité du cas 1, cas 2 ou cas 3, on remarque que la valeur de la demi-bande
de la matrice de rigidité est changée. Autrement, elle diminue ou reste identique. Pour le
cas 1, la valeur ancienne est égale à 6 et la nouvelle devient 6. Pour le cas 2, elle était 8
et devient égale à 4. Pour le cas 3, elle était 20 et devient 8. Pour notre cas de
numérotation, on a déjà donné sa nouvelle valeur, après l'introduction des conditions
aux limites par la formule (62) qui sera valable uniquement pour le cas 2.

Dans le cas général, on propose l’algorithme suivant permettant la détermination de
la valeur de la demi bande après introduction des conditions aux limites pour tout type
de numérotation des nœuds et cela par:

IBAP = 0
DO 70 I = 1, NN
IF (ICL (I) . EQ . 0) GOTO 70
       NPCNN = I
       NDCNN = NPCNN + IBAV - 1
       IF (NDCNN . GT . NN) NDCNN = NN
       IBB = 0
       DO 60 J = NPCNN, NDCNN
       IF (ICL (J) . EQ . 1) IBB = IBB + 1
60 CONTINUE
       IF (IBB . GT . IBAP) IBAP = IBB
70 CONTINUE

On procède maintenant aux étapes nécessaires à suivre pour préparer la matrice de
rigidité globale et le vecteur force finaux du domaine complet avec la proposition de
l'algorithme de calcul. On rappelle que l'élément fini utilisé est celui du triangulaire à
trois nœuds.

La matrice de rigidité du cas 2 est obtenue après avoir assembler les 30 éléments finis. Si les conditions aux limites sont non nulles, il faut modifier le vecteur { F } en ajoutant les effets des valeurs de frontière.

4. RESOLUTION NUMERIQUE DU PROBLEME

Cette partie est consacrée à la résolution du problème numériquement. On a réalisé
un programme en langage Fortran dans ce contexte. Le travail a été limité pour des
matrices [ K ] symétriques et bandes. Ces deux propriétés nous ont conduit à stocker la
partie supérieure (ou inférieure) de la matrice en éliminant les zéros qui se trouvent à
l'extérieur de la bande, dans un vecteur uni colonne (ou uni-ligne), qui donne
l'équivalence des cases pour la matrice [ K ] dans le vecteur { vk } qui représente la
matrice de rigidité fictive.

La méthode de calcul est simple et consiste à donner à chaque fois qu'on tombe dans
une case de la matrice K (i, j), sa position dans le vecteur { vk } à partir d'une certaine
formule de récurrence qui varie avec la façon de stockage. Si une fois, on tombe dans
une case K (i, j) au dessus de la diagonale (dans la partie supérieure), mais à
l'extérieur de la bande, le calcul sera ignoré puisque cette case est réellement nulle et
elle ne présente aucun changement de résultat.
Si une fois, on tombe dans une case K(i, j) au dessous de la diagonale (dans la partie inférieure), il suffit uniquement d'utiliser la propriété K(i, j) = K(j, i) (puisque dans le programme, on a choisi le stockage de la partie supérieure) et le traitement se fait comme une case K(i, j) avec les trois cas précédents.

## 4.1 Méthodes de stockage sous forme de vecteur

Considérons une matrice [K] d'ordre (N × N) symétrique et bande avec la largeur de la demi-bande est B. Le problème est de savoir stocker les cases de cette matrice dans un vecteur {$vk$}. Il est parfois plus difficile de faire une équivalence entre une case de la matrice [K] et la case de vecteur {$vk$}, ce qui a supposé plusieurs formes de stockage sous la forme de vecteur.

Ces formes de stockage sont différentes suivant la manière de stockage et la formule de récurrence. On propose et on utilise dans notre travail la forme de stockage en colonne dans le principe est présenté comme suit. On prend une matrice [K] d'ordre (7×7) symétrique et sa bande $B=4$ et cela pour fixer les idées [8], [16], [18], [20], [21], [28-32], [34].

La méthode de stockage choisie est représentée par la direction des flèches, c'est-à-dire, colonne par colonne, comme le montre la figure 19.

<figure><img src="image_19.png" /><figcaption>Fig. 19: Stockage en colonne</figcaption></figure>

On peut remarquer pour cet exemple, que la dimension nécessaire du vecteur {$vk$} est égale à $NN = 22$. Maintenant si on veut généraliser à une matrice [K] d'ordre (N × N) symétrique avec largeur de la demi-bande égale à B, on aura:

$$NN = \frac{B \times (2N - B + 1)}{2} \qquad (65)$$

On peut vérifier la validité de cette formule pour l'exemple précédent:

$$\left\{ \begin{array}{l} N = 7 \\ B = 4 \end{array} \rightarrow NN = \frac{4(2 \times 7 - 4 + 1)}{2} = 22 \right. \qquad (66)$$

On procède maintenant à la détermination de la formule de stockage des éléments de la matrice [K] (les éléments de la bande supérieure uniquement) dans le vecteur {$vk$}, c'est-à-dire, l'équivalence d'une case de la matrice [K] dans le vecteur {$vk$}. Pour trouver la formule de récurrence, on divise la matrice [K] en deux blocs.

Pour une case K(i, j), telle que $j \le B$, sa position dans le vecteur {$vk$} est $vk(L)$ avec:

$$L = 1+2+3+\cdots+(j-1)+i = \frac{(j-1)j}{2} + i \qquad (67)$$
On peut remarquer que le nombre $(j-1)j$ est un nombre pair. Pour une case K(i, j), telle que $j \le B$, sa position dans le vecteur $\{vk\}$ est donnée par:

$$L = i - j + \frac{B \times (1 - B + 2j)}{2} \quad \text{si} \quad j > B \quad (68)$$

Pour la formule (67), le compteur '$i$' varie $1 \le i \le j$.

Pour la formule (68), le compteur '$i$' varie $j(B-1) \le i \le j$.

La fonction $L(i, j)$ est continue au point $j = B$, c'est-à-dire, pour $j = B$, on peut appliquer n'importe quelles formules (67) ou (68).

La méthode de stockage des matrices bandes sous forme de vecteur nous fait gagner en espace mémoire l'ordre suivant. Pour une matrice $[K]$ d'ordre $(NN)$, symétrique et bande, on a:

$$NDCE = N^2 - NN = N^2 - \frac{B(2N - B + 1)}{2} \quad (69)$$

Au lieu de stocker $N^2$ cases, on stocke uniquement $NN$ cases. Pour voir un ordre sur les valeurs de NDCE et $NN$, on propose les exemples suivants:

$$N = 100 \quad N^2$$

<table><thead><tr><th>B</th><th>N<sup>2</sup></th><th>NN</th><th>NDCE</th><th>Cases gagnées</th></tr></thead><tbody><tr><td>100</td><td>10000</td><td>5050</td><td>4950</td><td>49.50 %</td></tr><tr><td>70</td><td>10000</td><td>4585</td><td>5415</td><td>54.15 %</td></tr><tr><td>35</td><td>10000</td><td>2905</td><td>7095</td><td>70.95 %</td></tr><tr><td>10</td><td>10000</td><td>955</td><td>9045</td><td>90.45 %</td></tr><tr><td>2</td><td>10000</td><td>199</td><td>9801</td><td>98.01 %</td></tr></tbody></table>

$$N = 5000 \quad N^2$$

<table><thead><tr><th>B</th><th>N<sup>2</sup></th><th>NN</th><th>NDCE</th><th>Cases gagnées</th></tr></thead><tbody><tr><td>480</td><td>2500000</td><td>2285040</td><td>22714960</td><td>90.85 %</td></tr><tr><td>270</td><td>2500000</td><td>1313685</td><td>23686315</td><td>94.74 %</td></tr><tr><td>150</td><td>2500000</td><td>738825</td><td>24261175</td><td>97.04 %</td></tr><tr><td>50</td><td>2500000</td><td>347585</td><td>24652415</td><td>98.60 %</td></tr></tbody></table>

Il est bien clair d'après ces exemples, que plus la largeur de la demi bande est petite, plus le nombre des cases gagnées augmente rapidement, ce qui veut dire que, durant la numérotation des nœuds du maillage généré, il est important de choisir celle qui donne la bande minimale, ce que l'on appelle par numérotation optimale. On peut démontrer que toute région plane ou surfacique peut être transformée ou considérée comme un carré ou rectangle.

On peut démontrer que le nombre des nœuds restant après l'entrée des conditions aux limites (si ces dernières sont de type Dirichlet), c'est-à-dire l'ordre de la matrice $[K]$ pour la résolution est donnée par $N = N_x \times N_y - 2$, $N_x - 2$, $N_y + 2$. Dans ce supposé que $N_x \le N_y$. Dans ce cas, la dimension nécessaire du vecteur $\{v_k\}$ pour le
cas le plus défavorable est donnée par $B_m=N_x+2$. Pour bien voir l'ordre de la dimension, ainsi que le pourcentage des cases gagnées, on considère l'exemple suivant pour $N_y=100$.

<table><thead><tr><th>N<sub>x</sub></th><th>B<sub>max</sub></th><th>N</th><th>NN</th><th>NDCE</th><th>Cases gagnées</th></tr></thead><tbody><tr><td>100</td><td>102</td><td>9604</td><td>974457</td><td>91262359</td><td>98.94 %</td></tr><tr><td>80</td><td>82</td><td>7644</td><td>623487</td><td>57807249</td><td>98.93 %</td></tr><tr><td>50</td><td>52</td><td>4704</td><td>243282</td><td>21884334</td><td>98.90 %</td></tr><tr><td>30</td><td>32</td><td>2744</td><td>87312</td><td>7442224</td><td>98.84 %</td></tr><tr><td>10</td><td>12</td><td>784</td><td>9342</td><td>605314</td><td>98.48 %</td></tr><tr><td>3</td><td>5</td><td>98</td><td>480</td><td>9124</td><td>95.00 %</td></tr></tbody></table>

On remarque bien que le nombre de cases gagnées dans les cas les plus défavorables dépasse 98 % de la méthode classique de stockage complet. Donc on a bien vu l'intérêt pratique de la technique de la matrice bande. Après avoir résolu le problème de stockage, on va procéder maintenant à la description qualitative des méthodes de résolution d'un système adapté pour des matrices bandes, stockées sous forme d'un vecteur.

## 4.2 Algorithme de résolution de système d'équations

On est amené à résoudre le problème après avoir rempli la matrice de rigidité [K] globale et le vecteur {F} afin de déterminer le vecteur {u} de solution à un système d'équations suivant [12, 13], [18], [20], [28, 29], [31, 32]:

$$[K] \times \{u\} = \{F\} \qquad (70)$$

On a choisi la méthode de résolution numérique directe de Khaletski [7], [25]. Pour une matrice symétrique et bande, l'algorithme de résolution du système (64) est présenté en [27]. Mettons la matrice [K] sous la forme d'un produit de deux matrices triangulaires inférieure $[L] = [L_{ij}]$ et de la matrice triangulaire supérieure $[H] = [H_{ij}]$ à diagonale unité, c'est-à-dire:

$$[K] = [L] \times [H] \qquad (71)$$

On aura par suite deux systèmes d'équations à matrice triangulaire:

$$[K] \times \{u\} = \{F\} \quad \rightarrow \quad [L] \times [H] = \{f\}$$

On pose:

$$[K] \times \{u\} = \{y\} \qquad (72)$$

$$[L] \times \{y\} = \{F\} \qquad (73)$$

Si la matrice [K] est symétrique, on peut démontrer que:

$$H_{ij} = L_{ji} / L_{ii} \qquad (74)$$

La matrice [K] possède la forme générale suivante:
$$
\left[
\begin{array}{ccccccc}
K_{11} & & & & & & \\
K_{21} & K_{22} & & & & & \\
K_{31} & K_{32} & K_{33} & & & & \\
\vdots & \vdots & & \ddots & & & \\
K_{B1} & K_{B2} & & & & & \\
0 & K_{B+1,2} & & & & & \\
\vdots & \vdots & K_{ii} & & & & \\
\vdots & \vdots & & \ddots & & & \\
0 & 0 & 0 & \cdots & K_{NN} & & \\
\end{array}
\right]
\quad (75)
$$

La matrice [ L ] qui vient de la décomposition de la matrice [ K ] possède la même
structure que [ K ], sauf qu’elle est triangulaire. Après l’élimination des opérations sur
les zéros (les éléments à l’extérieur de la bande), les éléments de la matrice [ L ] se
calculent par les relations suivantes [24]:

$$
\left\{
\begin{aligned}
L_{ij} &= K_{ij} - \sum_{r=1}^{r=j-1} L_{ir} \times \frac{L_{jr}}{L_{rr}} && \begin{cases} i=1,2,\dots,B \\ j=2,3,\dots,i \end{cases} \\
L_{ij} &= K_{ij} - \sum_{r=i-B+1}^{r=j-1} L_{ir} \times \frac{L_{jr}}{L_{rr}} && \begin{cases} i=B+1,\dots,N \\ j=i-B+1,\dots,i \end{cases}
\end{aligned}
\right.
\qquad (76)
$$

Egalement, les vecteurs { y } et { u } se calculent par l’algorithme suivant [24]:

$$
\left\{
\begin{array}{ll}
y_1 = F_1 / L_{11} & \\
y_i = \frac{1}{L_{ii}} \times \left( F_i - \sum_{r=1}^{r=i-1} L_{ir} \times y_r \right) & i=2,3,\dots,B \\
y_i = \frac{1}{L_{ii}} \times \left( F_i - \sum_{r=i-B+1}^{r=i-1} L_{ir} \times y_r \right) & i=B+1,\dots,N
\end{array}
\right.
\qquad (77)
$$

$$
\left\{
\begin{array}{ll}
u_N = y_N & \\
u_i = y_i - \frac{1}{L_{ii}} \times \sum_{r=i+1}^{r=N} L_{ir} \times u_r & i=N-1, \dots, N-B+1 \\
u_i = y_i - \frac{1}{L_{ii}} \times \sum_{r=i+1}^{r=N} L_{ir} \times u_r & i=N-B, \dots, 3, 2, 1
\end{array}
\right.
\qquad (78)
$$

Les éléments de la diagonale de la matrice [ K ] sont tous positifs et non nuls. De
même elle st définie positive.

4.3 Calcul des paramètres intéressants

Après la détermination de la solution du vecteur { u } dans chaque nœud du
maillage, notre intérêt s'oriente vers le calcul de quelques paramètres intéressants de
l'application de l'équation de Poisson. La discipline mécanique des fluides, par
exemple, à l'intérêt sur le calcul de l'écoulement dans des conduites, tels que l'air pour
pressuriser des cabines des avions, alimenter les réservoirs de carburant des engins,
conduire l'huile pour le graissage des arbres des moteurs etc. Donc, on a besoin de
connaître quelques résultats qui reflètent l'écoulement, tels que, le coefficient de perte de charge, vitesse maximale, débit du fluide traversant la conduite et vitesse moyenne.

La discipline 'Résistance des matériaux', donne le résultat de la torsion des poutres des sections uniquement circulaires. Si la section de la poutre est non circulaire, le calcul de la contrainte de torsion et l'angle de torsion seront données par l'élasticité linéaire, après faire la résolution numérique de l'équation de Poisson où la variable u de l'équation est appelée dans ce cas, flux de contrainte. Donc on a l'intérêt de faire le calcul de la formule de torsion sous la forme d'une certaine constante et la distribution de contrainte de torsion dans la surface de la section de la poutre. Le calcul de la valeur moyenne d'une certaine variable continue, et en particulier, vitesse moyenne de l'écoulement traversant une conduite a besoin de la valeur de l'aire de la surface de la conduite. Donc notre intérêt s'oriente encore vers le calcul de l'aire de section des conduites. Les positions des nœuds de maillage doivent être données par rapport à un repère passant par le centre de gravité de la section, ce qui veut dire, l'intérêt de calcul de la position du centre de gravité de cette section.

### 4.3.1 Aire et position du centre de gravité de la surface de la section

La surface du domaine a été subdivisée en petits éléments finis de forme triangulaires, ou quadrilatères de forme simples. Alors, l'aire totale de cette section sera approximer par la somme de toutes les aires des surfaces des éléments élémentaires simples. On peut écrire, par conséquent:

$$A = \sum_{i=1}^{i=NE} A_i \qquad (79)$$

Si l'élément fini est un quadrilatère (ou rectangle), la surface de cet élément peut être considéré comme l'union de deux triangles. Pour une section subdivisée en petits éléments de géométrie simples, le centre de gravité de cette surface composée est donnée par:

$$x_G = \frac{\sum_{i=1}^{i=NE} (x_G)_i A_i}{A} \quad ; \quad y_G = \frac{\sum_{i=1}^{i=NE} (y_G)_i A_i}{A} \qquad (80)$$

Si l'élément fini est un triangle, sa surface est donnée par la formule (13) et la position de son centre gravité d est donnée par les formules (21) et (22).

### 4.3.2 Intégrale de surface de la solution u du problème

Le calcul des paramètres physiques des phénomènes de la nature utilisent le résultat de certaine intégrale de surface de la solution du vecteur { u } pour déterminer par exemple le débit volumique du fluide à travers une section de forme quelconque, le moment de torsion M<sub>t</sub> exercé sur une poutre de section non circulaire, la vitesse moyenne u<sub>moy</sub>, le courant électrique à travers des fils. Ces paramètres sont donnés par le résultat suivant [9].

Par conséquent, l'intégrale suivante a un grand intérêt physique. Pour cela, on va déterminer mathématiquement l'approximation de résultat de cette intégrale en utilisant le principe de la discrétisation et les fonctions d'interpolations. Soit:
$$H = \iint_A u(x,y) \, dx \times dy \qquad (81)$$

Avec $H$, représentant le débit volumique par unité de profondeur pour le cas du mouvement des fluides et par exemple le moment de torsion pour le cas des problèmes de l'élasticité 2D, et $u(x,y)$ qui représente la distribution de vitesse dans ce cas de mouvement de fluide, la distribution de flux de contrainte pour le cas de problème de torsion des poutres de section non circulaire. La vitesse moyenne d'un fluide à travers une section droite est donnée par:

$$H_{\text{moy}} = \frac{1}{A} \times \iint_A u(x,y) \, dx \times dy \qquad (82)$$

On peut considérer l'intégrale comme la somme des intégrales dans les petits éléments finis du maillage. La fonction $u(x,y)$ dans l'élément fini est donnée par la formule (11) pour un élément triangulaire à trois nœuds.

Pour un élément triangulaire à trois nœuds, les fonctions d'interpolations sont données par la formule (12). Remplaçant la relation (12) dans la relation (81), en utilisant les relations (21), (22), (23), (24) et (25):

$$H = \frac{1}{3} \times \sum_{i=1}^{i=\text{NE}} A_i \times (u_1 + u_2 + u_3)_i \qquad (83)$$

La formule (83) donne le résultat de l'intégrale de surface de la solution $\{ u \}$ pour une discrétisation en éléments finis triangulaires à trois nœuds.

Parfois il est très intéressant pour certains phénomènes physiques, par exemple, le problème de l'étude de la torsion des poutres de section non circulaire, de calculer la dérivée de la solution $u$ dans tous les points du maillage considéré. On parle de la dérivée par rapport à $x$ et la dérivée par rapport à $y$.

Ces dérivées représentent la distribution des contraintes $\tau_{xz}$ et $\tau_{yz}$. Dans ce cas, il faut déterminer la contrainte maximale ainsi que sa position dans le domaine.

### 4.3.3 Position et valeur maximale de la solution $u$ du problème

Il est parfois utile de déterminer la valeur et la position de la solution maximale de $\{ u \}$ du problème. Par exemple, si le vecteur $\{ u \}$ du problème représente la distribution de contrainte d'une sollicitation élastique, il est très intéressant de déterminer sa valeur et sa position afin de localiser la section et le point dangereux qui provoque la rupture, ou, pour faire un dimensionnement d'une structure travaillante sous une contrainte inférieure à une certaine contrainte admissible de rupture, déterminée par les essais pratiques de la résistance des matériaux.

Si par exemple la solution $\{ u \}$ du problème représente la distribution de la vitesse de l'écoulement, il est utile de déterminer la position et la valeur maximale de cette vitesse afin de déterminer la vitesse débitante. Donc nous somme en présence d'un vecteur contient un nombre important de valeur et le problème consiste à déterminer la valeur maximale et la position dans le domaine physique. Il est possible de voir plusieurs points dans le maillage possédants la même valeur maximale de la solution $\{ u \}$.
# 5. RESULTATS EST COMMENTAIRES

Après une étude sur la résolution numérique de l'équation de Poisson et les techniques numériques employées pour résoudre ce problème par la méthode des éléments finis, on arrive au stade de la présentation des résultats obtenus, qu'on a divisés en deux parties. Pour chaque partie, on a réalisé un programme à part en langage Fortran. La première partie consiste à l'obtention et la visualisation des maillages obtenus de type 'H' dans des configurations simples et complexes. La deuxième partie consiste à la présentation des résultats obtenus sur la résolution de l'équation de Poisson pour des configurations choisies.

## 5.1 Résultats sur la génération de maillage

Les résultats donnés par ce programme sont choisis pour des configurations simples, telles que le cercle et pour d'autres complexes, telles que la configuration complexe. D'autres maillages sont aussi présentés pour des configurations touchant le domaine de l'aéronautique. Le premier exemple choisi, c'est celui du domaine circulaire, où on fait le changement des paramètres nécessaires et on voit pour chaque cas la forme du maillage obtenu. L'utilisation d'un tel maillage pour résoudre l'équation de Poisson sera effectuée sur la base de la variation du gradient de la solution u du problème.

La figure 20 présente un exemple sur la forme du maillage, lorsqu'on applique la procédure de filtrage des nœuds au maillage obtenu par notre technique. L'exemple choisi est pour $N_x = 35$ et $N_y = 30$. Le facteur de relaxation $\omega=0.3$ et le nombre d'itérations est choisi à 250. Il est à noter ici que cette option a été considéré pour tous les points internés du maillage. On note encore qu'on peut, à partir de cette option de filtrage de considérer plusieurs formes de maillage en jouant sur les paramètres, nombre d'itérations, la manière de filtrer les nœuds et le nombre et les numéros des nœuds à filtrer. Dans cet exemple, on a 1050 nœuds et 1972 triangles.

<figure><img src="image_20.png" /><figcaption>Fig. 20: Maillage triangulaire filtré dans un domaine circulaire</figcaption></figure>

Le deuxième domaine présenté dans ce modeste travail est le cas de la surface d'un profil d'aile. L'exemple choisi ici est le profil NACA 0012 symétrique. Ce domaine a un grand intérêt pour les chercheurs intéressés par l'aérodynamique et l'élasticité des structures. Par exemple, si on peut chercher et étudier le problème de la torsion des pales d'hélicoptères, des pales d'hélices des moteurs, les aubes du compresseur et autres cas, qu'il faut générer ce type de maillage.
On veut par exemple déterminer l'aire et les caractéristiques géométriques (position du centre de gravité, moments et produits d'inerties) d'un tel profil d'aile qui nous servira pour une étude de résistance et de rigidité. On s'intéresse encore à ce type de maillage.

<figure><img src="image_21.png" /><figcaption>Fig. 21: Maillage triangulaire filtré dans un domaine pour une configuration d'un profil d'aile</figcaption></figure>

La figure 21 illustre un exemple de maillage avec utilisation de la procédure de filtrage des nœuds. Il est obtenu après avoir utiliser la procédure de filtrage des nœuds, à l'ensemble des nœuds de maillage avec $ω=0.4$ et le nombre d'itérations égal à 1000. Il manque un peu de raffinement au niveau du bord d'attaque du profil.

On peut affiner encore le maillage, mais le dessin ne sera pas clair. On peut jouer encore sur les valeurs de $ω$, nombre d'itérations de la procédure de filtrage, positions des nœuds A, B, C et D de frontière, le nombre et les numéros des nœuds à filtrer, région du profil à filtrer, donc on peut obtenir des maillages assez efficaces (sous réserves) qui seront aptes à résoudre l'équation de Poisson avec une précision voulue.

La figure 22 représente le même maillage que celui du profil. Le but de ce maillage, est de voir l'efficacité du programme réalisé, pour générer des maillages dans des configurations très complexes. L'exemple, avec la procédure de filtrage des nœuds et des cellules triangulaires, est donc représenté en figure 25.

<figure><img src="image_22.png" /><figcaption>Fig. 22: Maillage triangulaire filtré dans une configuration complexe</figcaption></figure>

Le but de la figure 23 est de réaliser des maillages dans des configurations très complexes, telles que la navette spatiale. On a effectué avec intérêt la projection de cet engin sur un plan. On peut voir exactement l'effet de la procédure de filtrage sur les dimensions des cellules, surtout au niveau de la frontière du domaine. Ce type de domaine est très intéressant, lorsqu'on veut calculer la distribution de température sur la surface de l'engin au moment de son entrée dans l'atmosphère.

Un autre exemple de maillage, lorsque la frontière du domaine possède des grandes courbures répétées, et où au voisinage de ces courbures et coins demandant une distribution et une condensation de nœuds rigoureux, sont présentés dans la figure 24.
L'application du maillage de la figure 24 sera pour des problèmes de la torsion des clés mécaniques ayant des sections comme celles montrées dans ces figures.

<figure><img src="image_23.png" /><figcaption>Fig. 23: Maillage triangulaire filtré dans une configuration complexe</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Fig. 24: Maillage triangulaire dans une configuration complexe</figcaption></figure>

Le maillage dans un domaine rectangulaire dont le coté horizontal est égal à $x=a$ et le coté verticale est égale à $y=b$ est présenté dans la figure 25.

<figure><img src="image_25.png" /><figcaption>Fig. 25: Maillage triangulaire dans un domaine rectangulaire</figcaption></figure>

## 5.2 Résultats sur la résolution numérique de l'équation de Poisson

Dans la deuxième partie, on va présenter quelques résultats obtenus lors de la résolution de l'équation de Poisson avec les conditions aux limites de type Dirichlet. On a choisi ici le deuxième membre égal à l'unité pour faire des simplifications.
Notons ici que la majorité des problèmes physiques peuvent être modélisés par l'équation de Poisson, dont le deuxième membre est égal à une constante.

## 5.2.1 Domaine circulaire

Le premier domaine à étudier est celui le cercle de rayon R donné, avec les conditions aux limites $u=0$ sur $\Gamma$. La solution exacte pour ce domaine circulaire est présentée dans [35].

### 5.2.1.1 Exemple sur la distribution de la solution

L'exemple choisi ici est le cercle de rayon $R=1$ placé à l'origine avec la discrétisation de $N_x=6$ et $N_y=5$. On obtient par conséquent 30 points dans le domaine dont 18 points sont sur la frontière. Les points A, B, C et D de la frontière sont choisis aux milieux de chaque quadrant. On utilise la procédure de filtrage des nœuds avec 100 itérations et $\omega=0.1$, on obtient les résultats suivants qui seront comparés avec la solution exacte du cercle. On note ici qu'on n'a pas présenté les points de la frontière, puisque la valeur en ces points est égale à zéro.

<table><thead><tr><th colspan="2">Position du point</th><th colspan="2">Solution u(x,y) de l'équation</th></tr><tr><th>x</th><th>y</th><th>x<sub>calculée</sub></th><th>x<sub>exacte</sub></th></tr></thead><tbody><tr><td>0.49759</td><td>0.39798</td><td>0.14877</td><td>0.14850</td></tr><tr><td>0.15770</td><td>0.46540</td><td>0.19118</td><td>1.8963</td></tr><tr><td>0.17337</td><td>0.50985</td><td>0.18320</td><td>0.17750</td></tr><tr><td>0.48652</td><td>0.54406</td><td>0.12881</td><td>0.11682</td></tr><tr><td>0.56135</td><td>0.03968</td><td>0.17199</td><td>0.17083</td></tr><tr><td>0.18185</td><td>0.01866</td><td>0.24115</td><td>0.24165</td></tr><tr><td>0.16941</td><td>0.02500</td><td>0.24146</td><td>0.24267</td></tr><tr><td>0.54088</td><td>0.05012</td><td>0.17375</td><td>0.17623</td></tr><tr><td>0.57097</td><td>0.50104</td><td>0.11707</td><td>0.10574</td></tr><tr><td>0.19408</td><td>0.49930</td><td>0.18086</td><td>0.17826</td></tr><tr><td>0.14929</td><td>0.46111</td><td>0.19074</td><td>0.19127</td></tr><tr><td>0.49092</td><td>0.39458</td><td>0.14929</td><td>0.15082</td></tr></tbody></table>

On remarque bien que pour cette faible discrétisation, la solution déterminée est presque proche de la solution exacte avec une précision relative environ de 3 %. Les valeurs de A, H, $u_{moy}$ et $u_{max}$ sont données pour cette discrétisation par:

<table><thead><tr><th></th><th>A</th><th>u<sub>max</sub></th><th>u<sub>moy</sub></th><th>H</th></tr></thead><tbody><tr><th>Calculée</th><td>3.07581</td><td>0.24145</td><td>0.12145</td><td>0.35437</td></tr><tr><th>Exacte</th><td>3.14159</td><td>0.25000</td><td>0.12500</td><td>0.39269</td></tr></tbody></table>

### 5.2.1.2 Effet de la discrétisation sur la convergence du problème

Le but de ce paragraphe est de voir l'effet de la discrétisation sur la convergence du problème pour le même domaine circulaire précédent. On va changer le nombre de nœuds $N_x$ et $N_y$, plus précisément, augmenter ces valeurs, et pour chaque cas on calcule les valeurs de A, H, $u_{moy}$ et position et valeur de $u_{max}$ et encore on choisit
un point interne à l'intérieur du domaine et on calcule la valeur de $u$ de la solution en ce point, on aura les résultats suivants.

On remarque bien à partir des résultats obtenus par le programme réalisé, la convergence des valeurs vers ceux de la solution exacte. On remarque aussi que la convergence est assez lente, puisque chaque fois on augmente les valeurs de $N_x$ et $N_y$, on remarque une petite variation aux valeurs obtenues.

On remarque encore que ce n'est pas uniquement le but d'avoir des solutions précises d'augmenter le nombre de nœuds, il faut que la taille des cellules doivent être très petites, c'est-à-dire, la surface, la longueur des cotés et les trois angles internes doivent être aussi à discuter. Par exemple, on peut voir ce qu'on a dit entre les résultats obtenus entre la discrétisation de $N_x = 100$ et $N_y = 400$ et la discrétisation de $N_x = 200$ et $N_y = 200$. Pour les deux discrétisations, on a presque le même nombre de nœuds internes et le système à résoudre est presque le même. Malgré ça, les résultats obtenus pour la deuxième sont meilleurs par rapport aux premiers résultats. Il est à noter que la position de valeur maximale de la solution est obtenue au centre de gravité G du cercle $(0, 0)$.

$N_x = 50$

<table><thead><tr><th>N<sub>y</sub></th><th>A</th><th>u<sub>max</sub></th><th>u<sub>moy</sub></th><th>H</th></tr></thead><tbody><tr><td>50</td><td>3.14105</td><td>0.24975</td><td>0.12489</td><td>0.39229</td></tr><tr><td>100</td><td>3.14125</td><td>0.24984</td><td>0.12493</td><td>0.39244</td></tr><tr><td>150</td><td>3.14129</td><td>0.24985</td><td>0.12493</td><td>0.39247</td></tr><tr><td>200</td><td>3.14130</td><td>0.24986</td><td>0.12494</td><td>0.39248</td></tr><tr><td>250</td><td>3.14131</td><td>0.24986</td><td>0.12494</td><td>0.39248</td></tr><tr><td>300</td><td>3.14131</td><td>0.24986</td><td>0.12494</td><td>0.39248</td></tr><tr><td>400</td><td>3.14131</td><td>0.24987</td><td>0.12494</td><td>0.39248</td></tr><tr><td>500</td><td>3.14132</td><td>0.24987</td><td>0.12494</td><td>0.39248</td></tr></tbody></table>

$N_x = 200$

<table><thead><tr><th>N<sub>y</sub></th><th>A</th><th>u<sub>max</sub></th><th>u<sub>moy</sub></th><th>H</th></tr></thead><tbody><tr><td>200</td><td>3.14156</td><td>0.24998</td><td>0.12499</td><td>0.39267</td></tr><tr><td>300</td><td>3.14156</td><td>0.24998</td><td>0.12499</td><td>0.39268</td></tr><tr><td>400</td><td>3.14157</td><td>0.24999</td><td>0.12499</td><td>0.39268</td></tr><tr><td>500</td><td>3.14157</td><td>0.24999</td><td>0.12499</td><td>0.39268</td></tr></tbody></table>

On choisit maintenant un point interne arbitraire et voir l'effet de la discrétisation sur la valeur de la solution $u$ en ce point. Le point choisi se trouve à la position $x=0.0$, $y=0.33333$.

La solution exacte est égale à $u=0.22222$. On aura les résultats suivants. On note ici que la valeur de l'erreur relative est donnée par la formule suivante:

$$ \text{Érreur relative (\%)} = \frac{|u_{\text{calculée}} - u_{\text{exacte}}|}{u_{\text{exacte}}} \times 100 \qquad (84) $$
<table><thead><tr><th>N<sub>x</sub></th><th>N<sub>y</sub></th><th>u<sub>calculée</sub></th><th>u<sub>exacte</sub></th><th>Érreur relative %</th></tr></thead><tbody><tr><td>5</td><td>7</td><td>0.220934</td><td>0.222222</td><td>0.57941</td></tr><tr><td>9</td><td>13</td><td>0.221483</td><td>0.222222</td><td>0.33262</td></tr><tr><td>13</td><td>19</td><td>0.221839</td><td>0.222222</td><td>0.17339</td></tr><tr><td>37</td><td>55</td><td>0.222168</td><td>0.222222</td><td>0.02403</td></tr><tr><td>73</td><td>109</td><td>0.222207</td><td>0.222222</td><td>0.00667</td></tr><tr><td>145</td><td>217</td><td>0.222217</td><td>0.222222</td><td>0.00201</td></tr><tr><td>289</td><td>217</td><td>0.222217</td><td>0.222222</td><td>0.00126</td></tr></tbody></table>

On remarque clairement la puissance de la méthode des éléments finis, malgré le type de l'élément qui est le plus simple, un triangle à trois nœuds. On note ici que le temps d'exécution par ordinateur est assez grand pour le dernier exemple avec 289 et 217 nœuds suivant x et y respectivement.

## 5.2.2 Domaine rectangulaire

Les résultats exacts concernant le domaine rectangulaire sont donnés dans la référence [35].

On a varié le nombre n de 1 jusqu'à la valeur 10 et concernant le calcul par la méthode des éléments finis, on a utilisé un maillage de 120 par 120 nœuds et sans procédure de filtrage. Après exécution des programmes, on a obtenu les résultats suivants.

<table><thead><tr><th>n</th><th>H / a<sup>4</sup> (calculée)</th><th>u<sub>max</sub> / a<sup>2</sup> (calculée)</th><th>u<sub>moy</sub> / a<sup>2</sup> (calculée)</th></tr></thead><tbody><tr><td>1</td><td>0.03513</td><td>0.07365</td><td>0.03513</td></tr><tr><td>2</td><td>0.11431</td><td>0.11385</td><td>0.05715</td></tr><tr><td>3</td><td>0.19742</td><td>0.12266</td><td>0.06580</td></tr><tr><td>4</td><td>0.28071</td><td>0.12450</td><td>0.07017</td></tr><tr><td>5</td><td>0.36399</td><td>0.12489</td><td>0.07279</td></tr><tr><td>6</td><td>0.44727</td><td>0.12497</td><td>0.07454</td></tr><tr><td>7</td><td>0.53053</td><td>0.12498</td><td>0.07579</td></tr><tr><td>8</td><td>0.61379</td><td>0.12499</td><td>0.07672</td></tr><tr><td>9</td><td>0.69704</td><td>0.12499</td><td>0.07744</td></tr><tr><td>10</td><td>0.78028</td><td>0.12499</td><td>0.07802</td></tr></tbody></table>

En comparant les valeurs de la deuxième et troisième colonne, on peut remarquer la convergence par la méthode des éléments finis. On note ici, qu'on n'a pas trouvé dans la littérature les solutions exactes sur $u_{max}$ et $u_{moy}$.

## 5.2.3 Domaine elliptique

Le troisième domaine à étudier est celui de l'ellipse avec le coté vertical égale à b et la longueur horizontale égale à a. La solution exacte pour ce type de domaine est donnée par la référence [35]. Avant de passer aux calculs numériques de ces paramètres en utilisant le programme réalisé, on propose un exemple sur la distribution de la solution $u(x,y)$ en quelques points avec une discrétisation de 250 par 250 nœuds et une dimension de l'ellipse égale à a=3 et b=1. On aura les résultats suivants comparés avec ceux de la solution exacte.

De la même manière, on va déterminer les paramètres H, $u_{max}$ et $u_{moy}$ en utilisant le programme réalisé pour une discrétisation de 120 par 120 nœuds et cela en variant la dimension de l'ellipse, on aura les résultats numériques donnés dans le tableau suivant.
<table><thead><tr><th colspan="2">Position du point</th><th colspan="2">Solution u(x,y) de l'équation</th></tr><tr><th>x</th><th>y</th><th>x<sub>calculée</sub></th><th>x<sub>exacte</sub></th></tr></thead><tbody><tr><td>2.22200</td><td>0.23640</td><td>1.973377</td><td>1.973400</td></tr><tr><td>1.97340</td><td>0.27008</td><td>1.881731</td><td>1.881750</td></tr><tr><td>3.23392</td><td>0.27404</td><td>0.973490</td><td>0.973482</td></tr><tr><td>0.64630</td><td>0.42264</td><td>2.822653</td><td>2.822711</td></tr><tr><td>0.61246</td><td>0.45934</td><td>2.744908</td><td>2.744960</td></tr><tr><td>1.99497</td><td>0.46381</td><td>2.094747</td><td>2.094771</td></tr><tr><td>0.07928</td><td>0.07928</td><td>2.769831</td><td>2.769887</td></tr><tr><td>0.02603</td><td>0.70441</td><td>2.721037</td><td>2.721091</td></tr><tr><td>2.97457</td><td>0.61565</td><td>1.166052</td><td>1.166050</td></tr><tr><td>2.12272</td><td>0.75053</td><td>1.888660</td><td>1.888673</td></tr><tr><td>1.12118</td><td>0.88222</td><td>2.404637</td><td>2.404670</td></tr><tr><td>2.77569</td><td>0.83270</td><td>1.271308</td><td>1.271308</td></tr><tr><td>2.27322</td><td>1.05854</td><td>1.591269</td><td>1.591272</td></tr><tr><td>0.83347</td><td>1.28825</td><td>2.223859</td><td>2.223887</td></tr><tr><td>1.22324</td><td>0.38123</td><td>1.000142</td><td>2.000160</td></tr><tr><td>3.39219</td><td>1.13384</td><td>0.397339</td><td>0.397344</td></tr><tr><td>3.29681</td><td>1.23446</td><td>0.435928</td><td>0.435932</td></tr><tr><td>1.02283</td><td>1.78437</td><td>1.672796</td><td>1.672809</td></tr><tr><td>3.17710</td><td>1.44781</td><td>0.392305</td><td>0.392308</td></tr><tr><td>2.34871</td><td>1.71606</td><td>0.944684</td><td>0.944677</td></tr><tr><td>2.67003</td><td>1.72317</td><td>0.646585</td><td>0.646575</td></tr><tr><td>2.49107</td><td>1.91519</td><td>0.589281</td><td>0.589271</td></tr><tr><td>0.23415</td><td>2.47135</td><td>0.915709</td><td>0.915707</td></tr><tr><td>2.83975</td><td>1.93518</td><td>0.230071</td><td>0.230062</td></tr><tr><td>3.04116</td><td>1.91419</td><td>0.042710</td><td>0.042709</td></tr><tr><td>2.41676</td><td>2.24680</td><td>0.213264</td><td>0.213261</td></tr><tr><td>1.38596</td><td>2.72735</td><td>0.153922</td><td>0.153927</td></tr></tbody></table>

<table><thead><tr><th>n</th><th>H / a<sup>4</sup> (exacte)</th><th>H / a<sup>4</sup> (calculée)</th><th>n</th><th>u<sub>max</sub> / a<sup>2</sup> (exacte)</th><th>u<sub>max</sub> / a<sup>2</sup> (calculée)</th></tr></thead><tbody><tr><td>1</td><td>0.39269</td><td>0.39262</td><td>1</td><td>0.25000</td><td>0.24995</td></tr><tr><td>2</td><td>1.25663</td><td>1.25641</td><td>2</td><td>0.40000</td><td>0.39993</td></tr><tr><td>3</td><td>2.12057</td><td>2.12019</td><td>3</td><td>0.45000</td><td>0.44992</td></tr><tr><td>4</td><td>2.95679</td><td>2.95626</td><td>4</td><td>0.47058</td><td>0.47050</td></tr><tr><td>5</td><td>3.77595</td><td>3.77527</td><td>5</td><td>0.48076</td><td>0.48068</td></tr><tr><td>6</td><td>4.58502</td><td>4.58420</td><td>6</td><td>0.48648</td><td>0.48640</td></tr><tr><td>7</td><td>5.38783</td><td>5.38686</td><td>7</td><td>0.49000</td><td>0.48991</td></tr><tr><td>8</td><td>6.18652</td><td>6.18541</td><td>8</td><td>0.49230</td><td>0.49222</td></tr><tr><td>9</td><td>6.98238</td><td>6.98113</td><td>9</td><td>0.49390</td><td>0.49381</td></tr><tr><td>10</td><td>7.77621</td><td>7.77483</td><td>10</td><td>0.49504</td><td>0.49496</td></tr></tbody></table>

On remarque bien à partir des résultats obtenus dans ce tableau et la comparaison avec ceux de l'exacte, la convergence des résultats mais pour une discrétisation assez élevée. En modifiant la forme du domaine, le programme donne des résultats satisfaisants. Donc on peut dire et conclure que le programme peut être applicable pour des sections complexes et on peut estimer la convergence des résultats obtenus pour ce type de domaine malgré l'absence des résultats exacts. Les graphes suivants illustrent la
variation des paramètres H / a⁴, uₘₐₓ / a² et uₘョy / a² en fonction de n. On remarque bien les résultats du cercle à partir du cas de l'ellipse pour n = 1.

<table><thead><tr><th>n</th><th>u<sub>moy</sub> / a<sup>2</sup> (exacte)</th><th>u<sub>moy</sub> / a<sup>2</sup> (calculée)</th></tr></thead><tbody><tr><td>1</td><td>0.12500</td><td>0.12498</td></tr><tr><td>2</td><td>0.20000</td><td>0.19997</td></tr><tr><td>3</td><td>0.22500</td><td>0.22496</td></tr><tr><td>4</td><td>0.23529</td><td>0.23525</td></tr><tr><td>5</td><td>0.24038</td><td>0.24034</td></tr><tr><td>6</td><td>0.24324</td><td>0.24320</td></tr><tr><td>7</td><td>0.24500</td><td>0.24496</td></tr><tr><td>8</td><td>0.24615</td><td>0.24611</td></tr><tr><td>9</td><td>0.24695</td><td>0.24691</td></tr><tr><td>10</td><td>0.24752</td><td>0.24748</td></tr></tbody></table>

### 5.2.4 Domaine d'un profil d'aile

Le quatrième domaine à étudier est celui le cas d'un profil d'aile. L'exemple pris ici est de la série NACA 00 't', qui à la forme générale suivante [1]. C'est un profil symétrique donné par rapport au repère xoy par l'équation (85) [1]:

$$y(x) = t \times (a_1 \sqrt{x} + a_2 x + a_3 x^2 + a_4 x^3 + a_5 x^4) \quad 0 \le x \le 1 \quad (85)$$

$$a_1 = 1.4779155; a_2 = -0.624424; a_3 = -1.7478095; a_4 = 1.384087; a_5 = -0.489769$$

Si le profil est le NACA 0012, alors, l'épaisseur maximale est égale à t = 12 %, ou égale à 0.12. Elle est située à une distance de x = 0.2991671707 à partir du bord d'attaque.

<figure><img src="image_26.png" /><figcaption>Fig. 26: Profil d'aile type NACA 00 't'</figcaption></figure>

Pour l'intrados, on ajoute à l'équation, le signe moins (-) pour raison de symétrie. Ici la corde du profil est égale à l'unité. Malheureusement, pour ce domaine, on n'a pas trouvé dans la littérature disponible les résultats exactes concernant la distribution de la solution u(x,y) dans le domaine, ainsi que les paramètres H, uₘₐₓ / c² et uₘョy / c², et les valeurs de l'aire de la section et la position du centre de gravité du
profil. On peut noter ici que puisque le profil est symétrique, il faut trouver la position de son centre de gravité sur l'axe des abscisses. Autrement dit, $y_G = 0.0$.

Pour une discrétisation de 120 × 120 nœuds, on va déterminer pour un profil NACA 0012, les valeurs des paramètres en fonction de l'épaisseur 't'. Le tableau suivant résume les résultats obtenus. On note ici que la position du centre de gravité 'G' ne dépend pas de l'épaisseur et elle est située à une distance $x_G = 0.417287503$; $y_G = 0.0$.

<table><thead><tr><th>t</th><th>A</th><th>H</th><th>u<sub>max</sub> / c<sup>2</sup></th><th>u<sub>moy</sub> / c<sup>2</sup></th></tr></thead><tbody><tr><td>0.10</td><td>0.0677</td><td>0.00003</td><td>0.0012</td><td>0.0005</td></tr><tr><td>0.11</td><td>0.0744</td><td>0.00005</td><td>0.0014</td><td>0.0006</td></tr><tr><td>0.12</td><td>0.0812</td><td>0.00006</td><td>0.0017</td><td>0.0008</td></tr><tr><td>0.13</td><td>0.0880</td><td>0.00008</td><td>0.0020</td><td>0.0009</td></tr><tr><td>0.14</td><td>0.0947</td><td>0.00010</td><td>0.0023</td><td>0.0010</td></tr><tr><td>0.15</td><td>0.1015</td><td>0.00012</td><td>0.0026</td><td>0.0012</td></tr><tr><td>0.16</td><td>0.1083</td><td>0.00015</td><td>0.0030</td><td>0.0014</td></tr><tr><td>0.17</td><td>0.1150</td><td>0.00017</td><td>0.0034</td><td>0.0015</td></tr><tr><td>0.18</td><td>0.1218</td><td>0.00022</td><td>0.0037</td><td>0.0017</td></tr><tr><td>0.19</td><td>0.1286</td><td>0.00025</td><td>0.0042</td><td>0.0019</td></tr><tr><td>0.20</td><td>0.1354</td><td>0.00029</td><td>0.0046</td><td>0.0021</td></tr><tr><td>0.21</td><td>0.1421</td><td>0.00033</td><td>0.0050</td><td>0.0023</td></tr><tr><td>0.22</td><td>0.1489</td><td>0.00038</td><td>0.0055</td><td>0.0025</td></tr><tr><td>0.23</td><td>0.1557</td><td>0.00043</td><td>0.0059</td><td>0.0027</td></tr><tr><td>0.24</td><td>0.1624</td><td>0.00049</td><td>0.0064</td><td>0.0030</td></tr><tr><td>0.25</td><td>0.1692</td><td>0.00055</td><td>0.0069</td><td>0.0032</td></tr></tbody></table>

## 5.2.5 Divers domaines

Dans ce paragraphe, on va présenter en aperçu les différentes valeurs numériques obtenues sur les paramètres des modèles mathématiques de l'équation de Poisson pour une discrétisation de 300 par 300 nœuds. Les domaines considérés sont les figures 33, 40, 44, 46, demi cercle de rayon unité et quart de cercle de rayon unité.

### 5.2.5.1 Domaine de la figure 22

Les équations de l'extrados et de l'intrados dans ce domaine sont données par:

$$ 
\begin{aligned}
y_{\text{extrados}}(x) &= \sin((3\pi/2)x) + 1 \\
y_{\text{int rados}}(x) &= \cos((3\pi/2)x) \quad 0 \le x \le 1
\end{aligned} 
\qquad (86) $$

Les résultats des paramètres sont donnés sur le tableau suivant:

<table><tbody><tr><td>Position du point 'G' par rapport au repère xoy</td><td>x<sub>G</sub> = 0.5, y<sub>G</sub> = 0.5</td></tr><tr><td>A &amp; H</td><td>1.424409 - 0.030237</td></tr><tr><td>u<sub>max</sub> &amp; u<sub>moy</sub></td><td>0.039485 - 0.021228</td></tr><tr><td>Position (x,y) de u<sub>max</sub> par rapport au repère xoy</td><td>x = 0.28213<br>y = 1.29070</td></tr></tbody></table>

### 5.2.5.2 Domaine de la figure 23

L'équation de la frontière est choisie maintenant en coordonnées polaires par:
$$r(\theta) = 3 - \sin(3\theta + \cos 3\theta) \qquad 0 \le \theta \le 2\pi \qquad (87)$$

Les résultats des paramètres sont donnés sur le tableau suivant:

<table><tr><td>Position du point 'G' par rapport au repère xoy</td><td>x<sub>G</sub> = 0.0, y<sub>G</sub> = 0.0</td></tr><tr><td>A & H</td><td>22.102867 - 11.735912</td></tr><tr><td>u<sub>max</sub> & u<sub>moy</sub></td><td>1.254700 - 0.530967</td></tr><tr><td>Position (x,y) de u<sub>max</sub> par rapport au repère xoy</td><td>x = 1.425566<br>y = -0.391335</td></tr></table>

### 5.2.5.3 Domaine de la figure 24

L'équation de la frontière de la figure 24 est choisie maintenant en coordonnées polaires par:

$$r(\theta) = \text{Log}(\cos(\sin 4\theta \times \cos(15/5 \cdot \theta)) + 2) \qquad 0 \le \theta \le 2\pi \qquad (88)$$

Les résultats des paramètres sont donnés sur le tableau suivant:

<table><tr><td>Position du point 'G' par rapport au repère xoy</td><td>x<sub>G</sub> = -0.00633<br>y<sub>G</sub> = -0.00271</td></tr><tr><td>A & H</td><td>3.511725 - 0.461543</td></tr><tr><td>u<sub>max</sub> & u<sub>moy</sub></td><td>0.269628 - 0.131429</td></tr><tr><td>Position (x,y) de u<sub>max</sub> par rapport au repère xoy</td><td>x = 0.047670<br>y = -0.043356</td></tr></table>

### 5.2.5.4 Domaine d'un demi-cercle

Le domaine d'un demi-cercle de rayon R est choisi placé suivant l'axe vertical. Les résultats numériques des paramètres sont donnés dans le tableau suivant:

<table><tr><td>Position du point 'G' par rapport au repère xoy</td><td>x<sub>G</sub> / R = 0.0<br>y<sub>G</sub> / R = 0.42440</td></tr><tr><td>A & H</td><td>1.57076 - 0.07434</td></tr><tr><td>u<sub>max</sub> & u<sub>moy</sub></td><td>0.09758 - 0.04733</td></tr><tr><td>Position (x,y) de u<sub>max</sub> par rapport au repère xoy</td><td>x = 0.11788<br>y = 0.03551</td></tr></table>

### 5.2.5.5 Domaine d'un quart de cercle

Le domaine d'un quart de cercle de rayon R est choisi placé dans le premier quadrant. Les résultats numériques des paramètres sont donnés sur le tableau suivant:

<table><tr><td>Position du point 'G' par rapport au repère xoy</td><td>x<sub>G</sub> = 0.42441<br>y<sub>G</sub> = 0.42441</td></tr><tr><td>A & H</td><td>0.78539 - 0.02057</td></tr><tr><td>u<sub>max</sub> & u<sub>moy</sub></td><td>0.05508 - 0.02619</td></tr><tr><td>Position (x,y) de u<sub>max</sub> par rapport au repère xoy</td><td>x = 0.53936<br>y = 0.38612</td></tr></table>
## 6. CONCLUSION

Le présent travail nous a permis de faire la résolution de l'équation de Poisson dans un domaine simplement connexe bidimensionnel quelconque par la méthode des éléments finis en utilisant un élément triangulaire à trois nœuds avec les conditions aux limites de type Dirichlet, afin de déterminer les résultats des phénomènes physiques qui peuvent être présentés mathématiquement par l'équation de Poisson. Les conditions aux limites type Neumann et de Cauchy sont laissées comme des travaux de futurs. Le maillage développé dans ce conteste est de type 'H' avec des cellules triangulaires ou quadrilatères. Divers domaines ont été présentés dans ce travail touchant plusieurs disciplines d'une façon implicite. L'avantage de cette méthode est qu'elle peut modéliser l'importe quels domaines simplement connexes. Cependant, la comparaison des résultats afin de décider l'exactitude du programme élaboré a été faite pour les domaines circulaires, rectangulaires et elliptiques, dont la solution exacte existe.

La résolution du système d'équations obtenus est faite avec la technique des matrices bandes où le stockage de la partie bande supérieure des matrices est considéré dans un vecteur, vu que la matrice de rigidité est symétrique.

L'algorithmme utilisé est celui de Khaletski conçue spécialement pour la résolution directe des systèmes d'équations. Le transfert des cases entre la matrice et le vecteur correspondant est fait par la manière colonne par colonne. Le programme numérique réalisé peut résoudre un système d'équations d'ordre très élevé qui peut arriver à des millions ou même supérieures.

Un seul type d'élément fini est développé dans ce conteste. C'est le triangle à trois nœuds. Cependant on peut exploiter le programme numérique réalisé pour d'autres types d'éléments finis. Par exemple, l'élément triangulaire ou quadrilatère à plusieurs nœuds. Des changements nécessaires doivent être faits au programme.

**Remerciements-** Les auteurs remercient Khaoula, Abdelghani et Ritadj Zebbiche, ainsi que Mouza Ouahiba à la saisie de ce document.

## NOMENCLATURE

$u(x, y)$ : Solution de l'équation différentielle

$\Psi_i(x, y)$ : Fonctions d'interpolations

$(x_j, y_j)$ : Coordonnées d'un nœud

$v$ : Fonction test

$\Omega$ : Domaine plan

$K_{ij}$ : Eléments de la matrice de rigidité

$(f_j)^{(e)}$ : Vecteur force pour de élément fini

$c_1, c_2, c_3$ : Constantes pour l'approximation de la fonction $u(x, y)$

$F_i$ : Eléments du vecteur force

$I_{m,n}$ : Intégrales dans un triangle

$(m = 1, 2), (n = 1, 2)$

$f$ : Deuxième membre de l'équation de Poisson

$u_j$ : Valeur de $u$ au point $(x_j, y_j)$

$\delta_{ij}$ : Symbole de Kronecker

$\Gamma$ : Frontière du domaine $\Omega$

$\Omega^e$ : Domaine de l'élément fini

$[K]$ : Matrice de rigidité

$(f)$ : Vecteur force totale du domaine

$\alpha_i, \beta_i, \gamma_i$ : Coefficients des fonctions d'interpolations d'un triangle $(i=1, 2, 3)$

$I_{100}, I_{10}, I_{01}$ : Moments statiques

$A_e, A_i$ : Aire de la surface d'un triangle

$(m = 1, \dots, NE)$
$(x_G)_i, (y_G)_i$: Coordonnées du centre de gravité du triangle numéro 'i'

IBAV: Valeur de la bande avant l'application des conditions aux limites

B: Demi-bande de la matrice de rigidité

NN: Nombre total de nœuds dans le maillage

$N_x, N_y$: Nombre de nœud suivant les axes x et y

$\eta^*$: Variable normalisée

P, Q: Paramètres pour le contrôle des

points de maillage (Fonct. de condensation)

s: Variable intermédiaire de la fonction de condensation

H: Intégrale de la solution u du problème

$\{vk\}$: Vecteur remplaçant la matrice de rigidité

A, B, C, D: Points spécifiés sur la frontière du domaine

(e): Élément fini

nouveau: Nouvelle position d'un nœud

$B_L$: Demi-bande donnée par l'élément numéro L

IBAP: Valeur de la demi-bande après l'application des conditions aux limites

NE: Nombre total des éléments finis

NDDLE: Nombre de degré de liberté éliminé

NDDLR: Nombre de degré de liberté restant

NDCE: Nombre de case éliminée

$\omega$: Facteur de relaxation utilisé dans

l'algorithme de filtrage des nœuds $0 < \omega < 1$

m: Nombre de nœuds reliés par une arête d'un triangle avec le nœud considéré

i, j, k: Numéros des nœuds - élément fini

$[L]$ : Matrice triangulaire déterminée par la décomposition de la matrice de rigidité $[K]$

T: Transposée d'un vecteur ou matrice

t: Epaisseur de profil d'aile

i, j: Compteur sur les nœuds

ancien: Ancienne position d'un nœud

## REFERENCES

[1] C.A.J. Fletcher, 'Computational Techniques for Fluid Dynamics', Volume II, Specific Techniques for Different Flow Categories, Springer-Verlag, 1988, Austria.

[2] J.F. Imbert, 'Analyse des Structures par Eléments Finis', 3ème Edition, Sup'Aéro, Cepadues-Editions, 1991, France.

[3] J.Jr. Anderson, 'Modern Compressible Flow: With Historical Perspective', McGraw Hill Book company, New York, 1982.

[4] W.F. Hughes and E.W. Gaylord, 'Basic Equations of Engineering Science', Schaums Outline serie, Mc Graw-Hill Book Company, 1964.

[5] O.A. Ladyzenkaja et N.N. Ural'ceva, 'Equations aux Dérivées Partielles de Type Elliptique', Edition Dunod, Paris, 1988.

[6] S. Miranda, 'Partial Differential Equations of Elliptic Type', 2nd Edition, Springer, Berlin, 1970.

[7] W.G. Habashi, 'A Finite Element Approach to Subsonic Aerodynamics', International Journal of Numerical Methods in Engineering, Vol. 14, N°5, pp. 665 - 679, 1978.

[8] H. Kardestuncer, 'Finite Element Handbook', McGraw-Hill Book Company, USA, 2010.

[9] J.N. Reddy, 'An Introduction to the Finite Element Method', Edition Mc-Graw Hill Book Company, 2007.

[10] A.K. Aziz, 'The Mathematical Foundations of the Finite Element Method with Application to Partial Differential Equations', Academic Press, New York, 1989.

[11] A. Khenane, 'Méthodes des Eléments Finis, Enoncé des Principes de Base', Edition OPU, Algérie, 2007.
[12] K.J. Bathe, *Finite Element Procedures in Engineering Analysis*, Prentice-Hall, Englewood Cliffs, New Jersey, 1982.

[13] B. Demidovitch et I. Marron, 'Eléments de Calcul Numérique', Edition Mir, Moscou, 1987.

[14] A. Raltson and A. Rabinowitz, *A First Course in Numerical Analysis*, McGraw Hill Book Company, 2005.

[15] B. Irons and S. Ahmed, 'Techniques of Finite Elements', John Wiley and Sons, New York, 1990.

[16] D. Gouri et G. Touzot, 'Une Présentation de la Méthode des Eléments Finis', Université Laval, Moline S.A. Editeur, Paris, 1999.

[17] H.R. Schwarz and H. Rutishauser, 'Numerical Analysis of Symmetric Matrices', English Edition, Prentice-Hall, Englewood Cliffs, New Jersey, 2006.

[18] G. Forsythe and C.B. Moler, 'Computer Solution of Linear Algebraic Systems', Prentice-Hall, Inc., Englewood Cliffs, N.J., 1993.

[19] H. Kardestuncer, 'Elementary Matrix Analysis of Structures', McGraw-Hill, New York, 1999.

[20] J.R. Bunch and B.N. Parlett, 'Direct Methods for Solving Symmetric Indefinite Systems of Linear Equations', SIAM Journal of Numerical Analysis, Vol. 8, N°4, pp. 639 - 655, 1982.

[21] G.W. Stewart, 'Introduction to Matrix Computations', Academic Press: New York, 1987.

[22] K.J Bathe and E.L. Wilson, 'Numerical Methods in Finite Element Analysis', Prentice Hall, 2001.

[23] J.F. Thompson, S. Bharat and N.P. Weatherill, 'Handbook of Grid Generation', CRC Press, USA, 2008.

[24] E. Sadek, 'A Scheme for the Automatic Generation of Triangular Finite Elements', International Journal for Numerical Methods in Engineering, Vol. 15, N°12, pp. 1813 - 1822, 1994.

[25] N.P. Weatherill, 'A Method for Generating Irregular Computational Grids in Multy Connected Planar Domain', International Journal of Numerical Methods in Fluids, Vol. 8, pp. 181 - 197, 1993.

[26] J.H. Bramble and M. Zlamal, 'Triangular Elements in the Finite Element Methods', Mathematics of Computation, Vol. 24, N°4, pp. 809 - 820, 1970.

[27] B.H.V. Topping, J. Muylle, P. Ivanyi, R. Putanowicz and B. Cheng, 'Finite Element Mesh Generation', Saxe Coburg Publications, Hardback, USA, 2004.

[28] R. Rosen, 'Matrix Bandwidth Minimization', Proceedings of the 23rd National Conference, Association for Computing Machinery, Brandon Systems Press, Princeton, N.J., 1989.

[29] G.G. Always and D.W. Martin, 'An Algorithm for Reducing the Bandwidth of a Matrix of Symmetric Configuration', The Computer Journal, Vol. 8, pp. 264 - 272, 1988.

[30] R.K. Brayton, F.G. Gustavson and R.A. Willoughby, 'Some Results on Sparse Matrices', of Computation, Vol. 24, N°115, pp. 937 - 954, 1970.

[31] K.Y. Cheng, 'Minimizing the Bandwidth of Sparse Symmetric Matrices', Computing, Vol. 11, pp. 103 - 110, 1985.

[32] R.G. Henry, 'Algorithm for Matrix Bandwidth Reduction', Proceedings of the American Society of Civil Engineering, pp. 203-214, 1989.

[33] A.H. Stroud and D.H. Secrest, 'Gaussian Quadrature Formulas', Prentice-Hall, Inc., Englewood Cliffs, N.J., 1966.

[34] E.L. Wilson, K.J. Bathe and W.P. Doherty, 'Direct Solution of Large Systems of Linear Equations', Computers and Structures, Vol. 4, N°2, pp. 363 - 372, 1984.

[35] M. Boun-djad et T. Zebbiche, 'Résolution de l'Équation de Poisson dans un Domaine simplement Connexe par la Méthode des Eléments Finis', Editions Universitaires Européennes, Avril 2012.