# Application d'une méthode d'éléments finis d'ordre un à la résolution numérique des équations de Navier-Stokes

François Thomasset

► To cite this version:

François Thomasset. Application d'une méthode d'éléments finis d'ordre un à la résolution numérique des équations de Navier-Stokes. [Rapport de recherche] IRIA-RR-150, IRIA. 1975, pp.70. hal-04716071

HAL Id: hal-04716071

https://inria.hal.science/hal-04716071v1

Submitted on 1 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
RIA

laboria

laboratoire de recherche
en informatique
et automatique

Application d'une méthode
d'éléments finis d'ordre un
à la résolution numérique
des équations
de Navier-Stokes

F. Thomasset

<figure><img src="image_2.png" /></figure>

Rapport de Recherche n°150
Décembre 1975

Institut de Recherche
d'Informatique
et d'Automatique

Domaine de Voluceau
Rocquencourt
B. P. 105 78150 - Le Chesnay
France
Tél.: 954 90 20
APPLICATION D'UNE MÉTHODE D'ÉLÉMENTS FINIS
D'ORDRE UN A LA RÉSOLUTION NUMÉRIQUE DES
ÉQUATIONS DE NAVIER-STOKES

F. Thomasset

présenté par R. TEMAM

RÉSUMÉ :

On propose dans cet article l'étude numérique des équations stationnaires de Navier-Stokes en deux dimensions d'espace, par éléments finis non conformes de degré 1.

L'équation de conservation de la masse div $\vec{u} = 0$ est traitée par pénalisation-dualité.

ABSTRACT :

We study in this paper the numerical solution of stationary (2 D) Navier-Stokes equations.
Non conforming finite elements of degree one are used, and the mass conservation equation is
treated by a penalty-duality method.
# 1. - INTRODUCTION - FORMULATION VARIATIONNELLE DU PROBLEME CONTINU

On présente dans ce travail la mise en oeuvre numérique des éléments finis linéaires non conformes, appliqués à l'approximation des écoulements visqueux, incompressibles, en dimension 2 d'espace.

Soit $\Omega \subset \mathbb{R}^2$ le domaine contenant le fluide, de frontière $\partial\Omega$.

Dans le cas stationnaire, l'écoulement d'un fluide de viscosité $v$ est donné par les équations de Navier-Stokes :

$$ (1.1) \quad \left\| \begin{array}{l} - v \Delta \vec{u} + \sum_{j=1}^{2} u_j \frac{\partial \vec{u}}{\partial x_j} + \vec{\nabla} p = \vec{f} \quad \text{dans } \Omega \\ \\ \vec{\nabla} \cdot \vec{u} = 0 \quad \text{dans } \Omega \\ \\ \vec{u} \mid_{\partial \Omega} = \vec{g} \\ \\ \text{Inconnues} \quad \begin{array}{l} \vec{u} = \text{champ des vitesses.} \\ p = \text{pression.} \end{array} \\ \\ \text{Données} \quad \begin{array}{l} v = \text{viscosité.} \\ \vec{f} = \text{forces volumiques.} \\ \vec{g} = \text{conditions aux limites.} \end{array} \end{array} \right. $$

La contrainte $\vec{\nabla} \cdot \vec{u} = 0$ est la condition d'incompressibilité.

Introduisons quelques notations :

$$ W = H^1(\Omega) = \left\{ v \in L^2(\Omega), \frac{\partial v}{\partial u_j} \in L^2(\Omega), j = 1, 2 \right\} $$

$$ W_O = H_O^1(\Omega) = \left\{ v \in W, v|_{\partial\Omega} = 0 \right\} $$

$$ |W| = |H^1(\Omega)| = W \times W $$

$$ |W_O| = |H_O^1(\Omega)| = W_O \times W_O $$
$$
W = \left\{ \begin{array}{ll} \vec{v} \in W_0 : & \vec{v} \cdot \vec{v} > 0 \\ & \text{p. dans } \Omega \end{array} \right\}
$$

$\nabla \vec{u}, \vec{v}, \vec{w} \in W$, on pose :

$$
a(\vec{u}, \vec{v}) = \sum_{i,j=1}^{2} \int_{\Omega} \frac{\partial u_i}{\partial x_i} \frac{\partial v_j}{\partial x_j} dx_1 dx_2
$$

$$
b(\vec{u}, \vec{v}, \vec{w}) = \sum_{i,j=1}^{2} \int_{\Omega} u_i \frac{\partial v_j}{\partial x_i} w_j dx_1 dx_2
$$

$$
(\vec{f}, \vec{v}) = \int_{\Omega} \vec{f} \cdot \vec{v} \, dx_1 \, dx_2
$$

ceci étant posé, $\vec{u}$ est solution faible du problème (1.2.) si $\vec{u}$ vérifie :

$$
\begin{align*}
& \text{(1.2.)} \quad \left|
\begin{array}{@{}l@{}}
  \begin{aligned}[t]
    \vec{u} &\in W \\
    v \, a(\vec{u}, \vec{v}) + b(\vec{u}, \vec{u}, \vec{v}) &=(\vec{f}, \vec{v}), \quad \forall \vec{v} \in V
  \end{aligned} \\
  \vec{v} \cdot \vec{u} = 0 \text{ dans } \Omega \\
  \vec{u}\big|_{\partial\Omega} = \vec{g}
\end{array}
\right.
\end{align*}
$$

ce qui équivaut à l'existence de $(\vec{u}, p)$ tel que :

$$
\text{(1.3.)} \quad \left|
\begin{array}{@{}l@{}}
  \begin{aligned}[t]
    \vec{u} &\in W, p \in L^2(\Omega) \\
    v \, a(\vec{u}, \vec{v}) + b(\vec{u}, \vec{u}, \vec{v}) &\quad - (p, \vec{v} \cdot \vec{v}) = (\vec{f}, \vec{v}) \quad \forall \vec{v} \in W_0 \\
    \vec{v} \cdot \vec{u} &\quad = 0 \text{ dans } \Omega
  \end{aligned} \\
  \vec{u}\big|_{\partial\Omega} = \vec{g}
\end{array}
\right.
$$

(La pression $p$ est le multiplicateur de Lagrange associé à la contrainte $\vec{v} \cdot \vec{u} = 0$).
Remarque :

Si $\vec{v} \cdot \vec{u} = 0$, et $\vec{v}|_{\partial\Omega} = 0$ ou $\vec{w}|_{\partial\Omega} = 0$, on a immédiatement par la formule de Green :

$$ (1.4.) \quad b(\vec{u}, \vec{v}, \vec{w}) + b(\vec{u}, \vec{w}, \vec{v}) = 0 $$

On peut donc remplacer, dans (1.2.) ou (1.3.), la forme trilinéaire $\hat{b}(\vec{u}, \vec{v}, \vec{w})$ par : $\hat{b}(\vec{u}, \vec{v}, \vec{w}) = \frac{1}{2} (b(\vec{u}, \vec{v}, \vec{w}) - b(\vec{u}, \vec{w}, \vec{v}))$.

La forme $\hat{b}$ vérifie (1.4.) même lorsque $\vec{v} \cdot \vec{u} \neq 0$.

Ce procédé, dû à TEMAM [2.4], améliore la convergence des algorithmes numériques.

## 2. - DEFINITION DE L'APPROXIMATION.

### 2.1. - Choix de l'approximation

Soit $\Omega_h$ un domaine polygônal approchant $\Omega_i$, il existe un ensemble fini $\mathcal{T}_h$ de triangles deux à deux disjoints recouvrant $\Omega_h$ :

$$ i) \quad U T = \bar{\Omega}_h \\ T \in \mathcal{T}_h $$

$$ ii) \quad T \in \mathcal{T}_h, T' \in \mathcal{T}_h, T \neq T' \quad \Rightarrow \quad T \cap T' = \emptyset $$

iii) Pour deux triangles distincts $T$ et $T' \in \mathcal{T}_h$, on a :

- ou bien $T \cap T' = \text{un côté}$;

- ou bien $T \cap T' = \text{un sommet}$;

- ou bien $T \cap T' = \emptyset$.
<figure><img src="image_3.png" /></figure>

Notations :

M = nombre de triangles.

$N_I$ = nombre de sommets intérieurs.

$N_F$ = nombre de sommets frontière.

$C_I$ = nombre de côtés intérieurs.

$C_F$ = nombre de côtés frontière = $N_F$.

$N_T$ = nombre de trous dans le domaine.

On a les relations classiques suivantes :

$$ (2.1.) \qquad M = \frac{1}{3}(C_F + 2C_I) = N_F + 2N_I - 2(N_T - 1) $$

(Cf. EWING étal. [6])

Nous allons maintenant décrire les espaces $W^h$, $W_o^h$, $W^h$, $W_o^h$, $W^h$, approximations de $W$, $W_o$, $W$, $W_o$, $W$ en justifiant notre choix.
Pour $\vec{U}_h \in W_o^h$, on veut :

i) sur chaque $T \in \mathcal{T}_h$, pour chaque composante $U_{hj} (j = 1, 2)$ de $\vec{U}_h$ :
$U_{hj}|_T$ est un **polynôme de degré 1** ;

ii) $\vec{U}_h|_{\partial\Omega_h} = 0$ en un certain sens ;

iii) à la frontière entre les éléments $T \in \mathcal{T}_h$, on doit imposer des conditions de raccord que nous allons préciser.

La démarche la plus naturelle serait d'imposer le raccord de $\vec{U}_h$ aux sommets des triangles ; on aurait alors $(\vec{U}_h \in W_o^h)$.

$\vec{U}_h$ de classe $C^0$ sur $\bar{\Omega}_h$
$\vec{U}_h \equiv 0$ partout sur $\partial\Omega_h$

et $\vec{U}_h \in W_o^h$ serait défini par $2 \times N_I$ degrés de libertés

Si l'on veut imposer la condition d'incompressibilité discrète :

$$ (2.2.) \quad \vec{\nabla} \cdot (\vec{U}_h)_T = 0 \quad \forall T \in \mathcal{T}_h, $$

On aura autant de contraintes que de triangles. Mais les relations (2.1.) montrent que l'on a alors plus de contraintes que de degrés de liberté :

$$ M - (2 \times N_I) = N_F - 2 \times (N_T - 1) > 0 $$
En général les relations (2.2.) ne seront pas linéairement dépendantes et l'on a donc, pour des éléments finis linéaires conformes :

$$
\mathcal{W}^h = \left\{
\begin{array}{@{}l@{}}
  \begin{array}{@{}l@{}}
    \vec{v}_h : \vec{v}_h|_T \text{ polynôme de degré } 1, \\
    \phantom{\vec{v}_h : \vec{v}_h|_T } \vec{v}_h \text{ continue sur } \bar{\Omega}_h,
  \end{array} \\
  \vec{v} \cdot \vec{v}_h|_T = 0 \ \forall T
\end{array}
\right\} = \{\vec{0},\}
$$

et l'espace des fonctions à divergence discrète nulle est réduit à $\vec{0}$ :
cette remarque est due à M. FORTIN (cf [7] , où sont étudiées des
triangulations pour lesquelles $\mathcal{W}^h \neq \{\vec{0}\}$ avec éléments finis linéaires
conformes).

Afin d'obtenir $\mathcal{W}^h \neq \{\vec{0}\}$, on renonce à une continuité sur tout $\Omega_h$ pour n'imposer le raccord qu'aux milieux des côtés, obtenant ainsi des éléments non conformes.

On aura alors $(2 \times C_I)$ degrés de liberté et $2 \times C_I = 3M - C_F > M$.

2.2. - *Définition précise de l'approximation*

On appelle noeuds de la triangulation les milieux des côtés des
triangles. Soit $\mathcal{W}^h$ l'ensemble des fonctions scalaires $\varphi_h$ :

- linéaires sur chaque $T \in \mathcal{T}_h$ ;

- se raccordant aux noeuds , i.e. si T et T' sont deux triangles
  adjacents, M leur noeud commun :

