RIA

laboria

Institut de Recherche
d'Informatique
et d'Automatique

Domaine de Voluceau
Rocquencourt
B.P. 5 78150 - Le Chesnay
France
Tél.: 954 90 20

laboratoire de recherche
en informatique
et automatique

ARCHIVES

RÉSOLUTION NUMÉRIQUE
DES ÉQUATIONS
DE NAVIER - STOKES
PAR DES ÉLÉMENTS FINIS
DE TYPE MIXTE

M. FORTIN

Rapport de Recherche N° 184

Août 1976
# RÉSOLUTION NUMÉRIQUE DES ÉQUATIONS DE NAVIER-STOKES
PAR DES ÉLÉMENTS FINIS DE TYPE MIXTE

M. Fortin

## Résumé :

On construit pour l'approximation des équations de Navier-Stokes un schéma tenant compte de la domination des termes convectifs aux grands nombres de Reynolds. On approche d'abord les équations d'Euler par une méthode d'éléments finis de discontinus, dans laquelle on utilise une dérivation décentrée dans le sens du courant. La viscosité est ensuite introduite par l'introduction d'une variable auxiliaire, approche déduite des méthodes mixtes pour l'approximation du problème biharmonique.

## Abstract :

A Finite Element approximation is built for the Navier-Stokes equations which takes care of the predominance of non linear convective terms at large Reynolds numbers. Inviscid Euler equations are first approximated by discontinuous finite elements, which permit to define an upwind differentiation formula. Viscosity is then introduced through an auxiliary equation, this approach being related to mixed finite element methods for plate bending problems.
RESOLUTION NUMERIQUE DES EQUATIONS DE NAVIER-STOKES
PAR DES ELEMENTS FINIS DE TYPE MIXTE

M. FORTIN

# 1. INTRODUCTION

La résolution numérique des équations de Navier-Stokes présente un grand nombre de difficultés dont les plus délicates sont liées au comportement des solutions lorsque le paramètre de viscosité devient petit par rapport à d'autres constantes de l'écoulement. De façon précise les équations de Navier-Stokes mises sous forme adimensionnelle s'écrivent : (Nous utilisons la convention de sommation habituelle de l'indice répété).

$$ (1.1) \qquad \frac{\partial u}{\partial t} - \frac{1}{\Re} \Delta u + u_i \frac{\partial u}{\partial x_i} + \text{grad } p = f $$

$$ (1.2) \qquad \text{div } u = \frac{\partial u_i}{\partial x_i} = 0 $$

Le paramètre Re, appelé nombre de Reynolds, peut généralement s'écrire

$$ (1.3) \qquad \mathrm{Re} = \frac{Ud}{\nu} $$

où U est vitesse de référence, d une longueur de référence et ν la viscosité cinétique. Dans beaucoup de situations physiques intéressantes, le nombre de Reynolds peut prendre des valeurs très élevées ($10^6$ par exemple). Cette situation se réflète dans les solutions par l'apparition de couches-limites au voisinage des parois, i. e. de zones à très fort
gradient de vitesse où la viscosité joue un rôle important, alors que dans
le reste de l'écoulement les effets visqueux sont négligeables. Ces couches
limites peuvent éventuellement décoller de la paroi ce qui ne simplifie,
bien sûr, pas les choses.

On assiste d'autre part aux grands nombres de Reynolds à l'apparition de
**solutions turbulentes** i.e. d'écoulements en apparence désordonnés et ins-
tables.

Dans chacun de ces cas, la taille des phénomènes physiques, qu'il importe de
simuler pour obtenir une description fidèle de l'écoulement, devient prohibi-
tivement petite par rapport à tout maillage réalisable, même sur les plus gros
ordinateurs existants ou même imaginables. Par ailleurs, il est clair que
toute méthode numérique qui ne prend pas explicitement en compte, la présence
éventuelle de ces phénomènes, ne peut que fournir des résultats fantaisistes.

