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

HAL
open science

# Calcul numérique de l'écoulement avec ou sans sillage autour d'un profil bidimensionnel symétrique et sans incidence

Jean-François Bourgat, Georges Duvaut

► To cite this version:

Jean-François Bourgat, Georges Duvaut. Calcul numérique de l'écoulement avec ou sans sillage autour d'un profil bidimensionnel symétrique et sans incidence. [Rapport de recherche] IRIA-RR-145, IRIA. 1975, pp.33. hal-04716076

HAL Id: hal-04716076

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

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
IRIA
laboria

Institut de Recherche
d'Informatique
et d'Automatique

Domaine de Voluceau
Rocquencourt
78150 - Le Chesnay
Tél. 954.90.20

laboratoire de recherche
en informatique
et automatique

F

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

Calcul numérique
de l'écoulement
avec ou sans sillage
autour d'un profil
bidimensionnel
symétrique et
sans incidence

Jean-François

George

J.F. Bourgat, G. Duvaut

Rapport de Recherche n°145
Décembre 1975
CALCUL NUMÉRIQUE DE L'ÉCOULEMENT AVEC OU SANS SILLAGE
AUTOUR D'UN PROFIL BIDIMENSIONNEL SYMÉTRIQUE ET SANS INCIDENCE

J.F. Bourgat*, G. Duvaut**

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

RÉSUMÉ :

Formulé dans le plan de l'hodographe, le problème de l'écoulement bidimensionnel d'un fluide
parfait autour d'un profil symétrique sans incidence conduit à la résolution d'une inéquation varia-
tionnelle elliptique. On met ainsi en évidence une famille à un paramètre de solutions avec ou sans
sillage. De plus, l'étude débouche sur une méthode rapide et simple de calcul numérique de la solu-
tion. Les résultats obtenus pour quelques profils types sont indiqués et comparés à des travaux
antérieurs.

ABSTRACT :

Formulated in the hodograph plane, the bidimensional flow of a perfect fluid past a symmetric body lead to an elliptic variational inequality. By this way we show the existence of a one-parameter family of solutions. Furthermore the study gives us a simple and effective method for calculating the solutions. The results obtained with some typical profiles are given and compared with previous works.

<footnote>* IPIA-LABORIA. Domaine de Voluceau. Rocquencourt. 78150 Le Chesnay</footnote>

<footnote>** Département de Mécanique. Université Paris VI. 9 quai St-Bernard
75005 Paris</footnote>
## INTRODUCTION :

Le problème étudié est celui de l'écoulement bidimensionnel avec ou sans sillage d'un fluide parfait incompressible (*) autour d'un profil convexe symétrique et sans incidence. L'étude théorique du problème est faite au n°1 et conduit à la résolution d'une inéquation variationnelle elliptique dans le plan de l'hodographe. Elle met en évidence une famille à un paramètre de solutions allant de l'écoulement sans sillage à l'écoulement avec sillage se prolongeant jusqu'à l'infini. La méthode utilisée est particulièrement bien adaptée au calcul numérique de la solution présentée au n°2. Ce calcul numérique, réalisé par une méthode de relaxation utilise la géométrie simple du plan de l'hodographe et fournit directement une représentation paramétrique du sillage dans le plan physique. Nous donnons au n°3 les résultats obtenus sur plusieurs profils types. Comparée aux études antérieures sur ce problème (cf. LIGTHILL [12] et SOUTHWELL and VAISEY [14]), nous pensons que la présente méthode numérique, basée sur une description théorique rigoureuse de la classe de solutions, constitue un progrès substantiel du fait de la simplicité et de la rapidité de mise en oeuvre et d'exécution du calcul.

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

<footnote>(*) La méthode s'adapte également au cas des fluides compres-sibles comme il est indiqué dans H.BREZIS et G.DUVAUT [4].</footnote>
# 1 - FORMULATION DU PROBLEME.

## 1.1 Problème physique.

Nous proposons d'étudier l'écoulement plan, stationnaire, irrotationnel d'un fluide parfait incompressible autour d'un profil P convexe, symétrique et sans incidence. On cherche un écoulement avec ou sans sillage comme dans H. Brézis et G. Duvaut [4]

<figure><img src="image_6.png" /><figcaption>Figure N°1</figcaption></figure>

Écoulement avec sillage autour d'un demi-profil.

On désigne par (S) la frontière libre du sillage qui présente la même symétrie que P par rapport à l'axe x'x colinéaire à la direction de l'écoulement à l'infini amont de composantes ($q_∞, o$) dans un système d'axes x'Åx, y'Åy où A est le bord d'attaque. Par raison de symétrie on ne considère que la région $y \geq o$. (fig.1).
La fonction de courant $\psi(x,y)$ vérifie :