$$
(\varphi_h|_{T})(M) = (\varphi_h|_{T'})(M)
$$

<figure><img src="image_4.png" /></figure>
$W_o^h$ est le sous-espace des $\varphi_h \in W^h$ s'annulant en chaque noeud de la frontière $\partial\Omega_h$.

$$W^h = W^h \times W^h, \quad W_o^h = W_o^h \times W_o^h$$

Posons maintenant :

$$X^h = \left\{ \begin{array}{l} \pi_h : \forall T \in \mathcal{T}_h, \pi_h(x_1, x_2) \equiv \pi_h(T) \\ \text{pour } (x_1, x_2) \in \mathcal{T} \end{array} \right\}$$

= espace vectoriel des fonctions scalaires constantes par triangle.

Pour $\varphi_h \in W^h$, on note $(\delta_{jh} \varphi_h)(j = 1, 2)$ la fonction de $X^h$ valant sur chaque $T$ :

$$\delta_{jh} \varphi_h(T) = \frac{\partial}{\partial x_j} (\varphi_h|_T)$$

$$\vec{\nabla}_h \varphi_h = \begin{bmatrix} \delta_{1h} \varphi_h \\ \delta_{2h} \varphi_h \end{bmatrix}$$

$$(\vec{\nabla}_h \in W^h), \vec{\nabla}_h \cdot \vec{\nabla}_h = \delta_{1h} v_{1h} + \delta_{2h} v_{2h}$$

$$(2.3.) \quad V^h = \{ \vec{v}_h \in W_o^h : \vec{v}_h \cdot \vec{v}_h = 0 \}$$

Considérons un ensemble de triangulations $\mathcal{T}_h(h \in \mathcal{H})$ de $\Omega$ vérifiant la propriété :

- Si $\theta_h$ est le plus petit des angles de la triangulation,
$\theta_h \ge \theta_o > 0$ $\theta_o$ indépendant de $h$ (on dit que l'on a une famille **régulière** de triangulations).

Pour chaque $\mathcal{C}_h$,

$h = \sup$ des côtés des triangles de $\mathcal{C}_h$.

Alors (pour $h \to 0$), $W^h$ (resp. $W_o^h$, $W^h$, $W_o^h$, $W^h$) permet de définir une approximation stable et convergente de $W$ (resp. $W_o$, $W$, $W_o$, $W$). (Cf. TEMAM, Cours du Maryland [14], CROUZEIX et RAVIART [5], JAMET et RAVIART [12], pour la définition des opérateurs $p_h$ et $r_h$ et les démonstrations de stabilité et de convergence).

**Remarque :** $W^h, W_o^h, W^h, W_o^h, W^h$ sont de dimension finie :

$$ \dim W^h = 2(C_I + C_F), \text{ etc...} $$

### 2.3. - **Problème discret**

$$ a_h(\vec{u}_h, \vec{v}_h) = \sum_{i,j=1}^{2} \int_{\Omega_h} (\delta_{ih} u_{jh})(\delta_{ih} v_{jh}) \, dx_1 \, dx_2 $$

$$ b_h(\vec{u}_h, \vec{v}_h, \vec{w}_h) = \sum_{i,j=1}^{2} \int_{\Omega_h} u_{ih}(\delta_{ih} v_{jh}) w_{jh} \, dx_1 \, dx_2 $$

$$ \hat{b}_h(\vec{u}_h, \vec{v}_h, \vec{w}_h) = \frac{1}{2} (b_h(\vec{u}_h, \vec{v}_h, \vec{w}_h) - b_h(\vec{u}_h, \vec{w}_h, \vec{v}_h)) \\ (\vec{u}_h, \vec{v}_h, \vec{w}_h \in W^h) $$

$$ (\hat{f}, \vec{v}_h)_h = \sum_{i=1}^{2} \int_{\Omega_h} f_i v_{ih} \, dx_1 \, dx_2 $$
L'approximation discrète de (1.2.) s'écrit :

$$
\begin{array}{l@{\qquad}l}
\begin{array}{@{}l@{}}
  \text{(2.4.)} \\
  \quad \vec{u}_h \in W^h \\
  \quad \lor a_h(\vec{u}_h, \vec{v}_h) + \hat{b}_h(\vec{u}_h, \vec{u}_h, \vec{v}_h) = (\vec{f}, \vec{v}_h)_h \quad \forall \vec{v}_h \in V^h \\
  \quad \vec{v}_h \cdot \vec{u}_h = 0 \\
  \quad \vec{u}(M) = \vec{g}_h(M) \quad \text{pour tout noeud frontière M}
\end{array}
&
\begin{array}{@{}l@{}}
  \vec{u}_h \in W^h \\
  \lor a_h(\vec{u}_h, \vec{v}_h) + \hat{b}_h(\vec{u}_h, \vec{u}_h, \vec{v}_h) = (\vec{f}, \vec{v}_h)_h \quad \forall \vec{v}_h \in W^h \\
  \lor a_h(\vec{u}_h, \vec{v}_h) + \hat{b}_h(\vec{u}_h, \vec{u}_h, \vec{v}_h) = (\vec{f}, \vec{v}_h)_h \quad \forall \vec{v}_h \in V^h \\
  \vec{v}_h \cdot \vec{u}_h = 0 \\
  \vec{u}(M) = \vec{g}_h(M) \quad \text{pour tout noeud frontière M}
\end{array}
\end{array}
$$

où $\vec{g}_h$ est une approximation de $\vec{g}$ (précisée plus loin).

On sait (Cf TEMAM [14] ) que (2.4.) possède une solution $\vec{u}_h$
unique si $v$ n'est pas trop petit, et que $\vec{u}_h \to \vec{u}$ solution de (1.2.) si l'on
prend une famille régulière de triangulations.

Si l'on connaît une base de $V^h$, on peut essayer de résoudre
directement le système d'équations non-linéaires (2.4.). (Cf. CROUZEIX
(non publié) pour l'explicitation d'une telle base), mais on n'utilise
pas ce point de vue ici.

Considérons la forme linéaire sur $|W_o^h|$ définie par :

$$
\ell_h(\vec{v}_h) = \nu a_h(\vec{u}_h, \vec{v}_h) + \hat{b}_h(\vec{u}_h, \vec{u}_h, \vec{v}_h) - (\vec{f}, \vec{v}_h)_h
$$

$\ell_h$ s'annule sur le sous-espace $V^h$ défini par :

$$
\vec{v} \cdot (\vec{v}_h |_{T}) = 0 \quad \forall T \in \mathcal{T}_h
$$
Il existe donc des multiplicateurs $\lambda_T(T \in \mathcal{C}_h)$ tels que :

$$ (2.5.) \qquad \ell_h(\vec{v}_h) = \sum_{T \in \mathcal{C}_h} \lambda_T(\vec{v} \cdot (\vec{v}_h|_T)) \forall \vec{v}_h \in W_o^h $$

Introduisons la fonction $p_h \in X^h$ :

$$ p_h (T) = \frac{\lambda_T}{\text{aire}(T)} $$

alors :

$$ \begin{align*} \sum_T \lambda_T (\vec{v} \cdot (\vec{u}_h^\top |_T)) &= \int_{\Omega_h} p_h (\vec{v}_h^\top \cdot \vec{v}_h^\top) \, dx_1 \, dx_2 \\ &= (p_h, \vec{v}_h^\top \cdot \vec{v}_h^\top) \end{align*} $$

et (2.5.) s'écrit :

$$ (2.6.) \qquad \left| \begin{array}{l} \vec{u}_h \in W_o^h \\ \lor a_h(\vec{u}_h, \vec{v}_h) + \hat{b}_h(\vec{u}_h, \vec{u}_h, \vec{v}_h) - (p_h, \vec{v}_h \cdot \vec{v}_h) = (\vec{f}, \vec{v}_h)_h \forall \vec{v}_h \in W_o^h \\ \vec{v}_h \cdot \vec{u}_h = 0 \\ \vec{u}_h(M) = \vec{g}_h(M) \text{ pour tout } M \text{ noeud frontière} \end{array} \right. $$

(2.6.) est l'analogue discret de (1.3.) ; $p_h$ est la pression discrète : on peut montrer que $p_h \to p$ dans $L^2(\Omega) / \mathbb{R}$ quand $\ell \to 0$.
### 3. - DESCRIPTION DES ALGORITHMES NUMERIQUES

Les algorithmes utilisés sont dérivés d'algorithmes de recherche de point-selle ; on a tout d'abord examiné le problème de Stolses, qui peut être considéré comme un vrai problème de point-selle.

#### 3.1. - Approximation du problème de Stolses homogène

$$ (3.1.) \quad \left| \begin{array}{l} - v \Delta \vec{u} + \vec{\nabla} \cdot \vec{p} = \vec{f} \\ \vec{\nabla} \cdot \vec{u} = 0 \\ \vec{u}|_{\partial\Omega} = \vec{0} \end{array} \right. $$

La formulation variationnelle de (3.1.) s'écrit :

$$ (3.2.) \quad \left| \begin{array}{l} \text{Trouver } \vec{u} \in W_o^h \, , \ p \in L^2(\Omega) \text{ tels que :} \\ \qquad v \ a(\vec{u}, \vec{v}) - (p, \vec{v} \cdot \vec{v}) = (\vec{f}, \vec{v}) \ \forall \vec{v} \in W_o \\ \qquad \vec{\nabla} \cdot \vec{u} = 0 \end{array} \right. $$

L'approximation de (3.2.) est :

$$ (3.3.) \quad \left| \begin{array}{l} u_h \in W_o^h \, , \ p_h \in X^h \\ \qquad v \ a_h(\vec{u}_h, \vec{v}_h) - (p_h, \vec{v}_h \cdot \vec{v}_h)_h = (\vec{f}_h, \vec{v}_h) \ \forall \vec{v}_h \in W_o^h \\ \qquad \vec{\nabla}_h \cdot \vec{u}_h = 0 \end{array} \right. $$

Introduisons le lagrangien :

$$ (3.4.) \quad \mathcal{L}(\vec{v}_h, q_h) = \frac{1}{2} a_h(\vec{v}_h, \vec{v}_h) - (q_h, \vec{v}_h \cdot \vec{v}_h) - (\vec{f}_h, \vec{v}_h)_h (\vec{v}_h \in W_o^h, q_h \in X^h). $$
(3.3.) est équivalent à la recherche de $\vec{u}_h \in W_o^h$, $p_h \in X^h$ :

$$ (3.5.) \quad \left| \begin{array}{l} \mathcal{L}(\vec{u}_h, q_h) \le \mathcal{L}(\vec{u}_h, p_h) \le \mathcal{L}(\vec{v}_h, p_h) \\ \forall \vec{v}_h \in W_o^h, q_h \in X^h \end{array} \right. $$

**Algorithme d'Uzawa :**

$$ \vec{u}_h^{(o)} \in W_o^h \ ; \ p_h^{(o)} \in X^h $$

$p_h^{(m)}$ étant connu, on définit $\vec{u}_h^{(m+1)}$ par :

$$ (3.6.) \quad \mathcal{L}(\vec{u}_h^{(m+1)}, p_h^{(m)}) = \text{Min} \ \mathcal{L}(\vec{v}_h, p_h^{(m)}) \ \vec{v}_h \in W_o^h $$

Puis on modifie la pression suivant l'équation :

$$ (3.7.) \quad p_h^{(m+1)} = p_h^{(m)} - \rho \vec{v}_h \cdot \vec{u}_h^{(m+1)} $$

N.B. : le problème de minimisation (3.7.) se ramène à l'équation :

$$ (3.8.) \quad \nu a_h(\vec{u}_h^{(m+1)}, \vec{v}_h) = (\vec{f}, \vec{v}_h)_h + (p_h, \vec{v}_h, \vec{v}_h)_h \quad \forall \vec{v}_h \in W_o^h $$

**Algorithme d'Arrow-Hurwicz**

$$ \vec{u}_h^{(o)} \in W_o^h, \ p_h^{(o)} \in X^h $$

- quand $p_h^{(m)}$ est connu, au lieu de minimiser $\mathcal{L}(\vec{v}_h, p_h^{(m)})$

sur tout l'espace, on fait un seul pas de gradient :
$$
(3.9.) \quad \left\langle (\vec{u}_h^{(m+1)} - \vec{u}_h^{(m)}) + \rho \frac{\partial \mathcal{L}}{\partial \vec{v}_h} (\vec{u}_h^{(m)}, p_h^{(m)}) \right\rangle, \vec{v}_h > = 0 \quad (\vec{v}_h \in W_o^h)
$$

<,> est un produit scalaire sur W_o_h; par exemple, avec

<\vec{u}_h, \vec{v}_h> = a_h(\vec{u}_h, \vec{v}_h), (3.9.) s'écrit :

$$
\begin{equation}
\begin{split}
& (3.10.) \quad a_h(\vec{u}_h^{(m+1)} - \vec{u}_h^{(m)}, \vec{v}_h) + \rho(v a_h(\vec{u}_h^{(m)}, \vec{v}) - (\vec{f}, \vec{v}_h)_h)_h \\
& \qquad \qquad \qquad \qquad \qquad \qquad - (p_h, \vec{v}_h \cdot \vec{v}_h)_h) = 0 \quad \forall \vec{v}_h \in W_o^h
\end{split}
\end{equation}
$$

(Cf. TEMAM [14] . Ce schéma est donc implicite).

Connaissant $\vec{u}_h^{(m+1)}$, on modifie $p_h$ :

$$
(3.11.) \quad \alpha(p_h^{(m+1)} - p_h^{(m)}) = - \rho \vec{v}_h \cdot \vec{u}_h^{(m+1)}
$$

Remarque : Il s'avère avantageux d'utiliser les algorithmes
précédents en ajoutant au lagrangien L (v_h, q_h) un terme
de pénalisation :

$$
(3.12) \quad \mathcal{L}_{\varepsilon} (\vec{v}_h, q_h) = \mathcal{L} (\vec{v}_h, q_h) + \frac{1}{\varepsilon} (\vec{v}_h \cdot \vec{v}_h^*, \vec{v}_h^* \cdot \vec{v}_h^*)
$$

C'est la méthode des lagrangiens augmentés de Rockafellar
et Hestenes (voir [13] et la bibliographie de cet article).

(Cf. GLOWINSKI - FORTIN-MAROCCO [10], pour diverses applications
de cette méthode aux problèmes de mécanique).

Au lieu de (3.8.) on a alors l'équation :

$$
\begin{equation}
\begin{split}
& v a_h (\vec{u}_h^{(m+1)}, \vec{v}_h) + \frac{1}{\epsilon} (\vec{v}_h \cdot \vec{u}_h^{(m+1)}, \vec{v}_h \cdot \vec{v}_h) \\
& \qquad = (\vec{f}, \vec{v}_h)_h + (p_h^{(m)}, \vec{v}_h \cdot \vec{v}_h).
\end{split}
\tag{3.13}
\end{equation}
$$
### 3.2. - Quelques algorithmes pour les équations de Navier - Stok Stationnaire.

Les algorithmes utilisés sont dérivés des méthodes de résolution pour Stolses décrites plus haut.

**Algorithme (I) (UZAWA ; Cf. FORTIN - TEMAM - PEYRET) [8]**

$$ (3.14) \qquad \begin{array}{l} \vec{u}_h^{(o)} \in W^h, \quad \vec{u}_h^{(o)}|_{\partial\Omega_h} = \vec{g}_h \quad \vec{u}_h^{(o)}, p_h^{(o)} \text{ quelconques} \\ \qquad \qquad \qquad \qquad \text{vérifiant ces conditions.} \\ p_h^{(o)} \in X_h^h \end{array} $$

Connaissant $p_h^{(m)}$, on définit $\vec{u}_h^{(m+1)}$ par :

$$ (3.15) \qquad \begin{aligned} & \vec{u}_h^{(m+1)} \in W^h, \quad \vec{u}_h^{(m+1)}|_{\partial\Omega_h} = \vec{g}_h \\ & va_h(\vec{u}_h^{(m+1)}, \vec{v}_h) + \frac{1}{\epsilon}(\vec{v}_h \cdot \vec{u}_h^{(m+1)}, \vec{v}_h \cdot \vec{v}_h) + \hat{b}_h(\vec{u}_h^{(m+1)}, \vec{u}_h^{(m+1)}, \vec{v}_h) \\ & = (\vec{f}, \vec{v}_h)_h + (p_h^{(m)}, \vec{v}_h \cdot \vec{v}_h) \quad \forall v_h \in W_o^h \\ & p_h^{(m+1)} = p_h^{(m)} - \rho \vec{v}_h \cdot \vec{u}_h^{(m+1)} \end{aligned} $$

N.B. : On expliquera au § 4 comment résoudre effectivement (3.15).

**Algorithme (II) (UZAWA) (Cf. CROUZEIX [4])**

On donne $\vec{u}_h^{(o)}, p_h^{(o)}$ vérifiant (3.14).

Connaissant $p_h^{(m)}$, on définit un vecteur intermédiaire

$\vec{W}_h^{(m+1)} \in W_o^h$ par :
$$
\begin{equation} \tag{3.16.}
\begin{aligned}
& \left\|
  \begin{array}{@{}l@{}}
    \nu a_h(\vec{w}_h^{(m+1)}, \vec{v}_h) + \frac{1}{\epsilon}(\vec{\nabla}_h \cdot \vec{w}_h^{(m+1)}, \vec{\nabla}_h \cdot \vec{v}_h)_h + \hat{b}_h(\vec{u}_h^{(m)}, \vec{w}_h^{(m+1)}, \vec{v}_h) \\
    \qquad = (\vec{f}, \vec{v}_h)_h + (p_h^{(m)}, \vec{\nabla}_h \cdot \vec{v}_h)_h, \forall \vec{v}_h \in W_o^h
  \end{array}
\right\| \\
& \qquad = p_h^{(m+1)} = p_h^{(m)} - \omega^{(m)} \rho (\vec{\nabla}_h \cdot \vec{u}_h^{(m+1)}) \\
& \qquad \qquad \qquad \qquad \qquad \vec{u}_h^{(m+1)} = (1 - \omega^{(m)}) \vec{u}_h^{(m)} + \vec{w}_h^{(m+1)}
\end{aligned}
\end{equation}
$$

$\omega^{(m)}$ est un paramètre de "relaxation" compris entre 0 et 2.

Algorithme (III) (ARROW-MURWICZ) (Cf. TEMAM [14]).

On part de $u_h^{(o)}$, $p_h^{(o)}$ vérifiant (3.14)

$$
\begin{equation}
\left|
\begin{array}{@{}l@{}}
  a_h(\vec{u}_h^{(m+1)} - \vec{u}_h^{(m)}, \vec{v}_h) + \rho \left( \begin{array}{@{}l@{}}
    \nu a_h(\vec{u}_h^{(m)}, \vec{v}_h) + \frac{1}{\epsilon}(\vec{\nabla}_h \cdot \vec{u}_h^{(m+1)}, \vec{\nabla}_h \cdot \vec{v}_h)_h \\
    \qquad + \hat{b}_h(\vec{u}_h^{(m)}, \vec{u}_h^{(m+1)}, \vec{v}_h) - (p_h^{(m)}, \vec{\nabla}_h \cdot \vec{v}_h)_h \\
    \qquad - (\vec{f}, \vec{v}_h)_h
  \end{array} \right) \\
  \qquad \forall \vec{v}_h \in W_o^h \\
  \\
  \alpha(p_h^{(m+1)} - p_h^{(m)}) = - \rho \vec{\nabla}_h \cdot \vec{u}_h^{(m+1)}
\end{array}
\right|
\tag{3.17.}
\end{equation}
$$

Les résultats de ces algorithmes pour divers types de problèmes sont donnés au § 5.
# 4. - REMARQUES PRATIQUES

## 4.1. - Triangulation

Pour générer la triangulation automatiquement (ou du moins avec le minimum de manipulations "manuelles") deux types d'algorithmes ont été utilisés :

### a) La méthode des divisions successives (1)

Partant d'une triangulation initiale, grossière, on divise chaque triangle en quatre (2) nouveaux triangles, et on itère l'opération jusqu'à l'obtention d'une triangulation satisfaisante.

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<footnote>(1) Cf. BOURGAT [3]</footnote>

<footnote>(2) ou plus généralement en $k^2$ nouveaux triangles, $k = 2, 3, ...$</footnote>
Les avantages de cette méthode sont les suivants :

- Comme tous les triangles générés sont semblables aux triangles de départ, on ne risque pas d'introduire de petits angles ;

- on peut prévoir exactement le nombre de triangles et de sommets générés ;

$\text{Si } N_I^{(n)} = \text{nombre de sommets intérieurs de la } (n^{\text{ième}}) \text{ triangulation,}$

$N_F^{(n)} = \text{nombre de sommets frontière de la } (n^{\text{ième}}) \text{ triangulation,}$

$M^{(o)} = \text{nombre de triangles de la } (n^{\text{ième}}) \text{ triangulation,}$

$N_I^{(o)} = \text{nombre de sommets intérieurs de la triangulation initiale,}$

$N_F^{(o)} = \text{nombre de sommets frontière de la triangulation initiale,}$

$M^{(o)} = \text{nombre de triangles de la triangulation initiale,}$

les relations (2.1.) donnent :

$$
\left|
\begin{array}{l}
N_I^{(n)} = \frac{1}{2} \left( (k^{2n} - 1) M^{(o)} - (k^n - 1) N_F^{(o)} \right) + N_I^{(o)} \\
\\
N_F^{(n)} = k^n N_F^{(o)} \\
\\
M^{(n)} = k^{2n} M^{(o)}
\end{array}
\right.
$$

- on adapte facilement la méthode au cas où $\partial\Omega$ est courbe (mais pas trop compliquée).
Soient P et P' deux sommets voisins sur $\partial\Omega$, dans la $n^{\text{ième}}$ triangulation, Q le milieu de PP' ; dans la $(n+1)^{\text{ième}}$ triangulation, au lieu de générer un sommet en Q, on projette Q en Q' sur $\partial\Omega$ et on génère un sommet en Q'.
<figure><img src="image_7.png" /></figure>

(Il arrive alors que l'on génère des angles un peu plus petits que dans la triangulation initiale).

* On génère très facilement par cette méthode des maillages réguliers (voir pages......) ; par contre la méthode est peu pratique lorsque l'on doit beaucoup raffiner la triangulation dans certaines parties du domaine seulement (par exemple en mécanique des fluides pour tenir compte des couches limites, ou de zones où l'écoulement est particulièrement actif).

b) La méthode de la frontière initiale (Cf. AUSSEMS [1] BORSENBERGER [2] GEORGE [9]).