Disons immédiatement que nous n'aborderons pas ici la question de la modéli-
sation des écoulements turbulents. Toutefois il est probable que les techni-
ques que nous décrirons seraient mieux adaptées que les techniques classiques
pour incorporer de tels modèles et nous tenterons de souligner au passage
les raisons qui nous conduisent à avancer cette hypothèse.

Pour tenter de décrire le fil directeur de notre démarche, nous prendrons
pour acquis, le fait non démontré, mais plausible, que les solutions
(laminaires, i.e. non turbulentes) des équations de Navier-Stokes, approxi-
ment aux grands nombres de Reynolds, les solutions des équations d'Euler
des fluides non-visqueux. Dans cette optique, une couche-limite dans une
solution des équations visqueuses correspond à la présence d'une discon-
tinuité dans la solution des équations l'Euler. Nous tenterons donc d'obtenir,
pour l'approximation des équations de Navier-Stokes des schémas, donnant dans
le cas limite v = 0, des schémas valables pour l'approximation des équations
d'Euler. Ce n'est pas le cas, signalons-le immédiatement, des schémas habi-
tuels par différences finies ou par éléments finis. Nous aurons donc besoin,
à la lumière de ce que nous avons décrit plus haut, de considérer des
**approximations discontinues**. Par ailleurs une approche basée sur des schémas
d'approximation mixte pour les équations linéaires de Stokes nous permettre d'introduire la viscosité comme un phénomène en quelque sorte indépendant. Nous obtiendrons ainsi des schémas numériques dans lesquels, les couches-limites seront modélisées par des discontinuités, la viscosité modulant d'importance de celles-ci par un phénomène de frottement.

## 2. APPROXIMATION DES EQUATIONS D'EULER

Nous rappellerons d'abord les formulations variationnelles classiques des équations de Navier-Stokes. Nous nous placerons pour simplifier l'exposé dans le cas des conditions aux limites de Dirichlet homogènes. Définissons donc dans Ω ouvert (borné) de $\mathbb{R}^n$ de frontière $\Gamma$.

$$ (2.1) \qquad V = \{v | v = (v_1, v_2, \dots, v_n) \mid v \in (H_o^1(\Omega))^n, \operatorname{div} v = 0\} $$

$$ (2.2) \qquad a(u,v) = \sum_{i,j=1}^{n} \int_{\Omega} \frac{\partial u_i}{\partial x_j} \frac{\partial v_j}{\partial x_i} dx $$

$$ (2.3) \qquad b(u,u,v) = \sum_{i,j=1}^{n} \int_{\Omega} u_j \frac{\partial u_i}{\partial x_j} v_i dx $$

On vérifie alors classiquement (cf. e.g. LIONS [1]), que la solution du problème variationnel

