# Calcul des écoulements transoniques par des méthodes d'éléments finis et de contrôle optimal

Roland Glowinski, Olivier Pironneau

► To cite this version:

Roland Glowinski, Olivier Pironneau. Calcul des écoulements transoniques par des méthodes d'éléments finis et de contrôle optimal. [Rapport de recherche] IRIA-RR-138, IRIA. 1975, pp.44. hal-04716108

HAL Id: hal-04716108

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

Submitted on 1 Oct 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
IRIA
laboria

laboratoire de recherche
en informatique
et automatique

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

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

Calcul des écoulements
transoniques par des
méthodes d'éléments finis
et de contrôle optimal

R. Glowinski, O. Pironneau

Institut de Recherche
d'Informatique
et d'Automatique

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

Rapport de Recherche n°138
Octobre 1975
<figure><img src="image_4.png" /></figure>

# CALCUL DES ÉCOULEMENTS TRANSONIQUES PAR DES MÉTHODES D'ÉLÉMENTS FINIS ET DE CONTROLE OPTIMAL

R. Glowinski, O. Pironneau

## RESUME.

Dans ce travail on montre comment on peut résoudre les équations qui gouvernent les écoulements adiabatiques compressibles irrotationnels aux vitesses transoniques en leur associant des problèmes de contrôle optimal de systèmes à paramètres distribués. Ces problèmes sont discrétisés par éléments finis et résolus par des méthodes de gradient. On obtient ainsi une généralisation des méthodes de quasi-linéarisation classiques d'intégration par éléments finis des équations transoniques. Mais cette fois les schémas convergent aussi dans les régions supersoniques du fluide, en bi et en tridimensionnel, même en présence de chocs.

## ABSTRACT.

In this study a new method for solving the transonic equations for inviscid compressible potential flows is obtained. This method associates to each transonic equation a problem of optimal control for a distributed parameter system. These optimal control problems are discretized by a finite elements method and solved by a gradient algorithm. The resulting algorithm is a generalization of the classical methods for solving the transonic equations by the finite elements method except that these new schemes converge also in the supersonic regions and in the presence of shocks, for two and three dimensional flows.
Les écoulements autour de profils d'ailes ou dans les tuyères aux vitesses voisines de la vitesse du son sont régies par des équations non linéaires assez complexes (cf. LANDAU-LIFCHTZ [1], COURANT-FRIEDRICHS [2], RICHMEYER-MORTON [3], par exemple). La présence de discontinuités comme les chocs, le caractère mixte elliptique-hyperbolique des équations et la géométrie souvent très compliquée dans le cas des ailes d'avion, rendent l'intégration numérique de ces équations extrêmement difficile. Il existe déjà quelques programmes tridimensionnels opérationnels pour des écoulements mixtes subsonique-supersonique (cf. JAMESON [4], par exemple) ; mais ce sont des schémas aux différences finies donc peu pratiques pour des géométries compliquées.

Parmi les problèmes pratiques à résoudre on peut citer le calcul de l'écoulement autour d'une aile réelle (donc en plusieurs morceaux : ailerons, becs ...) et à géométrie variable. Les méthodes d'éléments finis se prètent bien à ce type de problème car on dispose de programmes de triangulation automatique. Malheureusement les méthodes en éléments finis comme celles de GELDER [5], DEVRIES-NORRIE [6] ne convergent pas dans les régions supersoniques du fluide.

On se propose dans cette étude d'utiliser les méthodes du contrôle optimal des systèmes gouvernés par des équations aux dérivées partielles (cf. LIONS [7]) pour traiter la non linéarité des équations transoniques. Le principe de base de la méthode est le suivant :

Soit à intégrer l'équation $\Delta\varphi=f(\varphi)$ dans un domaine $\Omega$ avec des conditions aux limites appropriées. Alors que les méthodes de quasi-linéarisation classiques consistent à calculer $\{\varphi^n\}$ tel que $\Delta\varphi^{n+1}=f(\varphi^n)$ (cf. PERRIAUX [8]); on va résoudre le problème :

$$ \min_{u} \left\{ \int_{\Omega} (u - f(\varphi))^{2} d\varphi \mid \Delta\varphi = u \right\} $$

par une méthode de gradient ou de gradient **conjugué**.
On se propose de développer cette méthode sur les trois formulations d'écoulements transoniques suivantes :

1. la formulation *divergentielle* : $\nabla \cdot \rho u = 0$ ; où la densité $\rho$ est reliée à la vitesse $u$ par une loi $\rho = G(|u|^2)$ et où $u$ dérive d'un potentiel $\varphi$.

2. la formulation *potentielle* : $\Delta\varphi = \frac{\nabla\varphi \cdot \nabla|\nabla\varphi|^2}{2c^2}$ où $c$ est la vitesse locale du son.

3. la formulation potentielle "linéarisée": $\alpha \frac{\partial \theta}{\partial x} + \frac{\partial^2 \theta}{\partial x^2} = \frac{\partial^2 \theta}{\partial y \partial x}$, où $\theta = \varphi - x$

Au §1 on rappelle les conditions d'application de ces trois équations. Puis au §2 on donne les formulations des problèmes de contrôle optimal associés à chaque équation ainsi que des conditions d'équivalence pour ces différents problèmes. Au §3 on obtient les conditions d'optimalité des problèmes de contrôle ce qui permet d'utiliser des méthodes de gradient pour leur résolution. L'approximation par éléments finis des algorithmes continus est exposé dans le §4.

Les résultats numériques pour les équations 1 et 3 sont exposés au § 5. Les essais numériques pour l'équation 2 seront exposés dans un article ultérieur. Pour la simplicité de l'exposé on s'est limité dans la partie théorique à des écoulements plans mais en fait ce qui suit vaut également pour des écoulements tridimensionnels. Les essais numériques ont été effectués sur des exemples bidimensionnels.
## Notation -

$\nabla \cdot \mathbf{u}$ : divergence u

$\nabla \times \mathbf{u}$ : rotationnel u

$\Delta \varphi$ : laplacien de $\varphi$

$\nabla \varphi$ : gradient $\varphi$

$H^m(\Omega)$ : espace de Sobolev d'ordre $m$ sur le domaine $\Omega$

$H^s(\Gamma)$ : espace de Sobolev d'ordre $s$ sur la frontière $\Gamma = \partial\Omega$

$|\cdot|_{L^\infty(\Omega)}$ : essup $|u(x)|$ $x \in \Omega$

## Plan -

### 0. Introduction

#### 1. Formulation du problème physique

1.1 Equations générales de l'écoulement

1.2 Cas des écoulements irrotationnels

1.3 Approximation des équations de l'écoulement irrotationnel

#### 2. Réduction des problèmes 1.16, 1.19; 1.20 à des problèmes de contrôle optimal.

2.1 Cas des problèmes mêlés de type Neuman Dirichlet

2.2 Cas des problèmes mêlés de type Cauchy

#### 3. Résolution des problèmes (2.3) (2.7) (2.9) par des méthodes de gradient

3.1 Résolution du problème (2.3)

3.2 Résolution du problème (2.7)

3.3.Résolution du problème (2.9)

#### 4. Approximation par éléments finis

4.1 Discrétisation et résolution du problème (2.9)

4.2.Discrétisation du problème (2.3.)

4.3 Discrétisation du problème (2.7)
5. Résolution numérique du problème (2.9)

5.1 Analyse des résultats relatifs à l'exemple 1

5.2 Analyse des résultats relatifs à l'exemple 2

5.3 Analyse des résultats relatifs à l'exemple 3

6. Remarques diverses et comparaison des méthodes

6.1 Cas des profils portants

6.2 Résolution numérique du problème (2.3)

6.3. Comparaison des méthodes

6.4 Conclusion.
# 1. FORMULATION DU PROBLEME PHYSIQUE

## 1.1 Equations générales de l'écoulement

On s'intéresse aux écoulements **adiabatiques stationnaires** de fluides **compressibles monophasiques** dans lesquels les **effets de viscosité sont négligés** (approximation des **fluides parfaits**).

Soit $\Omega$ le domaine occupé par le fluide du type ci-dessus.
On notera $\Gamma = \partial\Omega$, la frontière de $\Omega$. Les équations d'Euler du mouvement stationnaire sont alors données par

$$1.1 \quad \nabla \cdot (\rho \mathbf{u}) = 0 \quad (\text{conservation de la matière})$$

$$1.2 \quad \rho u \cdot \nabla u = -\nabla p \quad (\text{conservation de la quantité de mouvement})$$