Le principe consiste à construire des triangles aussi voisins que possible de triangles équilatéraux (qui donnent théoriquement les meilleurs résultats pour la convergence) et s'appuyant sur une discrétisation de la frontière $\partial\Omega$, fournie par l'utilisateur. Celui-ci décrit donc la fontière $\partial\Omega_h$ comme une suite de segments qui doit former une boucle fermée : si le domaine n'est pas simplement connexe, il faut introduire
des "coupures" réalisant une fausse frontière intérieure qui sera "effacée" une fois la triangulation générée.

Le polygône frontière étant donné, le programme se charge de le remplir de triangles.

La création d'un triangle peut se faire de deux manières :

- en créant un sommet à l'intérieur ;

<figure><img src="image_8.png" /></figure>

(il convient de contrôler qu'il se trouve effectivement à l'intérieur !).

- en joignant deux sommets consécutifs de la frontière :

<figure><img src="image_9.png" /></figure>
(Cf. BORSENBERGER [2] pour les détails de la stratégie suivie).

Lorsqu'un triangle est créé, on reprend le processus pour créer un nouveau triangle, avec le domaine précédent diminué du dernier triangle créé ; et ainsi de suite jusqu'à la réduction du domaine à 1 triangle.

Le principal avantage de la méthode est la possibilité de traiter simplement les frontières courbes, avec raffinement éventuel dans certaines parties du domaine.

(Voir triangulations pages.....)

Les inconvénients sont les suivants :

- il est pratiquement impossible de connaître exactement à l'avance le nombre de triangles que l'on va obtenir ;

- à moins de travailler avec un écran de visualisation, on ne maîtrise pas la génération des triangles.

## 4.2. - Les tableaux de nombres définissant la triangulation

On pose :

$N_I = N_I + N_F$ = nombre total de sommets.

$S(1,i), S(2,i)$ = coordonnées du ième sommet.

$\text{Nu}(1,j), \text{Nu}(2,j), \text{Nu}(3,j)$ = numéros des sommets du jème triangle.

$N_m$ = nombre de noeuds intérieurs. (Les noeuds sont les milieux des côtés).

$N_{lm}$ = nombre total de noeuds.

$S_m(1,i), S_m(2,i)$ = coordonnées du ième noeud.

$\text{Nu}_m(1,j), \text{Nu}_m(2,j), \text{Nu}_m(3,j)$ = numéros des noeuds du jème triangle.
Pour repérer les sommets-frontière, on leur donne habituellement
un numéro compris entre N + 1 et N<sub>1</sub> ; de même les noeuds frontières
ont un numéro compris entre N<sub>m+1</sub> et N<sub>1m</sub>.

**Remarque** : pour traiter des conditions aux limites variées
(Dirichlet homogène ou non homogène, Neumann) sur les vitesses,
on devra introduire un tableau supplémentaire, associant à
chaque noeud-frontière un code de condition aux limites. ■

On formera aussi le tableau des voisins :

TAB(1,i), TAB(2,i), TAB(3,i), TAB(4,i) = numéros des noeuds voisins du noeud i.

<figure><img src="image_10.png" /></figure>

Lorsque le noeud i n'a que 2 voisins (noeud frontière)

$$
\mathrm{TAB}(3,i) = \mathrm{TAB}(4,i) = 0 -
$$

4.3. - <u>Les fonctions de base.</u>

Explicitons d'abord la base de $W_i^h$.

A chaque noeud $M_i$ de la triangulation $(i = 1, \dots, N_{1m})$ est
associée la fonction $w_{(i)}$ définie par :

$$
(4.1.) \quad \begin{cases}
- w_{(i)} \text{ vaut } 1 \text{ en } M_i ; \\
- w_{(i)} \text{ vaut } 0 \text{ en tous les autres noeuds } M_j \neq M_i.
\end{cases}
$$
Le support de $w_{(i)}$ est formé de 2 triangles adjacents à $M_i$ (1 seul triangle si $M_i \in \partial\Omega_h$)

<figure><img src="image_11.png" /></figure>

On note :

$$ w_{(1,i)} = \begin{pmatrix} w_{(i)} \\ 0 \end{pmatrix}, \quad w_{(2,i)} = \begin{pmatrix} 0 \\ w_{(i)} \end{pmatrix}. $$

Les $w_{(i)}$ ($i = 1, ..., N_{1m}$) forment une base de $W_h^h$.

Les $w_{(i)}$ ($i = 1, ..., N_m$) forment une base de $W_o^h$.

Les $\vec{w}_{(1,i)}, \vec{w}_{(2,i)}$ ($i = 1, ..., N_{1m}$) forment une base de $W_h^h$.

Les $\vec{w}_{(1,i)}, \vec{w}_{(2,i)}$ ($i = 1, ..., N_m$) forment une base de $W_o^h$.

On peut maintenant réécrire les algorithmes du paragraphe 3.

Alors (3.13) (par exemple) s'écrit :

(STOKES)

$$ (4.2.) \quad \vec{U}_h^{(m+1)} = \sum_{i=1}^{N_m} \begin{pmatrix} X_i^{(m+1)} \\ Y_i^{(m+1)} \end{pmatrix} w_{(i)} $$

$(1 \le i \le N_m)$
$$
(4.3.) \text{ a)} \quad \left( \bigvee_{j=1}^{N_m} \sum_{j=(m-1)}^{m} X_j^{(m \to 1)} \alpha_{ij} + \frac{1}{\epsilon} \sum_{j=1}^{N} \left[ \begin{aligned} & (X_j^{(m+1)} (\delta_{1h} w(j), \delta_{1h} w(i)))_h \\ & + Y_j^{(m+1)} (\delta_{2h} w(j), \delta_{2h} w(i))_h \end{aligned} \right] \right) \\
= (f_1, w(i))_h + (p_h, \delta_{1h} w(i))
$$

$$
(4.3.) \text{ b)} \quad \left( \bigvee_{j=1}^{N} \sum_{j=(m-1)}^{m} Y_j^{(m+1)} \alpha_{ij} + \frac{1}{\epsilon} \sum_j \left[ \begin{aligned} & X_j^{(m+1)} (\delta_{1h} w(j), \delta_{2h} w(i)) \\ & + Y_j^{(m+1)} (\delta_{2h} w(j), \delta_{2h} w(i)) \end{aligned} \right] \right) \\
= (f_2, w(i))_h + (p_h^{(m)}, \delta_{2h} w(i))
$$

$$
où \alpha_{ij} = \int_{\Omega_h} \overrightarrow{\nabla_h w}(j) \cdot \overrightarrow{\nabla_h w}(i) \, dx_1 \, dx_2
$$

(4.3.) est un système de (2 × N_m) équations linéaires.

On notera que, si $\frac{1}{\epsilon} = 0$ (pas de pénalisation, dualité pure)

les équations (4.3.)<sub>a</sub> et (4.3.)<sub>b</sub> sont <u>découplées</u>, et que l'on
a en fait un système de N<sub>m</sub> équations à résoudre.

<u>Remarques</u> : si les noeuds M<sub>i</sub> et M<sub>j</sub> ne sont pas voisins :

$\alpha_{ij} = 0 = (\delta_{1h} w(j), \delta_{1h} w(i)) = \dots$ (conséquence du caractère
<u>local</u> des fonctions de base).

La matrice du système est <u>creuse</u>.
(NAVIER-STOKES-STATIONNAIRE)

ALGORITHME (I)

$$
(4.4.) \quad \vec{U}_h^{(m+1)} = \sum_{j=1}^{N} \sum_{m=1}^{M} \begin{pmatrix} X_j^{(m+1)} \\ Y_j^{(m+1)} \end{pmatrix} w_{(j)} + \sum_{j=N_{m+1}}^{N} \sum_{k=N_{m+1}}^{1m} \vec{g}_h(u_j) w_{(j)}
$$

On pose pour simplifier les notations :

$$
\sum_{j=N_{m+1}}^{N_{1m}} g_h(u_j) w_{(j)} = \vec{G}_h = \begin{pmatrix} G_{h1} \\ G_{h2} \end{pmatrix} \in \mathcal{I} W^h
$$

On a :

$$
\vec{u}_h^{(m+1)} - \vec{G}_h \in \mathcal{I} W_o^h
$$

(3.15.) s'écrit :