$$ (2.4) \qquad (u',v) + \mu a(u,v) + b(u,u,v) = (f,v) \quad \forall v \in V, \text{ p.p. en t.} $$

$$ (2.5) \qquad u(0) = u_o, $$

est solution du problème (1.1)-(1.2) avec la condition initiale (2.5) et les conditions de Dirichlet.

$$ (2.6) \qquad u|_{\Gamma} = 0 $$

Si n désigne la normale extérieure à $\Gamma$ et $\tau$ la vecteur unitaire tangent, la condition (2.6) se décompose évidemment en

$$ (2.7) \qquad u \cdot n|_{\Gamma} = 0 $$

$$ (2.8) \qquad u \cdot \tau|_{\Gamma} = 0 $$
La condition (2.8) est la condition d'adhérence à la paroi dans un fluide visqueux. L'approximation naturelle par éléments finis de ce problème (cf. CROUZEIX-RAVIART [1], FORTIN [1]) est donc de construire un espace approximant V.

Nous nous intéressons ici au comportement de la solution dans le cas limite $\mu = 0$, considérons d'abord le cas continu. Si on note

$$ (2.9) \qquad H = \{v | v \in (L^2(\Omega))^n, \operatorname{div} v = 0, u \cdot n|_T = 0\} $$

on montre facilement que dans certaines conditions ($f=0$ par exemple), la solution de (2.4)-(2.5) reste fornée dans $L^2(0,T;H)$. On peut donc espérer (mais le fait n'est pas démontré à ma connaissance), que la solution des équations de Navier-Stokes converge quand $\mu \to 0$ ver une solution des équations d'Euler correspondantes soit

$$ (2.10) \qquad \frac{\partial u}{\partial t} + u_i \frac{\partial u}{\partial x_i} + \operatorname{grad} p = f $$

$$ (2.11) \qquad \operatorname{div} u = 0 $$

$$ (2.12) \qquad u \cdot n|_T = 0, u(0) = u_o $$

On obtiendrait ainsi une solution des équations d'Euler plus faible que celles dont on sait démontrer l'existence (e.g. TEMAM [1]). La difficulté est de passer à la limite dans les termes non-linéaires, les arguments habituels de connaissance n'étant plus valables. Remarquons qu'une telle convergence n'est possible, si elle a lieu, que grâce à la densité de V dans H de sorte que la condition (2.8) peut à la limite disparaître. Dans le cas d'une approximation par éléments finis, ce n'est pas le cas, à moins de raffiner indéfiniment le maillage. Pour une triangulation donnée on ne peut donc pas espérer que les approximations classiques puissent donner, à $\mu = 0$, une approximation consistante des équations d'Euler, la raison principale tenant aux différences de condition aux limites. Nous allons donc tenter la démarche inverse, à savoir construire des approximations de l'équation d'Euler pour ensuite les modifier de façon à y inclure le cas échéant les effets visqueux.
Nous nous placerons dans le cas bidimensionnel et nous chercherons d'abord à construire des approximations de l'espace H. Rappelons d'abord qu'on peut écrire pour toute fonction $v \in H$,

$$ (2.13) \qquad v = \operatorname{rot} \psi = \left\{ \frac{\partial \psi}{\partial x_2} , - \frac{\partial \psi}{\partial x_1} \right\} , \psi \in H_o^1(\Omega) $$

Le résultat est classique pour les fonctions de V avec $\psi \in H_o^2(\Omega)$ et on déduit immédiatement (2.13) par densité.

L'approximation des fonctions de H peut donc se faire d'une façon très simple :

Soit $W_n$ une approximation de $H_o^1(\Omega)$ par des éléments finis conformes de type Lagrange, alors l'espace $H_h$ défini par

$$ (2.14) \qquad H_h = \{ v_h \mid v_h = \operatorname{rot} \psi_h, \psi_h \in W_h \} $$

est une approximation interne de H.

Nous supposerons par la suite que nous sommes dans le cas d'une telle approximation. On vérifiera aisément que les éléments de $H_h$ vérifient à l'interface de deux triangles des propriétés de continuité portant sur la *trace normale*. En effet si $\Gamma_K$ désigne la frontière du triangle K, on a $u \cdot n_K = \frac{\partial \psi}{\partial \tau_K}$ sur $\Gamma_K$ et la dérivée tangentielle est continue à l'interface par la conformité de l'approximation $W_h$ ; par contre on n'a pas de condition de raccord sur la trace tangentielle. En particulier les fonctions de $H_h$ vérifient à la frontière (2.7) mais non (2.8). Les espaces ainsi construits apparaissent donc comme les espaces minimaux pour l'approximation (interne) consistante des équations d'Euler. Notons au passage que ces espaces ont été considérés dans un cadre plus général dans RAVIART-THOMAS [1], pour l'approximation de $H(\operatorname{div};\Omega)$ où

$$ (2.15) \qquad H(\operatorname{div};\Omega) = \{ v \mid v \in (L^2(\Omega))^n, \operatorname{div} v = 0 \}. $$

Nous allons donc chercher maintenant à construire un principe variationnel permettant de résoudre de façon approchée les équations d'Euler.
Remarquons qu'il n'est plus possible d'utiliser directement (2.4).
En effet les espaces que nous avons construits ne sont pas assez
réguliers pour pouvoir définir les termes de la forme $u_i \frac{\partial u_j}{\partial x_i}$.
Nous allons utiliser, pour contourner cette difficulté, une idée
introduite par LESAINT [1], (cf. également LESAINT-RAVIART [1])
pour l'approximation des équations hyperboliques du premier ordre
par des éléments finis discontinus. La technique ainsi utilisée a
été reprise dans FORTIN [2] dans un cadre général d'approximations
externes par éléments finis.

Notons d'abord que l'on peut en utilisant la condition de divergence nulle réécrire la formulation variationnelle (2.4) sous une forme légèrement différente (dite conservative). (Nous faisons ici $\mu = 0$)

$$
(2.16) \quad (u',v) + \sum_{i,j=1}^{n} \int_{\Omega} \frac{\partial}{\partial x_j} (u_i u_j) v_i dx = (f,v) \quad \forall v.
$$

Soit donc $T_h$ une triangulation de $\Omega$, sur laquelle on a construit un espace $H_h$. Nous allons chercher à affaiblir la formulation (2.16) de façon à la rendre compatible avec notre approximation. Pour ce remarquons qu'on a en intégrant par parties sur un triangle $K$ de $T_h$ (nous utilisons la convention de sommation)

$$
(2.17) \quad \int_K \frac{\partial}{\partial x_j} (u_i u_j) v_i dx = - \int_K (u_i u_j) \frac{\partial v_i}{\partial x_j} + \int_{\partial K} u \cdot n (u_i v_i) d\Gamma
$$

L'idée essentielle du schéma de Lesaint est de distinguer dans la
frontière de K une partie $\partial K_-$ sur laquelle le flux est entrant
$(u \cdot n < 0)$ et une partie $\partial K_+$ sur laquelle il est sortant $(u \cdot n > 0)$.

Considérant le caractère hyperbolique du système, il est naturel
en raisonnant sur les caractéristiques de considérer que dans K,
l'écoulement n'est influencé par les triangles voisins qu'à travers
$\partial K_-$.

Remarquons que $u \cdot n$ est continu sur $\partial K$ et que la définition de $\partial K_+$ et de $\partial K_-$ est consistante d'un triangle à l'autre. Notant alors
sur l'interface S, entre le triangle K et un triangle voisin $K_1$

$$ (2.18) \quad u^+|_S = \begin{cases} \text{trace de } u \text{ dans } K \text{ sur } \partial K_+ \\ \text{trace de } u \text{ dans } K_1 \text{ sur } \partial K_- \end{cases} $$

Nous prendrons comme formulation variationnelle semi-discrete de l'équation d'Euler.

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

$$ (2.19) \quad (u',v) - \sum_K \left( \int_K (u_i u_j) \frac{\partial v_i}{\partial x_j} dx + \int_K u \cdot n (u_i^+ v_j^-) d\Gamma \right) = \int_{\Omega} f_i v_i dx \quad \forall v \in H_h^1 $$

En intégrant de nouveau par parties le terme non-linéaire on obtiendrait la forme équivalente

$$ (2.20) \quad (u',v) + \sum_K \left( \int_K \frac{\partial}{\partial x_j} (u_i u_j) v_i dx + \int_{\partial K_-} u \cdot n (u_i^{\text{ext}} - u_i^{\text{int}}) v_i d\Gamma \right) = \int_{\Omega} f_i v_i dx $$

où $u^{\text{ext}}$ est évidemment la valeur de $u$ dans le triangle adjacent et $u^{\text{int}}$ la valeur de $u$ dans le triangle K. Cette formulation consiste essentiellement à introduire des dérivations discrètes décentrées dans le sens du courant, ce qui est classique pour l'approximation des systèmes hyperboliques. On peut considérer que le terme de saut dans (2.20) exprime la présence d'une masse de Dirac aux interfaces, les fonctions $u$ étant discontinues. Lesaint a montré la convergence de telles approximations dans le cas des systèmes de Friedrichs linéaires. Le cas non-linéaire traité ici est évidemment beaucoup plus difficile à manipuler. Nous allons nous contenter ici de montrer que dans ce schéma, les termes de saut sont dissipatifs.

Nous allons d'abord, d'après une idée de Raviart, considérer une version légèrement généralisée de ce schéma permettant de contrôler l'aspect décentré (et comme nous le verrons, dissipatif) de ce schéma. Sur une interface S nous définissons toujours comme en (2.18), $u^+$
comme étant la trace de u dans le triangle pour lequel S correspond à
∂K+, et u^- la valeur dans le triangle
adjacent. La situation est résumée sur
la figure où les flèches représentent
le sens du courant (u·n s'annulle au
point A).

Nous considérons alors le schéma semi-
discret

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

$$
(2.21) \quad (u',v) - \sum_K \left( \int_K (u_i u_j) \frac{\partial v_i}{\partial x_j} dx + \int_{\partial K} u \cdot n \left( (1-\alpha) u_i^- + \alpha u_i^+ \right) v_i d\Gamma \right) = (f,v)
$$

On vérifiera que pour $\alpha = \frac{1}{2}$, on obtient des dérivations discrètes
"centrées" et que pour $\alpha < \frac{1}{2}$ les dérivations sont décentrées à contre-
courant. Le cas $\alpha=1$ correspond à (2.19) ou (2.20).

Faisons $v=u$ dans (2.21). On obtient alors

$$
(2.22) \quad \left\{
\begin{aligned}
& \frac{1}{2} \frac{d}{dt} |u|_L^2 + \sum_K \left( -\frac{1}{2} \int_{\partial K} u \cdot n (u_1^+ u_1 + (1-\alpha) u_1^- u_1) \right. \\
& \qquad \left. + \alpha u_2^+ u_2 + (1-\alpha) u_2^- u_2 \right] d\Gamma
\end{aligned}
\right\} = (f, u)
$$

On a noté dans (2.22), $u_1$ et $u_2$ les traces de u intérieures au triangle
K. Sur la partie sortante $\partial K_+$ on $u_i = u_i^+$ d'où la décomposition des inté-
grales de bord

$$
(2.23) \quad \left\{
\begin{aligned}
& (\alpha - \frac{1}{2}) \int_{\partial K_+} u \cdot n \left[ (u_1^+)^2 + (u_2^+)^2 \right] d\Gamma + (1-\alpha) \int_{\partial K_+} \left[ u \cdot n \cdot u_1^+ u_1^- + u_2^+ u_2^- \right] d\Gamma \\
& + (1-\alpha - \frac{1}{2}) \int_{\partial K_-} (u \cdot n) \left[ (u_1^-)^2 + (u_2^-)^2 \right] d\Gamma + \alpha \int_{\partial K_-} u \cdot n \cdot u_1^+ u_1^- + u_2^+ u_2^- d\Gamma
\end{aligned}
\right.
$$

Sommons sur tous les triangles et considérons une interface S. Notant
que dans deux triangles adjacents la normale est de signe opposé et
que les parties ∂K' et ∂K+(resp ∂K_ et ∂K'_+) se correspondent mutuel-
lement, on a si on retient la normale n orientée de sorte que u·n soit
positive :
$$
(2.24) \quad \left( \frac{2\alpha-1}{2} \right) \int_S u \cdot n \left( |u_1^+ - u_1^-|^2 + |u_2^+ - u_2^-|^2 \right) d\Gamma
$$

Reportant dans (2.22) on a

$$
(2.25) \quad \frac{1}{2} \frac{d}{dt} |u|^2 + \frac{(2\alpha-1)}{2} \sum_S \int_S (|u_1^+ - u_1^-|^2 + |u_2^+ - u_2^-|^2) u \cdot n \, d\Gamma = (f, n)
$$

Dans le cas f=0 par exemple, on aura dissipation (décroissance de u
en norme L²) pour α > ½. Pour α = ½ le schéma est conservatif et pour
α < ½ il sera a priori instable. Remarquons un fait important ; la
dissipation n'intervient que dans le cas u·n ≠ 0. Il sera donc possible
d'améliorer le schéma en ajustant interactivement le maillage de façon
à ce que les interfaces présentant de fortes discontinuité soient
orientées dans le sens du courant. Par ailleurs on montrerait facilement
que la dissipation n'est fonction que du saut de la composante tangen-
tielle de u, la composante normale étant continue.

En conclusion sous réserve d'une _discrétisation en temps_ convenable,
on peut par les techniques que nous venons de décrire obtenir pour les
équations d'Euler des schémas stables et consistants. L'étude précise
de l'erreur est pour le moment un problème ouvert. En ce qui concerne
le temps, nous aurions pu effectuer le raisonnement qui précède sur
des éléments finis discontinus espace-temps (cf. LESAINT-RAVIART [1],
JAMET [1].)

Nous aurions ainsi obtenu des schémas de Ruge-Kutta implicites
(CROUZEIX [1]), possédant de très bonnes qualités de stabilité mais
difficiles à manipuler. Nous effectuons actuellement la mise en oeuvre
par des schémas à plusieurs pas avec prédicteur et correcteur. Remar-
quons au passage que le caractère non-linéaire des équations implique
que tous schéma implicite devra comporter une part quelconque d'itéra-
tion. Par ailleurs, si on repère les fonctions de $H_h$ par les degrés
de liberté de la fonction de courant associée dans $v_h$, on remarquera
que la résolution du schéma explicite le plus simple,

$$
(2.26) \quad (u^{n+1}, v) = (u^n, v) - \Delta t \quad \tilde{b}(u^n, u^n, v) \quad \forall v \in H_h
$$
(b(u^n, u^n, v) représente les termes non-linéaires traités comme en (2.21))
nécessité la résolution d'un problème de Dirichlet en ψ. En effet
si u^{n+1} = rot ψ^{n+1} et v = rot ψ, (2.26) est de la forme

$$
(2.27) \quad \int_{\Omega} \operatorname{rot} \psi^{n+1} \operatorname{rot} \psi \,dx = (F, \psi) \quad \forall \psi \in W_h^\bullet.
$$

■

Nous allons maintant voir comment on peut greffer à un tel schéma
des termes de viscosité en utilisant des méthodes d'éléments finis
mixtes.

3. DISCRETISATION DES EQUATIONS DE NAVIER-STOKES

Le principe directeur pour passer du cas des équations d'Euler à
celui des fluides visqueux, sera de traiter les termes du second
ordre au moyen d'un schéma mixte compatible avec la discrétisation
utilisée pour les équations d'Euler. Pour ce faire la première
 étape consiste essentiellement à se ramener à un système du premier
ordre par l'introduction d'une inconnue supplémentaire soit :

$$
(3.1) \qquad \frac{\partial u_i}{\partial t} - \mu \frac{\partial p_{ij}}{\partial x_j} + u_j \frac{\partial u_i}{\partial x_j} + \text{grad } p = f
$$

$$
(3.2) \qquad \qquad \qquad \qquad \qquad p_{ij} = \frac{\partial u_i}{\partial x_j}
$$

qu'il s'agit maintenant d'écrire dans une forme variationnelle conve-
nable. Nous multiplions donc (3.1) par une fonction test à divergence
nulle ce qui donne

$$
(3.2) \quad (u',v) - \mu \left( \frac{\partial p_{ij}}{\partial x_j} , v_i \right) + b(u,u,v) = (f,v)
$$

équivalent de 2.4. La discrétisation brutale de (3.2) demanderait des
p_{ij} discrets dont la dérivée soit dans L^2(Ω). Tout comme nous avons
affaibli la formulation des termes non linéaires, nous allons modifier
le terme ( \frac{\partial p_{ij}}{\partial x_j} , v_i ) de façon à utiliser des espaces discrets possédant
moins de régularité. Nous utiliserons pour y arriver une formulation
mixte introduite dans le cas des problèmes du quatrième ordre par
HERMANN [1] et étudiée au niveau des propriétés de convergence dans
JOHNSON [1] et BREZZI-RAVIART [1]. Si on note toujours $n$ et $\tau$ les vecteurs unitaires normal et tangent à $\partial K$ et si l'on pose

$$ (3.3) \qquad p_{nn} = p_{ij} n_i n_j $$

$$ (3.4) \qquad p_{n\tau} = p_{ij} n_i \tau_j $$

une formulation possible sera alors

$$ (3.5) \quad \left\{ \begin{aligned} & \left( u_i^!, v_i \right) + \sum_K \left\{ \int_K \frac{\partial p_{ij}}{\partial x_j} v_i \, dx - \int_{\partial K} v \cdot n \, p_{nn} \, d\Gamma \right\} \\ & + \sum_K \left\{ - \int_K u_i u_j \frac{\partial v_i}{\partial x_j} \, dx + \int_{\partial K} u \cdot n \left[ \alpha u_i^+ + (1-\alpha) u_i^- \right] v_i \, d\Gamma \right\} \\ & = \int_{\Omega} f_i v_i , \quad \forall v \in H_h \end{aligned} \right. $$

$$ (3.6) \quad \int_{\Omega} p_{ij} q_{ij} \, dx = \sum_{K} \left\{ - \int_{K} u_i \frac{\partial q_{ij}}{\partial x_j} + \int_{\partial K} u \cdot n \, q_{nn} \, d\Gamma \right\}, \quad \forall q \in Q_h $$

ou encore

$$ (3.7) \quad \int_{\Omega} p_{ij} q_{ij} \, dx = \sum_{K} \left\{ \int_{K} \frac{\partial u_i}{\partial x_j} q_{ij} - \int_{\partial K} u \cdot \tau q_{n\tau} \, d\Gamma \right\}, \quad \forall q \in Q_h $$

où l'on a

$$ (3.8) \quad Q_h = \{ p_{ij} \mid \sum_i p_{ii} = 0, \frac{\partial p_{ij}}{\partial x_j} \in P_k(K) . \text{ "p}_{n\tau} \text{ continu}\} . $$

**Remarque 3.1** : Les éléments de $Q_h$ sont donc des tenseurs de trace nulle, non symétriques, dont la restriction à tout triangle K est pour chaque composante un polynôme de degré k. On impose de plus aux interfaces la continuité de la trace $p_{n\tau}$. On vérifiera (cf. BREZZI-RAVIART [1]) que l'on peut effectivement construire de tels espaces en choisissant les degrés de liberté de la façon suivante :
Soit K un triangle de $T_h$, les éléments de $\phi_h$ sont déterminés par

i) la donnée sur chacun des côtés de K de $k+1$ valeurs de $p_{n\tau}$ ou encore des moments d'ordre $\le k$ de $p_{n\tau}$

ii) les moments d'ordre $\le k-1$ de chacune des composantes. ■

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

En ce qui concerne la précision d'un tel schéma on peut donner une réponse complète dans le cas du problème linéarisé de Stokes stationnaire. Dans ce cas, si les fonctions à divergence nulle de $H_h$ sont aussi de degré $k$ sur chaque triangle et, par conséquent, à des fonctions de courant de degré $k+1$, on montre que l'erreur sur $u$ et sur $p$ est en norme $L^2(\Omega)$, d'ordre $h^{k+1}$. L'erreur $L^2$ en $u$ correspond évidemment à une erreur sur $\psi$ en norme $H_o^1$. En regard des résultats habituels sur les méthodes classiques d'éléments finis, on a donc atteint un résultat optimal. Le cas linéaire évolutif ne présente pour sa part aucune difficulté spéciale par rapport aux cas classiques (équation de la chaleur par exemple).

**Remarque 3.2.** : Les équations (3.6) ou (3.7) apparaissent comme un système linéaire auxiliaire à résoudre parallèlement à (3.5). Si on utilise comme on l'a suggéré au §2 un schéma en temps du type prédicteur correcteur, on aura donc pour chaque itération deux systèmes linéaires à résoudre. Toutefois ce sont des systèmes symétriques et définis positifs dont on peut effectuer la factorisation une fois pour toutes au début du calcul. Seule l'expérience pourra nous dire si la précision de ces schémas est suffisante pour justifier leur lourdeur apparente.

**Remarque 3.3.** : Il semblerait plus naturel d'utiliser dans les formulations mixtes de type (3.5)-(3.6) le tenseur des contraintes $\sigma_{ij}$ comme variable duale. En particulier une telle formulation présentait de grands avantages pour le traitement des problèmes à frontières libres. La difficulté est alors de construire des espaces de tenseurs
**symétriques.** On vérifie en particulier que la construction de la remarque 3.1. n'est plus possible. Par contre si on impose seulement que les moments d'ordre $\le k-1$ de $\sigma_{12}$ et de $\sigma_{21}$ soient égaux dans chaque triangle K, on obtient des espaces pour lesquels on peut démontrer des majorations d'erreur en $h^k$ (et non plus en $h^{k+1}$).

**Remarque 3.4** : Il semble aussi que dans le cas de formulations hybrides, il soit possible de construire des approximations utilisant le tenseur des contraintes sans perdre l'ordre de précision optimal. Ce procédé est plus ou moins équivalent à utiliser des approximations non-conformes pour les $\sigma_{ij}$. L'étude détaillée et a fortiori la mise en oeuvre reste à faire.

**Remarque 3.5** : La simulation des écoulements turbulents se fait généralement par l'introduction de termes supplémentaires (non linéaires) de dissipation ou encore par l'utilisation d'une loi de comportement modélisant les effets turbulents. L'introduction de tels modèles peut se faire dans le schéma proposé en modifiant de façon appropriée le second membre de (3.6) ou de (3.7). En particulier le schéma peut s'appliquer à tout fluide non newtonien.

En conclusion, disons que l'utilisation des éléments finis mixtes, semble ouvrir pour l'approximation des équations de l'hydrodynamique des perspectives nouvelles qui sont encore largement à explorer. En particulier, la possibilité d'obtenir des schémas valables à la limite pour les équations d'Euler devrait permettre d'effectuer des calculs valables pour des nombres de Reynolds plus élévés que par les méthodes classiques. On doit toutefois se garder d'un trop grand optimisme tant que les résultats numériques ne seront pas venus confirmer des espoirs encore fragiles.
REFERENCES

BONNEROT R., JAMET P. [1] - *A second order finite element method for the one dimensional Stefan problem*.

Int. J. Num. Meth. Eng., 8, (1974), 811-820.

BREZZI F., RAVIART P.A. [1] - *A paraître*.

CROUZEIX M. [1] - *Sur l'approximation des équations différentielles opérationnelles linéaires par des méthodes de Runge-Kutta*.

Thèse, Paris VI, 1975.

CROUZEIX M., RAVIART P.A. [1] - *Conforming and non-conforming finite element methods for solving the stationary Stokes equations*.

RAIRO.

FORTIN M. [1] - *Calcul numérique des écoulements des fluides de Bingham et des fluides newtoniens incompressibles par la méthode des éléments finis*.

Thèse Paris VI, 1972.

FORTIN M. [2] - *Approximation des équations de l'hydrodynamique par des méthodes d'éléments finis*.

Séminaire de Mathématiques Supérieures, Montréal 1975, à paraître.

HERMANN L.R. [1] - *Finite element bending analysis for plates*.

J. Eng. Mech. Div. ASCE, EM5, 93 (1967) 49-83.

JOHNSON C. [1] - *On the convergence of some mixed finite element methods in plate bending problems*.

Num. Math. 21 (1973) 43-62.

LESAINT P. [1] - *Sur la résolution des systèmes hyperboliques du premier ordre par des méthodes d'éléments finis*.

Thèse, Paris VI, 1975.

LESAINT P., RAVIART P.A. [1] - *On a finite element method for solving the Newton transport equation*.

Math. Aspect of Finite El. in Partial Diff. Equations. Academic Press, 1974.

RAVIART P.A., THOMAS J.M. [1] - *A mixed Finite Element Method for 2nd order elliptic problems*.

A paraître.

TEMAM R. [1]