$$
(1) \quad \Delta\psi = 0 \quad \text{sur l'extérieur D du profil et du sillage}
$$

$$
(2) \quad \psi = 0 \quad \text{sur } x'A, P, S,
$$

$$
(3) \quad |\operatorname{grad} \psi| = q_S \quad \text{sur } S,
$$

où $q_S$ est une constante positive donnée avec $q_S \le q_\infty$.

$$
(4) \quad \operatorname{grad} \psi \rightarrow (0, q_{\infty}) \quad \text{quand } x^2 + y^2 \rightarrow +\infty, (x,y) \in D.
$$

On suppose que la ligne S a une ordonnée décroissante quand $x$ croit et se raccorde au profil en aval du maitre couple. Le problème ainsi posé va être étudié par passage au plan de l'hodographe.

1.2 Plan de l'hodographe (fig. 2).

Nous effectuons la transformation (x,y) → (θ,q) où :

$$
(5) \quad q e^{i\theta} = \frac{\partial \psi}{\partial y} - i \frac{\partial \psi}{\partial x} = \frac{\partial \varphi}{\partial x} + i \frac{\partial \varphi}{\partial y}
$$

si $\varphi(x,y)$ est la fonction potentiel.

Nous avons alors :

$$
(6) \quad \begin{cases} d\varphi = q \cos \theta dx + q \sin \theta dy \\ d\psi = -q \sin \theta dx + q \cos \theta dy \end{cases}
$$

d'où il vient que :

$$
(7) \quad dx + i dy = \frac{e^{i\theta}}{q} (d\varphi + i d\psi)
$$
<figure><img src="image_7.png" /><figcaption>Figure N°2</figcaption></figure>

Représentation de l'écoulement dans le plan de l'hodographe.

Il est commode d'introduire la variable $\sigma$ par

$$ (8) \qquad q = q_{\infty} e^{-\sigma} $$

Considérant alors $\varphi$ et $\psi$ comme des fonctions de $\theta$ et $\sigma$ et exigeant que le second membre soit une différentielle totale, on obtient

$$ (9) \qquad \frac{\partial \varphi}{\partial \sigma} = \frac{\partial \psi}{\partial \theta} ; -\frac{\partial \psi}{\partial \sigma} = \frac{\partial \varphi}{\partial \theta} $$

D'où il résulte en éliminant $\varphi$

$$ (11) \qquad \frac{\partial^2 \psi}{\partial \theta^2} + \frac{\partial^2 \varphi}{\partial \sigma^2} \equiv \Delta \psi = 0 \quad \text{sur } \mathcal{H} = T(D) $$
où T désigne la transformation faisant passer du plan physique au plan $(\theta, \sigma)$. Soit $\Gamma = T(P)$ l'image du profil, d'équation $\sigma = 1(\theta)$. Les conditions aux limites sur $\psi$ se présentent de la manière suivante :

$$
(12) \quad \psi(o, \sigma) = o \quad \text{pour } \sigma \ge o
$$

$$
(13) \quad \psi(\theta, \sigma) = 0 \quad \text{sur } \Gamma
$$

$$
(14) \quad \frac{\partial \psi}{\partial \sigma} = - \frac{q_{\infty} e^{-1(\theta)}}{1+1'^2} R(\theta), \quad \frac{\partial \psi}{\partial \theta} = -1' \frac{\partial \psi}{\partial \sigma} \quad \text{sur } \Gamma
$$

(15) $\psi(\theta, \sigma_S) = 0$, $\theta_S \le \theta \le 0$, où $\sigma_S = -\log \frac{q_S}{q_\infty}$, et où $\theta_S$ est l'angle polaire de la tangente commune au sillage et au profil au point de raccordement. On a noté $R(\theta)$ le rayon de courbure algébrique (négatif) du profil au point où sa tangente fait l'angle orienté $\theta$ avec $\vec{Ox}$, les valeurs extrêmes de $\theta$ étant $\theta_o$ et $\theta_l$ avec $\theta_o > 0$, $\theta_l < 0$. Insistons sur le fait que $q_S$ est donné alors que $\theta_S$ et $\Gamma$ sont à priori inconnus.

1.3 Changement de fonction inconnue.

Suivant les idées de C. Baiocchi [1], nous introduisons la
fonction u(θ,σ) par

$$
(16) \quad u(\theta, \sigma) = \frac{1}{q_{\infty}} e^{-\int_{1(\theta)}^{\sigma} e^{\tau} \psi(\theta, \tau) d\tau}
$$

Cette fonction vérifie :

$$
(17) \quad \Delta u = - R(\theta) e^{-\sigma} \quad \text{dans } \mathcal{B}
$$

$$
(18) \quad u = 0, \operatorname{grad} u = 0 \quad \text{sur } \Gamma
$$

$$
(19) \quad u(o, \sigma) = H e^{-\sigma} \quad \text{pour } \sigma > o
$$

$$
(10) \quad u + \frac{\partial u}{\partial \sigma} = 0 \quad \text{sur} \quad \theta_S \le \theta \le o, \quad \sigma = \sigma_S \\
\text{où H désigne une constante.}
$$
On a alors :

Propriété 1.
-----------

La constante H est égale à la demi-hauteur du profil.

Démonstration.
----------

Considérons la fonction $u_1(\theta, \sigma)$ définie par

$$
(21) \qquad u_1(\theta, \sigma) = u(\theta, \sigma) e^{\sigma}
$$

qui satisfait :

$$
\[
\begin{array}{ll}
(22) \quad \left\{
\begin{array}{@{}l@{\quad}l@{}}
\Delta u_1 - 2 \frac{\partial u_1}{\partial \sigma} + u_1 = -R & \text{dans } \mathcal{B} \\[1em]
u_1 = H & \text{sur } \theta = 0, \sigma \ge 0 \\[1em]
\frac{\partial u_1}{\partial \sigma} = 0 & \text{sur } \sigma = \sigma_s, \theta_s \le \theta \le 0 \\[1em]
u_1 \ge 0 & \text{dans } \mathcal{B} \\[1em]
u_1 = 0, \text{ grad } u_1 = 0 & \text{sur } \Gamma
\end{array}
\right.
\end{array}
\]
$$

Posons :

$$
(23) \quad w(0) = \lim_{\sigma \to +\infty} u_1(\theta, \sigma)
$$

qui doit satisfaire :

$$
\[
(24) \quad \left\{
\begin{array}{@{}l@{\quad}l@{}}
& \displaystyle \frac{d^2 w}{d\theta^2} + w = -R(\theta) \\[\medskipamount]
& w(0) = H \\[\medskipamount]
& \displaystyle \frac{dw}{d\theta}(\theta_0) = w(\theta_0) = 0 \quad (\text{car } \Gamma \text{ est asymptote à } \theta = \theta_0)
\end{array}
\right.
\]
$$

($\Gamma$ est asymptote à $\theta = \theta_o$)
On vérifie alors par résolution explicite de (24) que H
doit être égale à la demi-hauteur du profil.

Propriété 2. (H. Brézis [3]).

Sous l'hypothèse :

$$
(25) \quad \rho = \min_{\theta \in [\theta_o, \theta_1]} |R(\theta)| > 0
\quad \text{On a}
$$

$$
(26) \quad u(\theta, \sigma) = 0 \quad \text{pour } \sigma < \sigma_o < 0
$$

où σ₀ est l'unique solution négative de

$$
(27) \quad (1 + \sigma_o) e^{-\sigma_o} = 1 - \frac{H}{\rho}
$$

Démonstration.

On considère la bande :

$$
(28) \quad G = \{ (\theta, \sigma) \mid \theta_1 < \theta < \theta_o, \sigma < 0 \}
$$

Dans G on a

$$
\[
(29) \quad \left\{
\begin{array}{ll}
u \ge 0 & \\
-\Delta u = R(\theta) e^{-\sigma} & \text{si } u > 0 \\[1em]
\text{grad } u \text{ continue dans } G & \\
u(\theta_1, \sigma) = u(\theta_o, \sigma) = 0 & \\
u(\theta_o, 0) \le H & \forall \theta
\end{array}
\right.
\]
$$

la dernière inégalité résultant du principe du maximum dans 𝒩,
où $\Delta u \ge 0$ : le maximum de u est alors atteint sur $\partial \beta$ ce qui entraine l'inégalité.

Définissons alors sur G une fonction $\varphi(\sigma)$, $\sigma \le 0$ telle que

$$ (30) \quad \left\{ \begin{array}{l} \varphi(0) = H \\[1em] \sigma \to \frac{d\varphi}{d\sigma}(\sigma) \text{ continue,} \\[1em] \varphi \ge 0 \\[1em] \frac{d^2\varphi}{d\sigma^2} = \rho e^{-\sigma}, \text{ si } \varphi > 0. \end{array} \right. $$

Le principe du maximum implique que

$$ (31) \qquad u \le \varphi $$

Par ailleurs il est facile de trouver explicitement $\varphi(\sigma)$ :

$$ (32) \qquad \frac{d^2\varphi}{d\sigma^2} = \rho e^{-\sigma} $$

implique

$$ (33) \qquad \varphi(\sigma) = \rho e^{-\sigma} + A\sigma + B $$

La condition $\varphi(0) = H$ entraine $B = H - \rho$. On choisit alors $A$ tel qu'il existe $\sigma_o < 0$ satisfaisant :

$$ (34) \qquad \varphi(\sigma_o) = 0, \quad \frac{d\varphi}{d\sigma}(\sigma_o) = 0 $$

c'est à dire :

$$ (35) \quad \left\{ \begin{array}{l} \rho e^{-\sigma_o} + A\sigma_o + H - \rho = 0 \\[1em] - \rho e^{-\sigma_o} + A = 0 \end{array} \right. $$

Le scalaire $\sigma_o$ est donc bien donné par (27) et (31) implique (26).
1.4 Cadre fonctionnel et résultats.

Introduisons l'ouvert Ω par

$$
(36) \quad \Omega = \{ (\theta, \sigma) \mid \theta_1 < \theta < \theta_o ; \sigma > \sigma_o ; \sigma < \sigma_s \text{ pour } \theta < \theta_o \},
$$

l'espace V par

$$
(37) \quad V = \{ w \mid w \in H^1(\Omega) \ ; w(\theta_o, \sigma) = w(\theta_1, \sigma) = w(\theta_o, \sigma_o) = o \}
$$

qui est un espace de Hilbert quand on le munit de la structure in-
duite par la norme

$$
(38) \quad \| w \| = \left( \int_{\Omega} (\operatorname{grad} w)^2 \cdot dx \right)^{1/2},
$$

et l'ensemble convexe fermé de V,

$$
(39) \quad K_H = \{ w \mid w \in V, w(o, \sigma) = H e^{-\sigma} \text{ pour } \sigma \ge o; w \ge o \}.
$$

Posons par ailleurs

$$
(40) \quad a(u,w) = \int_{\Omega} \operatorname{grad} u \cdot \operatorname{grad} w \, dx + \int_{\Gamma_{\partial O}} u w \, d\theta
$$

où

$$
(41) \quad \Gamma_o = \{ (\theta, \sigma) \mid \sigma = \sigma_\xi, \theta_1 \le \theta \le \theta_o \}.
$$

Il est alors clair que si u satisfait (17)-(20), on a :

$$
(42) \quad u \in K_H
$$

Formons alors, ∀ v ∈ K_H

$$
(43) \quad \int_{\mathcal{B}} -\Delta u \cdot (v - u) \, d\theta \, d\sigma = \int_{\mathcal{B}} R(\theta) e^{-\sigma} (v - u) \, d\theta \, d\sigma.
$$
Mais nous avons :

$$ \int_{\mathcal{B}} R(\theta) e^{-\sigma} (v-u) d\theta d\sigma \geq \int_{\Omega} R(\theta) e^{-\sigma} (v-u) d\theta d\sigma $$

car $R(\theta) < 0$ et $u = 0$ sur $\Omega - \mathcal{B}$. Par ailleurs ;

$$ \begin{aligned} \int_{\mathcal{B}} -\Delta u \cdot (v-u) d\theta \, l\sigma &= \int_{\partial \mathcal{B}} - \frac{\partial u}{\partial n} (v-u) d\Gamma + \int_{\mathcal{B}} \text{grad } u \cdot \text{grad}(v-u) d\theta d\sigma \\ &= \int_{\Omega} \text{grad } u \cdot \text{grad}(-u) d\theta \, d\sigma + \int_{\Omega} u \cdot (v-u) d\Gamma, \end{aligned} $$

de sorte que (43) fournit :

$$ (44) \quad a(u,v-u) \geq \int_{\Omega} R(\theta) e^{-\sigma} (v-u) d\theta \, d\sigma, \quad \forall v \in K_H $$

### Lemme 1.

Si l'application

$\theta \to R(\theta)$ est de carré intégrable sur $(\theta_1, \theta_o)$, alors

$$ (45) \quad v \to \int_{\Omega} R(\theta) e^{-\sigma} v \, d\theta \, l\sigma \text{ est une forme linéaire continue sur V.} $$

### Démonstration.

Sous l'hypothèse (45) l'application $(\theta, \sigma) \to R(\theta)e^{-\sigma}$ est de carré intégrable sur $\Omega$.

Par conséquent

$$ \left| \int_{\Omega} R(\theta) e^{-\sigma} v(\theta, \sigma) d\theta \, d\sigma \right| \le C \left( \int_{\Omega} v^2 d\theta \, d\sigma \right)^{1/2} $$

Par ailleurs nous avons :

$$ (46) \quad v(\theta, \sigma) = \int_{\theta_o}^{\theta} \frac{\partial v}{\partial \theta} \cdot \theta' \, , \sigma \, d\theta' $$
d'où

$$
(47) \quad \begin{cases} |v(\theta, \sigma)| \le |\theta_1| \dots |\theta_n|^{1/2} \left( \int_0^\theta \left( \frac{\partial v}{\partial \theta} \right)^2 d\theta \right)^{1/2} & \text{si } \sigma_0 < \sigma \le \sigma_s \\ |v(\theta, \sigma)| \le |\theta_0|^{1/2} \left( \int_0^\theta \left( \frac{\partial v}{\partial \theta} \right)^2 d\theta \right)^{1/2} & \text{si } \sigma > \sigma_s \end{cases}
$$

Par conséquent

$$
(48) \quad \begin{aligned}[t] \int_{\Omega} v^2(\theta, \sigma) d\theta \, d\sigma &= \int_{\sigma_0}^{\sigma} \int_{\theta_1}^{\theta} v^2 d\theta \, d\sigma + \int_{\sigma_0}^{+\infty} \int_{\theta_0}^{\theta} v^2 d\theta \, d\sigma \\ &\le \int_{\sigma_0}^{\sigma_s} (\theta_1 - \theta_0)^2 \int_{\theta}^{\theta_1} \left(\frac{\partial v}{\partial \theta}\right)^2 d\theta \, d\sigma + \epsilon_0^2 \int_{\sigma_s}^{+\infty} \int_{\theta_0}^{\theta} \left(\frac{\partial v}{\partial \theta}\right)^2 d\theta \, d\sigma \\ &\le (\theta_1 - \theta_0)^2 \int_{\Omega} \left(\frac{\partial v}{\partial \theta}\right)^2 d\theta \, d\sigma \le (\theta_1 - \theta_0)^2 \int_{\Omega} \|\text{grad } v\|^2 d\theta \, d\sigma \end{aligned}
$$

ce qui achève la démonstration du lemme ...

On a alors le

Théorème 1.

Sous les hypothèses (25) et (45), il existe un unique élé-
ment u satisfaisant (42) et (44).

Démonstration.

La forme bilinéaire $a(u,v)$ est évidemment coercive. Par ailleurs elle est continue car, d'après le lemme 1, la norme sur $V$ est équivalente à la norme induite par $H^1(\Omega)$. Le second membre de (44) étant une forme linéaire continue sur $V$, nous sommes alors dans le cas d'application d'un théorème d'existence classique pour les inéquations variationnelles (cf. J.L. Lions et G. Stampacchia [6]).

Théorème 2.

Le sillage est donné par la représentation paramétrique :
$$
(49) \quad \begin{cases}
x = x_o(\theta) + \frac{q_{\infty}}{q_S} [u(\theta, \sigma_S) \sin \theta + \frac{\partial u}{\partial \theta}(\theta, \sigma_S) \cos \theta] \\
y = y_o(\theta) + \frac{q_{S0}}{q_S} [-u(\theta, \sigma_S) \sin \theta + \frac{\partial u}{\partial \theta}(\theta, \sigma_S) \cos \theta]
\end{cases}
$$

où $x_o(\theta)$, $y_o(\theta)$ sont les coordonnées du point du profil où la tan-
gente fait l'angle $\theta$ avec Ox et où $\theta$ varie de $\theta_S$ à 0.

Si $q_S < q_\infty$, le sillage se referme à distance finie.

Si $q_S = q_\infty$, le sillage se prolonge à l'infini aval et sa largeur tend vers zéro.

Démonstration.

Le bord du sillage est l'image dans le plan physique du segment de Γ₀ d'abscisse $\epsilon \in [\epsilon_S, \infty)$, et d'ordonnée $\sigma = \sigma_S$. Nous avons alors, d'après (7),

$$
(50) \quad dx + i dy = \frac{e^{i\theta}}{q} d\varphi = \frac{e^{i\theta}}{q} \left( \frac{\partial \varphi}{\partial \theta} d\theta + \frac{\partial \varphi}{\partial \sigma} d\sigma \right).
$$

Mais il résulte de (8) et (10) que :

$$
\frac{\partial \varphi}{\partial \theta} = - \frac{\partial \varphi}{\partial \sigma} , \quad \frac{\partial \varphi}{\partial \sigma} = \frac{\partial \psi}{\partial \theta}
$$

d'où

$$
\mathrm{d}x + i \mathrm{d}y = \frac{e^{i\theta}}{q_S} \left( -\frac{\partial\psi}{\partial\sigma} \mathrm{d}\theta \right)
$$

puisque d σ = o.

Par ailleurs (16) fournit

$$
(51) \quad \frac{\partial u}{\partial \sigma} = -u + \frac{1}{q_{\infty}} \psi(\theta, \sigma)
$$

d'où

$$
\frac{\partial \psi}{\partial \sigma} = q_{\infty} \left( \frac{\partial^2 u}{\partial \sigma^2} + \frac{\partial u}{\partial \sigma} \right) = q_{\infty} (-R(\theta) e^{-\sigma} - \frac{\partial^2 u}{\partial \theta^2} + \frac{\partial u}{\partial \sigma})
$$
Sur $\Gamma_o$, cette égalité donne :

$$ \frac{\partial \psi}{\partial \sigma} = - q_{\infty} (R(\theta) e^{-\sigma S} + \frac{\partial^2 u}{\partial \theta^2} + u) $$

et par conséquent

$$ (52) \quad dx + idy = \frac{q_{\infty}}{q_S} \left( \frac{\partial^2 u}{\partial \theta^2} + u \right) e^{i\theta} d\theta + R(\theta) e^{i\theta} d\theta. $$

Si nous désignons par

$$ x_o = x_o(\theta), \quad y_o = y_o(\theta) $$

la représentation paramétrique du profil, nous avons

$$ (53) \quad dx_o + idy_o = e^{i\theta} ds = e^{i\theta} \frac{ds}{d\theta} d\theta = e^{i\theta} R(\theta) d\theta $$

où s représente l'abscisse curbiligne sur le profil. Nous avons donc

$$ (54) \quad dx + i dy = dx_o + i dy_o + \frac{q_{\infty}}{q_S} \frac{d}{d\theta} \left( \frac{\partial u}{\partial \theta} (\theta, \sigma_S) e^{i\theta} - i u(\theta, \sigma_S) e^{i\theta} \right) $$

soit (49) après intégration.

Les résultats généraux de régularités des solutions des inéquations variationnelles (cf. H. Brézis [ 2 ]) montrent que u est de classe $C^{1,\alpha}$ en dehors du segment $\theta = 0, \sigma_S \ge \sigma \ge 0$. Au point (0,0), u est seulement höldérien et aux points du segment $\theta = 0, \sigma_S \ge \sigma > 0$ les restrictions de u à $\theta > 0$ et $\theta < 0$ sont de classe $C^{1,\alpha}$. Il résulte de ces résultats que si $\sigma_S > 0$, soit $q_S < q_{\infty}$ le point $(x+iy)(\theta)$ tend vers un point à distance finie sur l'axe $\vec{0}x$ lorsque $\theta$ tend vers zéro.

Si $\sigma_S = 0$, c'est à dire $q_S = q_{\infty}$, $y(\theta)$ tend vers zéro lorsque $\theta$ tend vers zéro, $x(\theta)$ tend vers $+\infty$, car $u(\theta, \sigma)$ est seulement höldérien en (0,0).
## 2 - APPROXIMATION.

### 2.1 Approximation de $\Omega$ par un domaine borné.

Considérons le problème $\Omega_1$ obtenu en bornant $\sigma$ inférieurement par $\sigma_1 < \sigma_0$ ($\sigma_0$ défini en 27) et supérieurement par $\sigma_2 > \sigma_s$.

<figure><img src="image_8.png" /><figcaption>Figure N°3 : le domaine $\Omega_1$</figcaption></figure>

La frontière de $\Omega_1$ est composée des segments $\Gamma_o, \Gamma_1, \dots, \Gamma_5$ comme l'indique la figure 3.

Nous noterons $\Gamma'_5$ le segment défini par

$$ \Gamma'_5 = \{(0, \sigma), 0 \le \sigma \le \sigma_2\} $$
La solution u du problème (44) vérifie

$$
(55) \quad u = 0 \text{ sur } \Gamma_2 \text{ d'après la propriété 2 du paragraphe 1.3}
$$

$$
(56) \quad u + \frac{\partial u}{\partial n} = \frac{1}{q_{\infty}} \psi \quad \text{sur } \Gamma_4 \text{ d'après (51)}
$$

Introduisons

$$
(57) \quad K_1 = \{ v \mid v \in H^1(\Omega_1), \begin{array}{l} v \ge 0, v = 0 \text{ sur } \Gamma_1, \Gamma_2, \Gamma_3, \\ v(o, \sigma) = H e^{-\sigma} \text{ sur } \Gamma_5' \end{array} \}
$$

et

$$
(58) \quad a_1(u,v) = \int_{\Omega_1} \operatorname{grad} u \cdot \operatorname{grad} v \, d\theta \, d\sigma + \int_{\Gamma_o} uv \, d\theta + \int_{\Gamma_4} uv \, d\theta
$$

Le problème :

Trouver $u_1 \in K_1$ tel que

$$
(59) \quad a_1(u_1, v-u_1) \ge \int_{\Omega_1} R(\theta) e^{-\sigma} (v-u_1) d\theta \, d\sigma \quad \forall v \in K_1
$$

admet une solution unique puisque $a_1$ est une forme bilinéaire, symé-
trique, continue, coercive sur $H^1(\Omega_1)$ et que $K_1$ est un convexe fer-
mé de $H^1(\Omega_1)$.

Mais d'après (55) et (56), $u|_{\Omega_1}$ est solution unique de l'i-
néquation variationnelle

Trouver $w \in K_1$ tel que

$$
(60) \quad a_1(w,v-w) \ge \int_{\Omega_1} R(\theta) e^{-\sigma} (v-w) d\theta \, d\sigma + \frac{1}{q_{\infty}} \int_{\Gamma_4} \phi(v-w) d\theta \quad \forall v \in K_1
$$

Or $\phi \to 0$ lorsque $\sigma \to +\infty$ ce qui, par suite de la continuité de la solution par rapport au second membre, entraine

$$
u_1 \to u|_{\Omega_1} \quad \text{quand} \quad \sigma_2 \to +\infty
$$
Nous résoudrons donc le problème (59) qui s'interprète par

$$
\[
\begin{array}{@{}l@{\qquad}l@{}}
(61) &
  \left\{
    \begin{array}{@{}l@{}}
      -\Delta u_1 = - R(\theta) e^{-\sigma} \text{ dans } \Omega_1 \cap \mathcal{B} \\
      u_1 = 0 \text{ dans } \Omega_1 - (\Omega_1 \cap \mathcal{B}) \\
      u_1 = 0, \text{ grad } u_1 = 0 \text{ sur } \Gamma \cap \Omega_1 \\
      u_1 = 0 \text{ sur } \Gamma_1, \Gamma_2, \Gamma_3 \\
      u_1 = He^{-\sigma} \text{ sur } \Gamma_5' \\
      u_1 + \frac{\partial u_1}{\partial n} = 0 \text{ sur } \Gamma_o \text{ et } \Gamma_4
    \end{array}
  \right\}
\end{array}
\]
$$

2.2 Approximation du problème (59) par éléments finis
__d'ordre 1.__

Pour obtenir un maillage assez fin au voisinage de Γ₀
nous choisisons d'approcher le problème (59) par éléments finis.
D'autre part la contrainte v ≥ o suggère de prendre des éléments finis
linéaires car il suffit alors de vérifier la contrainte aux noeuds.

On trouvera dans ZIENKIEWICZ [16], PAVIART [13],
CIARLET-RAVIART [8], GLOWINSKI-LIONS-TREMOLIERES [10], des études
approfondies sur les éléments finis.

a) *Formulation du problème approché.*

On se donne une triangulation $\mathcal{L}_h$ de $\Omega_1$, avec h égal au
plus grand côté de tous les triangles de $\mathcal{L}_h$, telle que

$$
(62) \quad \bigcup_{T \in \mathcal{L}_h} T = \bar{\Omega}_1
$$
$$
(63) \quad \left\{
\begin{array}{l}
\text{Si } T_1 \text{ et } T_2 \in \mathcal{T}_h \text{ on a, ou } T_1 \cap T_2 = \emptyset \text{ ou} \\
\quad T_1 \text{ et } T_2 \text{ ont un côté commun ou} \\
\quad T_1 \text{ et } T_2 \text{ ont un sommet commun}
\end{array}
\right.
$$

On définit

$$
(64) \quad \omega_h = \{ P/P \text{ sommet de } T \in \mathcal{T}_h \}
$$

$$
(65) \quad \gamma_{ih} = \{ P/P \in \omega_h, P \in \Gamma_i \}, i = 0, \dots, 5
$$

$$
(66) \quad \gamma'_{5h} = \{ P/P \in \omega_h, P \in \Gamma_5' \}
$$

$$
(67) \quad V_h = \{ v_h / v_h \in \mathcal{C}^0(\bar{\Omega}_1) \cap H^1(\bar{\Omega}_1) : v_h \text{ linéaire sur } T, \forall T \in \mathcal{T}_h \}
$$

$$
(68) \quad K_h = \{ v_h \in V_h / v_h \ge 0, v_h(P) = He^{-c(P)} \text{ pour } P \in \gamma_5'h, \\
\phantom{(68) \quad K_h = \{ v_h \in V_h / v_h \ge 0, v_h(P) = He^{-c(P)} \text{ pour } } v_h(P) = 0 \text{ pour } P \in \gamma_{1h} \cup \gamma_{2h} \cup \gamma_{3h} \}
$$

Si l'on désigne par $w_h^P$, $P \in \omega_h$, les fonctions de base définies par

$$
(69) \quad w_h^P \in V_h, w_h^P(P) = 1, w_h^P(M) = 0 \text{ pour } M \neq P
$$

on a pour tout $v_h \in V_h$

$$
(70) \quad v_h = \sum_{P \in \omega_h} v_h(P) w_h^P
$$

On vérifie aisément que $V_h$ est un sous espace de $H^1(\Omega_1)$
donc un espace de Hilbert pour la norme induite et que $K_h$ est un convexe
fermé de $V_h$.

Par suite le problème approché relatif au problème (59)
que nous définissons par

$$
(71) \quad \left\{
\begin{array}{l}
\text{Trouver } u_h \in K_h \text{ tel que} \\
\quad a_1(u_h, v_h - u_h) \ge \int_{\Omega_1} R(\theta) e^{-\sigma} (v_h - u_h) d\theta d\sigma \quad \forall v_h \in K_h
\end{array}
\right.
$$

admet une solution unique.
b) <u>Convergence des solutions approchées</u> $u_h$ <u>lorsque</u> $h \to 0$.

On considère une famille $\mathcal{T}_h$ de triangulations telles que
(72) $$ \{
$$ Les angles de tous les triangles $T \in \mathcal{T}_h$ sont bornés inférieurement, uniformément en $h$, par $\theta_o > 0$.

On a alors le

Théorème.

Si les triangulations $\mathcal{T}_h$ vérifient (62), (63), (72), la
solution $u_h$ du problème (71) converge, lorsque $h \to 0$, vers la solution
$u_1$ du problème (59) dans $H^1(\Omega_1)$.

Démonstration.

Il suffit de vérifier (cf. GLOWINSKI-LIONS-TREMOLIERES [10])
que

i) les $v_h$ définissent une approximation interne de $H^1(\Omega_1)$,
c'est à dire que l'on peut trouver un espace $V$ dense dans $H^1(\Omega_1)$ tel
que $\forall v \in V$ on puisse construire $v_h \in V_h$ avec $v_h \to v$ dans $H^1(\Omega_1)$.

Ceci est démontré dans CIARLET-RAVIART [ 8 ] et
GLOWINSKI-LIONS-TREMOLIERES [ 10 ] avec $V = \mathcal{B}(\bar{\Omega})$ et $v_h = \sum_{P \in \omega_h} v(P) w_h^P$.

ii) les ensembles $K_h$ approchent $K$ au sens suivant :

- $\mathcal{B}(\bar{\Omega}) \cap K$ est dense dans $K$

- $\forall v \in \mathcal{B}(\bar{\Omega}) \cap K$ on peut trouver $v_h \in K_h$ avec $v_h \to v$ dans $H^1(\Omega_1)$.

Il suffit de prendre $v_h = \sum_{P \in \omega_h} v(P) w_h^P$ et l'on a $v_h \to v$ dans $H^1(\Omega_1)$
d'après i) et $v_h \in K_h$.
## 2.3 Résolution du problème approché (71)

Pour résoudre le problème (71) nous utilisons l'algorithme de surrelaxation ou sous-relaxation pour contraintes locales étudié dans GLOWINSKI [9] et CEA-GLOWINSKI [7].

### a) Description de l'algorithme.

Soit J la fonctionnelle définie par

$$ (75) \qquad J(v) = a(v,v) - 2(f,v), \quad v \in \mathbb{R}^N, $$

$$ (76) \qquad \text{où } a \text{ est une forme bilinéaire, continue, symétrique et coercive sur } \mathbb{R}^N \text{ et } f \in \mathbb{R}^N. $$

Soit $K \subset \mathbb{R}^N$ un convexe de type local, i.e.

$$ (77) \qquad K = \bigcap_{i=1}^{N} K_i \quad \text{où } K_i = \{ v_i \in \mathbb{R} \mid a_i \le v_i \le b_i \} $$

$a_i$ et $b_i$ pouvant éventuellement devenir infinis.

Posons $a_{ij} = a(w_i, w_j)$ et $f_i = (f, w_i)$, où $w_i$, $i=1, \dots, N$, est la base canonique de $\mathbb{R}^N$.

La fonctionnelle J admet un minimum unique sur le convexe fermé K que l'on recherche par l'algorithme suivant :

on se donne $u^o \in K$ et $\omega > 0$

on suppose $u^n$ connu et l'on détermine $u_i^{n+1}$, $i=1, \dots, N$

en trois étapes.

$$ u_i^{n+1/3} = - \frac{1}{a_{ii}} \left( \sum_{j=1}^{i-1} a_{ij} u_j^{n+1} + \sum_{j=i+1}^{N} a_{ij} u_j^N - f_i \right) $$

$$ u_i^{n+2/3} = \omega u_i^{n+1/3} + (1-\omega) u_i^n $$

$$ u_i^{n+1} = P_i(u_i^{n+2/3}) $$

où $P_i$ désigne la projection sur le convexe $K_i$.
On démontre dans GLOWINSKI [9] que si $0 < \omega < 2$, on a,
lorsque $n \to +\infty$

$u^n \to u$ dans $\mathbb{R}^N$

b) Application au problème (71)

Posons :

$$ \omega_{hl} = \{ P \in \omega_h / P \in \gamma_{1h} \cup \gamma_{2h} \cup \gamma_{3h} \cup \gamma'_{5h} \} $$

$N$ = nombre d'éléments de $\omega_{hl}$

$$ v = \sum_{P \in \omega_h} v(P) \cdot w_h^P $$

$$ (78) \qquad J(v) = a_1(v,v) - 2 \int_1^\Omega R(\theta) e^{-\sigma} v d\theta d\sigma + 2 \sum_{M \in \gamma'_{5h}} \mathrm{He}^{-\sigma(M)} a_1(w_h^M, v) $$

$$ (79) \qquad K = \bigcup_{P \in \omega_{hl}} K_p \quad \text{avec} \quad K_p = \{ v(P) \in \mathbb{R} | v(P) \ge 0 \} $$

J et K vérifiant les hypothèses (76) et (77), nous pou-
vons appliquer l'algorithme décrit en a).
3 - APPLICATIONS NUMERIQUES.

3.1 Choix des paramètres.

Nous avons traité différents exemples en faisant varier
le profil P autour duquel se fait l'écoulement ainsi que la vitesse
q_s sur le bord du sillage.

Les profils sont formés d'arcs de cercle raccordant des
angles θ₀ et θ₁ donnés.

- Nous considérons d'abord trois demi-profil s symétriques
(i.e. θ₀ = -θ₁), formés chacun d'un arc de cercle de flèche H=1, et
nous prenons pour valeurs de θ₀ : θ₀ = π/6, π/4, π/2.

<figure><img src="image_9.png" /><figcaption>Figure n°4 demi-profil symétrique</figcaption></figure>

Le rayon de courbure algébrique est R(θ) = -1/(1-cos θ₀).

- Nous traitons ensuite le cas d'un demi-profil non symétrique formé d'un quart de cercle de rayon 1 se raccordant (raccord

<figure><img src="image_10.png" /><figcaption>Figure N°5 demi-profil non symétrique</figcaption></figure>
de classe $C^1$) à un arc de cercle tel que : $\theta_o = \pi/2$, $\theta_l = -\pi/4$.

Nous obtenons dans ce cas

$$H=1$$
$$R(\theta) = \begin{cases} -1 & \text{pour } 0 \le \theta \le \pi/2 \\ -2-\sqrt{2} & \text{pour } -\pi/4 \le \theta < 0 \end{cases}$$

- Les vitesses sur le bord du sillage sont choisies à partir de $\sigma_s$. Nous prendrons successivement

$$\sigma_s = 0, 0.2, 0.4$$

ce qui correspond approximativement à

$$q_s = q_\infty, \quad q_s = \frac{5}{6} q_\infty, \quad q_s = \frac{2}{3} q_\infty$$

- On vérifie dans tous les cas que la valeur de $\sigma_o$ donnée par la formule (27) est telle que $\sigma_o \ge -1$ ; on pourra donc prendre pour $\Omega_l$ le domaine $\Omega$ borné dans la direction $\sigma$ par

$$\sigma_1 = -1, \sigma_2 = 4$$

## 3.2 Triangulation.

Les triangulations $\mathcal{L}_h$ du domaine $\Omega_l$ sont obtenues à partir d'une triangulation initiale que l'on subdivise successivement en joignant les milieux des côtés des triangles.

On vérifie que h, plus grand côté de tous les triangles converge vers o et que la propriété (72) est satisfaite.

La figure 6 représente la triangulation initiale choisie. La densité des triangles a été augmentée dans la région de $\Gamma_o$ car la détermination du sillage dépend des valeurs de la solution et de sa dérivée sur $\Gamma_o$. D'autre part, d'après les résultats de régularité du paragraphe 1.4, on sait que la solution u est de classe $C^1$ partout sauf sur le segment $\{\theta = 0, \sigma_s \ge \sigma \ge 0\}$ où elle est continue mais non dérivable (cf. fig. 13), aussi a-t-on appuyé la triangulation sur ce segment pour approcher au mieux la solution dans cette région.
<figure><img src="image_11.png" /><figcaption>Figure N°6
Triangulation initiale</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure N°7 : Triangulation pour $\theta_o = \frac{\pi}{2}$, $\theta_1 = -\frac{\pi}{4}$, $\sigma_s=0,4$</figcaption></figure>
Les essais numériques sont faits sur une triangulation
obtenue après 4 divisions de la triangulation initiale, (cf. figure 7)
ce qui nous donne

1088 triangles

581 noeuds (éléments de $\omega_h$)

533 variables (éléments de $\omega_{h1}$)

### 3.3 Résultats de l'algorithme.

Nous utilisons l'algorithme de surrelaxation décrit au paragraphe 1.3 après avoir, au préalable, fait une recherche du paramètre optimal de relaxation $\omega_{opt}$, qui optimise la rapidité de convergence de l'algorithme. Il s'agit en fait de rechercher le paramètre optimal du problème sans contraintes. C'est justifié dans GLOWINSKI [9]. Pour une description de l'algorithme qui recherche $\omega_{opt}$, cf. VARGA [12] chapitre 9.

Nous choisissons pour test d'arrêt des itérations

$$ \left( \frac{\sum_{i=1}^{N} |u_i^{n+1} - u_i^n|}{\sum_{i=1}^{N} |u_i^n|} \right) < 10^{-5} $$

Nous obtenons alors pour $\theta_o = -\theta_1 = \pi/6$ et $\sigma_s = 0$ :

$$ \omega_{opt} = 1.78 $$

nombre d'itérations = 49

Temps d'exécution sur CII-IRIS 80 : 1 mn

$$ \text{et pour } \theta_o = -\theta_1 = \frac{\pi}{2}, \sigma_s = 0.4 $$

$$ \omega_{opt} = 1.68 $$

nombre d'itérations = 32

Temps d'exécution 40 s
Tous les autres cas donnent des résultats situés entre
ces bornes.

3.4 Détermination de la ligne Γ :

La ligne Γ, image du profil par la transformation de
l'hodographe, est la ligne qui sépare la région de Ω₁ où u > 0 et la
région où u = 0. Elle permet de connaître les vitesses sur la ligne
d'eau du profil et surtout l'angle θs qui détermine le point où le
sillage se raccorde au profil.

On trace la ligne Γh en remarquant que, si les trois som-
mets T₁, T₂, T₃ d'un triangle sont tels que u_h(T₁)=u_h(T₂)=0 et
u_h(T₃) > 0, alors le segment T₁T₂ appartient à la ligne Γh.

On obtient ainsi une ligne brisée qui approche Γ et que
l'on pourra lisser pour une meilleure présentation des résultats
(cf. figure 8).

On trouvera sur les figures 9, ..., 12 les lignes $\Gamma_h$
obtenues pour les différents cas traités.

Remarque 2.2.

Sur les figures 9,..., 12 on voit que le point du profil
amont où la vitesse est égale à q∞ (i.e. c = 0) est indépendant de σs
pour un profil donné. Par contre le point du profil aval où la vitesse, après avoir augmenté, repasse par q∞, se déplace vers l'aval
quand σs augmente.

3.5 Tracé du sillage.

Le sillage est défini par la formule (49).

On trace le sillage en joignant linéairement les points
(x(θ), y(θ)) pour (θ, σs) ∈ γoh.

La dérivée $\frac{\partial u}{\partial u}$ pourra être approchée de façon simple car
la solution u est de classe $C^1, \alpha$ (cf. résultats de régularité du pa-
ragraphe 1.4). On trouvera sur la figure 13 une coupe de la solution
<figure><img src="image_13.png" /><figcaption>Figure n°8 - Lignes isovaleurs de la fonction u pour $\theta_o = \pi/2$, $\theta_1 = -\pi/4$ et lissaage de la ligne $\Gamma_h$ (isovalleur o)</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Figure N°13 : Coupe de u suivant la direction $\Gamma_o$ pour $\theta_o = -\theta_1 = \pi/4$</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>Figure N°9 : Tracé des lignes Γ pour θ₀ = -θ₁ = π/6</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Figure N°10 : Tracé des lignes l' pour θ₀ = -θ₁ = π/4</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Figure N°11 : Tracé des lignes Γ pour $\theta_o = -\theta_1 = \frac{\pi}{2}$.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Figure N°12 : Tracé des lignes Γ pour $\theta_o = \frac{\pi}{2}$ et $\theta_1 = -\frac{\pi}{4}$.</figcaption></figure>
selon le segment $\Gamma_o$.

Nous prendrons donc, puisque les segments de la frontière $\Gamma_o$ sont égaux

$$ \frac{\partial u}{\partial \theta}(M) = \frac{u(M_2) - u(M_1)}{\theta(M_2) - \theta(M_1)} \quad \underset{M_1}{\text{---}} \text{---} \underset{M_2}{\text{---}} $$

et si M est une extrémité de $\Gamma_o$.

$$ \frac{\partial u}{\partial \theta}(M) = \frac{3u(M) - 4u(M_1) - u(M_2)}{\theta(M) - \theta(M_2)} \quad \underset{M_2}{\text{---}} \text{---} \underset{M_1}{\text{---}} \text{---} \underset{M}{\text{---}} $$

On trouvera sur les figures 14,..., 17 les tracés des sillages pour les différents cas traités.

### Remarque 2.3.

On peut observer sur les figures 10, 12, 15, 17 dans les cas $\theta_o = -\theta_1 = \pi/4$ et $\theta_o = \pi/2$, $\theta_1 = -\pi/4$ que l'angle $\theta_s$ est le même et que les sillages sont identiques pour une même valeur du paramètre $\theta_s$. Le sillage ne dépend donc que de la forme aval du profil.

### Remarque 2.4

Outre les travaux théoriques de J.LERAY [1] qui mettent en évidence une famille différente de solutions avec sillage en utilisant des méthodes de fonctions analytiques, ce problème a été étudié par LIGHTHILL [12] et SOUTHWELL and VAISEY [14] dans le cadre de la cavitation. Des résultats comparables à ceux décrits ici sont obtenus avec des coefficients de cavitation K négatifs ($K = (q^2/q_\infty^2 - 1)$). Ils ont été fournis par des méthodes de transformation conforme ou par des méthodes numériques. La méthode étudiée ici présente l'avantage sur ces dernières de reposer sur une étude théorique rigoureuse, de s'adapter à tous les cas de profils symétriques sans incidence et, utilisant la
géométrie simple dans le plan de l'hodographe et indépendante du profil, de conduire à un calcul particulièrement court (de l'ordre de 40 secondes d'ordinateur CII IRIS 80 par profil).

## CONCLUSION

Au vu des résultats précédents, la méthode qui consiste à ramener la détermination du sillage à la résolution d'une inéquation variationnelle, fournit un moyen commode et efficace de résolution du problème de l'écoulement d'un fluide autour d'un profil symétrique.
<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /><figcaption>Figure N°14</figcaption></figure>

Sillages pour $\theta_o = -\theta_1 = \pi/6$

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

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

<figure><img src="image_23.png" /><figcaption>Figure N°15</figcaption></figure>

Sillages pour $\theta_o = -\theta_1 = \pi/4$

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

<figure><img src="image_26.png" /><figcaption>Figure N°16</figcaption></figure>

Sillages pour $\theta_o = -\theta_1 = \pi/2$

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

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

<figure><img src="image_29.png" /><figcaption>Figure N°17</figcaption></figure>

Sillages pour $\theta_o = \pi/2$ et $\theta_1 = -\pi/4$

<figure><img src="image_30.png" /></figure>
BIBLIOGRAPHIE

[1] C. BAIOCCHI : C.R. Acad. Sc. Paris, 273, 1971, p. 1215.

[2] H. BREZIS : Thèse, J. Math. Pures et Appl., 51, 1972, p. 1-168.

[3] H. BREZIS : Communication personnelle.

[4] H. BREZIS-G. DUVAUT : C.R. Acad. Sc. Paris, 276, 1973, p. 875.

[5] H. BREZIS-G. STAMPACCHIA : C.R. Acad. Sc. Paris, 276, 1973, p. 129.

[6] J.L. LIONS-G. STAMPACCHIA : Comm. Pure and Appl. Math., 20, 1967, p. 493-519.

[7] J. CEA-R. GLOWINSKI : Sur des méthodes d'optimisation par relaxation. RAIRO, R-3, Déc. 1973. pp. 5-32.

[8] P. CIARLET-P.A. RAVIART : General Lagrange and Hermite interpolation in $R^n$ with application to finite element method. Arch. Rat. Mech. Anal., 46, 1972, p. 177-199.

[9] R. GLOWINSKI : La méthode de relaxation, Rendiconti di Matematica, 14, Univ. de Rome, 1971.

[10] R. GLOWINSKI-J.L. LIONS-R. TREMOLIERES : Méthodes numériques de résolution des problèmes d'inéquation en Mécanique et en Physique. Dunod-Paris (à paraître).

[11] J.LERAY : Comm. Math. Helv., 8, 1935, pages 149 - 180 et 250 - 263.

[12] LIGHTHILL M.J. : 1949 Aero. Res. Coun. Rep. and Mem n°2328

[13] P.A. RAVIART : Cours de 3ème cycle-uer 48-Université PARIS VI.
[14] SOUTHWELL R.V. and VAISEY, 1946 Phil. Trans.Roy, Soc.A.240,
117.

[15] R.S. VARGA : Matrix Iterative Analysis. Prentice Hall. 1962.

[16] O.C. ZIENKIEWICZ : The finite element method. Mac Graw Hill
Publishing Co Ltd.. London, 1971.