$$
(4.5.) \text{ a)} \quad
\left[
\begin{array}{ll}
\displaystyle \sum_{j=1}^{N_m} X_j^{(m+1)} \quad \alpha_{ij} + \frac{1}{\epsilon} \sum_{j=1}^{N_m} \left[ \begin{array}{@{}l@{}} X_j^{(m+1)} \left( \delta_{2h^W(j)} , \delta_{1h^W(i)} \right) \\ \quad + Y_j^{(m+1)} \left( \delta_{2h^W(j)} , \delta_{1h^W(i)} \right) \end{array} \right] & \\
\displaystyle + \sum_{j=1}^{N_{1m}} X_j^{(m+1)} B_{ij} (X_1^{(M+1)}, \dots, X_{N_{1m}}^{(m+1)} ; Y_1^{(m+1)}, \dots, Y_{N_{1m}}^{(m+1)}) & \\
\multicolumn{2}{c}{j \neq i} \\
\multicolumn{2}{c}{= - \sum a_h (G_h, w_{(1,i)}) - \frac{1}{\epsilon} (\vec{\nabla}_h \cdot \vec{G}_h, \vec{\nabla}_h \cdot \vec{w}_{(1,i)}')} \\
\multicolumn{2}{c}{\quad + (f_1, w_{(i)}) + (p_h^{(m)}, \delta_{1h^W(i)}) \quad (i=1, \dots, N_m)}
\end{array}
\right]
$$
$$
\text{(4.5.) b)} \quad \left|
\begin{array}{@{}l@{}}
  \displaystyle \sum_{j=1}^{N_m} Y_j^{(m+1)} \quad \alpha_{ij} + \frac{1}{\epsilon} \sum_{j=1}^{N_m} \left[ \begin{aligned}[t]
  & X_j^{(m+1)} (\delta_{2h^w(j)}^{\delta_{2h^w(i)}}) \\
  & + Y_j^{(m+1)} (\delta_{2h^w(j)}^{\delta_{2h^w(i)}}) \end{aligned} \right] \\[2ex]
  + \sum_{\substack{j=1 \\ (j \neq i)}}^{N_{1m}} Y_j^{(m+1)} B_{ij} (X_1^{(m+1)}, \dots, X_{N_{1m}}^{(m+1)} ; Y_1^{(m+1)}, \dots, Y_{N_{1m}}^{(m+1)}) \\[1.5ex]
  = - \nu a_h (\vec{G}_h, \vec{w}_{(2,i)}) - \frac{1}{\epsilon} (\vec{\nabla}_h \cdot \vec{G}_h, \vec{\nabla}_h \cdot \vec{w}_{(2,i)}) \\[1.5ex]
  \qquad + (f_2, w_{(i)}) + (p_h^{(m)}, \delta_{2h^w(i)}) \quad (i = 1, \dots, N_m)
\end{array}
\right|
$$

où l'on a posé :

$$
\begin{equation}
\begin{aligned}
B_{ij}(X^{(m+1)}; Y^{(m+1)}) &= \hat{b}_h (\vec{u}_h^{(m+1)}, \vec{w}_{1,j}^{(m+1)}, \vec{w}_{(1,i)}^{(m+1)}) \\
&= \hat{b}_h (\vec{u}_h^{(m+1)}, \vec{w}_{(2,j)}^{(m+1)}, \vec{w}_{(2,i)}^{(m+1)})
\end{aligned}
\tag{4.6.}
\end{equation}
$$

(on a $B_{ij} = 0$ si $i = j$)

On peut expliciter le calcul de $B_{ij}$ :

- Si les noeuds $M_i$ et $M_j$ ne sont pas voisins, $B_{ij} = 0$ ;

- sinon, soit $r_{ij}$ l'aire du triangle $T_{ij}$ qui leur est commun,

$$
(4.7.) \quad B_{ij}(X; Y) = \frac{r_{ij}}{6} \left[ \begin{array}{@{}l@{}}
  X_i (\delta_{hx^w} (j) (T_{ij})) + Y_i (\delta_{hy^w} (j) (T_{ij})) \\
  - X_j (\delta_{hx^w} (i) (T_{ij})) - Y_j (\delta_{hy^w} (i) (T_{ij}))
\end{array} \right]
$$
(4.5.) est un système de 2 × N<sub>m</sub> équations non-linéaires que l'on
résout par l'algorithme suivant, du type sous-relaxation (suggéré
par M. FORTIN [7]).

(On pose $F_{1,i}$, $F_{2,i}$ = seconds membres dans (4.5.)<sub>a</sub>, (4.5.)<sub>b</sub>).

On se donne : MAXREL = nombre maximum d'itérations.

$\epsilon_{rel}$ = test d'arrêt,

$\omega_{rel}$ = paramètre de relaxation, compris entre 0 et 1.
(Calcul de $X^{(m+1)}$, $Y^{(m+1)}$ dans l'algorithme (I))

On pose :

$$X = X^{(m)}, Y = Y^{(m)}$$

<figure><img src="image_12.png" /></figure>
<figure><img src="image_13.png" /></figure>

ALGORITHME (II)

$$ \rightarrow_{U_h}^{(m)} = \sum_{j=1}^{N_{lm}} \begin{pmatrix} X_j^{(m)} \\ Y_j^{(m)} \end{pmatrix} \quad w(j), \quad \begin{pmatrix} X_j^{(m)} \\ Y_j^{(m)} \end{pmatrix} = \rightarrow_{g_h}^{(m_j)}, \quad j > N_m $$

$$ \rightarrow_{W_h}^{(m)} = \sum_{j=1}^{N_m} \begin{pmatrix} \tilde{X}_j^{(m)} \\ \tilde{Y}_j^{(m)} \end{pmatrix} \quad w(j) $$

On réécrit (3.16.) sous la forme :

$$ \begin{align} & \sum_{j=1}^{N_m} \tilde{X}_j^{(m+1)} \quad (v\alpha_{ij} + \frac{1}{\epsilon} (\delta_{1h}w(j), \delta_{2h}w(i))_h) \nonumber \\ & + \frac{1}{\epsilon} \sum_{j=1}^{N} \tilde{Y}_j^{(m+1)} \quad (\delta_{2h}w(j), \delta_{1h}w(i))_h \tag{4.8.) a} \\ & + \sum_{\substack{j=1 \\ j \neq i}}^{N_m} \tilde{X}_j^{(m+1)} \quad B_{ij}(X^{(m)}; Y^{(m)}) \nonumber \\ & = (f_1, w(i))_h + (p_h^{(m)}, \delta_{1h}w(i))_h, \quad (i=1, \dots, N_m) \nonumber \end{align} $$
$$
\begin{align*}
& \sum_{j=1}^{N} \int_{Y_j}^{\sim} (m+1) (v_{\alpha_{ij}} + \frac{1}{\epsilon} (\delta_{2h^W(j)}, \delta_{2h^W(i)})_h) \\
& \qquad + \frac{1}{\epsilon} \sum_{j=1}^{N_m} \int_{X_j}^{\sim} (m+1) (\delta_{2h^W(j)}, \delta_{2h^W(i)})_h \\
& \qquad (4.8.) \quad \text{b} \qquad + \sum_{\substack{j=1 \\ j \neq i}}^{N_m} \int_{Y_j}^{\sim} (m+1) B_{ij}(X^{(m)}; Y^{(m)}) \\
& \qquad \qquad \qquad \qquad \qquad \qquad \qquad j \neq i
\end{align*}
$$

$$
\begin{align*}
& \sum_{j=1}^{N_m} \int_{Y_j}^{\sim} (m+1) B_{ij}(X^{(m)}; Y^{(m)}) \\
& \qquad j \neq i \\
&= (f_2, w_{(i)})_h^{(m)} + (p_h^{(m)}, \delta_{2h^W(i)})_h^{(m)}, \quad (i = 1, \dots, N_m)
\end{align*}
$$

(4.8.) a et (4.8.)$_b$ sont deux systèmes couplés de $N_m$ équations linéaires : on résout (4.8.) par surrelaxation.

Remarque : Si $\frac{1}{\epsilon} = 0$ (pas de pénalisation), (4.8.)$_a$ et (4.8.)$_b$ sont découplés et ont la même matrice ; on peut les résoudre parallèlement.

ALGORITHM (III)

On obtient un système d'équations linéaires analogues à (4.8.)
que l'on ne détaillera ; on résout par surrelaxation.

N.B. : Pour tester la convergence, on posera :

$$
(4.9.) \quad C_h^{(m)} = \max_{T \in \mathcal{T}_h} | \vec{\nabla} \cdot (\vec{U}_h^{(m)} |_T) |
$$

4.4. - Propriétés des fonctions de base et de la matrice αᵢⱼ

Soit T un triangle de $\mathcal{C}_h$, $A_1, A_2, A_3$ ses sommets et $M_1, M_2, M_3$
les milieux des côtés.
<figure><img src="image_14.png" /></figure>

Soient $\mu_1^T, \mu_2^T, \mu_3^T$ les fonctions de base sur T, et $\lambda_1^T, \lambda_2^T, \lambda_3^T$
les coordonnées barycentriques : $\mu_i(M_j) = \delta_{ij}, \lambda_i(A_j) = \delta_{ij}$ ;
$i, j = 1, 2, 3$.

On a les relations :

$$
\begin{align*}
& \lambda_1^T + \lambda_2^T + \lambda_3^T = 2 = \mu_1^T + \mu_2^T + \mu_3^T \\
(4.10.) \quad & \mu_1^T = \lambda_1^T + \lambda_2^T - \lambda_3^T \\
& \mu_2^T = \lambda_2^T + \lambda_3^T - \lambda_1^T \\
& \mu_3^T = \lambda_3^T + \lambda_1^T - \lambda_2^T
\end{align*}
$$
<figure><img src="image_15.png" /></figure>

On vérifie que les gradients des fonctions de base sont portés par les hauteurs du triangle :

$$ (4.11.) \quad \left| \begin{array}{l} \vec{\nabla}_{\mu_1}^T = \frac{2}{|N_3 M_2|^2} \cdot \vec{N_3 M_1}, \\ \vec{\nabla}_{\mu_2}^T = \frac{2}{|N_1 M_2|^2} \cdot \vec{N_1 M_2}, \\ \vec{\nabla}_{\mu_3}^T = \frac{2}{|N_2 M_3|^2} \cdot \vec{N_2 M_3} \end{array} \right. $$

d'où l'on déduit, à l'aide de (4.9.):

$$ (4.12.) \quad \left| \begin{array}{l} \vec{\nabla}_{\mu_1}^T \cdot \vec{\nabla}_{\mu_2}^T \le 0, \quad \vec{\nabla}_{\mu_2}^T \cdot \vec{\nabla}_{\mu_3}^T \le 0 \\ |\vec{\nabla}_{\mu_1}^T|^2 = |\vec{\nabla}_{\mu_1}^T \cdot \vec{\nabla}_{\mu_2}^T| + |\vec{\nabla}_{\mu_1}^T \cdot \vec{\nabla}_{\mu_3}^T| \end{array} \right. $$
Par définition de $\alpha_{ij}$ :

$$ \alpha_{ij} = \sum_{T \in \mathcal{T}_h} \text{aire}(T) \stackrel{\vee}{V}_{w(i)}(T) \cdot \stackrel{\vee}{V}_{w(j)}(T) $$

(en fait la sommation s'étend seulement aux triangles voisins de i et j).

On peut donc résumer les propriétés de $(\alpha_{ij})_{i,j = 1,\dots, N_m}$

**Proposition :**

i) la matrice $A = (\alpha_{ij})_{i,j = 1,\dots, N_m}$ est symétrique définie positive.

ii) $\alpha_{ij} \le 0$ pour $i \ne j$.

iii) A est à dominance diagonale forte :

$$ \alpha_{ii} \ge \sum_{j \ne i} |\alpha_{ji}| $$

(avec inégalité stricte pour au moins un i).

La matrice A n'est pas toujours irréductible ; cette propriété est liée à la présence d'angles droits dans le réseau i si j et k sont deux noeuds voisins, $\alpha_{jk} = 0$ si et seulement si les côtés correspondants sont perpendiculaires.
Par exemple, pour la triangulation suivante :

<figure><img src="image_16.png" /></figure>

A est irréductible ; les noeuds M, N, P, Q ne sont pas reliés entre eux. Par contre la matrice A associée à la triangulation Fig. 4.4. est irréductible et 2 cyclique (la matrice de Jacobi est cyclique d'ordre 2). On connaît alors un algorithme de calcul du paramètre de surrelaxation optimal associé à A (cf. VARGA [15]).
<figure><img src="image_17.png" /></figure>

Enfin on vérifie la propriété suivante :

$$ (4.13.) \quad \int_{T} \left( \mu_{i}^{T} \cdot \mu_{j}^{T} \right) dx_{1} dx_{2} = \begin{cases} 0 & \text{si } i \neq j \\ \frac{1}{3} & \text{aire (T) si } i = j \end{cases} $$

(voir paragraphe suivant).

### 4.5. - Calcul numérique des fonctions de base sur T

<figure><img src="image_18.png" /></figure>
On introduit le triangle de référence $\hat{T}$, et l'application affine $\gamma$ qui applique $\hat{T}$ sur $T$, qui est de la forme :

$$ M = \gamma (\hat{M}) = \begin{pmatrix} G_{1,1} & G_{2,1} \\ G_{1,2} & G_{2,2} \end{pmatrix} \cdot \hat{M} + \begin{pmatrix} G_{(3,1)} \\ G_{3,2} \end{pmatrix} $$

$$ G = \begin{pmatrix} G_{1,1} & G_{2,1} \\ G_{1,2} & G_{2,2} \end{pmatrix} \qquad g = \begin{pmatrix} G_{3,1} \\ G_{3,2} \end{pmatrix} $$

On calcule $G$ en utilisant les équations :

$$ A_i = G \cdot \hat{A}_i + g, \quad i = 1, 2, 3. $$

On obtient :

$$ G = \begin{pmatrix} (A_1 - A_3)_1 & (A_2 - A_3)_1 \\ (A_1 - A_3)_2 & (A_2 - A_3)_2 \end{pmatrix} \cdot \hat{C} $$

où la matrice $\hat{C}$ ne dépend pas de $T$ :

$$ \hat{C} = \frac{1}{(\hat{A}_1 - \hat{A}_3)_1 (\hat{A}_2 - \hat{A}_3)_2 - (\hat{A}_1 - \hat{A}_3)_2 (\hat{A}_2 - \hat{A}_3)_1} \begin{pmatrix} (\hat{A}_2 - \hat{A}_3)_2 - (\hat{A}_2 - \hat{A}_3)_1 \\ -(\hat{A}_1 - \hat{A}_3)_2 & (\hat{A}_1 - \hat{A}_3)_1 \end{pmatrix} $$
$$ g = A_1 - G \cdot \hat{M}_1 $$
Connaissant $\gamma$, on calcule facilement :

$$
(4.14) \quad \text{aire } (T) = \frac{1}{2} \cdot \text{dét } G
$$

et $\gamma^{-1}$ :

$$
(4.15) \quad \hat{M} = \gamma^{-1} (M) = G' \cdot M + g' -
$$

On a les relations suivantes :

$$
\left|
\begin{array}{l@{\quad}l@{\quad}l}
\lambda_1^T (M) = \hat{x}_2 (\hat{M}) & (M = \gamma (\hat{M})) & \\
\lambda_2^T (M) = 1 - \hat{x}_1 (\hat{M}) & \hat{x}_2 (\hat{M}) & \\
\lambda_3^T (M) = \hat{x}_1 (\hat{M}) & & \\
\multicolumn{1}{c@{\quad}}{\displaystyle\int_{\hat{T}}^{\hat{T}}} & \hat{x}_1^{i1} \quad \hat{x}_2^{i2} \quad d\hat{x}_1 \quad d\hat{x}_2 = \frac{i1! \ i2!}{(i1 + i2 + 2)!} &
\end{array}
\right.
$$

On déduit de (4.16) et (4.9.)

$$
(4.17) \quad \left|
\begin{array}{l}
\displaystyle \int_T \mu_i^T dx_1 dx_2 = \frac{\text{aire } (T)}{3} \\[1.5em]
\displaystyle \int_T \mu_i^T \mu_j^T dx_1 dx_2 = \delta_{ij} \frac{\text{aire } (T)}{3} , \quad i, j = 1, 2, 3.
\end{array}
\right.
$$

On calcule $\nabla \mu_i^T$ ($i=1, 2, 3$) à l'aide de $G' = G^{-1}$,

$$
G' = \begin{pmatrix} G'_{11} & G'_{21} \\ G'_{12} & G'_{22} \end{pmatrix}
$$
$$ \frac{\partial \mu_1^T}{\partial x_j} = -2 G'_{j1} $$

$$ \frac{\partial \mu_2^T}{\partial x_j} = -2 G'_{j2} $$

$$ \frac{\partial \mu_3^T}{\partial x_j} = 2(G'_{j1} + G'_{j2}) $$

## 5. - RESULTATS NUMERIQUES

### 5.1. - Problème de la cavité avec paroi glissante
(Cf. FORTIN - TEMAM - PEYRET [8] ).

<figure><img src="image_19.png" /></figure>

$$ \vec{f} = 0 $$

$\Omega$ est le carré $[0, 1] \times [0, 1]$, les conditions aux limites sont :

$$ \vec{u} = \begin{pmatrix} 1 \\ 0 \end{pmatrix} \text{ sur } x_2 = 1 $$

$$ \vec{u} = \begin{pmatrix} 0 \\ 0 \end{pmatrix} \text{ sur } x_2 < 1 $$

Le nombre de Reynolds est égal à $\frac{1}{\nu}$
TRIANGULATION 1.1.

<figure><img src="image_20.png" /></figure>

NOMBRE DE TRIANGLES = 128

NOMBRE DE SOMMETS = 81

NOMBRE DE SOMMETS INTERIEURS = 49

NOMBRE DE POINTS MILIEUX = 208

NOMBRE DE POINTS MILIEUX INTERIEURS = 176.
<figure><img src="image_21.png" /><figcaption>TRIANGULATION 1.2.</figcaption></figure>

NOMBRE DE TRIANGLES = 512

NOMBRE DE SOMMETS = 289

NOMBRE DE SOMMETS INTERIEURS = 225

NOMBRE DE POINTS MILIEUX = 800

NOMBRE DE POINTS MILIEUX INTERIEURS = 736.
<figure><img src="image_22.png" /><figcaption>TRIANGULATION 1.3.</figcaption></figure>

NOMBRE DE TRIANGLES = 164

NOMBRE DE SOMMETS = 100

NOMBRE DE SOMMETS INTERIEURS = 66

NOMBRE DE POINTS MILIEUX = 263

NOMBRE DE POINTS MILIEUX INTERIEURS = 229
5.1.1. - Algorithme (I)

$$ i) \nu = 10^{-2} $$

