RIA
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

ARCHIVES

Sur le calcul
des carenes optimales
en écoulement laminaire visqueux

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

Rapport de Recherche n° 127
Juin 1975
SUR LE CALCUL DES CARENES OPTIMALES
EN ÉCOULEMENT LAMINAIRE VISQUEUX

R. Glowinski, O. Pironneau

On se propose de calculer le profil d'aire donnée qui a la plus petite résistance de trainée dans un écoulement uniforme incompressible visqueux régi par les équations bidimensionnelles stationnaires de Navier-Stokes.
On utilise une méthode exposée par les auteurs dans un travail antérieur.

Le profil obtenu est long et mince, d'avant tangent à un dièdre d'angle 90° et d'arrière effilé ; sa résistance de trainée ne diffère que de très peu de celle d'un profil infiniment mince de même longueur. Compte-tenu de la complexité numérique du problème considéré, la précision obtenue n'est pas très élevé (de l'ordre de 5 %) mais ce travail illustre la méthode exposée antérieurement et montre qu'une meilleure précision peut être obtenue à condition d'y mettre le prix. Par ailleurs l'équation adjointe de l'équation de Navier-Stokes est résolue numériquement ; cette équation et sa solution étant amenées à jouer un rôle important dans les problèmes d'optimum désign" de la mécanique des fluides.

ABSTRACT :

An approximation of the profile of given area with smallest drag in laminar flow is obtained (Reynolds number between 10³ and 10⁵). It was shown previously by the Authors that the skin friction on such profiles had to satisfy certain optimality conditions ; the method used is based on those results. It was found that the optimum profile is long and thin (thickness to cord ratio about 10%) the front end shaped like a wedge of angle 90° and the rear end is a cusp. The drag is very close to the drag of the flat plate of equal length.
# PLAN

## 1. INTRODUCTION

## 2. RAPPELS DES RESULTATS THEORIQUES

2.1. Position du problème

2.2. Rappels des résultats théoriques. (I) Résultats d'existence.

2.3. Rappels des résultats théoriques. (II) Conditions nécessaires d'optimalité.

2.4. Rappels des résultats théoriques. (III) Propriétés du système adjoint.

2.5. Analyse de la singularité au bord de fuite.

## 3. DESCRIPTION DE L'ALGORITHMÉ DANS LE CAS CONTINU.

## 4. APPROXIMATION DU PROBLEME D'OPTIMISATION.

4.1. Remarques préliminaires.

4.2. Discrétisation du profil.

4.3. Discrétisation de l'écoulement potentiel.

4.4. Discrétisation des équations de Prandtl.

4.5. Résolution itérative.

4.6. Discrétisation du système adjoint.

4.7. Solution itérative.

4.8. Description de la méthode de descente pour le problème approché.

## 5. DESCRIPTION ET ANALYSE DES RESULTATS NUMÉRIQUES.

5.1. Description de l'exemple traité.

5.2. Valeur numérique des paramètres.

5.3. Description des résultats.

5.4. Analyse des résultats.

## 6. REMARQUES DIVERSES ET CONCLUSION.

6.1. Calcul des profils optimaux pour d'autres volumes et d'autres nombres de Reynolds.

6.2. Conclusion.
SUR LE CALCUL DES CARENES OPTIMALES
___
EN ECOULEMENT LAMINAIRE VISQUEUX

# 1. INTRODUCTION :

Le tracé optimal d'ailes, de carenes, de carrosseries, sont autant de
problèmes importants soumis quotidiennement aux ingénieurs de l'aérodynamique,
de l'hydrodynamique etc... Ces problèmes sont souvent assez difficiles à
modéliser mathématiquement car de nombreux paramètres "extérieurs" comme
le coût de l'étude, les difficultés de réalisation, la résistance des
matériaux entrent en ligne de compte. Cette situation est d'ailleurs heureuse
car elle permet à l'ingénieur de faire valoir son intuition.

Dans ce travail, nous approchons des problèmes d'optimum désign de la
mécanique des fluides d'une manière "technocratique": on cherche quelle est
la forme optimale de l'objet à dessiner pour un seul critère, à savoir sa
résistance de trainée.

L'objet de l'étude est d'abord de résoudre un problème abstrait que pose
la mécanique des fluides, mais également d'obtenir des renseignements qualitatifs
et de montrer quels sont les paramètres importants dans ces problèmes.
L'ingénieur y trouvera peut être un support pour développer son intuition à
propos de problèmes plus complexes.

Jusqu'à présent les renseignements que l'on possédait sur ce type de
problèmes sont assez succincts (on renvoie par exemple à LANDAU-LIFCHITZ
[1, n°46], pour des nombres de Reynolds pas trop élevés, i.e. < 10⁵) ; il faut
si possible que l'écoulement reste laminaire et que la couche limite ne décolle
pas. Pour ce type d'écoulement PIRONNEAU [2] a donné des conditions que doivent
satisfaire les profils optimaux. Ces conditions peuvent aussi - théoriquement -
permettre de calculer ces profils. Dans [2] on suppose les écoulements
incompressibles et gouvernés par les équations stationnaires de Navier-Stokes
(écoulement laminaire) ce qui correspond à des nombres de Reynolds moyens
et des profils aérodynamiques. Les conditions mentionnées ci-dessus avaient
été établies en utilisant les méthodes du contrôle optimal développées
par LIONS [3] pour les systèmes gouvernés par des équations aux dérivées
partielles, malheureusement ces conditions font intervenir l'état adjoint
du système de Navier-Stokes, c'est-à-dire une fonction qui, à notre connaissance
n'a pas d'interprétation physique.

L'objet de l'article est d'illustrer la méthode développée théoriquement
dans [2] et de donner un calcul de cette nouvelle quantité qu'est l'état
adjoint. Il avait été établi dans [2] que la variation de trainée d'un
profil en fonction d'une petite déformation de celui-ci est reliée à la
friction . La méthode itérative utilisée pour calculer les profils
symétriques bidimensionnels, d'aire donnée, qui ont la résistance de trainée
minimale est essentiellement la suivante :

On se donne l'enveloppe extérieure Ĉ_o de la couche-limite d'un profil C°;
on calcule ensuite la vitesse de glissement potentiel sur Ĉ_o, puis du/dn
et dw/du par intégration des équations de Prandtl et du système adjoint dans
la couche-limite, sillage compris. Les résultats permettent de savoir dans
quel sens modifier Ĉ_o pour obtenir un profil de résistance de trainée minimale
et ainsi de suite. La vitesse de glissement potentiel sur Ĉ_o est calculée
par la méthode des singularités, adaptée par LUU [4], les équations de Prandtl
sont discrétisées suivant un schéma aux différences finies du type de ceux
utilisées par KELLER [5] et les équations adjointes sont intégrées par un
schéma de discrétisation adjoint de celui utilisé pour les équations de Prandtl.
L'algorithme obtenu est comme on peut s'en douter assez complexe, chaque
itération prenant, dans l'exemple traité, 30 secondes d'IBM 370 168. Le
profil obtenu satisfait les conditions nécessaires d'optimalité à 10⁻² près,
il est long, mince, d'avant pointu à 90° (1) et d'arrière effilé en point de
rebroussement. Pour améliorer la précision il faudrait posséder des moyens de
calcul plus puissants ce qui permettrait d'utiliser des approximations
(1) Ce résultat était acquis dans [2] ; (on avait aussi démontré un résultat
identique dans PIRONNEAU [6] pour les écoulements de Stokes) il semble donc
que l'avant du profil dépende très peu du nombre de Reynolds : 90° pour les
écoulements bidimensionnels ; entre 110° et 120° pour les écoulements tridimen-
sionnels axisymétriques.
plus "élaborées" et traiter la singularité au bord de fuite d'une manière plus sophistiquée. Quoiqu'il en soit, cette étude montre que la méthode exposée dans [2] est fiable et qu'un laboratoire de mécanique des fluides, désireux de résoudre numériquement des problèmes de même type, peut s'en inspirer.