$$1.3 \quad \rho u \cdot \nabla E = - p \nabla \cdot \mathbf{u} \quad (\text{conservation de l'énergie}).$$

où dans (1.1) - (1.3)

* $\rho$: densité du fluide

* $u$: vitesse de l'écoulement

* $p$: pression dans le fluide

* $E$: énergie interne par unité de volume

Les fonctions ci-dessus sont liées par **l'équation d'état**

$$(1.4) \quad E = F(\rho, p)$$

par exemple dans le cas d'un gaz parfait en écoulement adiabatique
on a

$$(1.5) \quad E = \frac{1}{\gamma-1} \frac{p}{\rho}$$

où $\gamma$ est le rapport des chaleurs spécifiques à pression et volume constants, soit $\gamma = C_p/C_v(\gamma)$

<footnote>(+) pour les gaz monoatomiques $\gamma = 5/3$ pour les gaz diatomiques $\gamma = 7/5$.</footnote>
Il convient d'ajouter aux équations ci-dessus des conditions aux limites "appropriées" ; par exemple

(1.6) $u, n = 0$ sur les parois des obstacles n étant le vecteur normal à l'obstacle.

On rencontrera également des conditions aux limites du type (1.7) $u$ et/ou $\rho u$ donnés sur une partie de $\Gamma$.

### Remarque 1.1

Les équations (1.2) (1.3) peuvent, à l'aide de (1.1), se réécrire :

$$ (1.8) \qquad \nabla \cdot (\mathbf{P} \cdot \mathbf{u} + (\mathbf{u} \cdot \mathbf{P})\mathbf{u}) = 0 $$

$$ (1.9) \qquad \nabla \cdot \mathbf{P} \mathbf{u} + \frac{\mathbf{P}}{\rho} + \frac{\mathbf{u}^2}{2} = 0 $$

ou $\mathbf{I}$ est le tenseur identité et $(\mathbf{u} ; \mathbf{u})$ le tenseur de composantes $u_i u_j$.

Les écoulements décrits par les équations ci-dessus peuvent posséder des **discontinuités internes**, comme les **chocs**, sur lesquelles il est nécessaire de respecter certaines conditions aux limites (condition de Rankine-Hugoniot). Soit plus précisement dans le cas des chocs, $n$ désignant la normale unitaire à la surface de discontinuité :

$$ (1.10) \qquad \rho_1 u_1 n = \rho_2 u_2 \cdot n $$

$$ (1.11) \qquad \rho_1 u_1 n (E_1 + \frac{(u_1 \cdot n)^2}{2} + p_1 / \rho_1) = \rho_2 u_2 \cdot n (E_2 + \frac{(u_2 \cdot n)^2}{2} + p_2 / \rho_2) $$

$$ (1.12) \qquad p_1 + \rho_1 (u_1 \cdot n)^2 = p_2 + \rho_2 (u_2 \cdot n)^2 $$

Par définition d'une onde de chocs les composantes tangentielles de $\mathbf{u}$ restent continues.

Il convient de signaler qu'au passage d'une onde de chocs il y a augmentation de l'entropie ($\dagger$).

<footnote>(†) La densité d'entropie par unité de volume $s$ est reliée à $E$ par la relation thermodynamique $dE = T ds - pd(\frac{1}{\rho})$ à laquelle il faut ajouter une autre équation d'état reliant la température $T$ à $p$ et $\rho$ ($p/\rho = R T$ pour un gaz parfait).</footnote>
On renvoie par exemple à LANDAU-LIFCHITZ [1], RITCHMYER-MORTON [3], ou COURANT-FRIEDRICHS [2] pour une discussion plus approfondie des phénomènes d'ondes de chocs en mécanique des fluides.

## 1.2 Cas des écoulements irrotationnels

Parmi les solutions du système (1.1)-(1.3) celles qui satisfont à

$$ (1.11) \quad \nabla \times u = 0 $$

jouent un rôle particulièrement important ; de tels écoulements sont dits "irrotationnels". De tels écoulements ne sont pas rares car on peut montrer à partir de (1.1)-(1.3) que si l'écoulement est irrotationnel dans une région $\Omega_i$ de $\Omega$ alors il est aussi irrotationnel en tout point de $\Omega$ traversé par une ligne de courant issue de $\Omega_i$. Il faut noter qu'en toute rigueur un écoulement irrotationnel perd cette propriété au passage d'une onde de choc (cf. [1, ch.12] toutefois il est montré dans [1, ch.9 et 12] que si le choc est de faible intensité il est raisonnable de supposer la conservation de la propriété d'irrotationnalité.

Si le fluide est un gaz parfait la forme particulière de F ($p,e$) et la condition (1.11) permettent d'obtenir des intégrales premières du système (1.1)-(1.3) de la façon suivante :

De (1.9),(1,4) on déduit à l'aide de (1.1),(1,2),(1.11) :

$$ (1.12) \quad \nabla \left( \frac{1}{2} u^2 + \frac{\gamma}{\gamma-1} \frac{p}{\rho} \right) = 0 $$

De (1.2) (1.11) on déduit

$$ \rho \nabla \left( \frac{|u|^2}{2} \right) = \nabla p $$

Soit, après élimination de $p$ à l'aide de (1.12)

$$ (1.13) \quad \frac{\rho}{(\gamma-1)} \nabla \left( \frac{|u|^2}{2} \right) - \left( K - \frac{|u|^2}{2} \right) \nabla p $$

où K est une constante provenant de l'intégration de (1.12).
L'équation (1.13) s'intègre à vue, on a l'habitude de
l'écrire sous la forme

$$
(1.14) \quad \rho = \rho_0 \left( 1 - \frac{\gamma-1}{\gamma+1} \frac{|u|^2}{c_*^2} \right)^{\frac{1}{\gamma-1}}
$$

La constante C* est connue sous le nom de vitesse critique
du son.

De (1.11) on déduit l'existence d'une **fonction potentielle** φ telle
que

$$
(1.15) \quad u = \nabla \varphi
$$

Compte-tenu de (1.14), (1.15) le système (1.1) (1.4) équivaut
donc à

$$
(1.16) \quad \nabla \cdot (\rho \nabla \varphi) = 0
$$

$$
(1.17) \quad \rho = G(|\nabla \varphi|^2)
$$

où G (|u|^2) est la fonction définie par (1.14).

Des conditions de saut (1.10)-(1.12) au passage d'une onde
de choc on gardera la condition :

(1.18) $\varphi$ continue

Cependant il faut noter que le système (1.16) (1.17) n'admet
**pas toujours de solution unique** ; parmi celles qui ont des chocs
on ne pourra donc retenir que celles qui ont un **saut d'entropie croissant**.
1.3 Approximation des équations de l'écoulement potentiel

On peut, au moins formellement, éliminer $\rho$ de (1.16), (1.17)

Dans le cas d'un gaz parfait on obtient :

$$ (1.19) \quad - \Delta \varphi + \frac{\nabla \varphi \cdot \nabla (|\nabla \varphi|^2)}{\gamma + 1 - (\gamma - 1) |\nabla \varphi|^2} = 0 $$

(on a changé $\varphi$ en $C_* \varphi$)

Equation à laquelle il fait ajouter les conditions aux limites convenables ainsi que la condition (1.17).

Au cas où dans $\Omega$ la vitesse de l'écoulement diffère assez peu en grandeur et en direction d'une vitesse fixe $u_0$, prise comme direction $Ox$, l'équation (1.19) est assez bien approchée par

$$ (1.20) \quad 2\alpha_* \frac{\partial \theta}{\partial x} \frac{\partial^2 \theta}{\partial x^2} - \frac{\partial^2 \theta}{\partial y^2} = 0 $$

ou $\alpha_* := \frac{\gamma+1}{2}$ et $\theta = \varphi - x$
## 2. REDUCTION DES PROBLEMES (1.16); (1.19), (1.20) A DES PROBLEMES DE CONTROLE OPTIMAL.

### 2.1 Cas des problèmes mêlés de type Neuman-Dirichlet

L'écoulement dans $\Omega$ étant irrotationnel, on va supposer que les conditions aux limites sur $\varphi$ sont du type suivant :

$$ (2.1) \qquad \varphi|_{\Gamma_1} = \varphi_1 $$

$$ (2.2) \qquad \frac{\partial \varphi}{\partial n} |_{\Gamma_2} = u_2 $$

$$ \text{avec } \Gamma_1 \cup \Gamma_2 = \Gamma_1 \cup \Gamma_2 = \emptyset, \Gamma_1 \neq \emptyset $$

On suppose $\Omega$ borné et $\Gamma$ régulière.

On considère le problème de contrôle optimal suivant :

$$ (2.3) \quad \min_{\varphi \in \Phi_m} \left\{ \int_{\Omega} |\nabla \varphi|^2 d\Omega \left| -\Delta \varphi = \nabla \cdot (G(|\nabla \varphi|^2) \nabla \varphi), \varphi|_{\Gamma} = 0 \right. \right\} $$

$$ \text{où } \Phi_m = \left\{ \varphi \in H^1(\Omega) \mid \varphi|_{\Gamma_1} = \varphi_1, \frac{\partial \varphi}{\partial n}|_{\Gamma_2} = u_2 ; |\nabla \varphi(\cdot)| \le m \quad \text{pp x } \varphi \in \Omega \right\} $$
et où $m$ est une constante donnée telle que $m < C_* \left(\frac{\gamma+1}{\gamma-1}\right)^{\frac{1}{2}}$

#### **Proposition 2.1**

On suppose que :

$$ (2.4) \quad \text{Il existe } K \text{ tel que } 0 \le G(|u|^2) \le K \text{ pour tout } u \in [0, n] $$

$$ (2.5) \qquad \varphi_1 \in H^1(\Gamma_1) ; \quad u_2 \in H^2(\Gamma_2) $$
$$
(2.6) \left\{
\begin{array}{@{}l@{}}
\text{Le problème (1.16), (1.17), (2.1), (2.2) où les égalités} \\
\text{sont prises respectivement au sens de } H^{-1}(\Omega), L^2(\Omega), H^2(\Gamma_1), H^{-1/2}(\Gamma_2) \\
\text{admet au moins une solution telle que } |\nabla \varphi(x)| \le m \quad \text{pp } x \in \Omega.
\end{array}
\right.
\quad (\dagger)
$$

Alors toute solution de (2.3) est solution de (1.16)-(2.2).
Réciproquement toute solution de (1.16)-(2.2) au sens ci-dessus telle
que $|\nabla \varphi(x)| \le m \quad \text{pp } x \in \Omega$ est solution de (2.3).

Démonstration
---

Pour tout $\varphi \in \Phi_m$, $\varepsilon$ défini dans (2.3) existe et appartient
à $H^1(\Omega)$. Or $G(|\nabla \varphi|^2) \in L^\infty(\Omega)$ donc, $G(|\nabla \varphi|^2)_{\nabla \varphi}$ étant dans $(L^2(\Omega))^2$
et $\Omega$ étant borné, $\cdot \in (G(|\nabla \varphi|^2)_{\nabla \varphi}) \in H^{-1}(\Omega)$.

Le problème (2.3) est donc bien posé (cf. Lions [7]).

Toute solution $\bar{\varphi}$ de (1.16).(2.2) est une solution de (2.3)
car pour ce $\bar{\varphi}$, $\nabla \varepsilon = 0$. D'autre part toute solution de (2.3) devant
être telle que $\nabla \varepsilon = 0$, on doit avoir $\nabla \cdot (G(|\nabla \varphi|^2) \nabla \varphi) = 0$ dans $H^{-1}(\Omega)$.
Car l'opérateur $\Delta$ est continu de $H^1(\Omega)$ dans $H^{-1}(\Omega)$.

Remarque 2.1
---

On peut aussi en relation avec (1.16) résoudre le problème

$$
\min_{\rho} \left\{ \int_{\Omega} [G(|\nabla \varphi|^2) - \rho]^2 d\Omega \middle| \nabla \cdot (\rho \nabla \varphi) = 0, \varphi|_{\Gamma_1} = \varphi_1, \frac{\partial \varphi}{\partial n}|_{\Gamma_2} = u_2 \right\}
$$

Ce problème est un problème de contrôle optimal où le contrôle
intervient dans les coefficients de l'équation aux dérivées partielles.
Il sera étudié dans un travail ultérieur.

(†) On rappelle que a = b au sens de V' dual de V si < a, v > = C b, v >
pour tout v ∈ V.
De façon analogue, on considère en relation avec (1.19)
le problème

$$
(2.7) \quad \min_{u \in \mathcal{U}} \left\{
\begin{array}{@{}l@{}}
  \displaystyle \int_{\Omega} |u - \nabla \varphi|^2 d\Omega \quad \left| \begin{array}{@{}l@{}}
  -\Delta \varphi = - \frac{-1}{\gamma+1} |u|^2 \nabla \cdot u - \frac{1}{\gamma+1} u \cdot \nabla (|u|^2) \\
  \qquad \varphi|_{\Gamma_1} = \varphi_1, \quad \frac{\partial \varphi}{\partial u}|_{\Gamma_2} = u_2
  \end{array} \right\} \\
  \text{où } u = \left\{ u \in (L^\infty(\Omega))^2 \right\} \quad (\dagger)
\end{array}
\right.
$$

$$
\text{où } u = \left\{ u \in (L^{\infty}(\Omega))^2 \right\} \quad (\dagger)
$$
et on a alors la

**Proposition 2.2**

On suppose que (2.5) est vérifiée que (1.19), (2.1), (2.2), prises
au sens de $H^{-1}(\Omega)$, $H^2(\Gamma_1)$, $H^{-2}(\Gamma_2)$ admet au moins une solution.

Alors toute solution de (2.7) est solution de (1.1)-(2.2).
Réciproquement toute solution de (1.19)-(2.2) tel que $\nabla \varphi \in (L^\infty(\Omega))$ est
solution de (2.7).

Démonstration.

Comme précédemment pour $u \in \mathcal{U}, \Delta\varphi \in H^{-1}(\Omega)$ donc (2.7) est
bien posé et a l'optimum $u = \nabla\varphi$. Il faut cependant noter qu'il peut
exister des suites minimisantes de (2.7) ayant pour limite $\bar{u} \notin (L^\infty(\Omega))^2$.

Dans ce cas le potentiel correspondant $\bar{\varphi}$ peut ne pas être une
solution de (1.19)-(2.2). Pour éviter cette situation il faudra, lors
de la résolution numérique de (2.7), imposer une contrainte du type
$|u(x)| < m \quad \text{pp dans } \Omega$.

Remarque 2.2

L'équation (1.19) peut aussi être comparée au problème

$$
(2.8) \quad \min_{v \in L^\infty(\Omega)} \left\{
\begin{array}{@{}l@{}}
  \displaystyle \int_{\Omega} \left[ v \left(1 - \frac{\gamma-1}{\gamma+1} |\nabla\varphi|^2 \right) + \frac{1}{\gamma+1} \nabla\varphi \cdot \nabla (|\nabla\varphi|^2) \right]^2 d\Omega \quad \left| \begin{array}{@{}l@{}}
  \Delta\varphi = v, \\
  \varphi|_{\Gamma_1} = \varphi_1, \quad \frac{\partial\varphi}{\partial n}|_{\Gamma_2} = u_2
  \end{array} \right\} \\
  \text{où } u = \nabla\varphi.
\end{array}
\right.
$$

<footnote>(†) On peut aussi ajouter une contrainte $\nabla \times u = 0$ ou $u = \nabla \Psi$</footnote>
Enfin, pour résoudre (1.20) nous allons considérer le problème

$$
(2.9) \min_{v \in L^2(\Omega)} \left\{ \int_{\Omega} \left( \alpha \left( \frac{\partial \varphi}{\partial x} \right)^2 + \frac{\partial \varphi}{\partial x} - v \right)^2 d\Omega \left| \begin{array}{l} \Delta \varphi = \frac{\partial v}{\partial x}, \quad \varphi|_{\Gamma_1} = \varphi_1, \quad \frac{\partial \varphi}{\partial n}|_{\Gamma_2} = u_2 \\ \frac{1}{H^2}(\Gamma_2) \text{ admet au moins une solution telle que } \frac{\partial \varphi}{\partial x} \in L^4(\Omega) \text{ alors toute solution de (2.9) est solution de } (1.20)-(2.2) \text{ et réciproquement.} \end{array} \right. \right\}
$$

**Proposition 2.3**

Si (1.20), (2.1), (2.2), prises au sens de $H^1(\Omega)$, $H^2(\Gamma_1)$, $\frac{1}{H^2}(\Gamma_2)$ admet au moins une solution telle que $\frac{\partial \varphi}{\partial x} \in L^4(\Omega)$ alors toute solution de (2.9) est solution de $(1.20)-(2.2)$ et réciproquement.

Démonstration :

Identique aux précédentes car (1.20) se récrit :

$$
\Delta\varphi = \frac{\partial}{\partial x}\left(\alpha\left(\frac{\partial\varphi}{\partial x}\right)^2 + \frac{\partial\varphi}{\partial x}\right) ; \text{ et d'autre part ou bien } \frac{\partial\varphi}{\partial x} \in L^4(\Omega) \\
\text{ou bien le critère vaut } +\infty
$$

2.2 Cas des problèmes mêlés de type Cauchy

On s'intéresse aux écoulements transoniques dans les tuyères
avec sorties des écoulements supersoniques. On se donne alors φ et ∂φ/∂n
à l'entrée de la tuyère et ∂φ/∂n sur les parois. Même avec ce type de
de condition aux limites on peut utiliser le problème (2.3) à condition
de prendre

$$
\Phi_m = \left\{ \varphi \in H^1(\Omega) \middle| \begin{array}{l} \varphi|_{\Gamma_1} = \varphi_1, \quad \frac{\partial \varphi}{\partial n}|_{\Gamma_1} = u_1, \quad \frac{\partial \varphi}{\partial n}|_{\Gamma_2} = u_2; \\ |\nabla \varphi(\Omega)| \le m \quad \text{pp dans } \Omega \end{array} \right\}
\quad \text{(voir ci-dessous pour les définitions de } \Gamma_1, \Gamma_2, u_1, u_2).
$$

Pour les problèmes (1.19) et (1.20) on peut aussi résoudre
(2.7) et (2.9) avec les conditions aux limites appropriées. Mais on
peut aussi penser à utiliser l'opérateur $\frac{\partial^2}{\partial x^2} - \frac{\partial^2}{\partial y^2}$, mieux adapté
aux conditions aux limites.
On suppose donc que $\Gamma = \Gamma_1 U\Gamma_2 U\Gamma_3$.

Ainsi pour résoudre (1.20) avec les conditions aux limites

$$
(2.10) \quad \varphi|_{\Gamma_1} = \varphi_1
$$

$$
(2.11) \quad \left. \frac{\partial \varphi}{\partial n} \right|_{\Gamma_1} = u_1, \quad \left. \frac{\partial \varphi}{\partial n} \right|_{\Gamma_2} = u_2
$$

On peut considérer le problème

$$
(2.12) \quad \lim_{v \to \nu} \left\{ \int_{\Omega} |v - \frac{\partial \varphi}{\partial x}|^2 d\Omega \left| \begin{array}{l} \frac{\partial^2 \varphi}{\partial x^2} - \frac{\partial^2 \varphi}{\partial y^2} = \frac{\partial v}{\partial x} ; \quad \varphi|_{\Gamma_1} = \varphi_1 \\ \\ \frac{\partial \varphi}{\partial n}|_{\Gamma_1} = u_1, \quad \frac{\partial \varphi}{\partial n}|_{\Gamma_2} = u_2 \end{array} \right. \right\}
$$

Malheureusement la théorie des équations différentielles
hyperboliques sur des ouverts non cylindriques ne permet pas de
développer une théorie générale appropriée à l'étude du problème
(2.12) ; en particulier pour certain v, l'équation différentielle
peut ne pas avoir de solution.

Nous avons cependant tenu à mentionner cette approche car les
Laboratoires Aéronautiques possèdent en général des sous-programmes
d'intégration numérique d'équations différentielles du type ci-dessus
(celle qui régit les écoulements supersoniques, par exemple) il leur
est donc assez facile de résoudre le problème (2.12) par les méthodes
du paragraphe suivant.
### 3. RESOLUTION DES PROBLEMES (2.3), (2.7), (2.9) PAR DES METHODES DE GRADIENT.

#### 3.1 Résolution du problème (2.3)

Le problème (2.3) est un problème d'optimisation de la forme

$$ (3.1) \quad \min \left\{ \mathcal{E}(\varphi) \middle| \varphi \in \Phi_m \right\} \\ \text{où} $$

$$ (3.2) \quad \mathcal{E}(\varphi) = \int_{\Omega} |\nabla \varepsilon(\varphi)|^2 d \Omega \\ \text{et où } \varepsilon(\varphi) \text{ est l'unique solution de} $$

$$ (3.3) \quad -\Delta\varepsilon = \nabla \cdot (G(|\nabla\varphi|^2 \nabla\varphi)) \quad \varepsilon|_{\Gamma} = 0 $$

On peut montrer (cf. Lions [7]) que dans certaines conditions $\mathcal{E}$ est une forme continûment différentiable sur $\Phi_m$ munie de la topologie de $H^1(\Omega)$. Le gradient de $\mathcal{E}$ se calcule par identification des termes linéaires du premier ordre en $\delta\varphi$ de la quantité $E(\varphi+\delta\varphi) - E(\varphi)$. Dans ce qui suit on suppose $G$ continûment différentiable dans $[0,m]$. Pour $\varphi, \varphi + \delta\varphi \in \Phi_m$ on a :

$$ (3.4) \quad \mathcal{E}(\varphi+\delta\varphi) - \mathcal{E}(\varphi) = \int_{\Omega} 2\nabla \cdot \left(\varepsilon(\varphi)\nabla\cdot\left(\varepsilon(\varphi+\delta\varphi) - \varepsilon(\varphi)\right)\right) d\Omega + \int_{\Omega} \left|\nabla\cdot\left(\varepsilon(\varphi+\delta\varphi) - \varepsilon(\varphi)\right)\right|^2 d\Omega $$

D'après (3.3) on note que $\delta\varepsilon = \varepsilon(\varphi+\delta\varphi) - \varepsilon(\varphi)$ satisfait à

$$ (3.5) \quad -\Delta\delta\varepsilon = \nabla \cdot \left[ G\left(|\nabla\varphi|^2\right) \nabla\cdot\left(\varepsilon\varphi + 2G'\left(|\nabla\varphi|^2\right)\nabla\varphi \cdot \nabla\delta\varphi\right) \nabla\varphi + o(\delta\varphi) \right] $$

Car $G$ est continûment différentiable. En multipliant (3.5) par $\varepsilon(\varphi)$, et en intégrant par partie, on obtient

$$ (3.6) \quad \int_{\Omega} \nabla \varepsilon \cdot \delta \varepsilon \, d\Omega = - \int_{\Omega} \left[ G\left(|\nabla \varphi|^2 \nabla \varepsilon \cdot \nabla \delta \varphi\right) + 2G'\left(|\nabla \varphi|^2\right) \left( \nabla \varphi \cdot \nabla \varepsilon \right) \left( \nabla \varphi \cdot \nabla \delta \varphi \right) \right] d\Omega \\ + o(\delta\varphi) $$
et d'après (3.4) et le fait que $o(\delta\epsilon) = o(\delta\varphi)$

$$ (3.7) \quad E(\varphi+\delta\varphi) - E(\varphi) = 2 \int_{\Omega} [G(|\nabla\varphi|^2)\nabla\epsilon \cdot \nabla\delta\varphi + 2G'(|\nabla\varphi|^2)(\nabla\varphi \cdot \nabla\epsilon)(\nabla\varphi \cdot \nabla\delta\varphi] \, d\Omega \\ \phantom{(3.7) \quad E(\varphi+\delta\varphi) - E(\varphi) = 2 \int_{\Omega} [} + o(\delta\varphi) $$

qui, après intégration par partie peut aussi se récrire

$$ (3.8) \quad E(\varphi+\delta\varphi) - E(\varphi) = 2 \int_{\Omega} \nabla \cdot [G(|\nabla\varphi|^2)\nabla\epsilon + 2G'(|\nabla\varphi|^2)(\nabla\varphi \cdot \nabla\epsilon)\nabla\varphi] \, \delta\varphi \, d\Omega \\ \qquad - \int_{\Gamma} [G(|\nabla\varphi|^2) \frac{\partial\epsilon}{\partial n} + 2G'(|\nabla\varphi|^2)(\nabla\varphi \cdot \nabla\epsilon) \frac{\partial\varphi}{\partial n}] \, \delta\varphi \, d\Gamma + o(\delta\varphi) $$

Donc le "gradient" de $E$ (défini par $E(\varphi+\delta\varphi) - E(\varphi) = \langle \nabla E(\varphi), \delta\varphi \rangle + o(\delta\varphi)$ est

$$ (3.9) \quad \nabla E(\varphi) = \nabla \cdot (G(|\nabla\varphi|^2)\nabla\epsilon + 2G'(|\nabla\varphi|^2)(\nabla\varphi \cdot \nabla\epsilon)\nabla\varphi) - \delta_{\Gamma}[G(|\nabla\varphi|^2) \frac{\partial\epsilon}{\partial n} \\ \qquad + 2G'(|\nabla\varphi|^2)(\nabla\varphi \cdot \nabla\epsilon) \frac{\partial\varphi}{\partial n}] $$

où $\delta_{\Gamma}$ est la fonction de Dirac au bord $\Gamma$ de $\Omega$.

Pour résoudre (2.3) on peut donc utiliser l'algorithme suivant

### Algorithme 3.1

$$ (3.10) \quad \text{Choisir } \varphi^0 \in \Phi_m, \text{ poser } k = 0 $$

$$ (3.11) \quad \text{Calculer } \varepsilon^k \text{ solution de } - \Delta \varepsilon^k = \nabla \cdot [G(|\nabla \varphi|^2) \nabla \varphi^k], \quad \varepsilon^k|_{\Gamma} = 0 $$

$$ (3.12) \quad \text{Calculer à l'aide de (3.7) où (3.8) la solution } d^k \text{ de} \\ \min \left\{ \langle \nabla E(\varphi^k), d \rangle \mid \int_{\Omega} |\varphi|^2 \, d\Omega \le 1, d|_{\Gamma_1} = 0, \frac{\partial d}{\partial n}|_{\Gamma_2} = 0 \right\} $$

$$ (3.13) \quad \text{Calculer } \lambda^k \text{ solution de} \\ \min_{\lambda} E(\varphi^k + \lambda^k d^k) - E(\varphi^k) $$

$$ (3.14) \quad \text{poser } \varphi^{k+1} = \varphi^k + \lambda^k d^k, k = k+1 \text{ et retourner en (3.11).} $$

### Remarque 3.1

On n'a pas incorporé la contrainte $|\nabla\varphi|_{L^\infty} \le m$ dans l'algorithme pour des raisons de simplicité de l'exposé.
L'algorithme ci dessus est une méthode de gradient (Frank) et Wolfe pour être exact) appliquée au problème (3.1) dans le sous espace $\Phi_{\infty}$. Si l'étape (3.13) peut être calculée exactement on aura intérêt à utiliser une méthode de gradient conjugué. On rappelle que dans le cas sans contrainte la méthode du gradient conjugué consiste à prendre :

$$ (3.15) \qquad \varphi^{k+1} = \varphi^k + \lambda h^k $$

avec

$$ (3.16) \qquad h^k = - \nabla E(\varphi^k) + \frac{h^{k-1}}{\| \nabla E(\varphi^{k-1}) \|^2} \langle \nabla E(\varphi^k), \nabla E(\varphi^k) - \nabla E(\varphi^{k-1}) \rangle $$

### Remarque 3.2

Il existe des méthodes spécifiques dont nous parlerons au § 4 pour résoudre (3.12) et (3.13).

### Proposition 3.1

On se place dans les conditions d'application de la proposition (2.1). On suppose de plus que $G$ est continûment différentiable dans le domaine de travail de l'algorithme. Alors tout point d'accumulation $\bar{\varphi}$ d'une suite $\varphi^k$ engendrée par l'algorithme (3.1) tel que $|\bar{\varphi}|_{L^\infty} < m$ ; $|\varphi^k|_{L^\infty} < m$ est un point stationnaire du problème (2.3). Si de plus $\varepsilon > 0$ tend vers zéro alors $\bar{\varphi}$ est une solution de (1.16)-(2.2).

### Démonstration

La démonstration de la convergence d'un algorithme de gradient est classique : elle repose sur la continuité différentiabilité de l'application $E(\cdot, \cdot)$, (cf. Yvon [9] pour les problèmes de contrôle distribué). Toutefois nous redonnons l'essentiel de l'argument de la démonstration. $E(\varphi^k + \lambda d^k)$ est continûment différentiable en $\lambda$ et admet pour dérivée à l'origine $\langle \nabla E(\varphi^k), d^k \rangle >$. Donc d'après (3.13) si $\langle \nabla E \varphi^k, d^k \rangle >$ ne tend pas vers zéro $\lambda^k$ est borné par $\bar{\lambda}$ donc d'après (3.13) et (3.8)
$$E(\varphi^{k+1}) - E(\varphi^k) \le E(\varphi^k + \bar{\lambda} d^k) - E(\varphi^k) \le \bar{\lambda} \langle \nabla E(\varphi^k), d^k \rangle \\ + o(\bar{\lambda} d^k)$$

comme d'après (3.12) $||d^k|| = 1$, il existe $\bar{\lambda}$ petit tel que

$$E(\varphi^{k+1}) - E(\varphi^k) \le - \frac{\bar{\lambda}}{2} \langle \nabla E(\varphi^k), d^k \rangle$$

Comme la suite $E(\varphi^k)$ est strictement décroissante ceci est incompatible avec $E(\bar{\varphi}) \ge 0$ donc $\lim_{k \to \infty} \langle \nabla E(\varphi^k), d^k \rangle = 0$. D'après (3.7) $\langle \nabla E(\varphi), d \rangle$ est continue en $\varphi$ donc $\langle \nabla E(\bar{\varphi}), d \rangle = 0$.

### Remarque 3.3

Tous les points stationnaires de (2.3) ne sont pas des solutions de (1.16)-(2.2). La proposition (2.1) peut donc sembler sans intérêt. Toutefois ce type de difficulté est classique en optimisation. L'expérience a montré que si $\tilde{\varphi}$ est une solution de (1.16)-(2.2) alors après un nombre fini d'essais on peut en changeant $\varphi^0$ obtenir en général $\tilde{\varphi}$ comme limite de $\{\varphi^k\}$

## 3.2 Résolution du problème (2.7)

Tout comme (2.3), le problème (2.7) est un problème d'optimisation en dimension infinie. Pour le résoudre par une méthode de descente il faut calculer les termes linéaires du premier ordre en $\delta u$ de la variation du critère $\delta E = E(u+\delta u) - E(u)$ où

$$(3.19) \quad E(u) = \int_{\Omega} |u - \nabla \varphi(u)|^2 d\Omega$$

$\varphi(u)$ étant la solution du système différentiel dans (2.7). Ce calcul se fait de la même façon que pour (2.3). Cependant l'élimination des termes en $\delta\varphi = \varphi(u+\delta u) - \varphi(u)$ dans $\delta E$ nécessite l'introduction d'une fonction intermédiaire $\psi$ solution d'un système différentiel adjoint de $\varphi$.
Plutôt que de donner une démonstration rigoureuse de l'évaluation
de $\delta E$ on se propose de faire un calcul formel calqué sur les méthodes
du calcul des variations souvent plus familières aux lecteurs. Il
convient de noter qu'une justification rigoureuse de ce qui suit est
très délicate et qu'il vaut mieux, d'un point de vue théorique, opérer
comme dans LIONS [7].

Soit $\Psi(x)$ le multiplicateur de Lagrange associé à l'équation
différentielle dans (2.7). L'hamiltonien du système s'écrit :

$$
\mathcal{H} = \int_{\Omega} |u - \nabla \varphi|^2 d \Omega + \int_{\Omega} \left[ (-\Delta \varphi + \frac{\gamma-1}{\gamma+1} |u|^2 \nabla \cdot u + \frac{1}{\gamma+1} u \cdot \nabla (|u|^2)) \right] \Psi d \Omega
$$

où encore, d'après la formule de Green

$$
\mathcal{H} = \int_{\Omega} |u - \nabla \varphi|^2 d\Omega + \int_{\Omega} \nabla \varphi \Psi \, d\Omega + \int_{\Omega} \left[ \frac{\gamma-1}{\gamma+1} |u|^2 \nabla \cdot u + \frac{1}{\gamma+1} u \cdot \nabla (|u|^2) \right] \Psi \, d\Omega \\
+ \int_{\Gamma} \frac{\partial \varphi}{\partial n} \Psi \, d\Gamma
$$

$\mathcal{H}(u,\varphi,\Psi)$ est stationnaire en $\varphi$ et $\Psi$, à l'optimum, donc

$\mathcal{H}(u,\varphi + \delta\varphi, \Psi)$ - $\mathcal{H}(u,\varphi, \Psi)$ est du 2ème ordre en $\delta\varphi$ ce qui entraîne,
compte tenu des conditions aux limites sur $\varphi$

$$
\int_{\Omega} 2(u - \nabla \varphi) \cdot \nabla \delta\varphi \, d\Omega + \int_{\Omega} \nabla \delta\varphi \nabla \Psi \, d\Omega + \int_{\Gamma} \frac{\partial \delta\varphi}{\partial n} \Psi \, d\Gamma = 0
$$

c'est à dire, $\Psi|_{\Gamma_1} = 0$ et

$$
(3.20) \quad \int_{\Omega} \nabla \Psi \nabla \omega \, d\Omega = \int_{\Omega} -2(u - \nabla \varphi) \omega \, d\Omega \quad \text{pour tout } \omega \in H^1(\Omega), \omega|_{\Gamma_1} = 0
$$
On reconnaît la formulation variationnelle de l'équation

$$
(3.21) \quad - \Delta \Psi = + 2 \nabla \cdot (u - \nabla \varphi), \quad \left. \frac{\partial \Psi}{\partial n} \right|_{\Gamma_2} = - 2 (u_n |_{\Gamma_2} - u_2), \quad \Psi|_{\Gamma_1} = 0
$$

L'équation (3.20) (ou 3.21) définie $\Psi$ de façon unique. Pour cette valeur de $\Psi$, et pour $\varphi$ solution du système de (2.7), l'accroissement du cout $\delta E$ est obtenu par

$$
\delta E = \mathcal{H}(u + \delta u, \varphi, \Psi) - \mathcal{H}(u, \varphi, \Psi)
$$

c'est-à-dire

$$
(3.22) \quad \delta E = \int_{\Omega} \left\{ 2(u - \varphi) \cdot \delta u + \Psi \left[ \frac{\gamma-1}{\gamma+1} (2u \cdot \delta u \nabla \cdot u + |u|^2 \nabla \cdot \delta u) - \frac{1}{\gamma+1} \delta u \nabla (|u|^2) \right. \right. \\
\left. \left. + \frac{2}{\gamma+1} u \nabla (u \cdot \delta u) \right] \right\} d \Omega \\
+ o(\delta u)
$$

Par intégration par partie on peut éliminer les termes contenant
des dérivées de $\delta u$ :

$$
(3.23) \quad \delta E = \int_{\Omega} \left[ 2(u - \nabla \varphi) + 2\Psi \frac{\gamma-1}{\gamma+1} (\nabla \cdot u)u + \frac{\Psi}{\gamma+1} |\nabla u|^2 - \frac{\gamma-1}{\gamma+1} \nabla (|u|^2\Psi) - \frac{2u}{\gamma+1} \nabla \cdot (\Psi u) \right] \\
\qquad \delta u \, d\Omega - o(\delta u) + \int_{\Gamma} \left( \frac{\gamma-1}{\gamma+1} \Psi |u|^2 n + \frac{2\Psi}{\gamma+1} u_n u \right) \delta u \, d\Gamma \triangleq \langle \nabla E(u), \delta u \rangle + o(\delta u)
$$

La méthode de gradient appliquée au problème (2.7) donne
l'algorithme suivant :

Algorithme (3.2)

$$
(3.24) \quad \text{Choisir } u^0 \in U \text{ poser } k = 0
$$

$$
(3.25) \quad \text{Calculer } \varphi^k \text{ solution de} \\
\qquad - \varphi^k = \frac{\gamma-1}{\gamma+1} |u^k|^2 \nabla \cdot u^k - \frac{1}{\gamma+1} u^k \cdot \nabla |u^k|^2 ; \varphi^k|_{\Gamma_1} = \varphi_1, \frac{\partial \varphi^k}{\partial n}|_{\Gamma_2} = u_2
$$
$$
(3.26) \quad \text{Calculer } \Psi_i^k \text{ solution de}
\qquad -\Delta \Psi^k = - 2 \nabla \cdot (u^k \nabla \phi^k), \quad \left. \frac{\partial \phi^k}{\partial n} \right|_{\Gamma_2} = 2(u^k \cdot n)|_{\Gamma_2} - u_2 \qquad \Psi^k|_{\Gamma_1} = 0
$$

$$
(3.27) \quad \text{Calculer } \nabla E(u^{\mathbf{k}}) \text{ par (3.22) ou (3.23)}
$$

$$
(3.28) \quad \text{Calculer } h^k \text{ solution de}
$$
$$
\min \left\{ \langle \nabla E(u^k), h \rangle \;\middle|\; \int_{\Omega} |h|^2 \, d\Omega \le 1 \right\}
$$

$$
(3.29) \quad \text{Calculer } \lambda^k \text{ solution de}
$$

$$
\min E(u^k + \lambda h^k) - E(u^k)
$$

poser $u^k = u^k + \lambda^k h^k$, $k = k+1$ et retourner en (3.25)

Remarque (3.4)

Pour conserver la propriété : $\nabla \phi^k \in (L^\infty(\Omega))^2$, il faut
ajouter dans (3.28) $|h(\cdot)|_{L^\infty} < K'$

On obtient alors des propriétés de convergence tout à fait
analogues à la proposition (3.1).

Remarque (3.5)

A titre indicatif nous donnons les formules correspondantes
au problème (2.8)

$$
(3.30) \quad \delta E = 2 \int_{\Omega} \left[ \left(1 - \frac{\gamma-1}{\gamma+1} |\nabla \varphi|^2\right) \left(v\left(1 - \frac{\gamma-1}{\gamma+1} |\nabla \varphi|^2\right) + \frac{1}{\gamma+1} \nabla \varphi \cdot \nabla (|\nabla \varphi|^2) + \Psi\right) \right] \delta v \, d\Omega + o(\delta v)
$$

où $\Psi$ est la solution de

$$
(3.31) \quad \int_{\Omega} \nabla \Psi \nabla \omega = \int_{\Omega} \left[ \frac{1}{\gamma+1} \nabla \omega \nabla |\nabla \varphi|^2 + \frac{2}{\gamma+1} \nabla \varphi \nabla (\varphi |\nabla \varphi|^2) - 2 \frac{\gamma-1}{\gamma+1} \nabla \varphi \nabla \omega \right] \\
\qquad \qquad \qquad \qquad \left[ \nu(1-\frac{\gamma-1}{\gamma+1} |\nabla \varphi|^2) + \frac{1}{\gamma+1} \nabla \varphi \cdot \nabla |\nabla \varphi|^2 \right] d\Omega
$$

pour tout $\omega \in H^1(\Omega)$ tel que $\omega|_{\Gamma_1} = 0$
### 3.3 Résolution du problème (2.9)

On pose, pour $v \in L^2(\Omega)$

$$ (3.32) \quad E(v) = \int_{\Omega} \left( \alpha \left( \frac{\partial \varphi}{\partial x} \right)^2 + \frac{\partial \varphi}{\partial x} - v \right)^2 d\Omega $$

où $\varphi$ est la solution de l'équation variationnelle

$$ (3.33) \quad \int_{\Omega} \nabla \varphi \nabla \omega \, d\Omega = \int_{\Omega} -\frac{\partial v}{\partial x} \omega \, d\Omega + \int_{\Gamma_2} u_2 \omega \, d\Gamma \quad \text{pour tout } \omega \in H^1(\Omega) \text{ tel} $$

$$ \text{que } \omega|_{\Gamma_1} = 0 $$

$$ \text{avec } \varphi|_{\Gamma_1} = \varphi_1 $$

Le problème (3.9) se récrit :

$$ \min \left\{ E(u) \mid u \in L^2(\Omega) \right\} $$

Evaluons, pour $\delta u \in L^2(\Omega)$ les termes du premier ordre de $\delta E = E(v+\delta v) - E(v)$.

Par définition, si $\frac{\partial \varphi}{\partial x}$ et $\frac{\partial \delta \varphi}{\partial x}$ sont dans $L^\infty(\Omega)$ :

$$ (3.34) \quad \begin{aligned} \frac{1}{2} \delta E = & \int_{\Omega} \left( \alpha \left( \frac{\partial \varphi}{\partial x} \right)^2 + \frac{\partial \varphi}{\partial x} - v \right) \left( 2\alpha \frac{\partial \varphi}{\partial x} \frac{\partial \delta \varphi}{\partial x} + \frac{\partial \delta \varphi}{\partial x} - \delta v \right) d\Omega \\ & + o\left( \delta v \right) + o \left( \frac{\partial \delta \varphi}{\partial x} \right) \\ & + o\left( \delta v \frac{\partial \delta \varphi}{\partial x} \right) \end{aligned} $$

De (3.33) on déduit que $\delta\varphi$ est la solution de :

$$ (3.35) \quad \int_{\Omega} \nabla \delta \varphi \nabla w \, d\Omega = - \int_{\Omega} \frac{\partial \delta v}{\partial x} w \, d\Omega \quad \text{pour tout } w \in H^1(\Omega) \text{ tel} $$

$$ \text{que } \omega|_{\Gamma_1} = 0 $$

$$ \text{avec } \delta\varphi|_{\Gamma_1} = 0 $$

De (3.35) avec $\omega = \delta\varphi$ et de l'inégalité de Poincaré on déduit qu'il existe $c, c'$ :

$$ \text{tel que } \left\| \delta\varphi \right\|_{H^1(\Omega)} \le c \left\| \frac{\partial \delta v}{\partial x} \right\|_{H^{-1}(\Omega)} \le c' \left\| \delta u \right\|_{L^2(\Omega)} $$
$$
\text{Donc } o\left(\frac{\partial \delta \varphi}{\partial x}\right) = o\left(\delta v \frac{\partial \delta \varphi}{\partial x}\right) = o\left(\delta v\right)
$$

Soit $\Psi$ la solution de l'équation varationnelle

$$
(3.36) \quad \int_{\Omega} \nabla \Psi \cdot \nabla \omega \, d\Omega = \int_{\Omega} \frac{\partial \omega}{\partial x} \left( 1 + 2\alpha \frac{\partial \varphi}{\partial x} \right) \left( \alpha \left( \frac{\partial \varphi}{\partial x} \right)^2 + \frac{\partial \varphi}{\partial x} - v \right) \, d\Omega = 0
$$

$$
\text{avec } \Psi|_{\Gamma_1} = 0
$$

En portant dans (3.36) $\omega = \delta\varphi$ et dans (3.35) $\omega = \Psi$ on obtient

$$
\int_{\Omega} \frac{\partial \delta \varphi}{\partial x} \left( 1 + 2\alpha \frac{\partial \varphi}{\partial x} \right) \left( \alpha \left( \frac{\partial \varphi}{\partial x} \right)^2 + \frac{\partial \varphi}{\partial x} - v \right) d\Omega = - \int_{\Omega} \frac{\partial \delta v}{\partial x} \Psi d\Omega
$$

Ce qui permet de déduire de (3.34) que

$$
(3.37) \quad \frac{1}{2} \delta E = \int_{\Omega} \left[ -\delta v \left( \alpha \left( \frac{\partial \varphi}{\partial x} \right)^2 + \frac{\partial \varphi}{\partial x} - v \right) + \Psi \frac{\partial \delta v}{\partial x} \right] d \Omega + o(\delta v)
$$

Algorithme 3.3

$$
\text{Soit } v = \left\{ v \in L^2(\Omega) \mid \frac{\partial \varphi}{\partial x}(v) \in L^\infty(\Omega) \right\}
$$

$$
(3.38) \quad \text{Choisir } v^0 \in V \text{ poser } k = 0
$$

$$
(3.39) \quad \text{Calculer } \varphi^k \text{ solution de (3.33) avec } v = v^k
$$

$$
(3.40) \quad \text{Calculer } \Psi^k \text{ solution de (3.36) avec } v = v^k, \varphi = \varphi^k
$$

$$
(3.41) \quad \text{Calculer } \delta u^k \text{ solution de}
$$

$$
\min_{v \in V} \left\{ \int_{\Omega} - \left[ \delta v \left( \alpha \left( \frac{\partial \varphi}{\partial x} \right)^k \right)^2 + \frac{\partial \varphi}{\partial x} - v^k + \Psi^k \frac{\partial \delta v}{\partial x} \right] d\Omega \mid \int_{\Omega} \delta v^2 d\Omega \le 1 \right\}
$$

$$
(3.42) \quad \text{Calculer } \lambda^k \text{ solution de}
$$

$$
\begin{gather*}
\min \left\{ E(v^k + \lambda \, \mathrm{d}v^k) - \cdot E(v^k) \right\} \\
\text{poser } v^k = v^k + \lambda^k \delta v^k \quad k = k + 1 \text{ et retourner en (3.39)}
\end{gather*}
$$
### Proposition 3.2

Tout point d'accumulation $\bar{\phi}$ d'une séquence $\{\varphi_i\}$ engendrée par l'algorithme tel que $\frac{\partial \bar{\phi}}{\partial x} \in L^\infty(\Omega)$ est un point stationnaire du problème (2.9). Si $E(v_i)$ tend vers zéro alors $\bar{\phi}$ est une solution de (1.20)-(2.2).

### Remarque 3.6

En un point stationnaire régulier $\bar{u}, \bar{\phi}$ de (2.9) on a d'après (3.37) :

$$ \alpha \left( \frac{\partial \bar{\phi}}{\partial x} \right)^2 + \frac{\partial \bar{\phi}}{\partial x} - \bar{v} - \frac{\partial \bar{\Psi}}{\partial x} = 0 $$

$$ \bar{\Psi}|_{\Gamma_2} = 0 $$

## 4. APPROXIMATION PAR ELEMENTS FINIS

### 4.1 Discrétisation et résolution du problème (2.9)

On se donne une triangulation $\mathcal{T}_h$ de $\Omega$ avec $h$ égal au plus grand côté de tous les triangles de $\mathcal{T}_h$, tel que, en supposant tous les triangles fermés :

$$ \bigcup_{T \in \mathcal{T}_h} T \subset \bar{\Omega} $$

Si $T_1, T_2 \in \mathcal{T}_h$, $T_1 \cap T_2$ est soit vide soit réduit à un côté commun soit réduit à un sommet commun.

On définit :

$$ \Omega_h = \bigcup_{T \in \mathcal{T}_h} T, \quad \Gamma_h = \partial \Omega_h $$

On choisit $\mathcal{T}_h$ pour que $\Gamma_{1h}, \Gamma_{2h}$ soient bien définies (c'est-à-dire que les points communs à $\Gamma_1$ et $\Gamma_2$ sont des noeuds de la triangulation).

On définit l'espace $V_h$, approximation de $L^2(\Omega)$ par
(4.1) $$ V_h = [ v_h | v_h \text{ constant sur } T, \text{ pour tout } T \in \mathcal{T}_h ] $$
On notera par $V_{oh}$ l'ensemble destiné à approximer
$$ \{ \omega \in H^1(\Omega) | \omega|_{\Gamma_1} = 0 \} : $$

$$
(4.2) \quad V_{oh} = \{ \omega_h \in C^0(\bar{\Omega}_h) \mid o_h \text{ linéaire sur } T, \forall T \in \mathcal{T}_h ; \omega_h|_{\Gamma_{1h}} = 0 \}
$$

Tout élément de $V_{oh}$ est complètement déterminé par les valeurs prisent sur

$$
W_h = \{ P \mid P \in \Omega_h \cup \Gamma_{2h}, P \text{ sommet de } \mathcal{T}_h \}
$$

En supposant les éléments de $W_h$ ordonnés par $i$ on notera
pour $w_h \in V_{oh}$,

$$
w_i = w_h (P_i), i = 1, \dots, N_h = \dim (V_{oh}), P_i \in \omega_h
$$

Pour chaque noeud $P_i$ de $W_h$ on définit $w_{hi} \in V_{oh}$ par

$$
w_{hi} \in V_{oh}
$$

$$
w_{hi}(P_i) = 1
$$

$$
w_{hi} (P) = 0 \quad \forall P \in W_h - \{P_i\}
$$

De façon évidente, $\{w_{h,i}\}_{P_i \in W_h}$ est une base de $V_{oh}$.

Le problème (2.9) est approché par le problème de contrôle discret ci-dessous. Pour la simplicité de l'exposé on a supposé $\varphi_1 = 0$.

$$
(4.3) \quad \min_{v_h \in V_h} \left\{ \int_{\Omega} \left[ \alpha\left(-\frac{\partial \varphi_h}{\partial x}\right)^2 + \frac{\partial \varphi_h}{\partial \bar{x}} - v_h \right]^2 d\Omega \middle| \begin{array}{@{}l@{}} \displaystyle \int_{\Omega} \nabla \varphi_h \cdot \nabla w_{h,i} \, d\Omega = \int_{\Omega} v_h \frac{\partial w_{hi}}{\partial \bar{x}} \, d\Omega \\[1em] \displaystyle + \int_{\Gamma_2} (u_2 - v_h) n_x w_{hi} \, d\Gamma \end{array} \right.
$$

$$
\text{pour tout } P_i \in W_h ; \varphi_h \in V_{oh}
$$

où $n_x$ est la première composante de la normale $n$ à $\Gamma_2$.
Il faut noter que même lorsque $h$ tend vers zéro, $\varphi_h$ n'approxime la solution de (3.33) que si $v_h$ tend vers une fonction dont on peut définir la trace sur $\Gamma_2$ (ce qui n'est pas le cas si $v$ est seulement dans $L^2(\Omega)$. Si on désire absolument l'équivalence avec (3.33), il faut prendre $v_h$ plus régulier et travailler avec des éléments du 2ème ordre.

L'équation variationnelle du problème (4.3) se résoud facilement en décomposant $\varphi_h$ sur les fonctions de base $w_{h_i}$.

$$ (4.4) \qquad \varphi_h = \sum_{i \le N_h} \alpha_i w_{h_i} $$

Les $\alpha_i$ sont alors calculés en résolvant le système linéaire

$$ (4.5) \quad \sum_{i \le N_h} \alpha_i \int_{\Omega} \nabla w_{h_i} \cdot \nabla w_{h_j} \, d\Omega = \int_{\Omega} v_h \frac{\partial w_{h_i}}{\partial x} \, d\Omega + \int_{\Gamma_2} (u_2 - v_h) n_x \, w_{h,j} \, d\Gamma \quad ; \ j = 1, \dots, N_h $$

l'équation adjointe discrète :

$$ (4.6) \quad \int_{\Omega} \nabla \psi_h \cdot \nabla w_h \, d\Omega = \int_{\Omega} \frac{\partial w_{h_i}}{\partial x} \left(1 + 2\alpha \frac{\partial \varphi_h}{\partial x}\right) \left(\alpha \left(\frac{\partial \varphi_h}{\partial x}\right)^2 + \frac{\partial \varphi_h}{\partial x} - v_h\right) \, d\Omega $$

se résoud de la même façon

on obtient donc l'algorithme suivant

### Algorithme 4.2

$$ (4.7) \quad \begin{gathered} \text{Choisir une triangulation } \mathcal{T}_h \text{ choisir } t > 0 \text{ petit.} \\ \text{Choisir } v_h^0 \in V_h, \text{ poser } k = 0 \end{gathered} $$
$$
(4.8) \quad \text{Calculer } \varphi_h^k \text{ solution de l'équation variationnelle (4.5) avec } v_h = v_h^k
$$

$$
(4.9) \quad \text{Calculer } \Psi_h^k \text{ solution de (4.6) avec } \varphi_h = \varphi_h^k, v_h = v_h^k
$$

$$
(4.10) \quad \text{Calculer } d = \alpha \left( \frac{\partial \varphi_h^k}{\partial x} \right)^2 + \frac{\partial \varphi_h^k}{\partial x} - v_h^k - \frac{\partial \Psi_h}{\partial x} + \langle \delta_{\Gamma_2} \cdot \Psi_h \rangle
$$

où $\langle \delta_{\Gamma_2} \cdot \Psi_h \rangle$ est constant par morceaux sur les triangles et vaut zéro
sur tous les triangles sauf ceux qui ont un côté sur $\Gamma_2$ où :

$$
\langle \delta_{\Gamma_2} \cdot \Psi_h \rangle |_{\Gamma} = \frac{1}{\text{aire}(T)} \int_{\Gamma_2} \cap T \Psi_h \ d\Gamma_2
$$

$$
(4.11) \quad \text{Calculer } \lambda \text{ solution de min}_\lambda E\left\{\langle v_h + \lambda d \rangle - E\langle v_h \rangle\right\}
$$

(on peut calculer $\lambda$ exactement car $E(v_h + \lambda d)$ est une fonction du 4ème degré en $\lambda$ dont on connait une racine).

$$
(4.12) \quad \text{Poser } u_h^{k+1} = v_h^k + \lambda d, \ k = k+1 \text{ retourner en (4.8) sauf si } ||d|| < \varepsilon.
$$

Remarque 4.2

On peut pour assurer la convergence de l'algorithme inclure
dans $v_h$ la contrainte $|v_h| < K$, K grand. A l'étape (4.12) il faut
alors poser $|u_h^{k+1}| = K$ sur tous les triangles où $|v_h^k + \lambda d|$ est plus
grand que K.

4.2 **_Discrétisation du problème_** (2.3)

Le problème (2.3) est approché par le problème de contrôle discret

$$
(4.13) \quad \min_{\varphi_h \in \Phi_h} \left\{ \frac{\int_{\Omega} |\nabla \varepsilon_h|^2 d\Omega}{\int_{\Omega_h} \nabla \varepsilon_h \nabla w_h \, d\Omega} \right\} = - \frac{\int_{\Omega_h} G(|\nabla \varphi_h|^2) \nabla \varphi_h \nabla w_h \, d\Omega}{\nabla w_h \in V_{oh} ; \varepsilon_h \in V_{oh}}
$$
où

$$
(4.14) \quad \Phi_h = \left\{ \varphi_h \in C^0(\bar{\Omega}_h) \mid \varphi_h \text{ linéaires par morceaux: surtout } T \in \mathcal{T}_h, \right.
$$

$$
\varphi_h |_{\Gamma_{1h}} = \varphi_1, \quad \frac{\partial \varphi_h}{\partial n} |_{\Gamma_{2h}} = u_2
$$

et ou $V_{oh}$ est défini par (4.2) avec $\Gamma_{1h} = \partial \Omega_h$.

L'équation variationnelle dans (4.13) est résolue comme précédem-
ment par inversion d'un système linéaire sur la base de $V_{oh}$. On notera
que l'évaluation de la dernière intégrale dans (4.13) ne pose pas de
problème car $\nabla \varphi_h$ est constant par morceaux.

On trouvera en annexe l'implémentation de la discrétisation
dans l'algorithme (3.1).

4.3 Discrétisation du problème (2.7)

On procède pour (2.7) comme pour (2.3) ; l'espace des contrôles
u est approximé par uh, ensemble des fonctions à valeur dans R² continues
linéaires par morceaux sur Ωh. Toutefois dans le cas ou on conserve la
contrainte ∇·u = 0 Il faut prendre certaines précautions quant au nombre
de variables indépendantes restantes qui doit rester au moins supérieur
au nombre de sommets si l'on veut que le problème discret soit une bonne
approximation du problème continu.

Le potentiel φ_h est calculé par l'équation variationnelle

$$
(4.15) \quad \int_{\Omega} \nabla \varphi_h \cdot \nabla w_{h,i} \, d\Omega = \int_{\Omega} \left( \frac{\gamma-1}{\gamma+1} |u_h|^2 \nabla \cdot u_h + \frac{1}{\gamma+1} u_h \cdot \nabla |u_h|^2 \right) w_{h,i} \, d\Omega + \int_{\Omega} u_h w_{h,i} \, d\Gamma
$$

pour tout $i \in W_h \cup \gamma_h$

où γh est l'ensemble des indices des sommets dans Ωh des triangles
des triangles de t_h qui ont un côté dans Γ2h.
Le potentiel adjoint $\Psi_h$ est aussi calculé par

$$ (4.16) \quad \int_{\Omega} \nabla \varphi_h \cdot \nabla w_{h_i} \, d\Omega = \int_{\Omega} 2(u_h \cdot \nabla \varphi_h) \cdot w_{h_i} \, d\Omega \quad \text{pour tout } i \in W_h \cup U_h $$

Enfin le sous-problème (3.28) est résolu à l'aide de (3.22) si la contrainte $rotu = 0$ est maintenu ; sinon il se résoud a une à l'aide de (3.23).

## 5. RESOLUTION NUMERIQUE DU PROBLEME (2.9)

On a testé la méthode sur 3 exemples :

1. une demi-tuyère divergente a des nombres de Mach de 0.83 et 1.27

2. un demi-profil avec bord d'attaque anguleux $M_\infty = 0.96$

3. une demi-tuyère convergente-divergente avec surface de raccordement anguleuse.

Dans les 3 cas on a pris une triangulation de 754 triangles et 420 noeuds. Les équations variationnelles dans (4.13) et (4.14) ont été résolues par la méthode de Cholesky en utilisant un sous-programme IBM-SSP. Chaque itération de l'algorithme (4.1) prend environ 0.4" sur un IBM 370/168. La place mémoire nécessaire pour cette triangulation est de 200 K-octet. Dans les 3 exemples on a choisi $v^0$ égal à la vitesse d'entrée du fluide.

### 5.1 Analyse des résultats relatifs à l'exemple 1.

Sur la fig. 1 on a représenté la géométrie de la tuyère, sa triangulation ainsi que la précision avec laquelle au bout de 30 itérations l'équation (1.20) est satisfaite. Sur les fig. 2, 3, 4, 5 on a tracé les isomachs au bout de 2,7, 14, 30 itérations, pour une entrée subsonique uniforme $M_\infty = 0.83$, le potentiel étant imposé constant à la sortie de la tuyère. On voit que 14 itérations sont plus que suffisantes. On notera que les conditions aux limites sur $\Gamma_2$ ne sont pas très bien satisfaites (sans que la précision de calcul de l'écoulement soit mauvais) Ceci est dû à l'utilisation de la forme intégrale (4.11) de l'équation (3.33).
Sur la fig.6 on a tracé la courbe de décroissance du critère $E(v_h^k)$ et du gradient $||d(v_h^k)||^2$ (d donné par (4.19)). On trouvera aussi sur cette figure les mêmes courbes relatives à l'exemple 2.

Ces résultats sont en parfait accord avec ceux obtenus par Perriaux [8]. Les temps de calculs sont aussi de même ordre.

On peut donc conclure de cet exemple que la méthode (4.1) a des performances tout à fait semblables aux méthodes de quasi-linéarisation usuelles (cf. Perriaux [8], Gelder [5]) dans le domaine d'application de celles-ci.

Pour tester la méthode sur un cas mixte supersonique-subsonique avec choc on a pris pour $\Gamma_1$ l'entrée et la sortie de la demi-tuyère. Le potentiel est donc pris constant à l'entrée et à la sortie de la tuyère. On a d'abord pris $\varphi = 0$ à l'entrée et $\varphi = 0.20$ à la sortie. On a obtenu les figures 6-9 sur lesquelles on peut constater que la méthode est très bien adaptée à ce type très complexe d'écoulement. L'équation (1.20) est satisfaite, au bout de 60 itérations avec une erreur sur chaque triangle $\epsilon_i = |\alpha(\frac{\partial\varphi}{\partial x})^2 + \frac{\partial\varphi}{\partial x} - u| \text{ inférieur à } 0.04$. On notera le comportement oscillatoire de la solution dans la zone supersonique. Ceci est très probablement dû à la propagation de l'erreur de position de la ligne de choc comme le prouve la fig.10 où le saut de potentiel a été pris égal à - 0.22. Par construction le choc ne peut en effet se trouver que sur une ligne de discrétisation il conviendrait donc de raffiner le maillage au voisinage du choc.

## 5.2 Analyse des résultats relatifs à l'exemple 2

L'objectif de l'exemple 2 est de simuler l'écoulement autour d'un profil d'aile en régime transonique. Sur la fig.11 on a représenté la géométrie du domaine, la triangulation et le degré de précision avec lequel l'équation (1.20) est satisfaite au bout de 60 itérations.
Sur les fig. 12 et 13 on a tracé les isomachs au bout de 30 et 60 itérations. On peut donner plusieurs arguments pour expliquer l'absence d'onde de choc à l'arrière du profil : l'influence de la paroi supérieure, l'avant pointu du profil, mais aussi le fait que la précision de résolution du problème n'est pas suffisante dans la partie supersonique de l'écoulement où les erreurs se propagent plus rapidement. Il faudrait mettre un poids dans le critère E pour donner plus d'importance à la partie supersonique.

### 5.3 Analyse des résultats relatifs à l'exemple 3.

L'objet de l'exemple 3 est de montrer que la méthode se comporte aussi bien sur des problèmes à géométrie compliquée en écoulement mixte. On a tracé les isomachs au bout de 60 itérations sur la fig. 15. On notera que les directions privilégiées de l'écoulement sont bien celles de la géométrie.
# 6 REMARQUES DIVERSES ET COMPARAISON DES METHODES.

## 6.1 Cas des profils portants

Lorsque l'écoulement se fait autour d'un profil non symétrique à bord de fuite, en point de rebroussement, une condition supplémentaire est nécessaire pour définir le potentiel. On utilise la condition de Kutta-Joukowsky :

$$ (6.1) \qquad \frac{\partial^2 \varphi}{\partial n^2} (F^+) = \frac{\partial^2 \varphi}{\partial n^2} (F^-) $$

ou $F^+$ et $F^-$ sont les deux points (confondus) extrêmes intrados et extrados du bord de fuite.

Dans le cas du problème (2.9) on a donc à résoudre à chaque itération

$$ (6.2) \qquad \Delta \varphi = \frac{\partial u}{\partial x} , \quad \varphi|_{\Gamma_1} = \varphi_1, \quad \frac{\partial \varphi}{\partial n}|_{\Gamma_2} = u_2 $$

avec (6.1)

Une technique standard consiste à décomposer $\varphi$ sous la forme
$$ \varphi = \varphi^1 + \lambda \varphi^2 $$
où $\varphi^1, \varphi^2$ sont solutions de systèmes différentiels convenables (cf. Perriaux [8]) et $\lambda$ est un paramètre qu'on ajuste de manière à satisfaire (6.1).

L'utilisation des algorithmes (3.2), (3.3) se fait sans autre modification. L'algorithme (3.1) doit tenir compte de (6.1) dans la définition de $\Phi$.

## 6.2 Résolution numérique de (2.3)

L'algorithme (2.1) a été programmé conformément aux indications données au § 4.2, et testé sur l'exemple 1 avec la même triangulation ; le mach d'entrée est spécifié uniforme et égal à 0.85 ; le potentiel est donné constant à la sortie de la tuyère.
Les résultats numériques sont très décevants : la méthode est très lente, il faut environ 1000 itérations pour obtenir un résultat a peu près satisfaisant. On a donné sur la fig. 16 un graphe de la décroissance du cout au long des itérations.

A postériori, on peut penser que la lenteur de la méthode provient du fait que la fonction $\varphi \to E(\varphi)$ est très plate et qu'il est donc très difficile de s'approcher de la solution optimale en $\varphi$ alors qu'on est tout de suite très près du $\varepsilon$ optimal correspondant.

La méthode se conduit donc à peu près comme une méthode de pénalité de l'équation d'état.

## 6.3 Comparaison des méthodes

Comme on a pu s'en rendre compte, la méthode (2.3) est la plus simple du point de vue théorique : c'est aussi la meilleure pour le passage d'un choc ; cependant nous la déconseillons car elle converge trop lentement.

Lorsque la géométrie du système le permet, la méthode (2.9) est la plus simple et la plus rapide des méthodes car on peut aussi accélérer la convergence en utilisant une méthode de gradient conjugué ; cependant dans le cas des profils portants il faut utiliser une méthode d'éléments finis du 2ème ordre on peut alors préférer traiter l'équation de l'écoulement potentiel incompressible complète et donc résoudre (2.7). Signalons que dans les cas sans chocs on peut préférer résoudre (2.8) à l'aide d'une méthode de gradient conjugué. Les résultats numériques pour ces deux méthodes seront présentés ultérieurement.
## 6.4 CONCLUSION

Le calcul des écoulements compressibles aux vitesses transoniques par des méthodes d'éléments finis est donc possible grâce à la formulation intermédiaire sous forme de problèmes de contrôle optimal. Ces nouvelles méthodes sont très compétitives avec les méthodes de différences finies déjà existantes et font mieux que les méthodes de quasi linéarisation en éléments finis. Il reste malgré tout un certain nombre de problèmes théoriques à résoudre, comme celui de la propagation des erreurs dans les zones supersoniques.

Il faut aussi maintenant tester les multiples possibilités de ces méthodes sur des problèmes réels en liaison avec les praticiens de la mécanique des fluides (ceux de l'industrie en particulier).
REFERENCES

[1] LANDAU, L. et F. LIFCHITZ
*Mécanique des fluides* - édition de Moscou (1953).

[2] COURANT, R. et K.O.FRIEDRICHS.
*Supersonic flow and shock waves*
interscience (1948).

[3] RICHMEYER, D. et MORTON, K.
*Difference methods for initial value problems*. Interscience (1957).

[4] JAMESON, A.
Iterative solution of transonic flows
over airfoils and wings including flows
at Mach 1 (Com. Pure Appl. math. 1974).

[5] GELDER, D.
Solution of the compressible flow equations
Int. J. Num. Meth. Eng., Vol 3, pp 35 - 43,
1971.

[6] DEVRIES, G. et D.H. NORRIE.
The application of the finite element
technique to potential flow problems.
J.Applied. Mech. transaction of ASME
798/Dec.1971.

[7] LIONS, J.L.
*Optimal control of systems governed by
partial differential equations*. Springer
(1968).

[8] PERRIAUX, J.
A.M.D. Company report January (1974).

[9] YVON, Thèse doct. Paris 1974.
## Algorithme A.1

(A.10) Choisir une triangulation $\mathcal{T}_h$ du type ci-dessus
Choisir $\varphi_h^o \in \Phi_h$ ; choisir $\beta \in ]0,1[$, $\varepsilon > 0$ petit. Poser $k = 0$.

(A.11) Calculer $\epsilon_h^k = \sum_{i<j} \beta_i w_{h_i} \beta_j$ est solution de
$$ \int_{\Omega} \beta_i w_{h_i} w_{h_j} \, d\Omega = - \int_{\Omega} G(|w_{h_i}|^2) w_{h_j} \, d\Omega \quad j = 1, \dots, N_h $$

(A.12) Calculer $d$ par (A.1) avec (on utilise (3.7))
$$ \alpha_i = \int_{\Omega} \left[ G(|\nabla \varphi_h^k|^2) \nabla \epsilon_h^k \nabla (w_{h_i} + \delta_{h_i}^{\bar{w}}) + 2G'(|\nabla \varphi_h^k|^2) (\nabla \varphi_h^k \cdot \epsilon_h^k) (\nabla \varphi_h^k, \nabla (w_{h_i} + \delta_{h_i}^{\bar{w}})) \right] d\Omega $$
$i \in W_h$ où $\delta_i = 1$ si $i \in \gamma_h$, 0 ailleurs.

(A.13) Calculer le plus petit entier positif $j$ tel que
$$ E(\varphi_h^k + \beta^j d) - E(\varphi_h^k) \le -\frac{1}{2} \beta^j ||d||^2 $$

(A.14) Si $||d||^2 < \varepsilon$ terminer sinon poser $\varphi_h^{k+1} = \varphi_h^k + \beta^j d$
$k = k + 1$ er retourner en (A.6)

## Remarque A.1

L'algorithme A.1 est maintenant une méthode de gradient sur les $\alpha_i$ et non sur $\varphi_h$.
# ANNEXE

Les problèmes (2.3) et (4.3) sont de même type puisqu'ils ne diffèrent que par l'ensemble des contrôles admissibles soit $\Phi_h$ au lieu de $\Phi_m$.

Pour calculer facilement la solution du programme linéaire dans (3.12) on choisit des fonctions de base supplémentaires continues linéaires par morceaux $\{\bar{w}_{hi}\}_{i \in \gamma_h}$ ou $\gamma_h$ est l'ensemble des indices des sommets dans $\Omega_h$ des triangles de $\mathcal{T}_h$ qui ont un côté dans $\Gamma_{2h}$

On peut choisir $\mathcal{T}_h$ de telle sorte que toute combinaison linéaire de la forme (†)

$$ (A.6) \qquad d = \sum_{i \in \bar{W}_h} \alpha_i w_{hi} + \sum_{i \in \gamma_h} \alpha_i \bar{w}_{hi} $$

satisfasse

$$ (A.7) \qquad \frac{\partial d}{\partial n} = 0 $$

La solution discrète du problème (3.12) est alors, a un facteur multiplicatif près, donnée par (A.6) avec

$$ (A.8) \qquad \alpha_i = - \langle \nabla E(\varphi^k), w_{hi} + \bar{w}_{hi} \rangle \quad i \in \gamma_h $$

$$ (A.9) \qquad \alpha_i = - \langle \nabla E(\varphi^k), w_{hi} \rangle \quad i \in W_h - \gamma_h $$

On obtient donc l'adaptation de l'algorithme (3.1) suivante

<footnote>(†) une méthode simple consiste à prendre une triangulation dont les triangles de la frontière $\Gamma_{2h}$ sont rectangles et à prendre $w_{hi} = 1$ au sommet rectangle et zéro aux autres sommets.</footnote>
<figure><img src="image_5.png" /><figcaption>Figure 1. Example 1: triangulation. $\epsilon_i = \left\langle \frac{\partial^k \phi}{\partial x^k} \right\rangle 2 \cdot \frac{\partial^k \psi}{\partial x^k} - v^k, k=1,2,\dots,k=14$.</figcaption></figure>

LABORIA- O.PIRONNEAU

.842

.198

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

NB D ITERATIONS 2

$E = 2 \times 10^{-4}$

$|d|^2 = 2 \times 10^{-4}$

<caption>Figure 2. Isomachs.</caption>

.189

LABORIA- O.PIRONNEAU
<figure><img src="image_7.png" /><figcaption>Figure 4. Isomachs.</figcaption></figure>

LABORIA- O.PIRONNEAU

<figure><img src="image_8.png" /><figcaption>Figure 3. Isomachs.</figcaption></figure>

LABORIA- O.PIRONNEAU
<figure><img src="image_9.png" /><figcaption>Figure 8. Isomachs.</figcaption></figure>

LABORIA- O.PIRONNEAU

<figure><img src="image_10.png" /><figcaption>Figure 7. Isomachs the shock is not visible because it was smoothed out by the isomach drawing subroutine.</figcaption></figure>

LABORIA- O.PIRONNEAU
<figure><img src="image_11.png" /><figcaption>Figure 11. Triangulation of example 2 and precision after 60 iterations.</figcaption></figure>

LABORIA- O.PIRONNEAU

<figure><img src="image_12.png" /><figcaption>Figure 8. Isomachs.</figcaption></figure>

LABORIA- O.PIRONNEAU
<figure><img src="image_13.png" /><figcaption>NB D ITERATIONS 60

E = .9x10⁻⁴
|d|² = 10⁻⁶</figcaption></figure>

Figure 13. Isomachs.

LABORIA- O.PIRONNEAU

<figure><img src="image_14.png" /><figcaption>NB D ITERATIONS 30

E = 10⁻⁴
|d|² = 2x10⁻⁶</figcaption></figure>

Figure 12. Isomachs.

LABORIA- O.PIRONNEAU
<figure><img src="image_15.png" /><figcaption>Figure 15. Isomachs.</figcaption></figure>

LABORIA- O.PIRONNEAU

<figure><img src="image_16.png" /><figcaption>Figure 14. example 3.</figcaption></figure>

LABORIA- O.PIRONNEAU
<figure><img src="image_17.png" /><figcaption>Figure 10. Same as for fig. 9 with a different $\varphi_2 - \varphi_1$.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Figure 9. Variation of the Mach along 3 lines of discretisation. The horizontal segments are the step function, given by the computer, from which the continuous line was extrapolated.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Figure 6. E and |df/dk| are drawn as functions of the index of iterations k for example 1 and 2.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Figure 10. Cost as a function of k for Problem (2.3)</figcaption></figure>