(Nombre de Reynolds = 100, Triangulation (1.1). (128 Triangles)).

$$ \cdot 1/\epsilon = 0 $$

<table><tr><td>ρ</td><td>.988 × 10<sup>-2</sup></td><td>.98 × 10<sup>-2</sup></td><td>.96 × 10<sup>-2</sup></td><td>.90 × 10<sup>-2</sup></td></tr><tr><td>m tel que C<sub>h</sub><sup>(m)</sup> ≤ 10<sup>-2</sup></td><td>>50</td><td>45</td><td>36</td><td>30</td></tr><tr><td>ρ</td><td>.89 × 10<sup>-2</sup></td><td>.87 × 10<sup>-2</sup></td><td>.84 × 10<sup>-2</sup></td><td>.8 × 10<sup>-2</sup></td><td>.7 × 10<sup>-2</sup></td></tr><tr><td>m tel que C<sub>h</sub><sup>(m)</sup> ≤ 10<sup>-2</sup></td><td>31</td><td>31</td><td>33</td><td>34</td><td>42</td></tr></table>

$$ \cdot 1/\epsilon = 5 $$

<table><tr><td>ρ</td><td>12.</td><td>10.</td><td>8.</td><td>7.</td><td>6.</td><td>5.</td><td>3.</td><td>1.</td></tr><tr><td>m tel que C<sub>h</sub><sup>(m)</sup> ≤ 10<sup>-2</sup></td><td>OVERFLOW</td><td>>100</td><td>10</td><td>4</td><td>3</td><td>3</td><td>4</td><td>8</td></tr><tr><td>m tel que C<sub>h</sub><sup>(m)</sup> ≤ 10<sup>-3</sup></td><td></td><td>>50</td><td>6</td><td>6</td><td>6</td><td>8</td><td>13</td><td>29</td></tr><tr><td>m tel que C<sub>h</sub><sup>(m)</sup> ≤ 10<sup>-4</sup></td><td></td><td>>50</td><td>28</td><td>28</td><td>28</td><td>30</td><td>32</td><td>>50</td></tr></table>

Temps de calcul (*) pour ρ = 6. : | 3-itérations : 1' | 6-itérations : 2'30"

<footnote>(*) Tous les temps de calcul sont donnés pour CII 10070</footnote>
. $1/\epsilon = 10.$

<table><thead><tr><th>ρ</th><th>20.</th><th>17.</th><th>16.</th><th>15.</th><th>14.</th><th>13.</th><th>12.</th><th>11.</th><th>10.</th><th>6.</th><th>3.</th></tr></thead><tbody><tr><td>m tel que C<sub>h</sub><sup>(m)</sup> ≤ 10<sup>-3</sup></td><td colspan="2">OVERFLOW</td><td>8</td><td>6</td><td>6</td><td>5</td><td>4</td><td>4</td><td>4</td><td>6</td><td>10</td></tr><tr><td>m tel que C<sub>h</sub><sup>(m)</sup> ≤ 10<sup>-4</sup></td><td>22</td><td>20</td><td>22</td><td>24</td><td>25</td><td>28</td><td>30</td><td>32</td><td>36</td><td>40</td><td></td></tr></tbody></table>

Temps de calcul : pour ρ = 16, 20 itérations : 6'20"

. $1/\epsilon = 15$

<table><thead><tr><th>ρ</th><th>25.</th><th>20.</th><th>16.</th><th>8.</th></tr></thead><tbody><tr><td>m tel que C<sub>h</sub><sup>(m)</sup> ≤ 10<sup>-2</sup></td><td>>50</td><td>5</td><td>3</td><td>3</td></tr><tr><td>m tel que C<sub>h</sub><sup>(m)</sup> ≤ 10<sup>-3</sup></td><td></td><td>12</td><td>4</td><td>6</td></tr><tr><td>m tel que C<sub>h</sub><sup>(m)</sup> ≤ 10<sup>-4</sup></td><td>>50</td><td>18</td><td>38</td><td></td></tr></tbody></table>

ii) $v = 5.10^{-3}$ Triangulation (1.1). (128 triangles) Reynolds = 200

. $1/\epsilon = 5$

<table><thead><tr><th>ρ</th><th>9.</th><th>8.</th><th>7.</th><th>6.95</th><th>6.9</th><th>5.0</th><th>4.7</th><th>4.0</th><th>3.5.</th><th>3.0</th></tr></thead><tbody><tr><td>m tel que C<sub>h</sub><sup>(m)</sup> ≤ 10<sup>-3</sup></td><td>18</td><td>10</td><td>8</td><td>6</td><td>6</td><td>6</td><td>8</td><td>11</td><td>14</td><td>18</td></tr></tbody></table>
### 5.1.2. - Algorithme (III)

i) $v = 1$ Triangulation (1.1)
(128 triangles)

. $1/\epsilon = 0$

<table><tr><td>α</td><td>5.</td><td>5.</td><td>5.</td><td>5.</td><td>5.</td><td>5.</td></tr><tr><td>ρ</td><td>1.35</td><td>1.3</td><td>1.25</td><td>1.20</td><td>1.</td><td>0,9</td></tr><tr><td>m tel que C<sub>h</sub><sup>(m)</sup> ≤ 10<sup>-2</sup></td><td>OVERFLOW</td><td>>50</td><td>33</td><td>34</td><td>40</td><td>>50</td></tr></table>

<table><tr><td>α</td><td>3</td><td>3.</td><td>3.</td><td>3.</td><td>3.</td><td>3.</td></tr><tr><td>ρ</td><td>1.2</td><td>1.1</td><td>1.0</td><td>0,9</td><td>0,85</td><td>0,75</td></tr><tr><td>m tel que C<sub>h</sub><sup>(m)</sup> ≤ 10<sup>-2</sup></td><td>OVERFLOW</td><td>23</td><td>24</td><td>26</td><td>27</td><td>30</td></tr></table>

<table><tr><td>α</td><td>2.</td><td>2.</td><td>2.</td><td>2.</td><td>2.</td><td>2.</td><td>2.</td></tr><tr><td>ρ</td><td>1,0</td><td>0,9</td><td>0,8</td><td>0,7</td><td>0,6</td><td>0,5</td><td>0,4</td></tr><tr><td>m tel que C<sub>h</sub><sup>(m)</sup> ≤ 10<sup>-2</sup></td><td>OVERFLOW</td><td>16</td><td>18</td><td>20</td><td>23</td><td>25</td><td>29</td></tr></table>

<table><tr><td>α</td><td>1.</td><td>1.</td><td>1.</td><td>1.</td><td>1.</td><td>1.</td><td>1.</td></tr><tr><td>ρ</td><td>0,8</td><td>0,75</td><td>0,70</td><td>0,65</td><td>0,60</td><td>0,50</td><td>0,40</td></tr><tr><td>m tel que C<sub>h</sub><sup>(m)</sup> ≤ 10<sup>-2</sup></td><td>OVERFLOW</td><td>10</td><td>11</td><td>13</td><td>15</td><td>17</td><td>22</td></tr></table>

<table><tr><td>α</td><td>0.5</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>ρ</td><td>0.6</td><td>0.5</td><td>0.45</td><td>0.40</td><td>0.35</td><td></td><td></td></tr><tr><td>m tel que C<sub>h</sub><sup>(m)</sup> ≤ 10<sup>-2</sup></td><td>OVERFLOW</td><td>19</td><td>21</td><td>24</td><td>29</td><td></td><td></td></tr><tr><td>m tel que C<sub>h</sub><sup>(m)</sup> ≤ 10<sup>-3</sup></td><td></td><td>>100</td><td>33</td><td>39</td><td>59</td><td></td><td></td></tr><tr><td>m tel que C<sub>h</sub><sup>(m)</sup> ≤ 10<sup>-4</sup></td><td></td><td></td><td>76</td><td>87</td><td>100</td><td></td><td></td></tr></table>

Pour $α = 0,5$, $ρ = 0,45$

76 itérations

Temps d'exécution : 2'22"

Pour $α = 5, 3, 2, 1$, il faut plus de 100 itérations pour obtenir $C_{h}^{(m)} \le 10^{-3}$
• 1/ε = 5

<table><tr><td>α</td><td>2.</td><td>2.</td><td>2.</td><td>2.</td><td>2.</td><td>2.</td></tr><tr><td>ρ</td><td>2.</td><td>1.8</td><td>1.7</td><td>1.6</td><td>1.55</td><td>1.50</td><td>1.00</td></tr><tr><td>m tel que C<sub>h</sub><sup>(m)</sup> ≤ 10<sup>-2</sup></td><td>OVER-FLOW</td><td>21</td><td>23</td><td>24</td><td>25</td><td>26</td><td>38</td></tr><tr><td>m tel que C<sub>h</sub><sup>(m)</sup> ≤ 10<sup>-3</sup></td><td></td><td>>100</td><td>>100</td><td>35</td><td>35</td><td>36</td><td>52</td></tr><tr><td>m tel que C<sub>h</sub><sup>(m)</sup> ≤ 10<sup>-4</sup></td><td></td><td></td><td>60</td><td></td><td>53</td><td>55</td><td>74</td></tr></table>

<table><tr><td>α</td><td>1.</td><td>1.</td><td>1.</td><td>1.</td><td>1.0</td><td>1.0</td><td>1.0</td></tr><tr><td>ρ</td><td>2.</td><td>1.57</td><td>1.55</td><td>1.53</td><td>1.50</td><td>1.0</td><td>0.75</td></tr><tr><td>m tel que C<sub>h</sub><sup>(m)</sup> ≤ 10<sup>-2</sup></td><td>OVER-FLOW</td><td>12</td><td>25</td><td>12</td><td>12</td><td>19</td><td>10</td></tr><tr><td>m tel que C<sub>h</sub><sup>(m)</sup> ≤ 10<sup>-3</sup></td><td></td><td>>100</td><td>35</td><td>19</td><td>18</td><td>26</td><td>>100</td></tr><tr><td>m tel que C<sub>h</sub><sup>(m)</sup> ≤ 10<sup>-4</sup></td><td></td><td></td><td>53</td><td>>70</td><td>70</td><td>43</td><td></td></tr></table>
• 1/ε = 5

<table><tr><td>α</td><td>0.5</td><td>0.5</td><td>0.5</td><td>0.5</td><td>0.5</td><td>0.5</td><td>0.5</td></tr><tr><td>ρ</td><td>1.4</td><td>1.3</td><td>1.2</td><td>1.1</td><td>0,90</td><td>0,80</td><td>0,70</td><td>0,65</td><td>0,55</td></tr><tr><td>m tel que C<sub>h</sub><sup>(m)</sup> ≤ 10<sup>-2</sup></td><td>OVER-<br/>FLOW</td><td>7</td><td>7</td><td>8</td><td>9</td><td>10</td><td>12</td><td>13</td><td>15</td></tr><tr><td>m tel que C<sub>h</sub><sup>(m)</sup> ≤ 10<sup>-3</sup></td><td></td><td>10</td><td>10</td><td>11</td><td>14</td><td>15</td><td>17</td><td>19</td><td>22</td></tr><tr><td>m tel que C<sub>h</sub><sup>(m)</sup> ≤ 10<sup>-4</sup></td><td></td><td>34</td><td>22</td><td>34</td><td>41</td><td>44</td><td>45</td><td>49</td><td>56</td></tr></table>

Temps d'exécution
pour m itérations.

m = 34 : m = 22 : m = 34 : m = 41 : m = 44 :
1'11" 1' 1'10" 1'30" 1'36"