D'autre part, il a été procédé au calcul numérique de l'état adjoint et il est possible qu'à la lumière de ce calcul les mécaniciens des fluides lui trouveront une interprétation physique. Ce problème est important car la généralisation de cette approche aux écoulements à couche-limite *turbulente* ne pose pas trop de problèmes théoriques mais on ne dispose pas de méthode de calcul de l'état adjoint analogue à celle utilisée pour calculer l'écoulement moyen.
## 2. POSITION DU PROBLEME ET RAPPEL DES RESULTATS THEORIQUES.

### 2.1. Position du problème.

On considère le problème d' "optimum design"

$$ (2.1) \quad \text{Min} \left\{ \frac{2}{\nu} \int_{\Omega \cap C} e_{1j} e_{ij} \, dx \right\} $$

où $\nu$ est un ensemble de domaines d'aire donnée a et de frontière "régulière" et où

$$ (2.2) \quad e_{ij} = \frac{\nu}{2} \left( \frac{\partial u_j}{\partial x_i} + \frac{\partial u_i}{\partial x_j} \right), \quad 1 \le i,j \le 2, $$

$u = (u_1, u_2)$ est solution de l'équation de Navier-Stokes stationnaire

$$ (2.3) \quad \left\{ \begin{array}{ll} \nu \nabla^2 u - u \cdot \nabla u = \nabla p & \text{dans } \Omega \cap C \\ \nabla \cdot u = \mathbf{0} & \text{dans } \Omega \cap C \\ u = u_o & \text{sur } \partial\Omega, \\ u = \mathbf{0} & \text{sur } \partial C, \end{array} \right. $$

p dénotant la *pression*, $\nu$ la *viscosité cinématique*.

Le problème d'optimisation ci-dessus n'est bien posé, en toute rigueur, que si le système (2.3) n'admet qu'une seule solution, ce qui est le cas si $C$ est "aérodynamique" et $\nu$ suffisamment grand.

D'un point de vue physique le problème correspond à la recherche du profil $C$ de * traînée minimale* en écoulement laminaire visqueux incompressible (cf. PIRONNEAU [2])

En effet si $\Omega = \mathbb{R}^2$ la *traînée* $F(C)$ de $C$ est telle que

$$ u_o \cdot F(C) = \rho \int_{\mathbb{R}^2 \setminus C} e_{ij} e_{ij} \, dx $$

où $\rho$ est la *densité* du fluide.
## 2.2. Rappels des résultats théoriques (I) Résultat d'existence.

On peut démontrer en procédant comme dans GARABEDIAN [7, ch. 5] que si $\{\partial C | C \in \mathcal{S}\}$ est un ensemble de courbes équilipschitziennes sans point double, telles que les systèmes (2.3.) correspondants n'admettent qu'une solution, alors le problème (2.1) admet une solution.

Le problème de l'unicité reste ouvert à notre connaissance. Le problème de la dépendance de la solution par rapport à la constante de Lipschitz est aussi ouvert.

## 2.3. Rappels des résultats théoriques (II) Conditions nécessaires d'optimalité.

On démontre dans PIRONNEAU [2, n°3] le

**Théorème 2.1 :**

On suppose les contours $\partial C$ de $\mathcal{S}$, $c^1$ par morceaux ; si $\partial C \in \mathcal{S}$ est paramétré par $\sigma$ avec

$$ \partial C = \{ \xi(\sigma) \mid \sigma \in [0,1] , \xi(0) = \xi(1) \} , \\ \text{si } \partial C' \text{ est "} \alpha\text{-voisin" de } \partial C \text{ au sens} $$

$$ \partial C' = \{ \xi'(\sigma) \mid \xi'(\sigma) = \xi(\sigma) + \alpha(\sigma)\eta(\sigma), \sigma \in [0,1], \alpha(0)=\alpha(1) \} $$

où $n(\sigma)$ est la normale extérieure à $C$ au point $\xi(\sigma)$ alors la variation $\delta E$ de l'énergie dissipée dans la fluide, soit

$$ (2.4) \quad E(C) = \frac{2}{V} \int_{\Omega-C} e_{ij} e_{ij} \, dx $$

est donnée par

$$ (2.5) \quad \delta E = \int_{\partial C} \alpha(\sigma) \frac{\partial u}{\partial n} \cdot \left( \frac{\partial u}{\partial n} + 2 \frac{\partial w}{\partial n} \right) d(\partial C) + O(\alpha), $$

avec $O(\alpha)$ tel que $\lim_{\lambda \to 0} \lambda^{-1} O(\alpha \lambda) = 0$ et $W$ solution du système adjoint de Navier-Stokes stationnaire.
$$
(2.6) \quad \left\{
\begin{array}{ll}
\displaystyle \nabla^2 w + u \cdot \nabla w - \nabla u \cdot w = -u \cdot \nabla u + \nabla q & \text{dans } \Omega \cap C \\[1em]
\nabla \cdot w = 0 & \text{dans } \Omega^c \cap C \\[1em]
w |_{\partial C} = 0 & \\[1em]
w |_{\partial \Omega} = 0.
\end{array}
\right.
$$

De ce théorème on déduit le

Corollaire 2.1 : *Sous les hypothèses du théorème 2.1 toute solution C du problème (2.1) doit vérifier*

$$
\frac{\partial u}{\partial u} \cdot \left( \frac{\partial u}{\partial u} + 2 \frac{\partial w}{\partial u} \right) = C^{\text{te}} \quad \text{p.p. \underline{ sur } C.}
$$

On rappelle qu'il a été démontré dans PIRONNEAU [2] que cette condition ne peut
être satisfaite à l'avant du profil que si celui-ci est tangent à un dièdre
d'angle 90°, orienté dans l'axe de l'écoulement. En effet un changement
de variable du type Falkher-Skan v.PIRONNEAU [8] montre que les équations
de u et w sont "localees" à l'avant du profil et que seul le demi-dièdre
d'angle 45° satisfait à la condition d'optimalité.

2.4. Rappels des résultats théoriques. (III) Propriétés du système adjoint (26).

On sait depuis Prandtl que dans certaines conditions ( C aérodynamique), la
solution de (2.3) peut être raisonnablement approchée de la façon suivante :

Si $\Omega = \mathbb{R}^2$ on peut décomposer $\Omega - C$ en deux régions : $\Omega - \hat{C}$ et $\hat{C} - C$
(voir figure 1). Dans $\Omega - \hat{C}$ l'écoulement $u$ dérive du potentiel $\varphi$

$$
(2.7) \quad u = (u_1, u_2) = \left( \frac{\partial \varphi}{\partial x}, \frac{\partial \varphi}{\partial y} \right)
$$

où $\varphi$ est solution de

$$
(2.8) \quad \nabla^2 \varphi = 0 \text{ dans } \Omega - \hat{C}, \quad \varphi_\infty = u_0 x, \quad \frac{\partial \varphi}{\partial u}\Big|_{\hat{C}} = 0,
$$

$\hat{C} - C$ étant un domaine "très mince" (épaisseur de l'ordre de $R \frac{1}{2}$) de chaque coté de $C$ et dans le sillage de $C$. On utilise pour écrire des équations de l'écoulement dans $\hat{C} - C$ le système des coordonnées normales tangentielles et ; en première approximation $u = (u_s, u_n)$ est dans $\hat{C} - C$ solution de
$$
(2.9) \quad \left\{
\begin{array}{l}
- \nu \frac{\partial^2 u_s}{\partial n} + u_s \frac{\partial u_s}{\partial s} + u_n \frac{\partial u_s}{\partial n} = U(s) \frac{dU(s)}{ds} \\[1em]
\qquad \frac{\partial u_s}{\partial s} + \frac{\partial u_n}{\partial n} = 0
\end{array}
\right.
$$

avec $U(s)=\frac{\partial\varphi}{\partial s}$ et $\hat{\partial C}$ et $u_s=u_n=0$ sur $\partial C$ (n=0) et $u(s)=U(s)$ sur $\hat{\partial}C$.
Dans le sillage la condition $u_s=0$ est remplacée par $\frac{\partial u_s}{\partial n}(s,0)=0$.

Dans PIRONNEAU [2, n°4] on démontre la

**Proposition 2.1** : *Dans* $\hat{C} - C$, $w = (w_s, w_n)$ *est en première approximation*
*solution de*

$$
(2.10) \quad \left\{
\begin{array}{@{}l@{}}
  \displaystyle \nu \frac{\partial^2 w_s}{\partial n^2} + u_s \frac{\partial w_s}{\partial s} + u_n \frac{\partial w_s}{\partial n} - \frac{\partial u_s}{\partial s} w_s = \frac{\partial q}{\partial s} - \nu \frac{\partial^2 u_s}{\partial n^2} \\[1em]
  \displaystyle - \frac{\partial u_s}{\partial n} w_s = \frac{\partial q}{\partial n} \\[1em]
  \displaystyle \frac{\partial w_s}{\partial s} + \frac{\partial w_n}{\partial n} = 0
\end{array}
\right.
$$

avec les conditions aux limites

$w_n = 0$ **sur** $\partial C$, $w_s = 0$ **sur** $\partial C$ (ou $\frac{\partial w_s}{\partial n} = 0$ **dans le sillage**)

$w_s = 0, q = 0$ **sur** $\partial C$, $w = 0$ **dans** $\Omega - \hat{C}$.

2.5. Analyse de la singularité au bord de fuite.

Si le bord de fuite du profil est arrondi, ouônique d'angle non
nul, la vitesse de glissement U tend vers zéro lorsque s tend vers s₁,
s₁ : abscisse curviligne du bord de fuite. Comme il est montré dans LANDAU-
LIFCHITZ [1], ceci provoque un décollement dans la couche limite à partir
d'un point d'abscisse curviligne s<sub>d</sub> < s<sub>1</sub> et en première approximation par
rapport à (s<sub>d</sub> - s) ¹/² on a

$$
(2.11) \quad u_n = \frac{A}{2} \frac{U_s sd, n}{\sqrt{sd-s}}
$$

$$
(2.12) \quad u_s = u_s sd, n) + A \frac{\partial u_s}{\partial n} (sd, n) \sqrt{sd-s}.
$$
Si on élimine $q$ dans (2.10) on obtient

$$ (2.13) \quad \nu \frac{\partial^3 w}{\partial n^3} + \frac{\partial^2 w}{\partial n^2} \frac{\partial^2 w}{\partial n^2} u_n + \frac{\partial^2 w}{\partial n \partial s} \frac{\partial u}{\partial n} \frac{\partial u}{\partial n} \frac{\partial w}{\partial n} + 2 \frac{\partial u}{\partial n} \frac{\partial w}{\partial n} \frac{\partial s}{\partial n} = - \nu \frac{\partial^3 u}{\partial n^3} $$

Supposons qu'au voisinage de $s = sd$ on ait

$$ (2.11) \quad ws = \sum_{i=1}^{+\infty} a_i(n) (s_d - s)^{i/2} + a_o(n) + \sum_{j=1}^{m} a_{-j}(n) (s_d - s)^{-j/2} ; \text{ de} $$

(2.11) - (2.14) on déduit par identification de coefficients que

$$ a_{-j} = 0, \quad j = 1, \dots, m $$

et

$$ (2.15) \quad \frac{A}{2} a''_o u_s + \frac{A}{2} u'_s a'_o - \frac{1}{2} a'_1 u_s - a_1 u'_s = 0 $$

Au voisinage de $n=0$, $u_s = u''_s (s_d, 0)^{n/2}$ ; il résulte donc de (2.15)

que

$$ A a''_o n + A a'_o - a'_1 n - 2a_1 = 0 $$

ce qui entraîne, compte tenu de la condition $w(s,o) = 0$

$$ w_s(s,n) = \alpha n \sqrt{s_d - s} + O(n) + O(\sqrt{s_d - s}) $$

et donc

$$ \lim_{s \to sd} \frac{\partial u_s}{\partial n} \left( \frac{\partial u_s}{\partial n} + 2 \frac{\partial w_s}{\partial n} \right) = 0. $$

Pour satisfaire la condition d'optimalité du corollaire 2.1 il faut donc que le bord de fuite soit un *point de rebroussement* du profil, seule condition pour éviter le décollement de la couche-limite ; il n'en reste pas moins que le changement de conditions aux limites au bord de fuite provoque une discontinuité de $u$ et $w$ qui rend fort difficile le calcul de $w$.

### 3. DESCRIPTION DE L'ALGORITHME D'OPTIMISATION DANS LE CAS CONTINU

Bien que l'ensemble des contrôles $S$ n'ait pas de structure d'espace vectoriel, on voit d'après (2.5) que la quantité

$$ \left. \frac{\partial u}{\partial n} \left( \frac{\partial u}{\partial n} + 2 \frac{\partial w}{\partial n} \right) \right|_{\partial C} $$

va jouer le rôle du gradient de $E$. Pour cette raison on posera

$$ (3.1) \quad \text{Grad } E = \left. \frac{\partial u}{\partial n} \left( \frac{\partial u}{\partial n} + 2 \frac{\partial u}{\partial n} \right) \right|_{\partial C} $$
soit

$$
(3.2) \quad K_l = \int_{\partial C} \mathrm{Grad} \, E(\sigma) \, d(\partial C),
$$

on remarquera, d'après (2.5) que

$$
\int_{\partial O} \alpha(\sigma) \, d(\partial C) = \therefore \delta E = \int_{\partial C} \alpha(\sigma) (\mathrm{Grad} \ E(\sigma) - K_l) \, d(\partial C) + O(\alpha).
$$

En conséquence si C n'est pas solution de (2.1) en portant sur la
normale la quantité

$$
(3.3) \quad \alpha(\sigma) = - \lambda (\mathrm{Grad} \ E(\sigma) - K_1), \quad \lambda \text{ "petit"},
$$

on obtient un profil amélioré de même aire puisque

$$
(3.4) \quad \delta \mathfrak{P} = - \lambda \int_{\partial C} (\mathrm{Grad} \, E(\sigma) - K_1)^2 d(\partial C) + O(\alpha) < 0.
$$

Ceci conduit donc à l'algorithme ci-dessous

$$
(3.5) \text{ choisir } C_o \text{ d'aire a, poser } i = o, \text{ choisir } \lambda \text{ "petit"},
$$

$$
(3.6) \text{ calculer } \frac{\partial u_i}{\partial n}, u_i \text{ solution de (2.3) avec } C = C_i,
$$

$$
(3.7) \text{ Calculer } \frac{\partial w_i}{\partial n}, w_i \text{ solution de (2.6) avec } C = C_i,
$$

$$
(3.8) \text{ Construire } C_{i+1} \text{ en portant sur la normale extérieure à } C_i \text{ la quantité}
$$

$$
- \lambda \left[ \frac{\partial u_i}{\partial n} \left( \frac{\partial u_i}{\partial n} + 2 \frac{\partial w_i}{\partial n} \right) - \int_{\partial C_i} \frac{\partial u_i}{\partial n} \left( \frac{\partial u_i}{\partial n} + 2 \frac{\partial w_i}{\partial n} \right) d (\partial C_i) \int_{x_i} d(\partial C_i) \right]
$$

$$
(3.9) \text{ Poser } i = i+1 \text{ et retourner en } (3.6).
$$

L'algorithme (3.5) - (3.9) n'est autre qu'une méthode de gradient
à pas constant. En fait pour assurer la convergence de l'algorithme, il faudrait
utiliser un pas variable (une méthode à deux lignes par exemple, cf.POLAK [9]),
mais les calculs correspondants sont assez coûteux si bien que pratiquement il
vaut mieux opérer à pas fixe quitte à changer $\lambda$ en $\frac{\lambda}{2}$ si $E(C_i+1) > E(C_i)$.

Par ailleurs à grand nombre de Reynolds (100 <R <100 000) on peut utiliser
les approximations de Prandtl pour calculer U ; il faut alors remplacer (2.3)

par la résolution de (2.8); (2.9). De même, dans (3.7), $\frac{\partial w_i}{\partial n}$ est calculé à l'aide (2.10).
## 4. APPROXIMATION DU PROBLEME (2.1)

### 4.1. Remarques préliminaires.

**Remarque 4.1** : Lorsque *le profil est aérodynamique*, plus précisement lorsque la couche-limite ne "décolle pas" ($\frac{\partial u_s}{\partial n} > 0$), la région $\hat{C} - C$ de $\Omega \setminus C$ se réduit à une mince couche d'épaisseur de l'ordre de $R^{-1/2}$ autour de C et dans le sillage de C. Comme $\hat{C}$ est difficile à calculer à partir de C on itère sur $\hat{C}$ et non sur C. On déduira C de $\hat{C}$ en portant sur la normale à $\partial\hat{C}$ l'épaisseur de la couche limite (qui est facile à calculer à partir de $u_s, u_n$).

**Remarque 4.2** : La contrainte "C aérodynamique" est difficile à traduire quantitativement, elle impose de grandes restrictions sur $C_o$.

**Remarque 4.3** : En fait il est commode d'utiliser, plutôt que $(w,q)$ le vecteur $(Q,r)$ défini par

$$Q = u + 2w, \quad r = 2q - \frac{1}{2} |u|^2 + p;$$

$$ (4.1) \qquad \left\{ \begin{array}{l} Q \text{ satisfait aussi à} \\[1em] \nu \quad \frac{\partial^2 Q_s}{\partial n^2} + u_s \frac{\partial Q_s}{\partial s} + u_n \frac{\partial Q_s}{\partial n} - \frac{\partial u_s}{\partial s} Q_s = \frac{\partial r}{\partial s} \\[1em] \qquad - \frac{\partial u_s}{\partial n} Q_s = \frac{\partial r}{\partial n} \\[1em] \qquad \frac{\partial Q_s}{\partial s} + \frac{\partial Q_n}{\partial n} = 0 \quad \text{dans } \hat{C} - C, \\[1em] \qquad Q_s = Q_n = 0 \quad \text{sur } \partial C, \quad Q_s = \frac{dU}{ds}, \ r = 0 \quad \text{sur } \partial \hat{C} \end{array} \right. $$

**Démonstration :**

De (2.3) et (2.6) on déduit :

$$\nu \nabla^2 Q + u \cdot \nabla (Q-u) - \nabla u \cdot (Q-u) + u \cdot \nabla u = \nabla (p + 2q),$$
soit

$$
(4.2) \quad \nu \nabla^2 Q + u \cdot \nabla Q - \nabla u \cdot Q = \nabla (p + 2 q) - \nabla u \cdot \nabla u = \nabla (p + 2 q - \frac{1}{2} |u|^2).
$$

Les équations (4.1) sont alors déduites de (4.2) par les mêmes
méthodes de perturbation singulière qui ont permis de déduire (2.10) de
(2.6).

4.2 Discrétisation du profil.

On approche l'ensemble $\mathcal{g}$ des profils admissibles (cf. n°2.3.) par un
sous ensemble $\mathcal{g}_h$, h étant un paramètre qui sera défini ci-dessous ; de façon
plus précise on introduit $N+1$ point distincts de $[0,1]$ soit
$\sigma_1 = 0, \sigma_2, \dots, \sigma_{N+1} = 1$, alors

$$
\begin{gather*}
\partial_c \hat{h} = \left\{ \xi_h(\sigma) \mid \sigma \in [0,1], \xi_h \in C^0[0,1] \times C^0[0,1], \xi_h \text{ affine sur } [\sigma_j, \sigma_{j+1}] \right\}_{j=1,\dots,N}, \\
\text{on notera } h = \max_{\substack{\alpha_j \\ j=1,\dots,N}} |\alpha_{j+1} - \sigma_j|.
\end{gather*}
$$

Il est commode d'introduire ensuite les points $\xi_j = \xi_h (\alpha_j)$, on a donc
$\xi_1 = \xi_{N+1}$ et l'arc $\xi_j, \xi_{j+1}$ de $\partial \hat{C}_h$ est un *segment de droite*. On approche
ensuite la normale au point $\xi_j$ par la perpendiculaire à $\overrightarrow{\xi_{j-1} \xi_{j+1}}$ passant
par $\xi_j$, cf. fig.1 (bien entendu d'autres choix sont possibles), et on pose

$$
s_j = \sum_{i=1}^{j+1} |\xi_i \xi_{i+1}|, \quad s_1 = 0.
$$

4.3. Discrétisation de l'écoulement potentiel.

Le calcul de l'écoulement potentiel dans $\Omega - \hat{C}$ ne constitue qu'une
étape intermédiaire pour le calcul de $\frac{\partial u}{\partial n}|_{\partial C}$. En fait il faut remarquer
que seul $\varphi$ (ou $\varphi$ est solution de 2.8) est nécessaire pour la suite des
calculs. Pour cela on a utilisé l'adaptation de TS. LUU $[4]$ de la méthode
des singularités discrètes. On rappelle que cette méthode est fondée sur
l'application de la formule de Green avec la fonction de Green du problème ;
soit dans notre cas :
Toute fonction harmonique à l'extérieur de $\hat{C}$ et nulle à l'infini (il faut donc en fait prendre $\varphi - u_x$ comme fonction inconnue) satisfait à

$$2 \varphi(z) = \int_{\partial \hat{C}} \left[ \frac{\partial \varphi}{\partial n} \log|z-\xi(\sigma)| - \varphi(\xi(\sigma)) \frac{\partial}{\partial n_{\sigma}} \log|z-\xi(\sigma)| \right] d(\partial \hat{C}).$$

Si $\partial \hat{C}$ est donné par (4.3), en tout point $\xi_j^C$ milieu de $\xi_j, \xi_{j+1}$,

$\frac{\partial \varphi}{\partial s}$ et $\frac{\partial \varphi}{\partial n}$ satisfait à

$$(4.3) \quad \frac{\partial \varphi}{\partial s} (\xi_j^C) - i \frac{\partial \varphi}{\partial u} (\xi_j^C) = \sum_{\rho=1}^{N} \frac{q_{\rho} + i\gamma_{\rho}}{2\Pi} e^{i(\beta_j - \beta_{\rho})} \log\left(\frac{\xi_j^C - \xi_{\rho}}{\xi_j^C - \xi_{\rho+1}}\right)$$

où $\beta_j$ est l'angle $(\vec{\partial}x, \xi_j, \xi_{j+1})$ et $q_{\rho}, \gamma_{\rho}$ sont déterminés en résolvant le système linéaire constitué par la partie imaginaire de (4.3) où

$\frac{\partial \varphi}{\partial n}(\xi_j^C)$ est remplacé par sa valeur (donnée par (2.8)). Pour les profils **symétriques** on prend $\gamma_\ell=0$, (1) ce qui donne

$$(4.4) \quad \sum_{\ell=1}^{N} \frac{q_{\ell}}{2\Pi} \left( e^{i(\beta_j - \beta_{\ell})} \log \frac{\xi_j^C - \xi_{\ell}}{\xi_j^C - \xi_{\ell+1}} \right) = u_0 \sin \beta_j, \quad j = 1, \dots, N.$$

On obtient alors

$$(4.5) \quad \frac{\partial \varphi}{\partial s} (\xi_j^C) = \sum_{\ell=1}^{N} \frac{q_{\ell}}{2\Pi} \text{Re} \left[ e^{i(\beta_j - \beta_{\ell})} \log \left( \frac{\xi_j^C - \xi_{\ell}}{\xi_j^C - \xi_{\ell+1}} \right) \right] + u_0 \cos \beta_j$$

### 4.4. **Discrétisation des équations de Prandtl.**

Dans (2.9) on utilise un système de variables réduites (cf. par exemple JONES et WATSON dans ROSENHEAD [10] ce qui permet de réécrire (2.9) sous la

$$(4.6) \quad \left\{ \begin{aligned} & -\frac{\partial^2 u_s}{\partial n^2} + u_s \frac{\partial u_s}{\partial s} + u_n \frac{\partial u_s}{\partial n} = U(s) \frac{dU(s)}{s} \quad \text{dans } ]0, +\infty[ \times ]0, +\infty[, \\ & \frac{\partial u_s}{\partial s} + \frac{\partial u_n}{\partial n} = 0 \quad \text{dans } ]0, +\infty[ \times ]0, +\infty[ , \end{aligned} \right.$$

<footnote>(1) Dans le cas général $\partial_{\ell}$ est choisi de manière à satisfaire la condition de **KUTTA-JOUKOWSLSY** au bord de fuite.</footnote>
avec les conditions aux limites

$$u(s) = U(s) \quad \sin = +\infty$$

$$u_s = u_n \quad \sin = 0, s < L$$

$$u_s = 0, \quad \frac{u_s}{\sin} = 0, s > L$$

où L est la longueur du bord supérieur de $\hat{\partial}C$, sillage non compris (on obtiendra la friction véritable en multipliant par $\frac{1}{W}$ la valeur obtenue à partir de (4.6), (4.7)).

Pour les profils non symétriques il faut résoudre (4.6) de chaque côté de C.

On peut montrer en utilisant la transformation de Falkner-Skan (cf., par exemple, ROSENHEAD [10] que le système (4.6) a un caractère hyperbolique en S et qu'il s'auto-démarre sur la ligne $s=0$ ; cependant pour des raisons de stabilité numérique il est préférable de partir de $s = s_o > 0$ et d'utiliser une condition aux limites supplémentaires, soit :

$$u(n, s_o) = u_d(n), \quad n \in ]0, +\infty[$$

où $u_d$ est calculé à l'aide des équations de Falkner-Skan ($u_d$ est tabulé dans ROSENHEAD [10] pour un bord d'attaque cônique d'angle 90°).

On approche ensuite $]0, +\infty[$ x]0, + $\infty$[ par ]0, +$\infty$[ x]0, H[, H "suffisament "grand" et on discrétise (4.6), (4.7) *par différences finies*. Soit M un entier positif (destiné à tendre vers $+\infty$) et $k=\frac{H}{M}$, k sera le pas de discrétisation dans la direction normale. Au point{ s_j, ik} les équations (4.6) sont approchées par

$$ (4.8) \left\{ \begin{array}{l} - \frac{u_{j+1}^2 - 2u_j u_{j-1} + u_{i-1}^2}{k^2} + \frac{1}{2} \frac{(u_i^{j+1})^2 - (u_i^j)^2}{s_{j+1} - s_j} + v_i^{j+1} \frac{u_{i+1}^j - u_{i-1}^j}{2k} = \frac{1}{2} \frac{(U^{j+1})^2 - (U^j)^2}{s_{j+1} - s_j}, \\[1em] v_i^{j+1} = v_{i-1}^{j+1} - k \frac{u_i^{j+1} - u_i^j}{s_{j+1} - s_j}; \quad i = 2, \dots, M-1, \quad j \ge j_o, \end{array} \right. $$
avec les conditions aux limites

$$
(4.9) \quad \left\{
\begin{array}{l}
v_1^j = 0, \quad u_M^j = U^j ; \quad u_1^j = 0 \text{ si } s_j \le L, \quad u_1^j = u_2^j \text{ si } s_j > L, \\
u_i^{jo} = u_i^i, \quad i = 1, \dots, M.
\end{array}
\right.
$$

4.5. **Solution itérative du système** (4.8), (4.9).

Etant donné la forme du système (4.8), (4.9) et en particulier le
fait que le vecteur $\{u_i^1\}_{1<i<M}$ est connu de par les conditions initiales, il
est naturel, pour résoudre (4.8), (4.9), d'utiliser une méthode *pas à pas*
dans laquelle $\{u_i^j\}_{j \ge 1}$, étant connu on détermine $\{u_i^{j+h}\}_{1 \le i \le M}$ par
un procédé *itératif* compte-tenu des conditions aux limites (4.9) et du caractère
non linéaire des équations (4.8).

Avec des notations évidentes, et pour {$u_i^j$} $1 \le i < M$ donné le système
(4.8), (4.9) est de la forme :

$$
(4.10) \quad \left\{
\begin{array}{@{}l@{}}
g_1(x,y) = 0 \\
g_2(x,y) = 0
\end{array}
\right.
$$

avec $x,y \in \mathbb{R}^M$, $x = (x_1, \dots, x_M)$, $y = (y_1, \dots, y_M)$, le vecteur $\{u_i^{j+1}\}$ jouant le rôle
de $x$ et $\{v_i^{j+1}\}$ celui de $y$.

Pour résoudre (4.10) on peut envisager une méthode du type
surrelaxation ; dans ce travail on a utilisé l'algorithme suivant

$$
(4.11) \quad (x^0, y^0) \text{ donné} \\
\quad (x^m, y^m) \text{ connu on détermine } (x^{m+1}, y^{m+1}) \text{ par}
$$

$$
(4.12) \quad g_1^i (x_1^{m+1}, \dots, \bar{x}_i^{m+1}, \dots ; y^m) = 0, \ 1 \le i \le M
$$

$$
(4.13) \quad x_i^{m+1} = x_i^m + \omega_1 (\bar{x}_i^{m+1} - x_i^m), \omega_1 > 0,
$$

$$
(4.14) \quad g_2(x^{m+1}, y^{m+1} = 0.
$$

Dans (4.12) on calcule $\bar{x}_i^{m+1}$ en résolvant tout simplement une équation
du second degré dont on ne retient que la plus grande racine ; par ailleurs le
calcul de $y^{m+1}$ équivaut à la résolution d'un système linéaire, triangulaire
didiagonal, donc trivial.
En ce qui concerne l'initialisation (4.11), pour résoudre (4.8), (4.9)
on utilise {$u_i^j$} pour le calcul de {$u_i^{j+1}$}, {$v_i^{j+1}$}.

Signalons que la *convergence* d'un tel algorithme semble difficile à établir d'un point de vue théorique (on trouvera dans CEA-GLOWINSKI [11] des résultats de convergence pour des algorithmes du type ci-dessus lorsque les opérateurs $g_1, g_2$ vérifient des propriétés convenables de *monotonie*.

### 4.6. *Discrétisation du système adjoint (4.1).*

On utilise le même système de coordonnées réduites que pour $u$, afin d'éliminer $v$ dans (4.1) et transformer $\hat{C} - C$ en un quart de plan.

Le système (4.1) devient (il est inutile par la suite de calculer $Q_n$):

$$ (4.15) \quad \left\{ \begin{array}{@{}l@{}} \displaystyle \frac{\partial^2 Q_s}{\partial n^2} + u_s \frac{\partial Q_s}{\partial s} + u_u \frac{\partial Q_s}{\partial n} + \frac{\partial u_s}{\partial n} Q_s = \frac{\partial r}{\partial s} \\[1em] \displaystyle - \frac{\partial w}{\partial n} Q_s = \frac{\partial r}{\partial n} \quad \text{dans } ]0, +\infty[ \times ]0, +\infty[ , \end{array} \right. $$

avec les conditions aux limites

$$ (4.16) \quad \left\{ \begin{array}{@{}l@{}} Q_s(n=0) = 0 \quad \text{si } s < L, \quad \frac{\partial Q_s}{\partial n}(n=0) = 0 \quad \text{si } s > L, \\[1em] Q_s(n=+\infty) = \frac{\partial \varphi}{\partial s}(s), \quad r(n=+\infty) = 0, \\[1em] Q_s(s=+\infty) = u_s(s=+\infty). \end{array} \right. $$

Une discrétisation de (4.15) *indépendante* de celle de (4.6) risque d'introduire des *instabilités* numériques graves. Pour discrétiser (4.15) on remarque que (4.15) est l'adjoint (au sens des opérateurs différentiels) de (4.6) ; ceci n'est pas évident à priori puisqu'à l'origine ce sont les équations "mères" (2.3) et (2.6) de (4.6), (4.15) qui sont adjointes l'une de l'autre. Pour obtenir l'adjoint du discrétisé (4.8) de (4.6), on introduit $\mu_i^j$ (soit $\frac{\partial u_s}{\partial n}$ en continu) d'où pour (4.8) la forme équivalente
$$
(4.17) \quad \left\{
\begin{aligned}
\frac{1}{k}(\mu_{i+1}^{j+1} - \mu_i^{j+1}) &= v_i^{j+1} \frac{\mu_{i+1}^{j+1} + \mu_i^{j+1}}{2} + \frac{1}{2} \frac{(u_i^{j+1})^2 - (u_i^j)^2 - (v_i^{j+1})^2 + (v_i^j)^2}{s_{j+1} - s_j}, \\
\frac{1}{k}(u_i^{j+1} - u_{i-1}^{j+1}) &= \mu_i^{j+1}, \\
\frac{1}{k}(v_i^{j+1} - v_{i-1}^{j+1}) &= \frac{u_i^{j+1} - u_i^j}{s_{j+1} - s_j}, \quad i=2, \dots, M-1, \ j > j_o\ ;
\end{aligned}
\right.
$$

1' *Hamiltonien* de (4.17) s'écrit alors

$$
(4.18) \quad \left\{
\begin{aligned}
\mathcal{H} &= \sum_{i,j} Q_i^j \left[ \frac{1}{k} (\mu_{i+1}^{j+1} - \mu_i^{j+1}) - v_i^{j+1} (\mu_{i+1}^{j+1} + \mu_i^{j+1}) \right] - \frac{(u_i^{j+1})^2 - (u_i^j)^2 + (v_i^{j+1})^2 - (v_i^j)^2}{2(s_{j+1} - s_j)} \\
&\quad + g_i^j \left[ \frac{1}{k} (u_i^{j+1} - u_{i-1}^{j+1}) - \mu_i^{j+1} \right] + r_i^j \left[ \frac{1}{k} (v_i^{j+1} - v_{i-1}^{j+1}) + \frac{(u_i^{j+1} - u_i^j)}{s_{j+1} - s_j} \right]
\end{aligned}
\right.
$$

En dérivant partiellement $\mathcal{H}$ par rapport à $\mu_i^{j+1}$ on obtient

$$
(4.19) \quad \frac{\partial \mathcal{H}}{\partial \mu_i^{j+1}} = - \frac{1}{k} Q_i^j + \frac{1}{k} Q_{i-1}^j - \frac{1}{2} Q_i^j v_i^{j+1} - \frac{1}{2} Q_{i-1}^j v_{i-1}^{j+1} - g_i^j
$$

de même

$$
(4.20) \quad \frac{\partial \mathcal{H}}{\partial v_i^{j+1}} = - \frac{1}{2} Q_i^j (\mu_{i+1}^{j+1} + \mu_i^{j+1}) + \frac{1}{k} (r_i^j - r_{i+1}^j)
$$

$$
(4.21) \quad \frac{\partial \mathcal{H}}{\partial u_i^{j+1}} = \frac{Q_i^j u_i^{j+1}}{s_{j+1} - s_j} + \frac{Q_i^{j+1} u_i^{j+1}}{s_{j+2} - s_j} + \frac{g_i^j - g_{i+1}^j}{k} + \frac{r_i^j}{s_{j+1} - s_j} - \frac{r_i^{j+1}}{s_{j+2} - s_{j+1}}
$$

On obtient alors les équations adjointes de (4.17) discrétisées de (4.15)
à partir de (4.19) - (4.21) par

$$
\frac{\partial \mathcal{H}}{\partial \mu_i^{j+1}} = 0, \quad \frac{\partial \mathcal{H}}{\partial v_i^{j+1}} = 0, \quad \frac{\partial \mathcal{H}}{\partial u_i^{j+1}} = 0,
$$

soit après élimination des $g_i^j$,
$$
(4.22) \quad \left\{
\begin{array}{l}
\displaystyle \frac{Q_i^j - 2Q_i^j + Q_{i+1}^j}{k^2} + \left( \frac{Q_i^{j+1}}{s_{j+2} - s_{j+1}} - \frac{Q_i^j}{s_{j+1} - s_j} \right) u_i^{j+1} + \frac{Q_{i+1}^j v_{i+1}^{j+1} - Q_{i-1}^j v_{i-1}^{j+1}}{2k} = 0 \\[1em]
\displaystyle \frac{r_i^{j+1}}{s_{j+2} - s_{j+1}} - \frac{r_i^j}{s_{j+1} - s_j},
\end{array}
\right.
$$

$$
(4.23) \quad - Q_i^j \frac{u_{i+1}^{j+1} - u_{i-1}^{j+1}}{2k} = \frac{r_{i+1}^j - r_i^j}{k},
$$

avec des conditions aux limites

$$
Q_M^j = U_M^j, r_M^j = 0; Q_I^j = 0 \text{ si } s_j \le L, Q_I^j = Q_2^j \text{ si } s_j > L,
$$

$$
Q_i^J = u_i^J \text{ (J entier positif destiné à tendre vers l'infini).}
$$

4.7. **Résolution itérative de** (4.22, (4.23).

Le système (4.22), (4.23) est voisin du sytème (4.8) ; c'est un
système parabolique linéaire discret à intégrer numériquement, *pas à pas*,
dans le sens des j décroissants, en partant de j = J.

Pour résoudre (4.22), (4.23) on a utilisé une méthode de *surrelaxation* analogue à celle du n°4.5 (avec un paramètre de surrelaxation $\omega_2$ différent de $\omega_1$) :

Connaissant {Qᵢⱼ₊₁} , {rᵢⱼ₊₁} on calcule {Qᵢⱼ} à l'aide de (4.22), puis
{vᵢⱼ} à l'aide de (4.23) ; de manière plus précise l'algorithme utilisé est
le suivant :

$$
(4.24) \quad \text{on initialise} \quad Q_i^{j,0} \text{ et } r_i^{j,0} \text{ à } Q_i^{j+1} \text{ et } r_i^{j+1}
$$

$$
(4.25) \quad \left\{
\begin{array}{l}
\text{on calcule } Q_i^{j,m+1} \text{ par } Q_i^{j,m+1} = Q_i^{j,m} + \omega_2 (\bar{Q}_i^{j,m} - Q_i^{j,m}) \\
\text{où } \bar{Q}_i^{j,m} \text{ est la solution de (4.22),}
\end{array}
\right.
$$

$$
(4.26) \quad \text{on calcule } r_{i+1}^{j,m} \text{ par (4.23),}
$$
puis on retourne en (4.25) avec $m=m+1$.

Le test d'arrêt utilisé est

$$ \text{Max}_{1 \le i \le M} \left| Q_{i}^{j,m+1} - Q_{i}^{j,m} \right| < \varepsilon, \varepsilon > 0 \quad \text{donné.} $$

### 4.8. Description de la méthode de descente pour le problème approché.

Compte-tenu des méthodes de discrétisations analysées ci-dessus, l'algorithme (3.5) - (3.9) est approché par

$$ (4.27) \quad \left\{ \begin{array}{l} \text{on choisit } N'' = 2N'+1, \text{ nombre de points de discrétisation des profils, puis } M,k \text{ et } J > N' \text{ pour définir le réseau de discrétisation des équations de la couche-limite ; on choisit enfin } j_o, \text{ indice de la section initiale pour les équations de la couche-limite.} \\ \text{on choisit les paramètres de surrelaxation } \omega_1, \omega_2. \end{array} \right. $$

$$ (4.28) \quad \text{On choisit les paramètres de surrelaxation } \omega_1, \omega_2. $$

$$ (4.29) \quad \left\{ \begin{array}{l} \text{On choisit } \varepsilon' \text{ test d'arrêt des méthodes de surrelaxation ci-dessus et } \varepsilon \text{ test d'arrêt de la méthode de descente ci-après.} \\ \text{On choisit } \varepsilon' \text{ test d'arrêt des méthodes de surrelaxation ci-dessus et } \varepsilon \text{ test d'arrêt de la méthode de descente ci-après.} \end{array} \right. $$

$$ (4.30) \quad \text{On choisit } \lambda, \text{ qui est le pas dans la direction de descente.} $$

$$ (4.31) \quad \text{On choisit } U_o, \text{ vitesse de l'écoulement à l'infini.} $$

$$ (4.32) \quad \left\{ \begin{array}{l} \text{On choisit } \hat{C}_o \text{ "long et mince", symétrique, de volume donné d'avant-cône d'angle } 90^\circ \text{ et d'arrière effilé, le bord de fuite étant un point de \underline{rebroussement} du profil. On pose } i_g=0. \end{array} \right. $$

$$ (4.33) \quad \left\{ \begin{array}{l} \text{On calcule } U^j = \frac{\partial \psi}{\partial s} (\xi_j^c) \text{ par (4.5), } q_l \text{ étant déterminé par la résolution du système linéaire (4.4). On détermine alors } \{s_j\} \text{ et } L = s_{N+1}. \end{array} \right. $$

$$ (4.34) \quad \left\{ \begin{array}{l} \text{On détermine } \{u_i^j\}, \{v_i^j\}, i=1, \dots, M, j=1, \dots, J \text{ en résolvant } \\ \text{(4.8), (4.9) par surrelaxation de paramètre } \omega_1. \end{array} \right. $$
$$
(4.35) \quad \left\{
\begin{array}{l}
\text{On détermine } Q_i^j, i=1, \dots, M, j=j_o, \dots, J \text{ en résolvant} \\
\text{(4.22), (4.23) par surrelaxation de paramètre } \omega_2.
\end{array}
\right.
$$

$$
(4.36) \quad \left\{
\begin{array}{@{}l@{\qquad}l@{}}
\text{On calcule, Grad } E_j & \displaystyle = \frac{u_2^j - u_1^j}{k} \times \frac{Q_2^j - Q_1^j}{k} \text{, puis} \\[1em]
K = \frac{1}{L} \sum_{j=1}^{N} (s_{j+1} - s_j) \text{ Grad } E_j \text{, } K_2 = \frac{N'}{L} \sum_{j=1}^{N'} (s_{j+1} - s_j) (\text{Grad } E_j - K_1)^2 \text{ et} \\[1em]
F = \frac{N'}{L} \sum_{j=1}^{N'} (s_{j+1} - s_j) \left( \frac{u_2^j - u_1^j}{k} \right) \text{ si l'on veut connaître la traînée.}
\end{array}
\right.
$$

$$
(4.37) \quad \text{Si } K_2 \le \epsilon' : \text{ stop.}
$$

$$
(4.38) \quad \left\{
\begin{array}{@{}l@{\qquad}l@{}}
\text{Si } K_2 > \epsilon' \text{ on porte, pour obtenir } C_{ig+1} \text{, la quantité } \lambda \text{ (Grad } E_j \text{-K)} & \\
\text{sur la normale à } \hat{C}_{ig} \text{ (cf. n°4.2). On pose } i_g = i_g + 1 \text{ et on} & \\
\text{retourne en (4.33).} &
\end{array}
\right.
$$

Remarque 4.4 : Dans l'algorithme (4.27) - (4.38), l'avant du profil est
systématiquement maintenu cônique d'angle 90°, car les équations de la couche-
limite dégénèrent à l'avant. De même pour éviter tout décollement de la couche-
limite on impose systématiquement à l'arrière des profils d'être éffilés.

5. DESCRIPTION ET ANALYSE DES RESULTATS NUMERIQUES.

5.1. Description de l'exemple traité.

On cherche à obtenir le profil symétrique, d'aire donnée $a$,
minimisant la traînée dans un écoulement compressible ; laminaire, uniforme
plan, et à grand nombre de Reynolds ($10^2 < \pi < 10^5$).

5.2. Valeurs numériques des paramètres : Initialisation de la méthode itérative.

Pour des raisons d'homogénéité dans les formules ont peut prendre

$$
u_o = 1, v = 1, \rho = 1.
$$

Le profil initial utilisé est une approximation du profil continu
$\hat{c}_o$ défini par

$$
\hat{c}_o = \begin{cases} \xi & | \xi \in \mathbb{R}^2, \xi_2 = 0.1(1 - \xi_1^2)(0.9\xi_1 + 1.1) \text{ si } -1 < \xi_1 < 0.9 \text{ et pour} \\ \xi_1 \in [0.9, 0.95] & \text{un demi-dièdre d'angle } 90^\circ \text{ se raccordant avec continuit\'e des d\'eriv\'ees premi\`eres au reste du profil} \end{cases}
$$
Le contour $\hat{\partial C_o}$ a été approché par le polygone dont les sommets sont les points $\hat{\partial C_o}$ d'abscisses

$$ \xi_1 = 1 - 2 \left( \frac{i-1}{N'} \right)^{\alpha_i}, \quad i = 1, \dots, N' $$

avec

$$ \alpha_i = 1.7 - 0.7 \left( \frac{i-1}{N'} \right), \quad i = 1, \dots, N' $$

et $N' = 47$.

l'aire de ce profil polygonal étant $a = 0.147$ et la longueur $L = 1.98$.

La résistance de trainée correspondante $\frac{\rho}{u_o}$ E ($C_o$) vaut 2.66 et on a $K_1 = 0.102, K_2 = 0.5 \times 10^{-2}$.

Les paramètres déterminent le réseau de discrétisation de la couche-limite soit $M = 250, k = 0.035, J = 58, J_o = 4$.

On a pris comme paramètre de surrelaxation $\omega_1 = 1.5$ (resp. $\omega_2 = 1.2$) pour l'algorithme (4.11)-(4.14) appliqué à la résolution de (4.8), (4.9) (resp. (4.22), (4.23)$^{(1)}$. Dans les tests d'arrêt pour les itérations de relaxation on a pris $\epsilon = 1.5 \times 10^{-5}$ ; quant au pas $\lambda$ dans la direction de descente il a été pris égal à 0.04.

### 5.3. Description des résultats.

Après 1 minute 30 secondes de calcul sur I.B.M. 370/168 (correspondant à 3 itérations !) on a obtenu le profil $\hat{C}_3$ et les courbes des figures 2,3, 4 ; la valeur de la fonctionnelle n'a pas sensiblement diminuée (on reviendra sur ce point au n°6) par contre Grad $E_j$ s'est uniformisé (surtout à l'avant, cf. fig.2 et 5) et on a $K_1 = 0.109, K_2 = 0.4 \times 10^{-2}$.

Un autre calcul a été effectué dans lequel $\hat{C}_o$ a été prolongé par un demi-plan rattaché au bord de fuite et situé sur l'axe de symétrie de $C_o$ (dans le sillage), ceci afin de réduire les effets dus aux discontinuités au bord de fuite du

<footnote>$^{(1)}$ i.e. l'algorithme (4.24) - (4.26).</footnote>
du profil ; les résultats sont représentés sur les figures 5,6,7.

### 5.4. Analyse des résultats.

Avec les paramètres choisis on peut estimer que la précision obtenue sur $u, \frac{\partial u}{\partial n}$ est de l'ordre de $10^{-3}$, ce qui donne une précision de $2 \times 10^{-2}$ sur la résistance de traînée. La précision sur $Q, \frac{\partial Q}{\partial n}$ est de l'ordre de $5 \times 10^{-2}$ (on peut s'en rendre compte par la figure 4 qui donne $w$, car la relation $Q=2w+u$ n'est pas préservée par la discrétisation), on obtient donc une précision de $5 \times 10^{-2}$ sur $K_1$.

Les précisions obtenues sont insuffisantes, il faudrait pouvoir les améliorer d'un facteur 10 environ pour pouvoir calculer un profil meilleur que $C_3$, en particulier pour améliorer l'arrière de $C_3$, car $\frac{\partial u}{\partial n}$ et $\frac{\partial Q}{\partial n}$ y sont petits. Il faut ajouter à cela les incertitudes résultant des discontinuités au bord de fuite (on note en comparant les fig. 2 et 5 l'influence du changement de conditions aux limites).

Néanmoins on constate les trois points suivants :

1. Si on part d'un profil initial plus long que $C_3$, l'algorithme engendre des profils plus courts.

2. $C_3$ est plus long que $C_0$ défini au n°5.2.

3. Grad E ($\sigma$) on peut donc estimer que $C_3$, compte-tenu de la précision sur $\frac{\partial u}{\partial n}$ et $\frac{\partial Q}{\partial n}$, ne peut être amélioré par les procédés de calcul décrits ci-dessus (on est dans le bruit numérique) et que $C_3$ est à 5% environ de l'optimum.

## 6. REMARQUES DIVERSES ET CONCLUSIONS.

### 6.1. Calcul des profils optimaux pour d'autres volumes et d'autres nombres de Reynolds.

Soit $\hat{C}_3^1$ ($z, \gamma$) la transformée de $\hat{C}_3$ par l'homothétie de centre $z$ et de rapport $\gamma$. On peut montrer qu'il existe $T(z,z')$ tel que $\hat{C}_3'(z,\gamma)$ soit obtenu par translation $T$ de $\hat{C}_3$ ($z,\gamma$). Donc les écoulements potentiels autour de $\hat{C}_3'(z',\alpha)$ et autour de $C_3$ ($z,\gamma$) sont les mêmes, et de plus la vitesse de glissement $U(s)$ est $\gamma$fois la vitesse de glissement autour de $\hat{C}_3$.
De l'homogénéité des équations de la couche-limite en fonction du nombre de Reynolds on déduit que la condition d'optimalité du corrolaire 2.1 est aussi satisfaite autour de $\hat{C}_3'$ (z,γ) pour tout z et tout γ. On en déduit que tout _homothétique de_ $\hat{C}_3$ _est optimal dans sa classe_ (¹), et la résistance de traînée est donnée par

$$ (6.1) \qquad F(\hat{C}_3') = 1.33 \left(\frac{Lu}{\nu}\right)^{-\frac{1}{2}} \rho L u_0^2 $$

où L est la longueur de $\hat{C}_3'$ ; le volume de $\hat{C}_3'$ est 0.035 L². La formule (6.1) est à comparer avec la formule V.12 89 de Rosenhead, loc. cit. ; on voit que la résistance de traînée de $\hat{C}_3'$ est très proche de celle de la plaque plane de même longueur L (coefficient 1.328 au lieu de 1.33).

Les profils $\hat{C}_3$ ne dépendent pas de $\Re = \frac{Lu}{\nu}$ _par contre_ $C_3$ _obtenu en retranchant l'épaisseur de la couche-limite sur la normale à $\hat{C}_3$, _dépend de_ $\Re$ , mais avec la précision obtenue $C_3$ ne diffère apprécia- blement de $\hat{C}_3$ que si $\Re<3000$ (sauf peut-être au bord de fuite). On remarque donc que quelque soit le volume et le nombre de Reynolds, l'avant de $C_3$ est tangent un dièdre d'angle 90° (ce résultat est également vrai pour les très petits nombres de Reynolds , cf. Pironneau [6]).

## 6.2. Conclusions.

Par insuffisance de moyens de calcul (temps disponible essentiellement) on n'a pas pu améliorer de façon "spectaculaire" le profil initial $C_0$ (obtenu par tâtonnements) ; on vérifie cependant que bien que ce profil soit déjà très voisin de l'optimum on a pu améliorer l'avant de façon sensible ; il a donc été montré que la méthode était viable et indiqué comment la mettre en oeuvre. On a enfin donné un calcul du système _adjoint_ des équations de Navier-Stockes stationnaires, système qui apparaît naturellement dès que l'on utilise les méthodes variationnelles du contrôle optimal en mécanique des fluides visqueux incompressibles. Cette étude présente donc en détail les difficultés techniques que rencontrerais un laboratoire de mécanique des fluides désireux de traiter avec précision le problème ci-dessus en écoulement bi ou tridimensionnel ou les autres problèmes d'optimum désign en écoulement laminaire mentionnés dans Pironneau [2].

<footnote>(¹) Aux erreurs d'approximation près, cf. n°5.4.</footnote>
BIBLIOGRAPHIE

[1] LANDAU L. et E.LIFCHITZ. Mécanique des fluides - édition de Moscou (1963)

[2] PIRONNEAU O. On optimum désign in fluid mechanics. J.Fluid, Mech.64, 97-110 (1974)

[3] LIONS J.L. Contrôle optimal des systèmes gouvernés par des équations aux dérivées partielles. Dunod (1968).

[4] LUU T.S. Sur la technique des singularités en hydro et aérodynamique. C.R. du colloque CNRS de Marseille (1971).

[5] KELLER H. Computational problems in boundary layer theory. Fourth international conference on numerical methods in fluid mechanics. Boulder (1974).

[6] PIRONNEAU O.On optimum profiles in Stokes flow J.Fluid, Mech. 59, 117-128 (1973).

[7] GARABEDIAN P.Partial differential Equations. Wiley(1964)

[8] PIRONNEAU O. Optimisation de structure : application à la mécanique des fluides. I.F.I.P. conférence. Paris (1974).

[9] POLAK E. Computational methods in optimization Academic press (1971)

[10] ROSENHEAD L. Laminar boundary layers. Oxford University press (1963).

[11] C.E.A. J et R.GLOWINSKI Sur les méthodes d'optimisation par relaxation. Revue Française d'automatique -( Déc.1973, R.3, 5.32)
<figure><img src="image_2.png" /><figcaption>fig. 1</figcaption></figure>
<figure><img src="image_3.png" /></figure>

figure 2 : On a tracé la vitesse de glissement potentiel U(s), la friction $\frac{\partial u_s}{\partial n}$, la friction adjointe $\frac{\partial Q_s}{\partial n}$ et le gradient $\frac{\partial u_s}{\partial n}, \frac{\partial Q_s}{\partial n}$ pour le profil $\hat{c}_3$ dessiné dans le bas de la figure.
<figure><img src="image_4.png" /><figcaption>figure 4</figcaption></figure>

Pour $C_3$ on a tracé $\omega_5$ (n,s) pour 12 valeurs de s correspondant aux 12 premiers points "." de la figure 2.
<figure><img src="image_5.png" /><figcaption>figure 3</figcaption></figure>

Pour $\hat{C}_3$ on a tracé $u_s(n,s)$ pour 13 valeurs de s correspondant au "." de la figure 2.
<figure><img src="image_6.png" /><figcaption>figure 5
On a tracé toutes les caractéristiques de l'écoulement autour de C₀ en supposant que C₀ est prolongé par une demi plaque plane à l'arrière.</figcaption></figure>
<figure><img src="image_7.png" /></figure>

On a tracé toutes les caractéristiques de l'écoulement autour de C₀ en supposant que C₀ est prolongé par une demi plaque plane à l'arrière.
<figure><img src="image_8.png" /><figcaption>figure 7</figcaption></figure>

On a tracé toutes les caractéristiques de l'écoulement autour de C₀ en supposant que C₀ est prolongé par une demi plaque plane à l'arrière.