<table><tr><td>α</td><td>0.1</td><td>0.1</td><td>0.1</td><td>0.1</td><td>0.1</td><td>0.1</td><td>0.1</td><td>0.1</td></tr><tr><td>ρ</td><td>0.38</td><td>0.36</td><td>0.32</td><td>0.31</td><td>0.30</td><td>0.29</td><td>0.20</td><td>0.10</td></tr><tr><td>m tel que C<sub>h</sub><sup>(m)</sup> ≤ 10<sup>-2</sup></td><td>OVER-<br/>FLOW</td><td>11</td><td>12</td><td>13</td><td>13</td><td>14</td><td>19</td><td>38</td></tr><tr><td>m tel que C<sub>h</sub><sup>(m)</sup> ≤ 10<sup>-3</sup></td><td></td><td>16</td><td>19</td><td>19</td><td>29</td><td>21</td><td>30</td><td>64</td></tr><tr><td>m tel que C<sub>h</sub><sup>(m)</sup> ≤ 10<sup>-4</sup></td><td></td><td>53</td><td>59</td><td>61</td><td>46</td><td>52</td><td>50</td><td>>100</td></tr></table>
$$ \text{ii) } v = 10^{-2} $$

Triangulation (1.1) (128 triangles)

$ \cdot 1/\epsilon = 0 $

<table><tr><td>α</td><td>0.1</td><td>0.1</td><td>0.5</td><td>0.5</td><td>0.5</td></tr><tr><td>ρ</td><td>0.4</td><td>0.48</td><td>0.9</td><td>1.0</td><td>OVERFLOW</td></tr><tr><td>c<sub>h</sub><sup>(500)</sup></td><td>1.6 × 10<sup>-2</sup></td><td>OVERFLOW</td><td>10<sup>-2</sup></td><td>0.5</td><td></td></tr><tr><td>δw<sub>h</sub><sup>(500)</sup></td><td>1.3 × 10<sup>-3</sup></td><td></td><td>8 × 10<sup>-4</sup></td><td>2.8 × 10<sup>-3</sup></td><td></td></tr></table>

<table><tr><td>α</td><td>1.</td><td>1.</td><td>2.</td><td>2.</td><td>2.</td><td>2.</td></tr><tr><td>ρ</td><td>1.40</td><td>1.45</td><td>1.7</td><td>1.8</td><td>1.9</td><td>2.0</td></tr><tr><td>c<sub>h</sub><sup>(500)</sup></td><td>6.5 × 10<sup>-2</sup></td><td>OVERFLOW</td><td>3 × 10<sup>-3</sup></td><td>2.8 × 10<sup>-3</sup></td><td>4.2 × 10<sup>-3</sup></td><td>OVERFLOW</td></tr><tr><td>δw<sub>h</sub><sup>(500)</sup></td><td>3.5 × 10<sup>-4</sup></td><td></td><td>2.3 × 10<sup>-4</sup></td><td>2.6 × 10<sup>-4</sup></td><td>3. × 10<sup>-4</sup></td><td></td></tr></table>

<table><tr><td>α</td><td>3.</td><td>3.</td><td>3.</td><td>3.</td><td>3.</td></tr><tr><td>ρ</td><td>1.8</td><td>2.2</td><td>2.3</td><td>2.4</td><td>2.5</td></tr><tr><td>c<sub>h</sub><sup>(500)</sup></td><td>2.1 × 10<sup>-3</sup></td><td>1.3 × 10<sup>-3</sup></td><td>1.4 × 10<sup>-3</sup></td><td>4.7 × 10<sup>-3</sup></td><td>OVERFLOW</td></tr><tr><td>δw<sub>h</sub><sup>(500)</sup></td><td>1.9 × 10<sup>-4</sup></td><td>1.3 × 10<sup>-4</sup></td><td>1.1 × 10<sup>-4</sup></td><td>2.7 × 10<sup>-4</sup></td><td></td></tr></table>

<table><tr><td>α</td><td>4.</td><td>4.</td><td>4.</td><td>4.</td><td>4.</td></tr><tr><td>ρ</td><td>2.4</td><td><u>2.6</u></td><td>2.7</td><td>2.8</td><td>3.</td></tr><tr><td>m tel que c<sub>h</sub><sup>(m)</sup> × 10<sup>-3</sup></td><td>478</td><td><u>470</u></td><td>497</td><td>&gt; 500</td><td>OVERFLOW</td></tr><tr><td>m tel que δw<sub>h</sub><sup>(m)</sup> × 10<sup>-4</sup></td><td>454</td><td><u>451</u></td><td>484</td><td>&gt; 500</td><td></td></tr></table>

<table><tr><td>α</td><td>5.</td><td>5.</td><td>5.</td><td>5.</td><td>5.</td></tr><tr><td>ρ</td><td>2.7</td><td><u>2.9</u></td><td>3.0</td><td>3.1</td><td>3.2</td></tr><tr><td>m tel que c<sub>h</sub><sup>(m)</sup> &lt; 10<sup>-3</sup></td><td>435</td><td><u>417</u></td><td>433</td><td>&gt; 500</td><td>OVERFLOW</td></tr><tr><td>m tel que δw<sub>h</sub><sup>(m)</sup> &lt; 10<sup>-4</sup></td><td>419</td><td><u>409</u></td><td><u>407</u></td><td>475</td><td></td></tr></table>
<table><thead><tr><th>α</th><th>6.</th><th>6.</th><th>6.</th><th>6.</th><th>6.</th></tr></thead><tbody><tr><th>ρ</th><td>3.1</td><td><u>3.2</u></td><td>3.3</td><td>3.4</td><td>3.5</td></tr><tr><th>m tel que C<sub>h</sub><sup>(m)</sup> &lt; 10<sup>-3</sup></th><td>388</td><td><u>379</u></td><td>394</td><td>473</td><td>OVERFLOW</td></tr><tr><th>m tel que δw<sub>h</sub><sup>(m)</sup> &lt; 10<sup>-4</sup></th><td>385</td><td><u>376</u></td><td>395</td><td>438</td><td></td></tr></tbody></table>

<table><thead><tr><th>α.</th><th>8.</th><th>8.</th><th>8.</th><th>8.</th><th>8.</th><th>8.</th><th>8.</th></tr></thead><tbody><tr><th>ρ'</th><td>3.2</td><td>3.4</td><td>3.6</td><td><u>3.7</u></td><td>3.8</td><td>3.9</td><td>4.</td></tr><tr><th>m tel que C<sub>h</sub><sup>(m)</sup> &lt; 10<sup>-3</sup></th><td>381</td><td>362</td><td>354</td><td><u>337</u></td><td>343</td><td>376</td><td>OVERFLOW</td></tr><tr><th>m tel que δw<sub>h</sub><sup>(m)</sup> &lt; 10<sup>-4</sup></th><td>365</td><td>360</td><td>340</td><td>338</td><td><u>329</u></td><td>340</td><td></td></tr><tr><th>m tel que C<sub>h</sub><sup>(m)</sup> &lt; 10<sup>-4</sup></th><td>&gt; 500</td><td>&gt; 500</td><td>489</td><td><u>480</u></td><td>492</td><td>&gt; 500</td><td></td></tr><tr><th>m tel que δw<sub>h</sub><sup>(m)</sup> &lt; 10<sup>-5</sup></th><td>&gt; 500</td><td>&gt; 500</td><td>492</td><td><u>478</u></td><td>486</td><td>492</td><td></td></tr></tbody></table>

<table><thead><tr><th>α</th><th>9.</th><th>9.</th><th>9.</th><th>9.</th><th>9.</th><th>9.</th></tr></thead><tbody><tr><th>ρ</th><td>3.8</td><td>3.9</td><td><u>4.0</u></td><td>4.1</td><td>4.2</td><td>4.3</td></tr><tr><th>m tel que C<sub>h</sub><sup>(m)</sup> &lt; 10<sup>-3</sup></th><td>333</td><td>330</td><td><u>321</u></td><td>341</td><td>439</td><td>OVERFLOW</td></tr><tr><th>m tel que δw<sub>h</sub><sup>(m)</sup> &lt; 10<sup>-4</sup></th><td>324</td><td>315</td><td><u>303</u></td><td>315</td><td>402</td><td></td></tr><tr><th>m tel que C<sub>h</sub><sup>(m)</sup> &lt; 10<sup>-4</sup></th><td>467</td><td>474</td><td><u>452</u></td><td>486</td><td>&gt; 500</td><td></td></tr><tr><th>m tel que δw<sub>h</sub><sup>(m)</sup> &lt; 10<sup>-5</sup></th><td>473</td><td>473</td><td>469</td><td><u>465</u></td><td>&gt; 500</td><td></td></tr></tbody></table>

<table><thead><tr><th>α</th><th>10.</th><th>10.</th><th>10.</th><th>10.</th><th>10.</th><th>10.</th><th>10.</th></tr></thead><tbody><tr><th>ρ</th><td>3.8</td><td>4.0</td><td>4.1</td><td>4.2</td><td>4.3</td><td>4.4</td><td>4.6</td></tr><tr><th>m tel que C<sub>h</sub><sup>(m)</sup> &lt; 10<sup>-3</sup></th><td>331</td><td>322</td><td></td><td>306</td><td>314</td><td>364</td><td>OVERFLOW</td></tr><tr><th>m tel que δw<sub>h</sub><sup>(m)</sup> &lt; 10<sup>-4</sup></th><td>311</td><td>307</td><td></td><td>304</td><td>300</td><td>330</td><td></td></tr><tr><th>m tel que C<sub>h</sub><sup>(m)</sup> &lt; 10<sup>-4</sup></th><td>479</td><td>449</td><td></td><td>434</td><td>441</td><td>500</td><td></td></tr><tr><th>m tel que δw<sub>h</sub><sup>(m)</sup> &lt; 10<sup>-5</sup></th><td>464</td><td>437</td><td></td><td>429</td><td>441</td><td>485</td><td></td></tr></tbody></table>
<table><thead><tr><th>α</th><th>11.</th><th>11.</th><th>11.</th><th>11.</th><th>11.</th><th>11.</th><th>11.</th></tr></thead><tbody><tr><th>ρ</th><td>4.20</td><td>4.36</td><td>4.40</td><td>4.44</td><td>4.50</td><td>4.60</td><td>5.00</td></tr><tr><th>m tel que C<sub>h</sub><sup>(m)</sup> &lt; 10<sup>-3</sup></th><td>299</td><td>297</td><td>288</td><td>287</td><td>302</td><td>330</td><td>OVERFLOW</td></tr><tr><th>m tel que δw<sub>h</sub><sup>(m)</sup> &lt; 10<sup>-4</sup></th><td>285</td><td>283</td><td>280</td><td>295</td><td>289</td><td>301</td><td></td></tr><tr><th>m tel que C<sub>h</sub><sup>(m)</sup> &lt; 10<sup>-4</sup></th><td>422</td><td>422</td><td>409</td><td>411</td><td>418</td><td>500</td><td></td></tr><tr><th>m tel que δw<sub>h</sub><sup>(m)</sup> &lt; 10<sup>-5</sup></th><td>421</td><td>413</td><td>417</td><td>406</td><td>412</td><td>421</td><td></td></tr></tbody></table>

On donne ci-dessous les lignes de courant de la solution obtenue avec la triangulation 1.1., $α = 11.$, $ρ = 4.4$, le nombre de Reynolds étant égal à 100.

Le tracé de ces courbes est obtenu de la manière suivante : on calcule, par une méthode d'éléments finis conformes, la fonction de courant $\psi$ vérifiant :

$$ \left| \begin{array}{l} -\Delta \psi = \operatorname{div} u \\ \\ \psi|_{\Gamma} = 0. \end{array} \right. $$

(La frontière étant convexe $\psi|_{\Gamma}$ est une constante que l'on peut prendre nulle).

Les lignes de courant sont les lignes de niveau de $\psi$ - comme $\psi$ est approchée par une fonction linéaire par morceaux, on obtient des lignes brisées, chaque segment de droite correspondant au passage d'une ligne sur le triangle.
<figure><img src="image_23.png" /></figure>

NAVIER-STOKES STATIONNAIRE - VISCOSITE = 0.9999E-02
$$ \text{iii)} \quad \nu = 10^{-2} \quad \begin{array}{l} \text{Triangulation (1.3)} \\ \text{(164 triangles).} \end{array} $$

$$ \therefore 1/\varepsilon = 5. $$

<table><tr><td>α</td><td>.1</td><td>.1</td><td>.1</td><td>.1</td><td>.1</td><td>.1</td><td>.1</td></tr><tr><td>ρ</td><td>.32</td><td>.31</td><td>.30</td><td>.29</td><td>.28</td><td>.27</td><td>.26</td><td>.25</td></tr><tr><td>m tel que $c_h^{(m)} \le 10^{-1}$</td><td>OVER-FLOW</td><td>33</td><td>22</td><td>23</td><td>24</td><td>25</td><td>26</td><td>27</td></tr><tr><td>m tel que $c_h^{(m)} \le 10^{-2}$</td><td></td><td>83</td><td>90</td><td>95</td><td>98</td><td>&gt;100</td><td>&gt;100</td><td>&gt;100</td></tr></table>

<table><tr><td>α</td><td>.05</td><td>.05</td><td>.05</td><td>.05</td><td>.05</td></tr><tr><td>ρ</td><td>.22</td><td>.20</td><td>.29</td><td>.18</td><td>.17</td></tr><tr><td>m tel que $c_h^{(m)} \le 10^{-2}$</td><td>OVER-FLOW</td><td>63</td><td>63</td><td>66</td><td>70</td></tr></table>
5.2.-

<figure><img src="image_24.png" /></figure>

$$ \vec{f} = \vec{0} $$

$$ \Omega = \{ (x_1, x_2) : (x_1^2 + x_2^2)^{1/2} < 5, ((x_1 - 1)^2 + x_2^2)^{1/2} > 2 \} $$

Conditions aux limites : le cercle extérieur est immobile, le cercle intérieur tourne à une vitesse $\omega$.

$$ \text{Sur } \Gamma_1 : (x_1^2 + x_2^2 = 25) : \vec{u} = 0 $$

$$ \text{Sur } \Gamma_2 : ((x_1 - 1)^2 + x_2^2 = 4) : \vec{u}(x_1, x_2) = \begin{pmatrix} -\omega x_2 \\ \omega(x_1 - 1) \end{pmatrix} $$

$$ \text{Le nombre de Reynolds vaut : } \frac{(2\omega) \times d}{\nu} = \frac{4\omega}{\nu} $$
<figure><img src="image_25.png" /></figure>

NOMBRE DE TRIANGLES = 126

NOMBRE DE SOMMETS = 46

NOMBRE DE SOMMETS INTERIEURS = 80

NOMBRE DE POINTS MILIEUX = 172

NOMBRE DE POINTS MILIEUX INTERIEURS = 206
<figure><img src="image_26.png" /></figure>

NOMBRE DE TRIANGLES = 322

NOMBRE DE SOMMETS = 188

NOMBRE DE SOMMETS INTERIEURS = 134

NOMBRE DE POINTS MILIEUX = 510

NOMBRE DE POINTS MILIEUX INTERIEURS = 456

TRIANGULATION 2.2.
### 5.2.1. - Algorithme (I)

$$ i) \quad \nu = 4 \times 10^{-2}, \omega = 1 \qquad \text{Triangulation (2.1).} $$
$$ \cdot \quad 1/\epsilon = 5 \qquad \text{(126 triangles)} $$

<table><thead><tr><th>ρ</th><th>1.0</th><th>9.5</th><th>9.0</th><th>8.5</th><th>8.0</th><th>7.</th><th>6.</th><th>5.</th><th>4.</th></tr></thead><tbody><tr><td>m tel que C<sub>h</sub><sup>(m)</sup> &lt; 10<sup>-2</sup></td><td>OVER-FLOW</td><td>20</td><td>12</td><td>10</td><td>10</td><td>10</td><td>11</td><td>12</td><td>13</td></tr><tr><td>m tel que C<sub>h</sub><sup>(m)</sup> &lt; 10<sup>-3</sup></td><td></td><td></td><td></td><td></td><td>22</td><td></td><td></td><td></td><td></td></tr><tr><td>m tel que C<sub>h</sub><sup>(m)</sup> &lt; 10<sup>-4</sup></td><td></td><td></td><td></td><td></td><td>129</td><td></td><td></td><td></td><td></td></tr></tbody></table>

Temps d'exécution : m = 45 :
pour m itérations : 9'12"

m = 129 :
30'

### 5.2.2. - Algorithme (II). Triangulation (2.1).

(126 triangles)

On a pris le paramètre de relaxation $\omega^{(m)}$ indépendant de m :

$$ i) \quad \nu = 1 \quad 5 vitesses de rotation \omega = 1 $$

Triangulation (2.1).

• $1/\epsilon = 0$

<table><thead><tr><th>ρ</th><th>10.</th><th>10.</th><th>8.</th><th>8.</th><th>8.</th></tr></thead><tbody><tr><td>ω<sup>(0)</sup></td><td>.1</td><td>.09</td><td>.125</td><td>.120</td><td>.110</td></tr><tr><td>m tel que C<sub>h</sub><sup>(m)</sup> &lt; 10<sup>-3</sup></td><td>OVERFLOW</td><td>&gt;70</td><td>&gt; 70</td><td>69</td><td>&gt; 70</td></tr></tbody></table>
<table><tr><td>ρ</td><td>2.6</td><td>2.6</td><td>2.6</td><td>2.6</td><td>2.6</td></tr><tr><td>ω<sup>(0)</sup></td><td>.38</td><td>.37</td><td>.36</td><td>.35</td><td>.30</td></tr><tr><td>m tel que C<sub>h</sub><sup>(m)</sup> &le; 10<sup>-3</sup></td><td>&gt;100</td><td>67</td><td>64</td><td>66</td><td>&gt;70</td></tr></table>

<table><tr><td>ρ</td><td>1.5</td><td>1.5</td><td>1.5</td><td>1.5</td><td>1.5</td></tr><tr><td>ω<sup>(0)</sup></td><td>.66</td><td>.64</td><td>.63</td><td>.62</td><td>.60</td></tr><tr><td>m tel que C<sub>h</sub><sup>(m)</sup> &le; 10<sup>-3</sup></td><td>&gt;100</td><td>67</td><td>63</td><td>64</td><td>66</td></tr><tr><td colspan="6">Temps d'exécution pour ρ = 1,5, ω<sup>(0)</sup> = 0.63, 63 itérations : 2'36".</td></tr></table>

<table><tr><td>ρ</td><td>1.</td><td>1.</td><td>1.</td><td>0.75</td><td>0.75</td><td>0.75</td><td>0.75</td><td>0.75</td></tr><tr><td>ω<sup>(0)</sup></td><td>0.96</td><td>0.94</td><td>0.89</td><td>1.25</td><td>1.24</td><td>1.23</td><td>1.17</td><td>1.13</td></tr><tr><td>m tel que C<sub>h</sub><sup>(m)</sup> &le; 10<sup>-3</sup></td><td>&gt;100</td><td>65</td><td>67</td><td>&gt;100</td><td>79</td><td>65</td><td>68</td><td>70</td></tr></table>

<table><tr><td>ρ</td><td>.65</td><td>.65</td><td>.65</td><td>.65</td><td>0.3</td><td>0.3</td></tr><tr><td>ω<sup>(0)</sup></td><td>1.35</td><td>1.34</td><td>1.32</td><td>1.20</td><td>1.6</td><td>1.5</td></tr><tr><td>m tel que C<sub>h</sub><sup>(m)</sup> &le; 10<sup>-3</sup></td><td>&gt;100</td><td>69</td><td>70</td><td>75</td><td>OVERFLOW</td><td>&gt;70</td></tr></table>

1/ε = 5

<table><tr><td>ρ</td><td>10.</td><td>10.</td><td>10.</td><td>10.</td><td>10.</td><td>10.</td><td>10.</td></tr><tr><td>ω<sup>(0)</sup></td><td>.90</td><td>.85</td><td>.80</td><td>.75</td><td>.70</td><td>.6</td><td>.4</td></tr><tr><td>m tel que C<sub>h</sub><sup>(m)</sup> &le; 10<sup>-3</sup></td><td>&gt;70</td><td>9</td><td>9</td><td>10</td><td>11</td><td>12</td><td>19</td></tr><tr><td>m tel que C<sub>h</sub><sup>(m)</sup> &le; 10<sup>-5</sup></td><td></td><td>&gt;70</td><td>48</td><td>51</td><td>55</td><td>66</td><td>&gt;70</td></tr></table>

Temps d'exécution pour ρ = 10., ω<sup>(0)</sup> = .8, 48 itérations : 1'50".
<table><tr><td>ρ</td><td>6.</td><td>6.</td><td>6.</td><td>6.</td><td>3.</td><td>3.</td><td>3.</td><td>3.</td></tr><tr><td>ω(0)</td><td>1.3</td><td>1.2</td><td>1.1</td><td>0.9</td><td>1.3</td><td>1.5</td><td>1.4</td><td>1.3</td><td>1.0</td></tr><tr><td>m tel que C<sub>h</sub>(m) &lt; 10<sup>-3</sup></td><td>25</td><td>10</td><td>11</td><td>14</td><td>32</td><td>16</td><td>18</td><td>19</td><td>24</td></tr><tr><td>m tel que C<sub>h</sub>(m) &lt; 10<sup>-5</sup></td><td>60</td><td>49</td><td>55</td><td>61</td><td>&gt;70</td><td>52</td><td>60</td><td>65</td><td>&gt;70</td></tr></table>

<table><tr><td>ρ</td><td>2.</td><td>2.</td><td>2.</td><td>2.</td><td>2.</td><td>2.</td><td>2.</td><td>2.</td></tr><tr><td>ω(0)</td><td>1.70</td><td>1.65</td><td>1.60</td><td>1.55</td><td>1.50</td><td>1.4</td><td>1.0</td><td>0.9</td></tr><tr><td>m tel que C<sub>h</sub>(m) &lt; 10<sup>-3</sup></td><td>&gt;70</td><td>39</td><td>24</td><td>23</td><td>25</td><td>26</td><td>37</td><td>41</td></tr><tr><td>m tel que C<sub>h</sub>(m) &lt; 10<sup>-5</sup></td><td>&gt;70</td><td>62</td><td>65</td><td>70</td><td>69</td><td>&gt;70</td><td>&gt;70</td><td></td></tr></table>

<table><tr><td>ρ</td><td>1.5</td><td>1.5</td><td>1.5</td><td>1.5</td><td>1.5</td><td>1.5</td><td>1.5</td><td>1.5</td></tr><tr><td>ω(0)</td><td>1.7</td><td>1.65</td><td>1.60</td><td>1.55</td><td>1.5</td><td>1.5</td><td>1.5</td><td>1.5</td></tr><tr><td>m tel que C<sub>h</sub>(m) &lt; 10<sup>-3</sup></td><td>&gt;70</td><td>36</td><td>30</td><td>31</td><td>33</td><td></td><td></td><td></td></tr><tr><td>m tel que C<sub>h</sub>(m) &lt; 10<sup>-5</sup></td><td>&gt;70</td><td>&gt;70</td><td>&gt;70</td><td>&gt;70</td><td>&gt;70</td><td></td><td></td><td></td></tr></table>
ii)

<table><tr><td>v = 0,1 ; vitesse de rotation ω = 1</td></tr></table>

Triangulation (2.1). (126 triangles)

$$ \cdot 1/\epsilon = 5 $$

<table><tr><td>ρ</td><td>10.</td></tr><tr><td>ω<sup>(0)</sup></td><td>1.</td><td>.7</td><td>.6</td><td>.5</td></tr><tr><td>m tel que C<sub>h</sub><sup>(m)</sup> &lt; 10<sup>-2</sup></td><td>&gt;30</td><td>4</td><td>3</td><td>4</td></tr><tr><td>m tel que C<sub>h</sub><sup>(m)</sup> &lt; 10<sup>-3</sup></td><td></td><td>19</td><td>22</td><td>26</td></tr></table>

<table><tr><td>ρ</td><td>1.</td></tr><tr><td>ω<sup>(0)</sup></td><td>1.4</td><td>1.3</td><td>1.2</td><td>1.1</td><td>1.0</td><td>.9</td><td>.7</td></tr><tr><td>m tel que C<sub>h</sub><sup>(m)</sup> &lt; 10<sup>-2</sup></td><td>&gt;30</td><td>17</td><td>17</td><td>18</td><td>19</td><td>21</td><td>27</td></tr><tr><td>m tel que C<sub>h</sub><sup>(m)</sup> &lt; 10<sup>-3</sup></td><td></td><td>&gt;30</td><td>30</td><td>&gt;30</td><td>&gt;30</td><td>&gt;30</td><td>&gt;30</td></tr></table>

### 5.2.3. - Algorithme (III)

<table><tr><td>v = 4 × 10<sup>-2</sup> , vitesse de rotation ω = 1/</td></tr></table>

Triangulation (2.1). (126 triangles)

$$ \cdot 1/\epsilon = 5 $$

<table><tr><td>α</td><td>.5</td><td>.5</td><td>.5</td><td>.5</td><td>.5</td><td>.5</td><td>.5</td><td>.5</td></tr><tr><td>ρ</td><td>2.8</td><td>2.5</td><td>2.4</td><td>2.3</td><td>2.0</td><td>1.5</td><td>1.0</td><td>.5</td></tr><tr><td>m tel que C<sub>h</sub><sup>(m)</sup> &lt; 10<sup>-2</sup></td><td>18</td><td>18</td><td>14</td><td>15</td><td>17</td><td>22</td><td>33</td><td>66</td></tr></table>
<table><tr><td>α</td><td>.1</td><td>.1</td><td>.1</td><td>.1</td><td>.1</td><td>.1</td><td>.1</td></tr><tr><td>ρ</td><td></td><td>.57</td><td>.56</td><td>.54</td><td>.52</td><td>.51</td><td></td></tr><tr><td>m tel que C<sub>h</sub><sup>(m)</sup> ≪ 10<sup>-2</sup></td><td></td><td>17</td><td>18</td><td>18</td><td>19</td><td>19</td><td></td></tr></table>

<table><tr><td>α</td><td>.05</td><td>.05</td><td>.05</td><td>.05</td></tr><tr><td>ρ</td><td>.5</td><td>.4</td><td>.3</td><td>.2</td></tr><tr><td>m tel que C<sub>h</sub><sup>(m)</sup> ≪ 10<sup>-2</sup></td><td>OVERFLOW</td><td>3</td><td>3</td><td></td></tr><tr><td>m tel que C<sub>h</sub><sup>(m)</sup> ≪ 10<sup>-3</sup></td><td></td><td>11</td><td>16</td><td></td></tr></table>

On donne ci-dessous les lignes de courant de la solution de l'algorithme (II) avec :

$$v = 1, \quad \frac{1}{\varepsilon} = 5., \quad \rho = 10. \quad \omega^{(0)} = .8.$$

La fonction de courant, dont les lignes de niveau sont les lignes de courant, vérifie encore, comme dans le cas simplement connexe :

$$-\Delta \psi = \operatorname{div} \vec{u} \quad \text{dans } \Omega$$

$$\vec{\nabla} \psi = \vec{o} = \vec{u} \quad \text{sur } \Gamma_1 \cup \Gamma_2.$$

On peut choisir $\psi|_{\Gamma_1} = 0$ ; alors $\psi|_{\Gamma_2}$ est une constante inconnue. On approche $\psi$ par éléments finis linéaires conformes.
<figure><img src="image_27.png" /><figcaption>ROTATION DU CERCLE INTERIEUR : + 1
TRIANGULATION (2.1)
ν = 1.</figcaption></figure>
### 5.3. -

La géométrie est la même qu'en 5.2., mais le cercle intérieur reste immobile tandis que le cercle extérieur tourne à une vitesse angulaire $\omega = +1$.

Le nombre de Reynolds vaut donc :

$$ \frac{5 \times 2}{\nu} = \frac{10}{\nu} $$

On a utilisé la triangulation (2.3) :

412 triangles,
243 sommets,
655 noeuds

(triangulation obtenue par subdivision d'une triangulation initiale de 103 triangles).

#### Algorithme (III).

i) $\nu = .5$
Nombre de Reynolds = 20
$\frac{1}{\epsilon} = 10.$ , $\alpha = .05$ , $\rho = .1$
m tel que $C_h^{(m)} \ll 10^{-5}$ : 796.

ii) $\nu = .1$
Nombre de Reynolds = 100.
$\frac{1}{\epsilon} = 10.$ , $\alpha = .05$ , $\rho = .1$
m tel que $C_h^{(m)} \ll 10^{-5}$ : 2975
<figure><img src="image_28.png" /></figure>

TRIANGULATION

(2.2)

Rotation du cercle extérieur = +1

v = .5
<figure><img src="image_29.png" /></figure>

TRIANGULATION (2.2)

Rotation du cercle extérieur = +1.

v = .1
<figure><img src="image_30.png" /><figcaption>TRIANGULATION 2.3</figcaption></figure>

<figure><img src="image_31.png" /></figure>
### 5.4. - Écoulement entre deux plans.

Ce calcul est destiné à tester la précision de la méthode, la solution exacte étant connue.

Les deux plans sont supposés infinis, et ont pour équations

$$y = 0 \text{ et } y = 1.$$

Si l'on suppose donné un même profil des vitesses paraboliques à l'entrée et à la sortie, la solution exacte est indépendante de x :

$$ \begin{array}{ll} \text{avec} & u(x,y) = \begin{cases} \frac{g}{2v}y(1-y) & 0 \\ 0 & -g \\ \nabla p = & 0 \end{cases} \\ & y \text{ vertical, } p \text{ transversal} \end{array} $$

<figure><img src="image_32.png" /></figure>

$$ \text{Nombre de Reynolds : } \frac{g}{2v^2} $$
Résultats numériques avec l'Algorithme (III) :

On a pris $g = 2$, $v = \frac{1}{10}$ (Nombre de Reynolds = 100).

On a utilisé la triangulation (1.1) (128 triangles) et les
paramètres suivants :

$$
\frac{1}{\epsilon} = 5., \quad \rho = .75, \quad \alpha = 2.
$$

Le nombre m d'itérations tel que $c_h^{(m)} < 10^{-5}$ est 423 ; la
différence entre deux vitesses itérées $u_h^{(m)}$ et $u_h^{(m-1)}$ est inférieure
à $10^{-5}$ à partir de m = 486.

Soit $\vec{u}_h$ la solution obtenue et $\vec{u}$ la solution exacte.

On pose :

$$
\begin{align*}
e_i &= |\vec{u}(x_i, y_i) - \vec{u}_h(x_i, y_i)| \\
&= \text{erreur au noeud de numéro } i
\end{align*}
$$

$$
\begin{align*}
e_T &= |\vec{u}(x_B, y_B) - \vec{u}_h(x_B; y_B)| \\
&= \text{erreur au barycentre B du triangle T}
\end{align*}
$$

On a :

$$
\text{Max}_{i} e_{i} = 0.45, \text{ soit } 9\% \text{ par rapport au maximum de } |\vec{u}|
$$

$$
\left( \frac{\sum_i |e_i|^2}{\sum_i |\vec{u}(x_i, y_i)|^2} \right)^{1/2} = .12 \quad \text{soit } 12\%
$$

$$
\text{Max}_{T \in \mathcal{T}_h} e_T = .2 \quad \text{soit} \quad 4\%
$$

$$
\frac{\sum_{T \in \mathcal{T}_h} |e_T|^2}{\sum_T |\vec{u}(x_B, y_B)|^2} = 0.75 \text{ soit } 7,5\% .
$$

On voit que l'approximation de la vitesse est meilleure au barycentre
des triangles qu'aux noeuds. Pour comparer la pression calculée $p_h$ à la
pression exacte $p = -gx$, il convient de les normaliser en posant :
$$ \Pi_h = p_h - \int_{\Omega} p_h(x,y)dx dy $$

$$ \Pi = p - \int_{\Omega} p(x,y)dx dy $$

$$ \text{de sorte que } \int_{\Omega} \Pi_h dx dy = \int_{\Omega} \Pi dx dy = 0 $$

On compare $\Pi$ et $\Pi_h$ aux barycentres des triangles.

On trouve :

$$ \frac{\sum_{T} \text{aire}(T) |\Pi_h(x_B, y_B) - \Pi(x_B, y_B)|}{\sum_{T} \text{aire}(T) |\Pi(x_B, y_B)|} = .63. $$

### 5.5. Equations de Stokes.

Afin de tester la précision de la méthode d'éléments finis linéaires non conformes, on a également résolu le problème de Stokes :

$$ \left\{ \begin{array}{l} -\Delta \vec{u} = \vec{f} - \nabla \vec{p} \quad \text{dans } \Omega \\ \text{div } \vec{u} = 0 \quad \text{dans } \Omega \\ \vec{u}|_{\Gamma} = 0 \end{array} \right. $$

avec $\Omega = ]0,1[ \times ]0,1[$, soit :

$$ \left\{ \begin{array}{l} a(\vec{u}, \vec{v}) = (\vec{f}, \vec{v}) + (p, \text{div } \vec{v}) \quad \forall \vec{v} \in H_o^1(\Omega) \\ \text{div } \vec{u} = 0 \\ \vec{u} \in V \end{array} \right. $$

L'algorithme d'Uzawa classique a été utilisé :

$$ \left\{ \begin{array}{l} a(\vec{u}^{n+1}, \vec{v}) = (\vec{f}, \vec{v}) + (p^n, \text{div } \vec{v}) \quad \forall \vec{v} \in H_o^1(\Omega) \\ \vec{u}^{n+1} \in H_o^1(\Omega) \\ p^{n+1} = p^n - \rho \text{ div } \vec{u}^{n+1} \end{array} \right. $$

$\vec{f}$ est choisi de telle manière que la solution exacte $\vec{u}$ soit connue (voir ci-dessous).
NOMBRE DE TRIANGLES : 512 :

$h = 1/16$

Nombre total de sommets : 289.

Nombre de sommets intérieurs : 225.

Nombre total de points milieux : 800.

Nombre de points milieux intérieurs : 736.

$$ \frac{1}{\epsilon} = 0, \rho = .85. $$

Second membre : correspond à une solution exacte avec pression constante

$$ \psi = 256 \times x_1^2(1-x)^2 \times \frac{1}{2}(1-x_2)^2 : \psi(\frac{1}{2}, \frac{1}{2}) = 1 $$

$$ f_1 = - \Delta(\frac{\partial \psi}{\partial 2}) $$

$$ f_2 = - \Delta(-\frac{\partial \psi}{\partial 1}). $$

60 itérations - 4 minutes
(C11 10070)

<figure><img src="image_33.png" /></figure>

STOKES STATIONNAIRE - VISCOSITE = 0.1000E01.

$$ \psi_h(\frac{1}{2}, \frac{1}{2}) = .996 $$

Ereur relative sur les vitesses : $1,2 \times 10^{-2}$.
## CONCLUSION.

Le principal mérite de la méthode d'éléments finis qui vient d'être décrite est la simplicité de mise en oeuvre.

Son avantage sur des méthodes de différences finies réside dans la possibilité de traiter aisément des domaines à frontières courbes.
BIBLIOGRAPHIE

[1] A. AUSSEMS : "Une mise en oeuvre frontale des éléments finis".
Thèse de Docteur Ingénieur.
Université scientifique et médicale de Grenoble.
Décembre 1972.

[2] BORSENBERGER : "Programme de triangulation automatique".
Rapport interne, IRIA, 1972.

[3] J.F. BOURGAT : "Analyse numérique du problème de la torsion élastoplastique".
Thèse de 3ème cycle, Université Paris VI, 1971.

[4] M. CROUZEIX : "Etude d'une méthode de linéarisation. Résolution numérique des équations de Stokes stationnaires - Application aux équations de NAVIER - STOKES stationnaires".
Cahier IRIA n° 12, (mai 1974), pp 139-244.

[5] M. CROUZEIX et P.A. RAVIART : "Conforming and Non-conforming Finite Element Methods for solving the stationary Stokes", publié
au L. A. 189.

[6] D.J.F. EWING, A.J. FAWKES, J.R. GRIFFITHS : " Rules governing the number of nodes and elements in a finite element mesh".
Int. J. Num. Meth. Engrg., 2 (1970), pp 597-600.

[7] M. FORTIN : Thèse d'état, Université Paris VI, 1972.

[8] M. FORTIN, R. TEMAM, R. PEYRET : "Calcul des écoulements d'un fluide visqueux incompressible", Rapport IRIA, 1970 (INF. 7013).
[9] J. ALAN GEORGE : "Computer implementation of the finite element method"
Computer Science Dept., Standford University".
Février 1971 (STAN - CS - 71 - 208).

[10] R. GLOWINSKI, M. FORTIN, A. MARROCCO.
A paraître.

[11] HESTENES : "Multiplier and gradient methods" J. Optimization, Theory
and Appl., 4 (1969) pp 303-320.

[12] P. JAMET et P.A. RAVIART : "Numerical Solution of the Stationary
Navier - Stokes Equations by Finite Elements methods".
Colloque IRIA, décembre 1973. pp 103-133.

[13] R.T. ROCKAFELLAR. "Augmented Lagrange Multiplier Functions and Duality
in non convex Programming"

[14] R. TEMAM : "On the theory and numerical Analysis of the Navier-Stokes Equations". A paraître.

[15] R. VARGA : "Matrix iterative analysis".
Prentice - Hall, 1962.
Calcul "local" de la fonction de courant $\psi$

Au lieu d'utiliser l'équation – $\Delta\psi = \text{rot } u$, on remonte à la propriété
qui définit $\psi$ :

$$
\vec{\nabla} \psi = \begin{vmatrix} -u_2 \\ u_1 \end{vmatrix}
$$

Si $\vec{u}$ est linéaire entre les points A et B :

$$
\begin{gather*}
(1) \quad \psi(B) = \psi(A) + u_1 \left(\frac{A+B}{2}\right) \cdot (AB)_2 - u_2 \left(\frac{A+B}{2}\right) \cdot (AB)_1 \\
\qquad \left( \frac{A+B}{2} = \text{milieu de AB} \right)
\end{gather*}
$$

Soit $\vec{u} \in W_h$ (sur chaque $T \in \mathcal{T}_h$, div $\vec{u}=0$). La formule (1) permet de
calculer la fonction de courant $\psi$ correspondante (quadratique par morceaux)
en deux étapes ; prenant $\psi = 0$ sur la frontière extérieure :

i) de proche en proche on calcule $\psi$ en chaque sommet de la triangulation par (1) ; grâce à la condition d'incompressibilité div $\vec{u} = 0$, la valeur obtenue en un sommet A ne dépend pas du chemin d'intégration suivi de la frontière extérieure jusqu'à A.

ii) on calcule ensuite $\psi$ en un milieu de côté B, en intégrant par
(1) de A à B, où A est le sommet opposé à B dans l'un des deux triangles adjacents à B.

<figure><img src="image_34.png" /></figure>
On note que l'on obtient en B une autre valeur de $\psi$ en intégrant de A' à B (voir figure).

$\psi$ est donc discontinue le long des côtés (1) ; pour obtenir des lignes de courant sans discontinuité on pourra prendre en B la moyenne des deux valeurs obtenues dans T et T'.

**Remarque** : En pratique, si la contrainte $\text{div}_h \vec{u} = 0$ est traitée par dualité ou pénalisation-dualité, on obtient pour la solution $\vec{u}$ une divergence petite mais non nulle. Dans l'étape i) on obtiendra donc en un sommet A différentes valeurs de $\psi$ en suivant différents chemins d'intégration depuis la frontière extérieure ; on obtient une approximation satisfaisante en prenant en A la moyenne de toutes les valeurs trouvées en ce sommet.

Il reste à tracer les lignes de courant, qui sont les équipotentielles de $\psi$ (quadratique sur chaque triangle).

Le procédé le plus simple est de diviser chaque triangle en 4 sous-triangles et d'approcher $\psi$ sur chaque sous-triangle par une fonction affine.

On donne ci-dessous les résultats de cette méthode pour les problèmes 5.1 et 5.2 :

• problème 5.1, 512 triangles, $v = 10^{-2}$

• problème 5.2 : 412 triangles, $v = 10^{-1}$

(cette figure remplace la figure p. 62)

<footnote>(1) Ceci n'est pas surprenant : on montre que $\psi$ est aprochée par l'élément non conforme de Morley.</footnote>
APPENDICE AU RAPPORT LABORIA N° 150

Calcul "local" de la fonction de courant $\psi$

Au lieu d'utiliser l'équation – $\Delta\psi = \text{rot } u$, on remonte à la propriété
qui définit $\psi$ :

$$ \vec{\nabla} \psi = \begin{vmatrix} -u_2 \\ u_1 \end{vmatrix} $$

Si $\vec{u}$ est linéaire entre les points A et B :

$$ (1) \quad \begin{gathered} \psi(B) = \psi(A) + u_1 \left(\frac{A+B}{2}\right) \cdot (AB)_2 - u_2 \left(\frac{A+B}{2}\right) \cdot (AB)_1 \\ \left( \frac{A+B}{2} \right) = \text{milieu de AB} \end{gathered} $$

Soit $\vec{u} \in W_h$ (sur chaque $T \in \mathcal{T}_h$, div $\vec{u}=0$). La formule (1) permet de
calculer la fonction de courant $\psi$ correspondante (quadratique par morceaux)
en deux étapes ; prenant $\psi = 0$ sur la frontière extérieure :

i) de proche en proche on calcule $\psi$ en chaque sommet de la triangu-
lation par (1) ; grâce à la condition d'incompressibilité div$_h \vec{u} = 0$, la valeur
obtenue en un sommet A ne dépend pas du chemin d'intégration suivi de la fron-
tière extérieure jusqu'à A.

ii) on calcule ensuite $\psi$ en un milieu de côté B, en intégrant par
(1) de A à B, où A est le sommet opposé à B dans l'un des deux triangles adja-
cents à B.

<figure><img src="image_35.png" /></figure>
On note que l'on obtient en B une autre valeur de $\psi$ en intégrant de A' à B (voir figure).

$\psi$ est donc discontinue le long des côtés $\binom{1}{0}$ ; pour obtenir des lignes de courant sans discontinuité on pourra prendre en B la moyenne des deux valeurs obtenues dans T et T'.

**Remarque** : En pratique, si la contrainte $\text{div}_h \vec{u} = 0$ est traitée par dualité ou pénalisation-dualité, on obtient pour la solution $\vec{u}$ une divergence petite mais non nulle. Dans l'étape i) on obtiendra donc en un sommet A différentes valeurs de $\psi$ en suivant différents chemins d'intégration depuis la frontière extérieure ; on obtient une approximation satisfaisante en prenant en A la moyenne de toutes les valeurs trouvées en ce sommet.

Il reste à tracer les lignes de courant, qui sont les équipotentielles de $\psi$ (quadratique sur chaque triangle).

Le procédé le plus simple est de diviser chaque triangle en 4 sous-triangles et d'approcher $\psi$ sur chaque sous-triangle par une fonction affine.

On donne ci-dessous les résultats de cette méthode pour les problèmes 5.1 et 5.2 :

• problème 5.1, 512 triangles, $v = 10^{-2}$

• problème 5.2 : 412 triangles, $v = 10^{-1}$

(cette figure remplace la figure p. 62)

<footnote>$^{(1)}$ Ceci n'est pas surprenant : on montre que $\psi$ est aprochée par l'élément non conforme de Morley.</footnote>