# Expériences numériques sur des problèmes non-linéaires résolus par éléments finis et lagrangien augmente

Americo Marrocco

► To cite this version:

Americo Marrocco. Expériences numériques sur des problèmes non-linéaires résolus par éléments finis et lagrangien augmente. [Rapport de recherche] IRIA-RR-309, IRIA. 1978, pp.104. hal-04716521

HAL Id: hal-04716521

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

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
Dw: 04185 RR 309

RIA

laboria

Institut de Recherche
d'Informatique
et d'Automatique

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

laboratoire de recherche
en informatique
et automatique

EXPÉRIENCES NUMÉRIQUES
SUR DES PROBLÈMES
NON-LINÉAIRES RESOLUS
PAR ÉLÉMENTS FINIS
ET LAGRANGIEN AUGMENTÉ

Americo MARROCCO

Rapport de Recherche N° 309

Mai 1978

1040
# EXPÉRIENCES NUMÉRIQUES SUR DES PROBLEMES NON-LINÉAIRES RÉSOLUS PAR ÉLÉMENTS FINIS ET LAGRANGIEN AUGMENTÉ

## Americo Marrocco

### Résumé :

On présente dans ce rapport les résultats numériques obtenus lors de la résolution de différents problèmes aux limites non linéaires par la méthode du lagrangien augmenté, l'approximation étant faite par éléments finis $P_1$. On trouvera un support théorique de la méthode dans M. FORTIN-R. GLOWINSKI [1]. Les problèmes étudiés sont les suivants :

- Problème du pseudo-laplacien

- Problème de magnétostatique bidimensionnelle

- Problème de l'écoulement d'un fluide de Bingham dans une conduite cylindrique

- Problème de la torsion élasto-plastique

- Problème de l'écoulement dans une tuyère

- Problème des surfaces minima

Deux algorithmes principaux sont mis en œuvre numériquement, et des variantes sont envisagées. Cette technique du Lagrangien augmenté permet de résoudre une variété très grande de problèmes non linéaires, avec le même logiciel de base constitué essentiellement par un programme de résolution d'un système linéaire à matrice symétrique définie positive. La non linéarité spécifique à chaque problème étudié se ramène à l'étude d'une équation (ou inéquation) non linéaire à 1 inconnue explicitée pour les problèmes présentés ici, et à des systèmes non linéaires à 2 ou 3 équations (ou inéquations), selon la dimension de $\Omega$, pour certains autres problèmes. Cette technique du lagrangien augmenté se généralise évidemment au cas tridimensionnel et des essais concluant ont déjà été faits pour le problème 2.

### Abstract :

Some non linear partial differential equations are solved numerically by finite elements and augmented lagrangian technique. Theoretical support of this can be found in M. FORTIN-R. GLOWINSKI [1]. In this paper numerical results are illustrated by many figures. The different problems we have studied, are the following :

- A class of non linear Dirichlet problems

- The two-dimensional magnetostatic problem

- Bingham flow in a cylindrical pipe

- The elasto-plastic torsion problem

- A compressible flow through a nozzle

- Minimal surface problem

Two algorithms are tested numerically. With this augmented lagrangian technique, different non linear problems can be solved with the same package which consists principally of a routine for solving a symmetric positive definite linear system. The non linearity, which is specific to the problem, leads to a single non linear equation (or inequation) for each finite element and for the problems presented here. It can occur non linear systems with 2 or 3 equations (or inequations) for other type of problems. This augmented lagrangian technique can be easily generalised for the three-dimensional case and some results are already obtained for the problem number 2.
I. - INTRODUCTION

Nous appliquons la méthode du Lagrangien augmenté pour résoudre différents problèmes non-linéaires tels que le pseudo-laplacien [1]. (III-1.22), les équations de Maxwell de la magnétostatique, l'écoulement d'un fluide visco-plastique dans une conduite cylindrique [1] (III. Exemple 1), le problème de la torsion élasto plastique d'un cylindre de Section Ω [1] (III. Exemple 2), le problème des surfaces minima, l'écoulement subsonique dans une tuyère.

On supposera dans la suite que le domaine Ω dans lequel on doit résoudre le problème considéré est un polygone de R².

Soit $\mathcal{T}_h$ une triangulation finie de Ω vérifiant

$$
(1.1) \quad \left\{
\begin{array}{ll}
T \subset \bar{\Omega} \quad \forall T \in \mathcal{T}_h \\
\bigcup_{T \in \mathcal{T}_h} T = \bar{\Omega}
\end{array}
\right.
$$

$$
(1.2) \quad \left\{
\begin{array}{@{}l@{\quad}l@{}}
\multicolumn{2}{@{}l@{}}{T, T' \in \mathcal{T}_h \text{ alors}} \\
\text{ou } T \cap T' = \emptyset \\
\text{ou } T \text{ et } T' \text{ ont un sommet commun} \\
\text{ou } T \text{ et } T' \text{ ont un côté commun}
\end{array}
\right.
$$

Nous reprenons maintenant essentiellement les notations du Chapitre III de [1].
L'espace V sera approché par un espace du type

$$
(1.3) \quad \mathcal{V}_h = \{ v_h | v_h \in \mathcal{C}^0(\bar{\Omega}) , v_h|_T \in P_1 \quad \forall T \in \mathcal{T}_h \}
$$

avec $P_l = \text{espace des polynômes de degré } \le l$.

L'espace $V_h$ que l'on retiendra pour chaque application sera un sous-espace
de $\mathcal{V}_h$ caractérisé par certaines conditions imposées à $v_h \in \mathcal{V}_h$, sur le bord de $\Omega$.

L'espace H sera approché par

$$
(1.4) \quad L_h = \{ z_h | z_h = \sum_{T \in \mathcal{T}_h} z_T \chi_T , z_T \in \mathbb{R}^2 \}
$$

$\chi_T$ étant la fonction caractéristique de $T \in \mathcal{T}_h$.
Vu l'approximation choisie pour V et H par $V_h$ et $L_h$ respectivement,
l'application B est aussi une application linéaire de $V_h$ dans $L_h$ (pour les
problèmes étudiés)

$$
(1.5) \quad B \in \mathcal{L}(V_h, L_h) .
$$

Rappelons que si le problème initial est

$$
(1.6) \quad \text{Min} \quad \{F(Bv_h) + G(v_h)\} \\
v_h \in V_h
$$

on peut se ramener à la recherche d'un point-selle du Lagrangien augmenté

$$
(1.7) \quad \mathcal{L}_r(v_h, q_h; \mu_h) = F(q_h) + G(v_h) + (\mu_h, Bv_h - q_h) + \frac{r}{2} |Bv_h - q_h|^2
$$

avec $\mathcal{L}_r : V_h \times L_h \times L_h \to \mathbb{R}$.

Nous pouvons aussi rappeler brièvement les algorithmes utilisés pour la
résolution numérique.

/ ALG 1 /

$$
(1.8) \quad \lambda_h^0 \in L_h, \text{ donné}
$$

$\lambda_h^n$ étant connu, on détermine $\{u_h^n, p_h^n\}$ puis $\lambda_h^{n+1}$ par

$$
(1.9) \quad \left\{
\begin{array}{l}
\mathcal{L}_r(u_h^n, p_h^n, \lambda_h^n) \le \mathcal{L}_r(v_h, q_h, \lambda_h^n) \quad \forall \{v_h, q_h\} \in V_h \times L_h \\
\\
\{u_h^n, p_h^n\} \in V_h \times L_h
\end{array}
\right.
$$

$$
(1.10) \quad \lambda_h^{n+1} = \lambda_h^n + \rho_n (B u_h^n - p_h^n) .
$$

On déduit immédiatement de (1.9) que le couple $\{u_h^n, p_h^n\}$ est caractérisé
par le système couplé
$$
(1.11) \quad \left\{
\begin{array}{l}
G(v_h) - G(u_h^n) + (\lambda_h^n, B(v_h - u_h^n)) + r(Bu_h^n - p_h^n, B(v_h - u_h^n)) \ge 0 \\
\forall v_h \in V_h, u_h^n \in V_h
\end{array}
\right.
$$

$$
(1.12) \quad \left\{
\begin{array}{l}
F(q_h) - F(p_h^n) - (\lambda_h^n, q_h - p_h^n) + r(p_h^n - Bu_h^n, q_h - p_h^n) \ge 0 \\
\forall q_h \in L_h, \quad p_h^n \in L_h.
\end{array}
\right.
$$

La structure même du système (1.11)-(1.12) conduit à l'algorithme suivant
noté ALG2.

ALG 2

$$
(1.13) \quad \{p_h^0, \lambda_h^1\} \in L_h \times L_h \text{ donné}
$$

`{p_h^{n-1}, λ_h^n} étant connu, on détermine successivement u_h^n, p_h^n, λ_h^{n+1} par`

$$
(1.14) \quad \left\{
\begin{array}{l}
G(v_h) - G(u_h^n) + (\lambda_h^n, B(v_h - u_h^n)) + r(Bu_h^n - p_h^{n-1}, B(v_h - u_h^n)) \ge 0 \\
\forall v_h \in V_h, \quad u_h^n \in V_h,
\end{array}
\right.
$$

$$
(1.15) \quad \left\{
\begin{array}{l}
F(q_h) - F(p_h^n) - (\lambda_h^n, q_h - p_h^n) + r(p_h^n - Bu_h^n, q_h - p_h^n) \ge 0 \\
\forall q_h \in L_h, \quad p_h^n \in L_h
\end{array}
\right.
$$

$$
(1.16) \quad \lambda_h^{n+1} = \lambda_h^n + \rho_n (Bu_h^n - p_h^n) .
$$

La différence avec (1.11)-(1.12) c'est qu'il faut se donner un $p_h^o$ initial,
et "remplacer" dans (1.11) $p_h^n$ par $p_h^{n-1}$, mais cette façon de faire va con-
duire a un algorithme qui sera plus efficace du point de vue temps calcul.
Remarquons que l'algorithme ALG 2 peut s'interpréter aussi par rapport à
ALG 1.

Pour résoudre le système couplé (1.11)-(1.12) on utilisé une méthode de
relaxation par blocs (le 1er bloc étant constitué par (1.11) et le deuxième
bloc par (1.12)), il faut faire un certain nombre d'itérations (par blocs)
pour résoudre le système couple (1.11)-(1.12), dans l'algorithme ALG 2, on
ne fait qu'une itération et on passe à l'étape suivante.
Suivant la méthode de résolution utilisée pour résoudre (1.14) qui dans nos applications est un **problème linéaire**, on peut envisager des variantes ( et éventuellement des extrapolations) des algorithmes ALG 1 et ALG 2.

Remarquons aussi que le problème (1.15) se ramène facilement à NT (si NT est le nombre de triangles de $\mathcal{T}_h$) équations (ou inéquations) non linéaires, équations (ou inéquations) à 1 variable.

## ALG 1 - 1er Schéma

Méthode de Gauss-Seidel sur le système (1.11)-(1.12).

$$ (1.17) \qquad \left\{ \begin{array}{l} \text{On choisit } \{p_h^o, \lambda_h^o\} = \{0,0\} \in L_h \times L_h \\[1em] \phantom{\text{On choisit }} u_h^o \in V_h. \end{array} \right. $$

On suppose connus $\lambda_h^n, p_h^n, u_h^n$ et

a) $\lambda_h^n$ étant donné on résout le système non linéaire (1.11)-(1.12) par la méthode itérative suivante

$$ (1.18) \qquad \begin{array}{l} \text{On pose } v^o = u_h^n \\ z^o = p_h^n \end{array} $$

pour $k=0,1,2,3,...$

On fait un balayage de type relaxation (on ne résout pas) sur

$$ (1.19) \qquad \left\{ \begin{array}{l} G(v_h) - G(v^{k+1}) + (\lambda_h^n, B(v_h - v^{k+1})) + r(Bv^{k+1} - z^k, B(v_h - v^{k+1})) \ge 0 \\ \forall v_h \in V_h, \quad v^{k+1} \in V_h \end{array} \right. $$

ensuite on résout les NT equations (ou inéquations) non linéaires

$$ (1.20) \qquad \left\{ \begin{array}{l} F(q_h) - F(z^{k+1}) - (\lambda_h^n, q_h - z^{k+1}) + r(z^{k+1} - Bv^{k+1}, q_h - z^{k+1}) \ge 0 \\ \forall q_h \in L_h, \quad z^{k+1} \in L_h. \end{array} \right. $$
On arrête ce processus en k lorsque par exemple

$$
(1.21) \quad \text{Rez} = \frac{\sum_i |v_i^{k+1} - v_i^k|}{\sum_i |v_i^{k+1}|} \le \epsilon
$$

Lorsque (1.21) est vérifié avec $\epsilon=10^{-5}$, on pose

$$
(1.22) \quad \left\{
\begin{array}{l}
u_h^{n+1} = v^{k+1} \\
p_h^{n+1} = z^{k+1}
\end{array}
\right.
$$

b) On fait ensuite

$$
(1.23) \quad \lambda_h^{n+1} = \lambda_h^n + \rho (B u_h^{n+1} - p_h^{n+1}) .
$$

On retourne à a) en incrémentant n si le test (1.24) n'est pas vérifié

$$
(1.24) \quad \text{Res} = \frac{\sum_i |u_{i,h}^{n+1} - u_{i,h}^n|}{\sum_i |u_{i,h}^{n+1}|} \le \epsilon
$$

(En général $=10^{-5}$ dans (1.24)).

ALG 1 - 2ème Schéma

C'est la méthode de résolution par relaxation par blocs, (conduisant comme
cas particulier limite à l'algorithme ALG 2).

$$
(1.25) \quad \text{On choisit } \{p_h^o, \lambda_h^o\} = \{0,0\} \in L_h \times L_h
$$
$$
u_h^o \in V_h .
$$

On suppose connus $\lambda_h^n, p_h^n, u_h^n$ et

a) $\lambda_h^n$ étant donné on résout le système non-linéaire (1.11)-(1.12) par la méthode itérative par blocs suivante
$$
\begin{array}{ll}
\text{On pose} & v^o = u_h^n \\
               & z^o = p_h^n
\end{array}
$$

1°) On résout le problème

$$
(1.26) \quad \left\{
\begin{array}{l}
G(v_h) - G(v^{k+1}) + (\lambda_h^n, B(v_h - v^{k+1})) + r(Bv^{k+1} - z^k, B(v_h - v^{k+1})) \ge 0 \\
\forall v_h \in V_h, \quad v^{k+1} \in V_h
\end{array}
\right.
$$

qui est en général (dans nos applications), un **problème linéaire** que l'on
peut résoudre par une méthode itérative (du type surrelaxation) **ou mieux**
**par une méthode directe**. La matrice du problème linéaire étant toujours
la même pour chaque n, il suffit donc de faire la factorisation une fois
pour toute, la résolution étant ensuite très simple à chaque étape.

2°) On résout ensuite les NT équations (ou inéquations) non-linéaires

$$
(1.27) \quad \left\{
\begin{array}{l}
F(q_h) - F(z^{k+1}) - (\lambda_h^n, q_h - z^{k+1}) + r(z^{k+1} - Bv^{k+1}, q_h - z^{k+1}) \ge 0 \\
\forall q_h \in L_h, \quad z^{k+1} \in L_h .
\end{array}
\right.
$$

On arrête le processus en k lorsque

$$
(1.28) \quad \text{Rez} = \frac{\sum_i |v_i^{k+1} - v_i^k|}{\sum_i |v_i^{k+1}|} \le \epsilon
$$

Si (1.28) n'est pas vérifié (avec $\epsilon=10^{-5}$) on retourne à 1°) sinon on pose

$$
u_h^{n+1} = v^{k+1}
$$

$$
p_h^{n+1} = z^{k+1}
$$

b) On fait ensuite

$$
(1.29) \quad \lambda_h^{n+1} = \lambda_h^n + \rho (B u_h^{n+1} - p_h^{n+1})
$$

et on retourne à a) en incrémentant n si (1.30) n'est pas vérifié.
$$
(1.30) \quad \operatorname{Res} = \frac{\sum_i |u_{i,h}^{n+1} - u_{i,h}^n|}{\sum_i |u_{i,h}^{n+1}|} \le \epsilon
$$

On tire facilement de ALG 1 - 2ème Schéma, un schéma de résolution ALG 2 ;
Explicitons

ALG 2 - Schéma 1

$$
\begin{array}{ll}
(1.31) & \text{On choisit } \{p_h^0, \lambda_h^0\} = \{0,0\} \in L_h \times L_h \\
& u_h^0 \in V_h
\end{array}
$$

On suppose connus $\lambda_h^n, p_h^n, u_h^n$ et

1°) On calcule $u_h^{n+1}$ solution de

$$
(1.32) \quad \left\{
\begin{array}{l}
G(v_h) - G(u_h^{n+1}) + (\lambda_h^n, B(v_h - u_h^{n+1})) + r(Bu_h^{n+1} - p_h^n, B(v_h - u_h^{n+1})) \ge 0 \\
\forall v_h \in V_h, \quad u_h^{n+1} \in V_h
\end{array}
\right.
$$

(1.32) est un problème linéaire dans nos applications, résolu par une
méthode itérative ou mieux par une méthode directe.

2°) On calcule $p_h^{n+1}$ en résolvant les NT équations (ou inéquations) non linéaires

$$
(1.33) \quad \left\{
\begin{array}{l}
F(q_h) - F(p_h^{n+1}) - (\lambda_h^n, q_h - p_h^{n+1}) + r(p_h^{n+1} - Bu_h^{n+1}, q_h - p_h^{n+1}) \ge 0 \\
\forall q_h \in L_h, \quad p_h^{n+1} \in L_h
\end{array}
\right.
$$

3°) On calcule ensuite $\lambda_h^{n+1}$ par

$$
(1.34) \quad \lambda_h^{n+1} = \lambda_h^n + \rho (Bu_h^{n+1} - p_h^{n+1})
$$

on retourne à 1°) en incrémentant n si (1.35) n'est pas vérifiée.
$$ (1.35) \quad \mathrm{Res} = \frac{\sum_i |u_{i,h}^{n+1} - u_{i,h}^n|}{\sum_i |u_{i,h}^{n+1}|} \le \epsilon $$

Si on utilise une méthode itérative pour résoudre le problème (1.32), on peut en combinant ALG 1 - 1er Schéma et ALG 2 - 1er Schéma, aboutir au schéma suivant

## ALG 2 - 2ème Schéma

$$ (1.36) \quad \begin{gathered} \text{On choisit } \{p_h^o, \lambda_h^o\} = \{0,0\} \in L_h \times L_h \\ u_h^o \in V_h. \end{gathered} $$

On suppose connu $\lambda_h^n, p_h^n, u_h^n$ et

1°) On calcule $u_h^{n+1}$ par un seul balayage de type relaxation sur le système

$$ (1.37) \quad \left\{ \begin{array}{l} G(v_h) - G(u_h^{n+1}) + (\lambda_h^n, B(v_h - u_h^{n+1})) + r(Bu_h^{n+1} - p_h^n, B(v_h - u_h^{n+1})) \ge 0 \\ \forall v_h \in V_h, \ u_h^{n+1} \in V_h \end{array} \right. $$

On poursuit ensuite comme dans (ALG 2 - 1er Schéma), le test d'arrêt étant toujours (1.35).

Nous allons maintenant donner des indications générales sur le comportement des divers schémas exposés ci-dessus. Nous parlerons tout d'abord de [ALG 2 - 2ème Schéma] - (1.36) - Remarquons tout d'abord que cet algorithme n'a pas de justification théorique et on ne sait pas a priori s'il va converger ou non. Nous avons tenté plusieurs essais sur le problème du pseudo-laplacien (Exemple 1) dans un cas très compliqué il est vrai car nous avons choisi $p=1.1$. Nous avons pris $\rho=r$ (paramètre de pénalisation) dans (1.34) ; signalons tout de suite que cette valeur de $\rho$ a été un très bon choix pour les autres algorithmes. Pour les essais numériques nous avons toujours pris une valeur fixe de $\rho : \rho_n = \rho$.
Nous ne sommes pas arrivés à faire converger cet algorithme ; on constate des oscillations de $R\epsilon(1.35)$ et parfois même une nette divergence suivant la valeur du paramètre $r$. Peut-être le choix de $\rho=r$ est-il trop fort ici, mais précisons dès maintenant que cet algorithme ne présente que peu d'intérêt, surtout si on possède une méthode directe pour résoudre les systèmes tel que (1.32).

Revenons maintenant à l'algorithme ALG 1 - 1er Schéma (1.17)-(1.24). Cet algorithme est un de ceux qui est utilisé dans [2] ; tel quel, cet algorithme se bloque, "converge" assez rapidement au bout de 4 ou 5 itérations en $n$ (ou itérations sur $\lambda$). Le phénomène qui se produit est le suivant : dès la première itération en $k$, le test (1.21) est vérifié et par là même le test (1.24) d'où blocage et arrêt des itérations, mais nous n'avons cependant pas atteint la solution. On peut remédier à cet inconvénient soit en éliminant le test (1.21) à la lère itération en $k$ soit en faisant porter le test sur le couple $(v^k, z^k)$.

Notons cependant que dans le cas du pseudo-laplacien (Exemple 1), cet algorithme, tel qu'il est présenté en (1.17)-(1.24) est valable car la solution obtenue après "blocage" n'est pas plus mauvaise (et parfois meilleure) que la solution réelle du problème discret. (Sur cet exemple nous pouvions comparer à la solution analytique du problème voir [2]).

Pour les équations de la magnétostatique (Exemple 2), il se produit le même phénomène de blocage mais la solution obtenue ainsi est assez loin de la solution réelle du problème discret. Est-ce dû à la "forme de la non linéarité" ? Donc cet algorithme (1.17)-(1.24) est dangereux si on ne remédie pas à la "fausse convergence". Si nous prenons maintenant ALG 1 - 2ème Schéma - (1.25)-(1.30), on constate évidemment le même phénomène de blocage que dans l'algorithme précédent.

De façon générale pour l'algorithme ALG 1 (1er ou 2ème schéma) on a constaté que le nombre d'itérations en $n$ diminuait lorsque le paramètre de pénalisation $r$ augmentait, par contre le nombre d'itérations en $k$, lui augmentait avec $r$. On a aussi constaté qu'il valait mieux prendre $r$
plutôt grand et limiter les itérations en k (sans nécessairement atteindre la convergence au sens (1.21) ou (1.28)). En extrapolant cette idée, on limite à 1 les itérations en k et alors l'algorithme ALG 1 - 2ème schéma devient ALG 2 - 1er schéma.

Cet algorithme ALG 2 - 1er Schéma est l'autre algorithme utilisé dans [2] pour le problème du pseudo-laplacien. Il n'a pas conduit à des surprises comme l'algorithme ALG 1 (1er et 2ème schémas) et lorsque la convergence était obtenue au sens du test (1.35), on avait effectivement la solution du problème de point-selle (et évidemment la solution u du problème d'optimisation initial). Les essais numériques présentés dans [2] étaient réalisés avec la méthode itérative de surrelaxation pour la résolution de (1.32) mais cet algorithme ALG 2 1er Schéma devient évidemment plus efficace si on utilise une méthode directe pour la résolution du système linéaire (1.32).

Les algorithmes utilisés pour les applications numériques sont surtout les algorithmes ALG 2 1er Schéma et ALG 1 2ème schéma (évitant ou n'évitant pas le blocage intempestif pour ce dernier). La résolution du problème (1.32) ou (1.26) qui est pour nous un système linéaire à matrice symétrique définie positive, se fait par la méthode de Cholesky.

Ce qui sera typique du problème étudié sera l'équation (ou l'inéquation) du type (1.12). On détaillera donc ce problème pour chaque application.
## II - EXEMPLE 1 - PROBLEME DU PSEUDO-LAPLACIEN

### II.1. Présentation du problème

Rappelons ici le problème déjà étudié dans [2] et signalé au Chap. 3 (6.2)- de [1].

Soit $\Omega$ un ouvert borné de $\mathbb{R}^N$, de frontière $\Gamma$ régulière et $p \in ]1, +\infty[$.
On note $V$ l'espace $W_0^1, p(\Omega)$, Banach reflexif pour la norme

$$ \|v\|_1 = \left( \int_{\Omega} |\nabla v|^p dx \right)^{\frac{1}{p}} $$

Soit $p'$ tel que $\frac{1}{p'} + \frac{1}{p} = 1$ ; on note $V'$ le dual $W^{-1}, p'(\Omega)$ de $V$ et $\|\cdot\|_1$ la norme duale de $\|\cdot\|_1$.

Soit $f \in V$, on considère le problème aux limites non linéaire

$$ (2.1) \qquad \left\{ \begin{array}{l} - \nabla (|\nabla u|^{p-2} \nabla u) = f \text{ dans } \Omega \\ u|_{\Gamma} = 0 . \end{array} \right. $$

Si on note $A$ l'opérateur de $V \to V'$ défini par

$$ (2.2) \qquad A(v) = - \nabla (|\nabla v|^{p-2} \nabla v) , $$

le gradient $J'$, de $J$, ou $J$ est l'application de $V \to \mathbb{R}$ définie par

$$ (2.3) \qquad J(v) = \frac{1}{p} \int_{\Omega} |\nabla v|^p dx - \langle f, v \rangle $$

s'écrit

$$ (2.4) \qquad J'(v) = A(v) - f $$

Dans (2.3) $\langle \cdot, \cdot \rangle$ dénote la forme bilinéaire de la dualité entre $V$ et $V'$.

Il y a donc équivalence à résoudre le problème (2.1) dans $V$ et résoudre le problème de minimisation.

$$ (2.5) \qquad \left\{ \begin{array}{l} J(u) \le J(v) \quad \forall v \in V \\ u \in V. \end{array} \right. $$
Comme J est strictement convexe, différentiable, avec $\lim_{||v||_1 \to \infty} J(v) \to +\infty$, il en résulte (cf. LIONS [3]) que (2.5) admet une solution et une seule caractérisée par

$$ (2.6) \quad \left\{ \begin{array}{l} \displaystyle \int_{\Omega} |\nabla u|^{p-2} \nabla u \cdot \nabla v \, dx = \langle f, v \rangle \quad \forall v \in V \\[1em] \qquad u \in V \end{array} \right. $$

La relation avec la formulation générale, c'est-à-dire la définition des applications B, F, G de (1.6) est donnée au Chap. 3 (6.2) de [1].

## II.2. Formulation discrète en Lagrangien augmenté.

Pour le problème du pseudo-laplacien, le lagrangien augmenté s'écrit :

$$ (2.7) \quad \mathcal{L}_r(v_h, q_h, \mu_h) = \frac{1}{p} \int_{\Omega} |q_h|^p dx - \langle f, v_h \rangle + \frac{r}{2} \int_{\Omega} |\nabla v_h - q_h|^2 dx + \int_{\Omega} \mu_h \cdot (\nabla v_h - q_h) dx $$
avec $\mathcal{L}_r : V_h \times L_h \times L_h \to \mathbb{R}$.

Explicitons maintenant le système (1.11)-(1.12), que l'on peut écrire de la façon suivante :

$$ (2.8) \quad \left\{ \begin{array}{l} \displaystyle \int_{\Omega} \nabla \hat{u}_h \cdot \nabla v_h \, dx = \frac{1}{r} \langle f, v_h \rangle + \int_{\Omega} (\hat{p}_h - \frac{1}{r} \mu_h) \cdot \nabla v_h \, dx \\[1em] \qquad \forall v_h \in V_h, \quad \hat{u}_h \in V_h \end{array} \right. $$

$$ (2.9) \quad \left\{ \begin{array}{l} \displaystyle \int_{\Omega} (|\hat{p}_h|^{p-2} \hat{p}_h + r \hat{p}_h) q_h \, dx = \int_{\Omega} (r \nabla \hat{u}_h + \mu_h) \cdot q_h \, dx \\[1em] \qquad \forall q_h \in L_h, \quad \hat{p}_h \in L_h \end{array} \right. $$

Si NT est le nombre de triangles de la triangulation $\mathcal{T}_h$, le système non linéaire (2.9) se décompose en NT systèmes, indépendants les uns des autres, de 2 équations à 2 inconnues admettant chacun une solution unique dans $\mathbb{R}^2$. De façon plus précise on détermine $\hat{p}_h$ en résolvant pour chaque triangle $T \in \mathcal{T}_h$ le système non linéaire suivant :
$$
(2.10) \quad \begin{cases} |\hat{p}_h|^{p-2} \hat{p}_{h,1} + r \hat{p}_{h,1} = r \frac{\partial \hat{u}_h}{\partial x_1} + \mu_{h,1} \\ |\hat{p}_h|^{p-2} \hat{p}_{h,2} + r \hat{p}_{h,2} = r \frac{\partial \hat{u}_h}{\partial x_2} + \mu_{h,2} \end{cases}
$$

L'indice 1 ou 2 indique le numéro de la composante ; il n'y a pas d'ambiguité dans l'écriture de (2.10), toutes les quantités intervenant étant constantes par triangle. Le système (2.10) peut se ramener à une seule équation à une seule inconnue $|\hat{p}_h|$, les composantes $\hat{p}_{h,1}$ et $\hat{p}_{h,2}$ étant alors facilement calculables à partir de (2.10).

On tire $\hat{p}_{h,1}^2$ de la première équation 2.10 et $\hat{p}_{h,2}^2$ de la deuxième. On écrit ensuite que

$$
|\hat{p}_h|^2 = p_{h,1}^2 + p_{h,2}^2
$$

soit

$$
(2.11) \quad |\hat{p}_h|^2 = \frac{\left( r \frac{\partial \hat{u}_h}{\partial x_1} + \mu_{h,1} \right)^2 + \left( r \frac{\partial \hat{u}_h}{\partial x_2} + \mu_{h,2} \right)^2}{(|\hat{p}_h|^{p-2} + r)^2}
$$

Si on pose $q = |\hat{p}_h|$ on tire facilement de (2.11) que $q$ doit vérifier
l'équation

$$
(2.12) \quad g(q) = q^{p-1} + rq = |r \nabla \hat{u}_h + \mu_h|
$$

L'équation (2.12) est du type $g(q) = Cte \ q \in \mathbb{R}^+$. Nous donnons sur les
figures II.1 à II.12 les variations de $g(q)$ et des dérivées de $g$ pour
différentes valeurs de $r$ et $p$. La Figure II.1 donne pour $p=1.1$, les dif-
férentes courbes $g(q)$ pour les valeurs suivantes de $r : 10^{11}$, $2 \cdot 10^{11}$, $4 \cdot 10^{11}$,
$5 \cdot 10^{11}$, $7 \cdot 10^{11}$ et $10^{12}$. La figure II.2 donne les variations de $g'(q)$ pour
les mêmes valeurs de $r$ (la courbe correspondante à $r=10^{12}$ n'est pas tracée
car en dehors du cadre). Sur la figure II.3 on représente $g''(q)$ (évidemment
indépendant de $r$). Sur ces trois premières figures nous avons représenté
la variation au voisinage de l'origine ($q$ variant de $0$ à $10^{-12}$). Sur les
figures II.4, II.5, II.6 nous présentons des courbes analogues, toujours
pour $p=1.1$, mais $q$ variant cette fois-ci entre $0$ et $1$. Les valeurs de $r$
sont cette fois : .1, .2, .5, .7, 1., 2., 5., 7., 10..
La figure II.7 donne pour $p=3$, les différentes courbes $g(q)$ pour $r$ prenant les valeurs $10^{-2}$, $5 \cdot 10^{-2}$, $10^{-1}$, $0.5$, $1.$, $5.$, $10$. La figure II.8 donne les variations de $g'(q)$ pour les mêmes valeurs de $r$. La figure II.9 donne $g''(q)$ qui est évidemment constante vu la valeur de $p$. Toujours en conservant les mêmes valeurs de $r$, on donne sur les figures II.10, II.11, II.12, les courbes analogues pour $p=10$.

La fonction $g(q)$ est convexe et croissante si $p \ge 2$ et concave et croissante pour $p < 2$. La méthode utilisée pour résoudre (2.12) est la méthode de Newton classique à une variable soit :

$$ (2.13) \qquad \xi^0 \in \mathbb{R}^+ \text{ donné} $$

$$ (2.14) \qquad \xi^{m+1} = \xi^m - \frac{g(\xi^m) - \text{const.}}{g'(\xi^m)} $$

$\xi^0$ étant choisi de telle façon que $(g(\xi^0)-\text{const.}) g''(\xi^0) \ge 0$. Ce type d'initialisation a une grande importance surtout pour les valeurs extrêmes de $p$ c'est-à-dire au voisinage de 1 et pour $p$ assez grand.

## II.3. Résultats numériques

Nous avons repris certains exemples déjà traités dans [2] en les complétant avec les nouveaux algorithmes (et le test évitant le blocage). Pour les applications nous nous sommes limités à $\Omega$, cercle découpé en 256 triangles avec 145 noeuds (voir figure II.13).

Nous avons pris $f = \text{cte} = 0.1$ pour les applications. Si nous tenons compte de la formule déterminée empiriquement dans [4], nous pouvons déterminer le choix optimal du paramètre de pénalisation pour une autre valeur de la constante $f$, connaissant celui pour $f=0.1$. Rappelons la formule.

Pour $p$ donné, si $r_o$ est la valeur optimale du paramètre de pénalisation pour la constante $C_o$, alors pour la constante $C_l$, la valeur optimale $r_l$ sera donnée par

$$ (2.15) \qquad \frac{1}{r_l} = \frac{1}{r_o} \cdot \left(\frac{C_l}{C_o}\right)^{\frac{2-p}{p-1}}. $$
Les essais numériques ont été faits pour trois valeurs particulières de p :
p = 1.1 , p = 3. , p = 10.

Sur les figures II.14 à II.22 on a représenté le nombre d'itérations en fonction du paramètre de pénalisation r.

Pour l'algorithme ALG1 ler schéma (1.17)-(1.24), une itération représente

* 1 balayage de type relaxation sur le système linéaire (2.8)

* 1 résolution des équations non linéaires (2.9)

Pour l'algorithme ALG2 ler schéma (1.31)-(1.35) et pour l'algorithme ALG1 2ème schéma (1.25)-(1.30) une itération représente

* 1 résolution du système linéaire (2.8). (Une fois la matrice du problème factorisée on n'a plus que 2 systèmes triangulaires à résoudre à chaque étape).

* 1 résolution des équations non linéaires (2.9).

Sur les figures relatives aux algorithmes ALG1 (ler et 2ème schémas) il y a éventuellement 3 courbes.

- courbe 1 : correspond à l'algorithme ALG1, ler où 2ème schéma, tel qu'il est présenté en (1.17)-(1.24) ou (1.25)-(1.30) avec blocage intempestif à l'itération 4 ou 5 en n.

- courbe 2 : correspond à l'algorithme ALG1 pour lequel on évite le blocage intempestif ; c'est-à-dire qu'on résoud effectivement le problème de point-selle.

- courbe 3 : représente le nombre d'itérations en n (dans le cas de résolution effective du problème de point-selle).

Sur les courbes relatives à l'algorithme ALG2 ler schéma, il peut éventuellement y avoir 4 courbes :

- courbe 1 : elle correspond à $\epsilon = 5.10^{-4}$ dans (1.35)

- courbe 2 : " $\epsilon = 10^{-4}$ "

- courbe 3 : " $\epsilon = 5.10^{-5}$ "

- courbe 4 : " $\epsilon = 10^{-5}$ "
Les figures II.14 à II.16 correspondent à $p=1.1$. La figure II.14 est relative à l'algorithme ALG1 1er schéma, tel qu'il est présenté en (1.17)-(1.24) et cette courbe a déjà été donnée dans [2]. La figure II.15 correspond à l'algorithme ALG 2 1er schéma. On y distingue 4 courbes correspondant à 4 tests d'arrêt différents dans (1.35). La valeur optimale du paramètre de pénalisation se situe aux environs de $4 \cdot 10^{11}$ pour le test d'arrêt $\epsilon = 10^{-5}$.

La figure II.16 correspond à l'algorithme ALG1 2ème schéma. La courbe 1 correspond évidemment à l'algorithme tel qu'il est présenté, et la courbe 2 à l'algorithme dans lequel on évite le phénomène de "blocage". La courbe 3 donnant le nombre d'itérations d'UZAWA, montre que ce nombre diminue lorsque le paramètre de pénalisation augmente, mais le nombre global des itérations (somme des itérations en k dans l'algorithme), augmente lui aussi, la valeur optimale de r se situe aux alentours de $6 \cdot 10^{11}$. Il fallait moins de 100 itérations à l'optimum pour l'algorithme ALG2 1er schéma, il en faut plus de 200 pour l'algorithme ALG1 2ème schéma.

Les figures II.17 à II.19 correspondent à $p=3$. Bien que le choix du paramètre de pénalisation semble moins critique que pour $p=1.1$, on constate un comportement analogue des algorithmes.

Les figures II.20 à II.22 donnent le comportement des algorithmes pour $p=10$.

Sur la figure II.23 nous donnons les équipotentielles de la solution u, les paramètres étant $p=1.1$, $r=10^{12}$, $f=c=0.1$. Cette solution est obtenue avec l'algorithme ALG1 2ème schéma tel qu'il a été décrit, le test d'arrêt dans (1.28) et (1.30) tant $\epsilon = 10^{-5}$.

Sur la figure II.24 on représente la solution obtenue avec ce même algorithme lorsqu'on évite le phénomène de blocage, mais comme on peut le voir la solution n'est pas très différente de celle de la figure précédente. Sur les figures II.25, II.26 on représente la solution u obtenue avec l'algorithme ALG2 1er schéma, avec $\epsilon = 5 \cdot 10^{-4}$ dans 1.35 pour la figure II.25 et $\epsilon = 10^{-5}$ pour la figure II.26. On constate que les solutions obtenues diffèrent très peu et que l'on pourrait se contenter de prendre $\epsilon = 5 \cdot 10^{-4}$ dans (1.35).
Les figures II.27 à II.30 sont relatives à p=3. Comme précédemment on donne les équipotentielles de la solution u pour les algorithmes ALG1 2ème schéma et ALG2 1er schéma. Sur la figure II.27 on donne u obtenu avec ALG1 2ème schéma, tel qu'il est présenté en (1.25)-(1.30) et sur la figure II.28 on donne u obtenu avec ce même algorithme, en évitant le phénomène de "blocage". Le test d'arrêt est dans les 2 cas $ε = 10^{-5}$.

Les figures II.29 et II.30 correspondent à l'algorithme ALG2 1er schéma avec respectivement $ε = 5 \cdot 10^{-4}$ et $ε = 10^{-5}$ dans (1.35). Dans tous les autres cas les solutions obtenues sont très voisines.

Les figures II.31 à II.34 sont relatives à p=10, les deux premières à l'algorithme ALG1, les deux suivantes à l'algorithme ALG2.

La figure II.31 correspond à l'algorithme ALG1, 2ème schéma tel qu'il est présenté en (1.25)-(1.30), la figure suivante représente la solution calculée lorsqu'on évite le phénomène de blocage exposé au paragraphe 1.

Sur la figure II.33, relative à l'algorithme ALG2 1er schéma, on représente la solution obtenue avec $ε = 5 \cdot 10^{-4}$ dans (1.35) alors que sur la figure II.34 on a $ε = 10^{-5}$.
<figure><img src="image_2.png" /><figcaption>Fig. II.1 - Famille de courbes g(q) pour le pseudo-laplacien p=1.1 au voisinage de q=0</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. II.2 - Famille des courbes dérivées d'ordre 1 pour p = 1.1.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. II.3 - Dérivé seconde de g(q) pour p=1.1</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Fig. II.4 - Famille de courbes g(q) pour le pseudo-laplacien p=1.1</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. II.5 - Famille des courbes dérivées d'ordre 1 pour p=1.1</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. II.6 - Dérivée seconde de g(q) pour p=1.1</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Fig. II.7 - Famille de courbes g(q) pour le pseudo-laplacien p=3</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. II.8 - Famille des courbes dérivées d'ordre 1 pour p=3</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. II.9 - Dérivée seconde de g(q) pour p=3</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Fig. II.10 - Famille de courbes g(q) pour le pseudo-laplacien p=10</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. II.11 - Famille des courbes dérivées d'ordre 1 pour p=10</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. II.12 - Dérivée seconde de g(q) pour p=10</figcaption></figure>
<figure><img src="image_14.png" /></figure>

NOMBRE DE NOEUDS 145

NOMBRE D'ÉLÉMENTS 256

LABORIA- A.MARROCCO

<caption>Fig. II.13</caption>

Maillage de Ω

Comportement de l'algorithme en fonction du paramètre de pénalisation r.

<figure><img src="image_15.png" /><figcaption>Fig. II.14</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. II.15</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. II.16</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Fig. II.17</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. II.18</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Fig. II.19</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Fig. II.20</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Fig. II.21</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Fig. II.22</figcaption></figure>
LIGNES DE NIVEAU DE LA SOLUTION

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

PAS ENTRE 2 LIGNES 0.52E-15
VALEURS EXTREMES DE LA FONCTION
0.79E-14
PARAMETRES
C 0.100E+0
R 0.100E+13
P 1.10
-1 0.10E-4

LABORIA- A.MARROCCO

<figure><img src="image_25.png" /><figcaption>Fig. II.23</figcaption></figure>

PAS ENTRE 2 LIGNES 0.52E-15
VALEURS EXTREMES DE LA FONCTION
0.78E-14
PARAMETRES
C 0.100E+0
R 0.500E+12
P 1.10
1 0.10E-4

LABORIA- A.MARROCCO

<figure><img src="image_26.png" /><figcaption>Fig. II.24</figcaption></figure>

PAS ENTRE 2 LIGNES 0.50E-15
VALEURS EXTREMES DE LA FONCTION
0.78E-14
PARAMETRES
C 0.100E+0
R 0.500E+12
P 1.10
0 0.50E-3

LABORIA- A.MARROCCO

<figure><img src="image_27.png" /><figcaption>Fig. II.25</figcaption></figure>

PAS ENTRE 2 LIGNES 0.52E-15
VALEURS EXTREMES DE LA FONCTION
0.78E-14
PARAMETRES
C 0.100E+0
R 0.500E+12
P 1.10
0 0.10E-4

LABORIA- A.MARROCCO

<caption>Fig. II.26</caption>
<figure><img src="image_28.png" /></figure>

PAS ENTRE 2 LIGNES 0.979E-2
VALEURS EXTREMES DE LA FONCTION
0 0.147E+0
PARAMETRES
C 0.100E+0
R 0.100E+1
P 3.00
-1 0.10E-4

LABORIA- A.MARROCCO

<figure><img src="image_29.png" /><figcaption>Fig. II.27</figcaption></figure>

PAS ENTRE 2 LIGNES 0.978E-2
VALEURS EXTREMES DE LA FONCTION
0 0.147E+0
PARAMETRES
C 0.100E+0
R 0.100E+1
P 3.00
1 0.10E-4

LABORIA- A.MARROCCO

<figure><img src="image_30.png" /><figcaption>Fig. II.28</figcaption></figure>

PAS ENTRE 2 LIGNES 0.979E-2
VALEURS EXTREMES DE LA FONCTION
0 0.147E+0
PARAMETRES
C 0.100E+0
R 0.100E-1
P 3.00
0 0.50E-3

LABORIA- A.MARROCCO

<figure><img src="image_31.png" /><figcaption>Fig. II.29</figcaption></figure>

PAS ENTRE 2 LIGNES 0.979E-2
VALEURS EXTREMES DE LA FONCTION
0 0.147E-0
PARAMETRES
C 0.100E-0
R 0.100E-1
P 3.00
0 0.10E-4

LABORIA- A.MARROCCO

<caption>Fig. II.30</caption>
<figure><img src="image_32.png" /><figcaption>Fig. II.31</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>Fig. II.32</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>Fig. II.33</figcaption></figure>

<figure><img src="image_35.png" /><figcaption>Fig. II.34</figcaption></figure>
### III. EXEMPLE 2 - PROBLEME DE MAGNETOSTATIQUE

Nous reprenons ici le problème déjà exposé dans [5] [6] en y ajoutant des résultats numériques. La machine étudiée est un alternateur tetrapolaire.

#### III.1 Présentation du problème.

Les équations de Maxwell de la magnétostatique sont les suivantes

$$ (3.1) \qquad \nabla \times \vec{H} = \vec{j} $$

$$ (3.2) \qquad \vec{B} = \mu \vec{H} $$

$$ (3.3) \qquad \nabla \cdot \vec{B} = 0 $$

où $\vec{H}$ est le champ magnétique, $\vec{j}$ le vecteur densité de courant, $\vec{B}$ l'induction magnétique, $\mu$ la perméabilité magnétique. En introduisant le potentiel vecteur $\vec{A}$ tel que

$$ (3.4) \qquad \vec{B} = \nabla \times \vec{A} $$

les relations (3.1)-(3.2)-(3.3) conduisent à l'équation unique

$$ (3.5) \qquad \nabla \times (\nu_r \nabla \times \vec{A}) = \mu_0 \vec{j} $$

$$ \begin{gathered} \text{où } \nu_r = \frac{1}{\mu_r} \text{ avec } \mu = \mu_o \mu_r \\ \mu_o = C^{\text{te}} = 4\pi 10^{-7} \text{ MKSA (perméabilité du vide)} \end{gathered} $$

et $\mu_r$ étant la perméabilité relative.

Pour les problèmes bidimensionnels, $\vec{j}$ et $\vec{A}$ ont une structure particulière

$$ \vec{A} = (0,0,A_3) , \quad \vec{j} = (0,0,j_3) $$

si bien que (3.5) s'écrit aussi

$$ (3.6) \qquad - \sum_{i=1}^{2} \frac{\partial}{\partial x_i} \left( \nu_r \frac{\partial A_3}{\partial x_i} \right) = \mu_o j_3 $$
ou bien

$$
(3.7) \qquad -\nabla \cdot (v_r \nabla A_3) = \mu_o j_3 .
$$

En analysant de plus près la relation (3.2) qui est linéaire dans l'air ou le
cuivre µ=µ₀ = 4π 10⁻⁷ MKSA, et non linéaire dans les parties ferro-magnétiques,
on peut voir que v_r = 1/µ_r peut s'exprimer naturellement comme une fonction
monotone de |v̂_B|^2. (On ne tient pas compte de l'effet d'hystérésis). Pour les
problèmes bi-dimensionnels cartésiens on a la relation

$$
(3.8) \qquad |\vec{B}|^2 = |\nabla \times \vec{A}|^2 = |\nabla A_3|^2
$$

Ainsi l'équation (3.6) s'écrit plus explicitement

$$
(3.9) \quad - \sum_{i=1}^{2} \left( \frac{\partial}{\partial x_i} v_r(x, |\nabla A_3|^2) \right) \frac{\partial A_3}{\partial x_i} = \mu_o j_3
$$

Finalement si on note par $\psi$ la fonction telle que

$$
(3.10) \qquad \begin{cases} \displaystyle \frac{\partial \psi}{\partial |\vec{B}|^2} (x, |\vec{B}|^2) = v_r(x, |\vec{B}|^2) \\ \displaystyle \psi(x, 0) = 0 \end{cases}
$$

Il est facile de voir que (3.9) est l'équation d'Euler, dans un espace fonc-
tionnel adéquat du problème de l'optimisation de l'énergie magnétique, la
fonctionnelle d'énergie étant donnée par

$$
(3.11) \quad \mathcal{F}(A) = \frac{1}{2} \int_{\Omega} \psi(x, |\nabla A|^2) dx - \mu_o \int_{\Omega} j \cdot A \, dx
$$

III.2. *Approximation de* $v_r$

Dans les parties ferromagnétiques, l'inverse de la caractéristique magnétique (relative), exprimée comme fonction de |v̂_B|^2, est approchée, au sens des moindres carrés, par une fonction de la famille suivante
$$
(3.12) \qquad \tilde{\epsilon}_{\alpha, C, T}(z) = \epsilon + (C-\epsilon) \frac{z^{\alpha}}{z^{\alpha+T}}
$$

$\tilde{v}$ est positive, bornée, croissante, les paramètres $\epsilon, \alpha, C, T$ étant positifs et $C > \epsilon$.

III.3 Formulation du problème

Soit H<sup>1</sup><sub>O</sub>(Ω) l'espace de Sobolev

$$
(3.13) \quad H_o^1(\Omega) = \{v | v \in L^2(\Omega), \frac{\partial v}{\partial x_i} \in L^2(\Omega), v|_T = 0\}
$$

avec la norme

$$
(3.14) \quad \|v\|^2_{H_o^1(\Omega)} = \int_{\Omega} |\nabla v|^2 dx \quad \text{qui en fait un espace de Hilbert}
$$

Le problème de magnétostatique se formule ainsi

$$
(3.15) \quad \left\{
\begin{array}{l}
\text{Trouver } A \in H_o^1(\Omega) \text{ tel que} \\
\mathfrak{F}(A) \le \mathfrak{F}(v) \quad \forall v \in H_o^1(\Omega)
\end{array}
\right.
$$

où $\mathfrak{F}$ est donné par (3.11)

$\mathfrak{F}$ étant continue sur $H_o^1(\Omega)$, strictement convexe et $\mathfrak{F}(v) \to \infty$ lorsque
$\|v\|_{H_o^1(\Omega)} \to \infty$, le problème d'optimisation (3.15) admet une solution unique,

aussi caractérisée par :

$$
(3.16) \quad \left\{
\begin{array}{@{}l@{}}
\displaystyle \int_{\Omega} v_r(x, |\nabla A|^2) \nabla A \cdot \nabla v \, dx - \mu_o \int_{\Omega} j \cdot v \, dx = 0 \quad \forall v \in H_o^1(\Omega) \\[\medskipamount]
\quad A \in H_o^1(\Omega).
\end{array}
\right.
$$

### III.4 Formulation du problème en lagrangien augmenté (discret)

Si on traduit l'écriture (1.7) pour le problème de magnétostatique on obtient
le lagrangien augmenté

$$
(3.17) \quad \left\{
\begin{array}{@{}l@{}}
\displaystyle \mathcal{L}_r(v,z,\mu) = \frac{1}{2} \int_{\Omega} \psi(x, |z|^2) \, dx - \mu_o \int_{\Omega} j \cdot v \, dx + \frac{r}{2} \int_{\Omega} |\nabla v - z|^2 \, dx \\
\displaystyle \qquad + \int_{\Omega} \mu \cdot (\nabla v - z) \, dx.
\end{array}
\right.
$$
$\mathcal{L}_r$ étant toujours une application de $V_h \times L_h \times L_h \to R$.

Dans (3.17) on peut remplacer le terme pénalisateur $\frac{r}{2} \int_{\Omega} |\nabla v - z|^2 dx$ par un
terme du type

$$
(3.18) \quad \frac{r}{2} \int_{\Omega} \eta |\nabla v - z|^2 dx
$$

où $\eta$ est une fonction de x, strictement positive. D'un point de vue pratique
il est tout à fait indiqué de la prendre constante par triangle. Nous don-
nerons plus loin une interprétation de cette fonction $\eta$. En prenant une
telle fonction $\eta$ on peut rendre la pénalisation plus homogène. La théorie
et les algorithmes indiqués au Chapitre 3 de [1] ne changent pas si on remplace
aussi le terme de dualité par

$$
(3.19) \quad \int_{\Omega} \eta \, \mu \cdot (\nabla v - z) \, dx .
$$

Le lagrangien augmenté que nous considérerons pour le problème de magnéto-
statique sera donc le suivant :

$$
(3.20) \quad \left\{
\begin{array}{@{}l@{}}
\displaystyle \mathcal{L}_r(v,z,\mu) = \frac{1}{2} \int \psi(x,|z|^2) \, dx - \mu_o \int_{\Omega} j \cdot v \, dx + \frac{r}{2} \int_{\Omega} \eta |\nabla v - z|^2 \, dx \\
\displaystyle \qquad + \int_{\Omega} \eta \mu \cdot (\nabla v - z) \, dx.
\end{array}
\right.
$$

Explicitons maintenant le système (1.11)-(1.12) que l'on peut écrire de la
façon suivante

$$
(3.21) \quad \left\{
\begin{array}{@{}l@{}}
\displaystyle \int_{\Omega} \eta \nabla \hat{u}_h \cdot \nabla v_h \, dx = \frac{1}{r} \mu_o \int_{\Omega} j \cdot v_h \, dx + \int_{\Omega} \eta (\hat{p}_h - \frac{1}{r} \mu_h) \cdot \nabla v_h \, dx \\
\quad \forall v_h \in V_h, \hat{u}_h \in V_h
\end{array}
\right.
$$

$$
(3.22) \quad \left\{
\begin{array}{@{}l@{}}
\displaystyle \int_{\Omega} [v_r(x, |\hat{p}_h|^2) \hat{p}_h + r\eta \hat{p}_h] q_h \, dx = \int_{\Omega} \eta (r \nabla \hat{u}_h + \mu_h) \cdot q_h \, dx \\
\quad \forall q_h \in L_h, \hat{p}_h \in L_h.
\end{array}
\right.
$$
Il est à noter que le "pas" (1.10) reste inchangé, et que la fonction η étant constante par triangle, on peut faire avec le problème (3.22) la même démarche qu'avec le problème (2.9) et l'équation à l'inconnue que l'on doit résoudre (sur chaque triangle T de h) est

$$ (3.23) \qquad g(q) = (v_r(x, |q|^2) + r\eta)q = \eta |r \nabla u_h + \mu_h| . $$

La dépendance de $v_r$ en x signifie que nous avons 2 types d'équations suivant que l'on se trouve dans l'air ou dans le fer. Dans l'air (ou le cuivre), $v_r$, exprimé comme fonction de $|q|^2$ est une constante, on obtient donc $q$ explicitement dans (3.23).

Dans le fer, $v_r$ est donné par une expression du type (3.12). Ainsi la fonction $g(q)$ est une fonction strictement croissante mais possédant un point d'inflexion qui heureusement ne dépend ni de $r$, ni de $\eta$. On pourra donc utiliser la méthode de Newton pour résoudre le problème associé à (3.23), à condition de bien initialiser le processus.

Nous donnons sur les figures III.1 et III.2 les variations de $g(q) - (3.23)$ pour différentes valeurs du produit $r\eta$. La figure III.1 correspond à la fonction $v_r$ retenue pour le stator, la figure III.2 pour le rotor.

La figure III.3 correspond au découpage de la courbe représentant $g(q)$ (3.23) pour pouvoir initialiser le processus de Newton en toute sécurité.

### III.5. Résultats numériques.

Pour les applications numériques nous nous sommes limités à l'alternateur tetrapolaire. Le domaine $\Omega$ est formé d'un quart de cercle ; celui-ci est découpé en 812 triangles et il y a 384 noeuds intérieurs (voir figure III.4).

Nous pouvons dire ici pourquoi nous avons été amenés à introduire la fonction $\eta$ dans la formulation du lagrangien augmenté (3.20), et par la même voir comment on peut interpréter cette fonction $\eta$.

Le problème de magnétostatique, pour les faibles courants d'excitation $j$, est pratiquement un problème linéaire à milieu non homogène. En effet, le courant étant "petit", l'induction créée est aussi petite et si on regarde la
fonction $v_r$ (inverse de la perméabilité magnétique relative), on voit que cette fonction est pratiquement constante pour les faibles valeurs de l'argument (induction) ; celà est évidemment dû à la valeur très grande du paramètre T dans (3.12). On a représenté sur la figure III.5 la courbe $v_r(|\vec{B}|^2)$ pour le rotor. (Il y a évidemment une courbe analogue pour le stator).

Donc pour les faibles valeurs du courant j on est pratiquement conduit à résoudre un problème linéaire du type suivant :

$$ (3.24) \qquad \left\{ \begin{array}{l} \displaystyle \int_{\Omega} v_r(x) \nabla A \cdot \nabla v \, dx - \mu_o \int_{\Omega} j \cdot v \, dx = 0 \quad \forall v \in H_o^1(\Omega) \\[1em] \qquad A \in H_o^1(\Omega) \end{array} \right. $$

la fonction $v_r$ étant constante par morceaux, égale à 1 dans l'air ou le cuivre, à $\epsilon_S$ dans le stator et $\epsilon_R$ dans le rotor. Maintenant si on revient à la formulation en lagrangien augmenté standard donné en (1.7), c'est-à-dire

$$ (3.25) \qquad \left\{ \begin{array}{l} \displaystyle \ell_r(v,z,\mu) = \frac{1}{2} \int_{\Omega} \psi(x,|z|^2) dx - \mu_o \int_{\Omega} j \cdot v \, dx + \frac{r}{2} \int_{\Omega} |\nabla v - z|^2 dx + \\ \displaystyle \qquad \int_{\Omega} \mu \cdot (\nabla v - z) dx \end{array} \right. $$

et supposons qu'on veuille résoudre le problème linéaire (3.24), on n'obtient pas le résultat en une seule itération (même en prenant r=1).

On rétablit pour ainsi dire la "propriété de linéarité" en introduisant le fonction $\eta$ valant 1 dans l'air, $\epsilon_S$ dans le stator, et $\epsilon_R$ dans le rotor. Avec le lagragien (3.20) et la fonction $\eta$ ainsi définie, on obtient la solution du problème de magnétostatique linéaire en **1 seule itération** avec l'algorithme ALG2 pour le schéma, pour la valeur r=1 du paramètre de pénalisation.

Il est donc naturel d'introduire cette fonction $\eta$, et l'on vient de voir qu'elle a la même signification que l'inverse de la perméabilité relative ($v_r$). Ainsi, par l'introduction de cette fonction $\eta$, **on peut** dans les parties ferromagnétiques estimer la solution en termes de répartition de l'induction,
la valeur de η correspondante étant calculée par (3.12). Estimer cette
solution ainsi, l'électrotechnicien sait le faire, car c'est ce qu'il
fait lorsqu'il fait des calculs à la main !

Vu l'approximation du problème par éléments finis $P_1$ conformes, on en déduit
évidemment que la perméabilité (ou l'inverse) est constante sur chaque tri-
angle $T \subset \mathcal{T}_h$. Cette remarque liée à l'interprétation de la fonction $\eta$ justi-
fie pleinement qu'on prenne cette fonction $\eta$ constante par triangle.

Que fait donc l'algorithme de pénalisation-dualité ? il corrige l'estimation
de la solution faite par l'électrotechnicien ; si cette estimation est bonne
il n'aura que très peu de travail à faire, il convergera donc très vite.

Nous n'avons utilisé pour les essais numériques que l'algorithme ALG2
ler schéma (1.31)-(1.35). Cet algorithme s'est déjà montré le plus efficace
pour l'exemple précédent. Nous avons à notre disposition un paramètre sup-
plémentaire si l'on peut dire, la fonction η qui sera définie comme étant
constante sur chaque triangle $T \in \mathcal{T}_h$.

Nous connaissons les valeurs minima de la résistivité magnétique relative
dans le stator $\varepsilon_S$ et dans le rotor $\varepsilon_R$, et nous nous servirons de ces valeurs
pour initialiser la fonction $\eta$. Nous avons choisi les initialisations suivantes, on se donne

$$
(3.26) \qquad k \in \mathbb{R} - \{0\}
$$

et on définira $\eta_k$ par :

$$
(3.27) \quad \left\{
\begin{array}{ll}
\eta_k = 1 & \text{si le triangle est dans l'air ou dans le cuivre} \\
\eta_k = k \varepsilon_R & \text{si le triangle est dans le rotor} \\
\eta_k = k \varepsilon_S & \text{si le triangle est dans le stator.}
\end{array}
\right.
$$

En gardant la même notation $\eta_k$, on conviendra que *si k est négatif* on fait
l'initialisation suivante :
On résout le problème linéaire

$$
(3.28) \quad \left\{
\begin{array}{@{}l@{}}
\displaystyle \int_{\Omega} \eta_{|k|} \nabla A_h \cdot \nabla v_h \, dx = \mu_o \int_{\Omega} j \cdot v_h \, dx \\
\forall v_h \in V_h, \ A_h \in V_h.
\end{array}
\right.
$$

Ensuite dans les parties ferromagnétiques, on calcule $|\nabla A_h|^2$ et on initialise
$\eta_k$ par

$$
(3.29) \quad \left\{
\begin{array}{ll}
\eta_k = 1 & \text{si le triangle est dans l'air ou le cuivre} \\
\eta_h = v_r (|\nabla A_h|^2) & \text{dans les parties ferromagnétiques.}
\end{array}
\right.
$$

Bien sûr, ces initialisations (3.27) ou (3.29) sont un peu artificielles et
il vaudrait peut-être mieux se fier ici au bon sens de l'électrotechnicien
qui définirait, lui, une répartition plus vraisemblable et plus proche de
la réalité des inductions magnétiques et on aurait par la même une meilleure
définition des perméabilités magnétiques estimées. Ceci étant dit, nous
présentons tout de même les résultats obtenus avec les initialisations
(3.27), ou (3.29).

Nous avons d'abord fait des essais avec j = 0.5 Å/mm² ; cette valeur n'est
donnée qu'à titre indicatif, disons simplement qu'avec cette valeur du cou-
rant d'excitation j, le problème se comporte comme un problème linéaire.
On le voit bien d'ailleurs sur la figure III.6, correspondant à k=1 dans
(3.27) ; pour r=1 il faut seulement 2 itérations (ce qui est le minimum,
à cause du test, la solution est bien évidemment obtenue dès la lère itéra-
tion). Nous représentons sur cette figure III.6 la variation du nombre
d'itérations en fonction du paramètre de pénalisation r. La numérotation
des courbes (sur les figures suivantes également) a la même signification
qu'en II.3. Les résultats de la figure III.7 sont obtenus pour j=0.5 et
k=10 dans (3.27). Il est évidemment inutile de faire d'autres essais, car
dans ce cas particulier on sait comment il faut choisir les paramètres :
k=1, r=1 !
Passons maintenant à des essais numériques avec un courant d'excitation un peu plus fort, conduisant à un problème non linéaire (mais "pas trop" si l'on peut dire). Dans les figures III.8 à III.16 nous donnons la variation du nombre d'itérations de l'algorithme ALG2 1er schéma, nécessaire à la vérification du test (1.35), en fonction du paramètre de pénalisation r, la densité de courant étant $j=2 A/mm^2$.

Pour la figure III.8, correspondant à $k=1$ dans (3.27), la valeur optimale de r est voisine de 10, et il faut 90 itérations pour réaliser le test (1.35) avec $\epsilon = 10^{-5}$.

Sur la figure III.9, correspondant à $k=5$ dans (3.27) on voit que la valeur optimale de r a diminué, ainsi que le nombre d'itérations. Le paramètre r optimal se situe vers 4, et il ne faut plus que 57 itérations, toujours pour obtenir (1.35) avec $\epsilon = 10^{-5}$. Pour la figure III.10, correspondant à $k=10$ (dans 3.27), on obtient $r_{opt} \approx 2$ et 55 itérations. Sur la figure III.11, correspondant à $k=20$ dans (3.27), on voit que la valeur optimale du paramètre de pénalisation r est voisine de 1, le nombre d'itérations se situant toujours aux environs de 50.

Si on prend $k=50$ (fig. III.12), on surestime franchement la valeur des relativités relatives $v_r$, mais l'algorithme a un comportement meilleur, pour $r=0.3$, il faut moins de 40 itérations. (Peut-être que le conditionnement du système linéaire (3.21) y est pour quelque chose).

Les figures III.13 à III.16 sont relatives aux résultats obtenus avec k négatif c'est-à-dire avec l'initialisation de $\eta$ par (3.28)-(3.29).

La figure III.13 correspond à $k=-1$. Le creux que l'on distingue près de $r=1$ est certainement dû au fait que le problème n'est pas trop non linéaire.

Pour les figures suivantes, k augmentant en valeur absolue, on constate que la valeur optimale de r augmente ainsi d'ailleurs que le nombre d'itérations ; lorsque k est assez grand en valeur absolue, on retombe évidemment sur une initialisation équivalente à $k=+1$.
Les figures III.17 à III.25 sont relatives aux résultats obtenus avec une densité de courant $j=7.5 \text{ A/mm}^2$ (marche normale de la machine) ; la non linéarité du problème est ici assez prononcée.

Sur les figures III.17 à III.20 nous prenons successivement $k=1, 10, 20, 50$ ; nous constatons que le choix optimal du maramètre de pénalisation $r$ diminue de la valeur 15 à la valeur 1., le nombre d'itérations nécessaires à la réalisation de (1.35) diminue aussi. Il faut à peine plus de 50 itérations pour $k=50$ et $r=1$.

Sur la figure suivante III.21 nous représentons les résultats obtenus lorsque la fonction $\eta$ est partout initialisée à 1 c'est-à-dire on utilise le lagrangien augmenté (3.25). Les performances sont tout de même moins bonnes que celles montrées sur la figure III.20.

Les figures III.22 à III.25 donnent les résultats obtenus avec $k$ négatif (initialisation de $\eta$ par (3.28)-(3.29)). Nous avons pris successivement $k=-1, -10, -20, -50$. D'après les résultats obtenus, la meilleure valeur de $k$ se situerait entre -20 et -10. Lorsqu'on prend $k$ assez grand en valeur absolue ($|k| = 50$), on retrouve les résultats pour $k=+1$.

Sur les figures III.26 à III.31 on a représenté les courbes d'itération en fonction du paramètre de pénalisation $r$, pour une densité de courant $j=20 \text{ A/mm}^2$ conduisant à de très fortes saturations dans les parties ferromagnétiques.

Nous prenons successivement $k=1, 10, 50$ sur les figures III.26 à III.28 et $k=-1, -10, -50$ sur les figures III.29 - III.31. Nous avons constaté sur les figures III.6 à III.31 l'influence du choix de la fonction $\eta$ sur le nombre d'itérations ; (le meilleur choix de $\eta$ serait évidemment celui correspondant à la solution du problème !). Insistons encore sur le fait que l'algorithme de pénalisation dualité, surtout dans un contexte intéractif, est d'utilisation très agréable pour l'electrotechnicien, car il lui permet de fournir (si celui-ci le désire), une estimation de la solution ; l'algorithme corrige ensuite les erreurs éventuelles de cette estimation.
Reprenons l'équation variationnelle (3.21)

$$ (3.30) \quad \left\{ \begin{array}{l} \displaystyle \int_{\Omega} \eta \nabla A_h \cdot \nabla v_h \, dx = \frac{1}{r} \mu_o \int_{\Omega} j \cdot v_h \, dx + \int_{\Omega} \eta (\hat{p}_h - \frac{1}{r} \mu_h) \cdot \nabla v_h \, dx \\[1em] \forall v_h \in V_h, A_h \in V_h. \end{array} \right. $$

Le second membre de (3.30) est composé de 2 termes, il y apparaît (à un facture multiplicatif près), le second membre du problème initial (3.16).

Nous calculons $\int_{\Omega} j \cdot v_h \, dx$ pour chaque élément $v_h$ appartenant à la base canonique de $V_h$. On interpole linéairement ces valeurs sur la triangulation, on obtient une certaine fonction que l'on représente sur la figure III.32.

On calculera et représentera de la même façon le second membre de (3.30), c'est ce qui sera noté "second membre final" sur les dessins.

Nous présentons des résultats pour la densité de courant $j=7.5 \text{ A/mm}^2$ (soit $0.75 \times 10^{-7}$ dans le système MKSA).

Sur la figure III.33, nous preprésentons les lignes d'induction dans l'alternateur, c'est-à-dire les équipotentialles de la fonction $A_h$ (solution du problème de magnétostatique).

Sur la figure III.34, nous représentons les équipotentialles de la fonction "second membre final", lorsque $k=50$ dans (3.27), le paramètre de pénalisation valant 1. Voilà en quelque sorte la fonction qu'il faut fournir comme second membre de (3.30) pour que la solution de ce problème linéaire, soit aussi la solution du problème non linéaire initial. La figure III.35 donne le "second membre final" lorsqu'on prend $k=20$ dans (3.27), le paramètre de pénalisation $r=3$.

Lorsque $k$ est négatif (convention d'écriture précisée en (3.28), (3.29)) on peut représenter les équipotentialles de la fonction $A_h$ obtenue en résolvant (3.28).
Ainsi la figure III.36, pour laquelle $k=-1$, représente les équipotentielles de la solution $A_h$ de (3.28). Dans ce cas particulier ($k=-1$), on peut dire que cette "solution initiale" est la solution que l'on obtiendrait si le problème de magnétostatique était un problème linéaire (quelle que soit la valeur de $j$).

Sur la figure II.37, on voit en fait le second membre de (3.30) qu'il faut choisir pour obtenir la vraie solution du problème de magnétostatique (problème non linéaire), donnée sur la figure III.33.

Sur la figure III.38 ($k=-10$) on a représenté les équipotentielles de $A_h$, solution de (3.28) ; la solution ici, n'est pas "trop éloignée" de la solution réelle du problème non linéaire (la valeur maximum de $A_h$ vaut $2.54 \times 10^{-2}$) ; mais il y a tout de même certaines différences, ainsi sur la figure III.33 donnant la solution du problème non linéaire, on distingue des "fuites" dans les encoches, fuites dues à la saturation magnétique du stator ; cela n'apparaît pas sur la figure III.38. La figure III.39 donne dans ce cas les équipotentielles de la fonction "second membre final" dans (3.30).
<figure><img src="image_36.png" /><figcaption>Fig. III.1 - Famille de courbes g(q) pour le problème de magnétostatique dépendant des caractéristiques magnétiques du stator.</figcaption></figure>

<figure><img src="image_37.png" /><figcaption>Fig. II.3</figcaption></figure>
<figure><img src="image_38.png" /><figcaption>Fig. III.4
Maillage de Ω</figcaption></figure>

<figure><img src="image_39.png" /><figcaption>Fig. III.5
Caractéristiques magnétiques du rotor.</figcaption></figure>
<figure><img src="image_40.png" /><figcaption>Fig. III.6</figcaption></figure>

<figure><img src="image_41.png" /><figcaption>Fig. III.7</figcaption></figure>

<figure><img src="image_42.png" /><figcaption>Fig. III.8</figcaption></figure>

<figure><img src="image_43.png" /><figcaption>Fig. III.9</figcaption></figure>

<figure><img src="image_44.png" /><figcaption>Fig. III.10</figcaption></figure>

<figure><img src="image_45.png" /><figcaption>Fig. III.11</figcaption></figure>
<figure><img src="image_46.png" /><figcaption>Fig. III.12</figcaption></figure>

<figure><img src="image_47.png" /><figcaption>Fig. III.13</figcaption></figure>

<figure><img src="image_48.png" /><figcaption>Fig. III.14</figcaption></figure>

<figure><img src="image_49.png" /><figcaption>Fig. III.15</figcaption></figure>

<figure><img src="image_50.png" /><figcaption>Fig. III.16</figcaption></figure>

<figure><img src="image_51.png" /><figcaption>Fig. III.17</figcaption></figure>
<figure><img src="image_52.png" /><figcaption>Fig. III.18</figcaption></figure>

<figure><img src="image_53.png" /><figcaption>Fig. III.19</figcaption></figure>

<figure><img src="image_54.png" /><figcaption>Fig. III.20</figcaption></figure>

<figure><img src="image_55.png" /><figcaption>Fig. III.21</figcaption></figure>

<figure><img src="image_56.png" /><figcaption>Fig. III.22</figcaption></figure>

<figure><img src="image_57.png" /><figcaption>Fig. III.23</figcaption></figure>
<figure><img src="image_58.png" /><figcaption>Fig. III.24</figcaption></figure>

<figure><img src="image_59.png" /><figcaption>Fig. III.25</figcaption></figure>

<figure><img src="image_60.png" /><figcaption>Fig. III.26</figcaption></figure>

<figure><img src="image_61.png" /><figcaption>Fig. III.27</figcaption></figure>

<figure><img src="image_62.png" /><figcaption>Fig. III.28</figcaption></figure>

<figure><img src="image_63.png" /><figcaption>Fig. III.29</figcaption></figure>
<figure><img src="image_64.png" /><figcaption>Fig. III.30</figcaption></figure>

<figure><img src="image_65.png" /><figcaption>Fig. III.31</figcaption></figure>

<figure><img src="image_66.png" /><figcaption>Fig. III.32
Lignes de niveau du second membre "courant"</figcaption></figure>

<figure><img src="image_67.png" /><figcaption>Fig. III.33
Lignes de niveau de la solution (ou lignes d'induction) correspondant à un courant d'excitation de $7.5 \times 10^6$ MKSA</figcaption></figure>
<figure><img src="image_68.png" /><figcaption>Fig. III.34 - Lignes de niveau du second membre de (3.30) à la convergence du processus</figcaption></figure>

<figure><img src="image_69.png" /><figcaption>Fig. III.35 - Lignes de niveau du second membre de (3.30) à la convergence du processus</figcaption></figure>

<figure><img src="image_70.png" /><figcaption>Fig. III.36</figcaption></figure>

Solution initiale du potentiel vecteur correspondant aussi à la solution "linéaire" du problème de magnétostatique (K=-1).
<figure><img src="image_71.png" /><figcaption>Fig. III.37</figcaption></figure>

Lignes de niveau du second membre de (3.30) à la convergence du processus

<figure><img src="image_72.png" /><figcaption>Fig. III.38</figcaption></figure>

Solution initiale du potentiel vecteur correspondant aussi à la solution "linéaire" du problème de magnétostatique ($K = -10$)

<figure><img src="image_73.png" /><figcaption>Fig. III.39</figcaption></figure>

Lignes de niveau du second membre de (3.30) à la convergence du processus.
## IV. EXEMPLE 3 - ECOULEMENT D'UN FLUIDE DE BINGHAM

On reprend ici l'exemple 1 du chapitre 3 de [1].

### IV.1. Présentation du problème

$\Omega$ étant un ouvert borné de $\mathbb{R}^2$, de frontière $\Gamma$, on définit sur $H_0^1(\Omega)$ (3.13), la fonctionnelle suivante

$$ (4.1) \qquad J_g(v) = \frac{v}{2} \int_{\Omega} |\nabla v|^2 dx + g \int_{\Omega} |\nabla v| dx - \int_{\Omega} f \cdot v dx. $$
$$ v > 0, g \ge 0, f \in L^2(\Omega). $$

On cherche à minimiser (4.1) sur $H_0^1(\Omega)$.

$J_g$ étant continue, strictement convexe sur $H_0^1(\Omega)$, avec $\lim_{g \to \infty} J_g(v) = +\infty$ lorsque $\|v\|_{H_0^1(\Omega)} \to \infty$, il y a existence et unicité de la solution du

problème d'optimisation de $J_g$ ; la solution $u_g$ étant caractérisée par l'inéquation variationnelle

$$ (4.2) \qquad \left\{ \begin{array}{l} \displaystyle \int_{\Omega} \nabla u \cdot \nabla (v-u) dx + g \int_{\Omega} |\nabla v| dx - g \int_{\Omega} |\nabla u| dx \ge \int_{\Omega} f \cdot (v-u) dx \\[1em] \forall v \in H_0^1(\Omega), \ u \in H_0^1(\Omega) \end{array} \right. $$

On démontre dans [7] que dans le cas particulier d'un écoulement laminaire stationnaire d'un fluide de Bingham, dans une conduite cylindrique de section $\Omega$, que la vitesse de l'écoulement $u$ est solution de (4.2), $f$ étant la chute linéique de pression (constante dans la pratique), $g$ étant une constante positive représentant un seuil de contraintes en deçà duquel le fluide se comporte comme un matériau rigide et au delà duquel il se comporte comme un fluide incompressible de viscosité $v$.

### IV.2. Formulation discrète en Lagrangien augmenté

Le lagrangien augmenté pour le problème de l'écoulement d'un fluide de Bingham dans une conduite cylindrique s'écrit (voir aussi chap. 3 (1.23) de [1]).
$$
(4.3) \quad \left\{
\begin{aligned}
\mathcal{L}_r(v_h, q_h, \mu_h) &= \frac{v}{2} \int_{\Omega} |q_h|^2 dx + g \int_{\Omega} |q_h| dx - \int_{\Omega} f \cdot v_h dx + \\
&\quad + \frac{r}{2} \int_{\Omega} |\nabla v_h - q_h|^2 dx + \int_{\Omega} \mu_h \cdot (\nabla v_h - q_h) dx
\end{aligned}
\right.
$$

avec $\mathcal{L}_r : V_h \times L_h \times L_h \to \mathbb{R}$.

Le système (1.11)-(1.12) s'écrit ici :

$$
(4.4) \quad \left\{
\begin{array}{@{}l@{\quad}l@{}}
\displaystyle \int_{\Omega} \nabla \hat{u}_h \cdot \nabla v_h \, dx = \frac{1}{r} \int_{\Omega} f \cdot v_h \, dx + \int_{\Omega} (\hat{p}_h - \frac{1}{r} \mu_h) \cdot \nabla v_h \, dx & \\
\forall v_h \in V_h, \quad \hat{u}_h \in V_h. &
\end{array}
\right.
$$

On ne peut cette fois écrire (1.12) sous forme d'équation variationnelle comme
dans (2.9) ou (3.22), la fonction F (1.6) n'étant pas différentiable, on a
donc l'inéquation variationnelle en la variable $q_h$ de (4.3)

$$
(4.5) \quad \left\{
\begin{array}{@{}l@{\quad}l@{}}
\displaystyle \int_{\Omega} (v+r) \hat{p}_h \cdot (q_h - \hat{p}_h) \, dx + g \int_{\Omega} (|q_h| - |\hat{p}_h|) \, dx \ge \int_{\Omega} (r \nabla \hat{u}_h + \mu_h) (q_h - \hat{p}_h) \, dx & \\
\forall q_h \in L_h, \quad \hat{p}_h \in L_h &
\end{array}
\right.
$$

qui est évidemment la condition d'optimalité du problème d'optimisation
suivant :

$$
(4.6) \quad \left\{
\begin{array}{@{}l@{\quad}l@{}}
\text{Trouver } y_h \in L_h \text{ tel que} & \\
\mathcal{G}(y_h) \le \mathcal{G}(z_h) \quad \forall z_h \in L_h &
\end{array}
\right.
$$

avec

$$
(4.7) \quad \mathcal{G}(z) = \frac{(v+r)}{2} \int_{\Omega} |z|^2 dx + g \int_{\Omega} |z| dx - \int_{\Omega} (r \nabla u_h + \mu_h) \cdot z \, dx.
$$

L'inéquation variationnelle (4.5) se ramène à une famille d'inéquations indé-
pendantes sur chaque triangle (comme déjà indiqué au chapitre 3).

Utilisant les notations du chapitre 3, on aboutit à une famille de problèmes
du type
$$
(4.8) \quad \left\{
\begin{array}{l}
\xi \in \mathbb{R}^2 \\
\inf_{\xi} \left( \frac{v+r}{2} |\xi|^2 + g|\xi| - (d(x), \xi) \right)
\end{array}
\right.
$$

où $|\cdot|$ et $(\cdot, \cdot)$ désignent la norme et le produit scalaire dans $\mathbb{R}^2$, et où on
a posé pour $x \in \Omega$

$$
(4.9) \quad d(x) = \mu_h(x) + r \nabla u_h(x)
$$

(4.9) est constant par triangle T, on notera plutôt $d_T$, si $\bar{\xi}_T$ est la solution
de (4.8) on a alors

$$
(4.10) \quad \left\{
\begin{array}{ll}
\bar{\xi}_T = 0 & \text{si } g \ge d_T \\
\bar{\xi}_T = \frac{1}{v+r} (d_T - g \frac{d_T}{|d_T|}) & \text{si } g < d_T
\end{array}
\right.
$$

Si on utilise les notations de [8] , en termes de "sous-différentiel", on
aboutit pour chaque triangle à une équation généralisée du type

$$
(4.11) \quad (v+r) \hat{p}_h - (r \nabla \hat{u}_h + \mu_h) \in \partial f(\hat{p}_h)
$$

ou

$$
(4.12) \quad \partial f(y) = \begin{cases} g \frac{y}{|y|} & \text{si } y \neq 0 \\ g B(0,1) & \text{si } y=0 \end{cases}
$$

$$
\text{avec } B(0,1) = \{ \xi \in \mathbb{R}^2 , |\xi| \le 1 \} .
$$

On obtient évidemment la même solution (4.10) que l'on peut écrire en
condensé sous la forme

$$
(4.13) \quad \hat{p}_h = \frac{r \nabla \hat{u}_h + \mu_h}{v+r} \quad \text{Max}(0;1-\frac{g}{|r \nabla \hat{u}_h + \mu_h|})
$$

(4.13) est la traduction "algébrique" du problème graphique suivant
(on suppose pour faire le dessin que l'on est non plus dans $\mathbb{R}^2$ mais dans $\mathbb{R}$).
<figure><img src="image_74.png" /><figcaption>Figure IV.1</figcaption></figure>

L'analogue des problèmes (2.12) ou (3.23) est ici que l'on doit calculer l'abcisse de l'intersection de la courbe (fig. IV.1) avec la courbe constante $a = r \nabla \hat{u}_h + \mu_h$. Il est facile de voir que les différents cas possibles peuvent s'exprimer par la seule formule (4.13).

Si on veut une expression faisant absolument intervenir le "module de q" comme dans (2.12) ou (3.23), on la tire immédiatement de (4.13).

### IV.3. Résultats numériques.

Pour les applications numériques nous avons pris $v=1$ et $f=C=constante$ dans (4.1). Les essais numériques ont été faits avec $\Omega = cercle de rayon 1$, découpé en 256 triangles et 145 noeuds (fig. II.13).

Dans le cas particulier d'un cercle et avec $f = cte$, nous avons la solution analytique du problème (cf. par exemple [9]), le seuil critique pour la plasticité est $g_c = \frac{C}{2}$

$$ (4.14) \quad \begin{cases} \text{si } g \ge g_c \text{ alors } u=0 \text{ est solution du problème} \\ \text{si } g < g_c \\ u(r) = \left(\frac{1-r}{2}\right) \left[ \frac{C}{2} (1+r) - 2g \right] & R' \le r \le R \\ u(r) = \left(\frac{1-R'}{2}\right) \left[ \frac{C}{2} (1+R') - 2g \right] & 0 \le r \le R' \\ \text{ou } R' = \frac{2g}{C} \end{cases} $$

(le r qui intervient dans (4.14), n'a évidemment rien à voir avec le paramètre de pénalisation)
On voit d'après (4.14) que c'est le rapport $\frac{g}{c}$ qui est intéressant. Pour les essais numériques nous nous sommes donc limités à la valeur C=20, et nous avons pris 3 valeurs de g : 2, 5, 8.

Dans le cas g=2, la zone rigide doit être un cercle centré en O et de rayon 0.2, pour g=5 on doit avoir un cercle de rayon 0.5 et pour g=8 un cercle de rayon 0.8. Les figures IV.2 à IV.7 donnent le comportement des algorithmes ALG2 1er schéma (1.31)-(1.35), et ALG1 2ème schéma (1.25)-(1.30) par rapport au paramètre de pénalisation r.

Sur la figure IV.2, on a représenté le nombre d'itérations nécessaires à la réalisation du test (1.35) (ALG2 1er schéma), lorsque le paramètre de plasticité g vaut 2. Pour cette valeur de g le problème n'est pas trop "non linéaire".

La figure IV.3 représente les résultats obtenus avec l'algorithme ALG1, 2ème schéma, toujours pour g=2. La numérotation des courbes a la même signification que dans les paragraphes précédents.

Les figures IV.4 et IV.5 donnent les résultats pour g=5 avec les algorithmes ALG2 1er schéma et ALG1 2ème schéma respectivement.

Les figures IV.6 et IV.7 sont relatives à g=8 ; et toujours respectivement aux algorithmes ALG2 1er schéma et ALG1 2ème schéma. Dans le cas de l'algorithme ALG1 2ème schéma, si on résout correctement le problème de point selle (courbe 2), il faut toujours plus de 400 itérations, c'est pourquoi cette courbe n'apparaît pas sur la figure.

Sur les figures IV.8 à IV.16, on représente les équipotentielles de la solution u, calculée par les différents algorithmes. La valeur des paramètres est apparente sur chaque figure. Le nombre 0.10 E-4 en bas à gauche indique que le test d'arrêt retenu pour la convergence des algorithmes (1.28), (1.30), (1.35) est avec $\epsilon = 10^{-5}$, le nombre précédent qui est soit 0, soit -1, soit 1 indique

(4.15) $$ \begin{cases} 0 & : \text{algorithme ALG2, 1er schéma} \\ -1 & : \text{algorithme ALG1, 2ème schéma tel qu'il est présenté en (1.25)-(1.30)} \\ 1 & : \text{algorithme ALG1, 2ème schéma évitant le blocage intempestif.} \end{cases} $$
Pour g=5, la solution obtenue par les trois méthodes précédentes est la même.
Pour g=2 et g=8, la solution obtenue pour l'algorithme ALG1, 2ème schéma,
avec phénomène de blocage possible, est légèrement différente des autres.
Remarquons aussi que dans le cas g=5, on doit obtenir pour la zone rigide
un cercle de centre O et de rayon 0.5, or la triangulation est particuliè-
rement bien adaptée puisqu'elle s'appuie sur ce cercle de rayon 0.5 (voir
fig. II.13).

Sur les figures IV.17 à IV.19, nous avons représenté par des hachures les zones rigides. Elles ont été obtenues par résolution du problème avec l'algorithme ALG2, 1er schéma et de la façon suivante

Si $|\hat{p}_h| \le 10^{-2}$ sur le triangle T considéré, alors ce triangle est dans la zone rigide.

Sur la figure IV.17, correspondant à g=2, on constate que l'on obtient un
carré pour la zone rigide, au lieu du cercle de rayon 0.2.

Celà est évidemment dû à la triangulation et il suffit de regarder la figure
(II.13) pour découvrir qu'il est difficile d'obtenir autre chose. Par contre
sur les figures IV.18, IV.19 les résultats obtenus sont plus conformes à la
réalité.

Le phénomène que l'on constate sur la figure IV.17 pour g=2, n'est pas lié
à la méthode de résolution mais plutôt à la méthode d'approximation ou plus
simplement à la triangulation ; ce phénomène est général pour la résolution
des inéquations variationnelles, et nous verrons qu'il est encore plus crucial
pour le problème de l'élasto-plasticité.
<figure><img src="image_75.png" /><figcaption>Fig. IV.2</figcaption></figure>

<figure><img src="image_76.png" /><figcaption>Fig. IV.3</figcaption></figure>

<figure><img src="image_77.png" /><figcaption>Fig. IV.4</figcaption></figure>

<figure><img src="image_78.png" /><figcaption>Fig. IV.5</figcaption></figure>

<figure><img src="image_79.png" /><figcaption>Fig. IV.6</figcaption></figure>

<figure><img src="image_80.png" /><figcaption>Fig. IV.7</figcaption></figure>
<figure><img src="image_81.png" /><figcaption>Fig. IV.8</figcaption></figure>

<figure><img src="image_82.png" /><figcaption>Fig. IV.9</figcaption></figure>

<figure><img src="image_83.png" /><figcaption>Fig. IV.10</figcaption></figure>

<figure><img src="image_84.png" /><figcaption>Fig. IV.11</figcaption></figure>
LIGNES DE NIVEAU DE LA SOLUTION

<figure><img src="image_85.png" /><figcaption>Fig. IV.12</figcaption></figure>

<figure><img src="image_86.png" /><figcaption>Fig. IV.13</figcaption></figure>

<figure><img src="image_87.png" /><figcaption>Fig. IV.14</figcaption></figure>

<figure><img src="image_88.png" /><figcaption>Fig. IV.15</figcaption></figure>
<figure><img src="image_89.png" /><figcaption>Fig. IV.16
Lignes de niveau de la solution</figcaption></figure>

<figure><img src="image_90.png" /><figcaption>Fig. IV.17
Visualisation des zones rigides</figcaption></figure>

<figure><img src="image_91.png" /><figcaption>Fig. IV.18
Visualisation des zones rigides</figcaption></figure>

<figure><img src="image_92.png" /><figcaption>Fig. IV.19
Visualisation des zones rigides</figcaption></figure>
## V. EXEMPLE 4 - PROBLEME DE LA TORSION ELASTO-PLASTIQUE

Nous résolvons ici numériquement le problème de la torsion élasto-plastique d'une barre de section Ω ; ce problème classique a déjà été traité par diverses méthodes dans [10].

### V.1. Présentation du problème.

Soit Ω ouvert borné de R², de frontière Γ régulière ; on considère l'espace de Sobolev $H_0^1(\Omega)$

$$ (5.1) \qquad H_0^1(\Omega) = \{v | v \in L^2(\Omega), \frac{\partial v}{\partial x_i} \in L^2(\Omega), v|_{\Gamma} = 0\} $$

On définit K convexe fermé de $H_0^1(\Omega)$ par

$$ (5.2) \qquad K = \{v | v \in H_0^1(\Omega), |\nabla v| \le 1 \text{ p.p.}\} $$

On considère la fonctionnelle suivante sur $H_0^1(\Omega)$

$$ (5.3) \qquad J(v) = \frac{1}{2} \int_{\Omega} |\nabla v|^2 dx - \int_{\Omega} f \cdot v dx $$

avec $f \in L^2(\Omega)$.

J est continue, strictement convexe, différentiable, et telle que

$$ \limite \quad J(v) = \infty \\ \|v\| \qquad H_0^1(\Omega) \to \infty $$

Le problème suivant

$$ (5.4) \qquad \begin{cases} \text{trouver } u \in K \text{ tel que} \\ J(u) \le J(v) \quad \forall v \in K \end{cases} $$

possède donc une solution unique caractérisée par l'inéquation variationnelle

$$ (5.5) \qquad \begin{cases} \int_{\Omega} \nabla u \cdot \nabla (v-u) dx - \int_{\Omega} f \cdot (v-u) dx \ge 0 \\ \forall v \in K, u \in K \end{cases} $$
La motivation physique du problème est la suivante : considérons une barre cylindrique de section $\Omega$, infiniment longue (on suppose pour l'instant simplement convexe). Supposons que la barre soit composée d'un matériau isotrope, élastique parfaitement plastique, (avec champ de plasticité donné par le critère de Von Mises). Partant d'un état initial sans contraintes, un moment de torsion croissant est appliqué à la barre. La torsion est caractérisée par $C$ qui est l'angle de torsion par unité de longueur. Alors pour tout $C$, il résulte du principe de Haar-Karman, que la détermination du champ de contraintes est équivalent, (dans un système d'unités physiques adéquat), à la résolution du problème variationnel suivant.

$$ (5.6) \quad \min_{v \in K} \left[ \frac{1}{2} \int_{\Omega} |\nabla v|^2 dx - C \int_{\Omega} v dx \right] $$

Le vecteur contraintes $\sigma$ dans une section est lié à u par $\sigma = \nabla u$, si bien que u est un potentiel des contraintes ; on pourra donc déduire $\sigma$ quand la solution u de (5.4) ou (5.6) est connue.

### Remarque 5.1

Dans le cas d'une section multi connexe (voir fig. V.1) le problème de torsion doit être modifié (cf. [10-11])

<figure><img src="image_93.png" /><figcaption>Figure VI.1</figcaption></figure>

Nous noterons $\Omega$ la réunion de la section $\Omega^*$ et des trous $w_i$, $i=1,2,...,p$. et $\Gamma$ sera la frontière $\partial\Omega$. Le problème de la torsion élasto-plastique dans le cas de section multiconnexe peut-être formulé ainsi

$$ (5.7) \quad \begin{cases} \text{trouver } u \in \tilde{K} \text{ tel que} \\ J(u) \le J(v) \quad \forall v \in \tilde{K} \end{cases} $$
avec $\tilde{K}$ défini par

$$
(5.8) \quad \tilde{K} = \{ v | v \in H_o^1(\Omega), |\nabla v| \le 1 \text{ p.p., } \nabla v = 0 \text{ sur } \omega_i, i=1,2,\dots,p \} .
$$

V.2. Formulation discrète en lagrangien augmenté

Le lagrangien augmenté pour le problème de la torsion élasto-plastique s'écrit

$$
(5.9) \quad \mathcal{L}_r(w_h, q_h, \mu_h) = \frac{1}{2} \int_{\Omega} |q_h|^2 dx - C \int_{\Omega} v_h dx + \frac{r}{2} \int_{\Omega} |\nabla v_h - q_h|^2 dx + \int_{\Omega} \mu_h \cdot (\nabla v_h - q_h) dx
$$
$\mathcal{L}_r$ étant cette fois-ci une application de

$$
V_h \times K_h \times L_h \rightarrow \mathbb{R}
$$

avec $K_h$ convexe fermé de $L_h$ défini par

$$
(5.10) \quad K_h = \{ q_h \mid q_h \in L_h, |q_h| \le 1 \text{ p.p.} \}
$$

dans le cas d'une section connexe et par

$$
(5.11) \quad K_h = \{ q_h \mid q_h \in L_h, |q_h| \le 1, |q_h| = 0 \text{ sur } \omega_i \}
$$

dans le cas d'une action multiconnexe.

L'an analogue de (1.11)-(1.12) pour le problème de l'élastoplasticité conduit à
une équation pour u_h et une inéquation pour p_h ; il y a différentiabilité du
Lagrangien mais on optimise sur un convexe, on aboutit donc à une inéquation
variationnelle.

$$
(5.12) \quad \left\{
\begin{array}{l}
\displaystyle \int_{\Omega} \nabla \hat{u}_h \cdot \nabla v_h \, dx = \frac{C}{r} \int_{\Omega} v_h \, dx + \int_{\Omega} (\hat{p}_h - \frac{1}{r} \mu_h) \nabla v_h \, dx \\
\forall v_h \in V_h , \quad \hat{u}_h \in V_h
\end{array}
\right.
$$

$$
(5.13) \quad \left\{
\begin{array}{l}
(1+r) \displaystyle \int_{\Omega} \hat{p}_h (q_h - \hat{p}_h) \, dx - \displaystyle \int_{\Omega} (r \nabla \hat{u}_h + \mu_h) (q_h - \hat{p}_h) \, dx \ge 0 \\
\forall q_h \in K_h , \quad \hat{p}_h \in K_h .
\end{array}
\right.
$$
Comme dans le cas de l'écoulement du fluide de Bingham, on obtient des inéquations indépendantes sur chaque triangle T de $C_h^\alpha$ que l'on sait résoudre expliquement (voir chapitre 3 - 1.34).

Si on résout en $|\hat{p}_h|$ sur le triangle T on obtient

$$ (5.14) \quad |\hat{p}_h| = \min \left( 1, \frac{|r \nabla \hat{u}_h + \mu_h|}{1+r} \right) $$

Si le triangle T est dans un trou (cas de section multiconnexe) on pose alors

$$ (5.15) \quad |\hat{p}_h| = 0 $$

### V.3. Résultats numériques

Pour les essais numériques nous avons pris $\Omega = ]0,1[ \times ]0,1[ \text{ découpé en 384 triangles et 209 noeuds (fig. V.2). Nous avons pris différentes valeurs de C dans (5.6).}$

Les figures V.3 et V.4 montrent le comportement des algorithmes ALG2 et ALG1 pour $C=1$ de l'angle de torsion. Pour cette valeur de C, le problème traité est linéaire, il n'y a pas de zones plastiques. Pour l'algorithme ALG2 le schéma si on prend $r=1$, on obtient la solution dès la première itération (comparer cette courbe avec la courbe de la figure III.6). Pour l'algorithme ALG1 la deuxième schéma, le nombre d'itérations d'UZAWA (itérations en n dans $(1.25)-(1.31)$), reste constamment égal à 3, mais le nombre d'itérations totales lui, augmente avec $r$ (les courbes 1 et 2 sont confondues).

Sur les figures suivantes V.5 à V.10, on constate que le problème devient plus "difficile" à résoudre au fur et à mesure que $C$ augmente. On donne successivement les résultats pour $C=5$, $7.5$, $10$. et à chaque fois pour l'algorithme ALG2 et ALG1 la deuxième schéma.

Le nombre d'itérations minimal passe de 20 à 70 pour ALG2 et de 90 à 195 pour ALG1.
Remarque 5.2

Nous avons résolu le problème de l'élasto-plasticité, cette fois-ci en se fixant les données physiques (nous avons pris C=10), et en raffinant la triangulation pour voir le comportement des algorithmes suivant le nombre de variables du problème. Nous pouvons voir sur les figures V.11, V.12 et V.13 les triangulations $\mathcal{T}_h$ successives que nous avons choisi pour ce test. Les triangulations $\mathcal{T}_h$ comprennent successivement 128, 512 et puis 2048 triangles. Les résultats sont relatifs à l'algorithme ALG2 pour le schéma (1.31)-(1.35), bien entendu la conclusion reste valable pour l'algorithme ALG1. D'après les résultats obtenus sur la figure V.14 il semble que le nombre d'itérations pour résoudre le problème ne dépend pas de la discrétisation. Cela présente un avantage certain du point de vue numérique : étant en présence d'un problème nouveau dont on en peut pas dire grand chose quant au comportement des algorithmes ALG1 et ALG2, il suffit de faire des tests numériques avec une triangulation grossière (petit nombre de variables), pour déterminer avec une bonne précision une valeur optimale des paramètres (en particulier du paramètre de pénalisation).

On donne aussi à titre indicatif, sur la figure V.15, le comportement de l'algorithme ALG2 pour le schéma, par rapport au paramètre de pénalisation $r$, lorsque $\Omega$ est toujours le carré, découpé en 384 triangles (de la figure V.2) mais comportant un trou en forme de losange en son milieu (demi-diagonale : 0.5). Nous avons cette fois-ci un problème d'élasto-plasticité sur un domaine $\Omega$ multi-connexe.

La valeur de l'angle de torsion est $C=5$ ; il faut donc comparer la figure V.15 à la figure V.5. Suivant l'allure des courbes et la position, on constate que la plasticité a augmenté.

Sur les figures suivantes V.16 à V.39 nous donnons des résultats concernant la solution $u$ (lignes de niveau), ainsi que la visualisation des zones plastiques pour diverses géométries et diverses triangulations.

Sur la figure V.16 on voit la solution $u$ que l'on obtient avec la triangulation de la figure V.2 (384 triangles), pour un angle de torsion $C=10$. La valeur du paramètre de pénalisation retenue est $r=4$, (la valeur n'a évidemment pas d'influence sur la solution finale, mais seulement sur le travail qu'il faut fournir pour obtenir cette solution !) et l'algorithme utilisé est ALG2 pour le schéma. Sur la figure suivante V.17 on a représenté les
zones plastiques, déterminées pratiquement de la façon suivante :

$$
\text{Pour chaque } T \in \mathcal{T}_h \text{ on regarde la valeur de } |\hat{p}_h| \text{ (solution de (5.13) lorsque l'algorithme a convergé)}
$$
$$
(5.76) \quad \left\{
\begin{array}{ll}
\bullet & \text{si } |\hat{p}_h| \ge 1-\rho \Rightarrow T \in \text{zone plastique} \\
\bullet & \text{si non à la zone élastique}
\end{array}
\right.
$$

ρ a été choisi égal à 10⁻² sur la plupart des figures, cette valeur apparaît
d'ailleurs sur les figures, ainsi que l'algorithme utilisé (même convention
qu'en 4.15). Les zones plastiques sont très mal dessinées et n'ont pas une
interprétation physique facile. Le résultat obtenu est lié à la triangula-
tion retenue (figure V.2). Si on prend une triangulation plus fine (figure
V.18, le domaine Ω est subdivisé en 1536 triangles), on peut voir (figure
V.20) que les zones plastiques sont mieux dessinées, mais le défaut dû à
la triangulation (mauvaise orientation de certains triangles), subsiste avec
un effet moindre il est vrai. Les lignes de niveau de la solution u (figure
V.19) obtenu avec la triangulation de la figure V.18 sont beaucoup plus
régulières que celle de la figure V.16. Si on choisit une triangulation
adaptée au problème, (par exemple la triangulation donnée sur la figure V.13),
les zones plastiques sont beaucoup mieux dessinées et déterminées avec préci-
sion (voir figure V.22), les lignes de niveau de la solution u sont aussi
très régulières (figure V.21).

Nous donnons maintenant quelques résultats obtenus pour des domaines non
simplement connexes. La figure V.23 donne le découpage du domaine Ω pouvant
laisser apparaître un trou dont la forme sera indiquée sur les figures suivantes.
Cette triangulation comprend 1536 triangles et 809 noeuds. La figure V.24 repré-
sente les lignes de niveau de la solution u obtenue dans le domaine multiconnexe,
pour un angle de torsion C=5. L'interprétation de cette solution en zones plas-
tiques et élastiques (par 5.16) est donnée sur le figure V.25.

Les figures V.26 et V.27 donnent des résultats analogues, avec le même angle
de torsion C=5, mais la forme du trou ayant changé dans la section Ω.
Les résultats des figures V.28 à V.33 sont toujours relatives à la triangulation de la figure V.23.

Si on augmente la valeur de C, les zones plastiques vont augmenter en "surface". La figure V.28 donne les lignes de niveau de la solution u lorsque l'angle de torsion C=10, le trou dans la section Ω ayant une forme rectangulaire, alors que sur la figure V.30, le trou a une forme différente. Les figures V.29 et V.31 donnent l'interprétation en zones plastiques et élastiques toujours par (5.16). Lorsqu'on augmente la valeur de l'angle de torsion C, les zones plastiques sont plus difficiles à déterminer du point de vue numérique, mais cela reste dû essentiellement à la triangulation suivant ou ne suivant pas la surface libre du problème continu.

Les figures V.32 et V.33 donnent les résultats relatifs à un angle de torsion C=15., la figure V.32 donnant les lignes de niveau de u et la figure V.33 l'interprétation en zones plastiques et élastiques.

La détermination pratique des zones plastiques et élastiques devient encore plus cruciale si on choisit un domaine tel que celui représenté sur la figure V.34. Le domaine Ω est un cercle, contenant deux trous circulaires. Le découpage est fait en 1024 triangles et 545 noeuds. La surface u (fonction potentielle) est moins sympathique que dans les cas précédents. Pourtant si on regarde les lignes de niveau de la solution u (fig. V.35 et V.38), elles ont l'air assez régulières.

Si on détermine les zones plastiques comme précédemment, c'est-à-dire, triangle par triangle en décidant que le triangle est plastique si $|\hat{p}_h| \ge 1-\rho$ (avec $\rho=10^{-2}$), on obtient comme résultat ce qui est représenté sur la figure V.36. On obtient une sorte de damier ce qui n'est pas physique du tout. L'angle de torsion C vaut 10 pour ce premier exemple. Si on est un peu moins rigoureux sur le test de plasticité (si on prend $\rho = 5.10^{-2}$), on obtient les zones données par la figure V.37 (les zones plastiques ainsi obtenues sont un peu plus connexes).

Les figures V.38 et V.39 donnent les résultats obtenus pour un angle de torsion plus faible C=2.5. La détermination des zones plastiques est faite avec $\rho = 5.10^{-2}$.
<figure><img src="image_94.png" /><figcaption>Fig. V.2 - Maillage de Ω (0)</figcaption></figure>

<figure><img src="image_95.png" /><figcaption>Fig. V.3 - Comportement de l'algorithme en fonction du paramètre de pénalisation r.</figcaption></figure>

<figure><img src="image_96.png" /><figcaption>Fig. V.4 - Comportement de l'algorithme en fonction du paramètre de pénalisation r.</figcaption></figure>

<figure><img src="image_97.png" /><figcaption>Fig. V.5 - Comportement de l'algorithme en fonction du paramètre de pénalisation r.</figcaption></figure>
<figure><img src="image_98.png" /><figcaption>Fig. V.6
Comportement de l'algorithme en fonction du paramètre de pénalisation r</figcaption></figure>

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

<figure><img src="image_100.png" /><figcaption>Fig. V.8
Comportement de l'algorithme en fonction du paramètre de pénalisation r</figcaption></figure>

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

<figure><img src="image_102.png" /><figcaption>Fig. V.9
Comportement de l'algorithme en fonction du paramètre de pénalisation r</figcaption></figure>
<figure><img src="image_103.png" /><figcaption>Fig. V.11
Maillage 1 de Ω</figcaption></figure>

<figure><img src="image_104.png" /><figcaption>Fig. V.12
Maillage 2 de Ω</figcaption></figure>

<figure><img src="image_105.png" /><figcaption>Fig. V.13
Maillage 3 de Ω</figcaption></figure>
<figure><img src="image_106.png" /><figcaption>Fig. V.14</figcaption></figure>

Courbes de comportement de l'algorithme en fonction du paramètre de pénalisation r pour les 3 maillages des fig. V.11, V.12, V.13.

<figure><img src="image_107.png" /><figcaption>Fig. V.15</figcaption></figure>

Comportement de l'algorithme en fonction du paramètre de pénalisation r (domaine multi-connexe).

<figure><img src="image_108.png" /><figcaption>Fig. V.16</figcaption></figure>

Lignes de niveau de la solution

<figure><img src="image_109.png" /><figcaption>Fig. V.17</figcaption></figure>

Visualisation des zones plastiques (test $10^{-2}$)
<figure><img src="image_110.png" /><figcaption>Fig. V.18</figcaption></figure>

Raffinement du maillage O de la fig. V.2

<figure><img src="image_111.png" /><figcaption>Fig. V.19</figcaption></figure>

Lignes de niveau de la solution

<figure><img src="image_112.png" /><figcaption>Fig. V.20</figcaption></figure>

Visualisation des zones plastiques
(test $10^{-2}$)

<figure><img src="image_113.png" /><figcaption>Fig. V.21</figcaption></figure>

Lignes de niveau de la solution obtenue
avec maillage figure V.13
<figure><img src="image_114.png" /><figcaption>Fig. V.22
Visualisation des zones plastiques
(test $10^{-2}$)</figcaption></figure>

<figure><img src="image_115.png" /><figcaption>Fig. V.23
Maillage de Ω pour étude de section multi-connexe</figcaption></figure>

<figure><img src="image_116.png" /><figcaption>Fig. V.24
Lignes de niveau de la solution</figcaption></figure>

<figure><img src="image_117.png" /><figcaption>Fig. V.25
Visualisation des zones plastiques
(test $10^{-2}$)</figcaption></figure>
<figure><img src="image_118.png" /><figcaption>Fig. V.26
Lignes de niveau de la solution</figcaption></figure>

<figure><img src="image_119.png" /><figcaption>Fig. 5.27
Visualisation des zones plastiques
(test $10^{-2}$)</figcaption></figure>

<figure><img src="image_120.png" /><figcaption>Fig. V.28
Lignes de niveau de la solution</figcaption></figure>

<figure><img src="image_121.png" /><figcaption>Fig. V.29
Visualisation des zones plastiques
(test $10^{-2}$)</figcaption></figure>
<figure><img src="image_122.png" /><figcaption>Fig. V.30
Lignes de niveau de la solution</figcaption></figure>

<figure><img src="image_123.png" /><figcaption>Fig. V.31
Visualisation des zones plastiques
(test $10^{-2}$)</figcaption></figure>

<figure><img src="image_124.png" /><figcaption>Fig. 5.32
Lignes de niveau de la solution</figcaption></figure>

<figure><img src="image_125.png" /><figcaption>Fig. 5.33
Visualisation des zones plastiques</figcaption></figure>
<figure><img src="image_126.png" /><figcaption>Fig. V.34
Maillage d'un cercle - étude section multiconnexe</figcaption></figure>

<figure><img src="image_127.png" /><figcaption>Fig. V.35
Lignes de niveau de la solution</figcaption></figure>

<figure><img src="image_128.png" /><figcaption>Fig. 5.36
Visualisation des zones plastiques
(test $10^{-2}$)</figcaption></figure>

<figure><img src="image_129.png" /><figcaption>Fig. 5.37
Visualisation des zones plastiques
(test $5 \cdot 10^{-2}$)</figcaption></figure>
<figure><img src="image_130.png" /></figure>

PAS ENTRE 2 LIGNES 0.261E-1
VALEURS EXTREMES DE LA FONCTION
0.418E+0
PARAMETRES
C 0.250E+1
R 0.700E+1

0 0.10E-4

LABORIA- A.MARROCCO

<caption>Fig. V. 38</caption>

Lignes de niveau de la solution

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

ZONES PLASTIQUES

PARAMETRES

C 0.250E+1
R 0.700E+1

0 0.50E-1

LABORIA- A.MARROCCO

<caption>Fig.V.39</caption>

Visualisation des zones plastiques
(test 5.10⁻²)
# VI - EXEMPLE 5 - PROBLEME DE L'ECOULEMENT DANS UNE TUYERE

## VI.1. Présentation du problème.

Si on se réfère à [12] [13] par exemple, l'équation aux dérivées partielles régissant un écoulement transsonique est donnée par

$$ (6.1) \qquad \nabla \cdot (\rho \nabla \phi) = 0 $$

Dans (6.1), $\phi$ est le potentiel des vitesses et la vitesse de l'écoulement $u$ est donnée par $u = \nabla\phi$. Dans (6.1) la densité volumique $\rho$ est une fonction du carré de la vitesse $u$ (donc de $|\nabla\phi|^2$ - comparer à (3.7) ou (3.9)), donnée explicitement par

$$ (6.2) \qquad \rho(x) = \rho_o(1 - kx)^{2.5} \quad , \rho_o \text{ et } k \text{ constantes} \quad , \rho_o = 1 \text{ et } k=6. $$

On va résoudre l'équation (6.1) dans une tuyère convergente-divergente (domaine $\Omega \subset \mathbb{R}^2$) avec des conditions aux limites appropriées, essentiellement une différence de potentiel entre l'entrée et la sortie de la tuyère et une condition de glissement sur les bords de la tuyère $u \cdot n = 0$.

Comme nous venons de le remarquer, formellement le problème de l'écoulement transsonique, possède la même formulation que le problème de magnétostatique (paragraphe III, exemple 2), on peut donc penser à introduire une fonction potentielle telle que $\psi$ en (3.10), on aurait donc une fonctionnelle $\mathfrak{F}$ sur un espace V (analogue à $H_o^1(\Omega)$), mais malheureusement cette fonctionnelle ne possède pas les mêmes propriétés que (3.11), en particulier elle n'est pas strictement convexe, et on ne pourra rien dire quant à l'unicité de la solution du problème. Il est en fait bien connu, que pour les écoulements subsoniques (qui correspondront à une différence de potentiel pas trop élevée), le problème est bien posé et la solution est unique, mais si on augmente la différence de potentiel l'écoulement devient transsonique et il faut imposer une condition supplémentaire (sur la variation de l'entropie) pour assurer l'unicité de la solution et obtenir en fait la solution qui corresponde au problème physique.
Dans les expériences numériques qui sont présentées on ne prend pas en compte la condition variation imposée de l'entropie ; on assurera l'unicité de la solution calculée en imposant la condition

$$ (6.3) \qquad |\nabla \phi| \le 1 $$

Ainsi dans le cas de l'écoulement subsonique ($|\nabla \phi| < 1$) on résoudra effective-ment le problème, alors que dans le cas où la contrainte (6.3) est satutée, on résout un autre problème qui peut ressembler à une sorte d'élasto-plasti-cité (voir § 5, exemple 4).

## VI.2. Résolution par lagrangien augmenté.

Si on procède comme dans le cas du problème de magnétostatique le système (1.11), (1.12) que l'on doit résoudre s'écrit

$$ (6.4) \qquad \begin{cases} \displaystyle \int_{\Omega} \nabla \phi_h \cdot \nabla v_h \, dx = \int_{\Omega} (\hat{p}_h - \frac{1}{r} \mu_h) \cdot \nabla v_h \, dx \\ \forall v_h \in V_h, \quad \phi_h \in V_h \end{cases} $$

$$ (6.5) \qquad \begin{cases} \displaystyle \int_{\Omega} [ \rho (|\hat{p}_h|^2) \hat{p}_h + r \hat{p}_h ] \cdot q_h \, dx = \int_{\Omega} (r \nabla \phi_h + \mu_h) \cdot q_h \, dx \\ \forall q_h \in L_h, \quad \hat{p}_h \in L_h. \end{cases} $$

Si nous explicitons le système (6.5) on peut se ramener pour chaque triangle $T \in \mathcal{T}_h$ à l'équation suivante

$$ (6.6) \qquad g(q) = (\rho(q^2)+r)q = |r \nabla \phi_h + \mu_h| $$

On donne sur la figure VI.1 la variation de la fonction $g(q)$ pour différentes valeurs de $r$ (les valeurs de $r$ retenues sont $10^{-5}, 10^{-2}, 10^{-1}, .25, .5, .75, 1, 1.5, 2, 3, 4, 5, 7, 10$).

Nous voyons que suivant la valeur de $r$, l'allure de la courbe $g(q)$ est dif-férente, (on pourra voir aussi les courbes dérivées $g'(q)$ sur la figure VI.2). Si $r$ est supérieur (approximativement) à 0.7, la fonction $g(q)$ ayant une dérivée positive, est croissante, si $r$ est plus petit que cette valeur, il risque d'y avoir 2 solutions pour l'équation $g(q) = \text{const.}$ dans l'intervalle considéré.
Remarques 6.1

. Les courbes g(q) et les dérivées ont été représentées pour q ∈ [0,2.44],
la raison est tout simplement que pour {valeur environ q > 2.45} ρ(q²) et
donc g(q) n'est plus défini.

. Les courbes g(q) (6.6) présentent un point d'inflexion situé vers q=1.7
(voir aussi figure VI.3 représentant g''(q)).

. la courbe ρ(q) (analogue de la courbe v_r pour le problème de magnétostatique
fig. III.5) est représenté sur la figure VI.4.

. Au vu de l'allure des courbes g(q) (figure VI.1), il est conseillé de prendre
certaines précautions dans l'utilisation de l'algorithme de Newton pour résoudre
l'équation (6.6). Nous avons adopté la méthode suivante pour résoudre (6.6).

1) On se donne a priori une valeur maximum pour q : q<sub>max</sub>

2) On calcule g(q<sub>max</sub>)

3) a) si $|r \nabla \phi_h + \mu_h| \ge g(q_{max}) \Rightarrow$ la solution de (6.6) est prise $q = q_{max}$.

(Ce n'est donc pas une équation que l'on résout mais une inéquation, et il
faudrait remplacer $L_h$ par un $K_h$ du type (5.10) dans (6.5) et écrire en fait
une inégalité variationnelle).

3) b) si $|r \nabla \phi_h + \mu_h| < g(q_{max})$ on résout (6.6) par la méthode de Newton classique.

Ainsi si on prend $q_{max}=1$, le problème traité aura toujours une solution unique, quelle que soit la différence de potentiel appliquée aux bornes de la tuyère, et aussi quelle que soit la valeur du paramètre r. Si la valeur maximum $q_{max}=1$ n'est pas atteinte (différence de potentiel assez faible), on résout vraiment le problème de l'écoulement subsonique dans la tuyère. Lorsque cette borne est effectivement atteinte et qu'il y a saturation on résout un autre problème ayant ou n'ayant pas d'interprétation physique.

Si on prend r assez grand, la courbe g(q) est croissante pour q ∈ [0, q&#773;], avec
q&#773; > 1 et on pourrait espérer résoudre le problème de l'écoulement transsonique
par cette méthode, mais hélas ce n'est pas le cas.
Dans l'exemple numérique traité, la vitesse limite de l'écoulement subsonique, $q=1$ est atteinte lorsqu'on applique une différence de potentiel d'au moins 0.85. Maintenant si on applique une différence de potentiel de 0.865 qui doit certainement entraîner un écoulement transsonique, si on prend comme borne des vitesse $q_{max} = 1.35$, et si on choisit r assez grand (pour pouvoir résoudre (6.6)), $r > 0.5$, après maintes oscillations, l'algorithme utilisé semble converger, mais la solution obtenue dépend apparemment de r (cela est-il dû au fait que théoriquement il y a une infinité de solutions, et on atteindrait l'une ou l'autre suivant le valeur de r ?).

Si par contre la saturation de la contrainte $q = q_{max}$ est assez forte, la solution obtenue semble être unique, ainsi par exemples pour une différence de potentiel de 0.9, et pour $q_{max} = 1.35$ la solution que l'on obtient ne semble pas dépendre de r. La solution obtenue numériquement peut présenter des chocs, mais comme la solution calculée est symétrique, elle ne peut évidemment pas être physique (il ne faut pas oublier qu'on ne prend pas en compte la condition sur l'entropie).

## VI.3 Résultats numériques.

La tuyère étant symétrique, on peut faire les calculs sur une demi-tuyère. Cette demi-tuyère est découpée en 490 triangles et il y a 300 noeuds (voir figure VI 5).

Sur la figure VI.6, nous représentons les courbes nombre d'itérations en fonction du paramètre de pénalisation $r$, pour l'algorithme ALG2 ler schéma (1.31), (1.35). Le problème est traité pour une différence de potentiel de 0.7 qui conduit à un écoulement subsonique : $|\nabla\phi| < 1$ partout.

On a constaté numériquement que si on prenait r trop petit alors l'algorithme ne convergeait plus, en ce sens que l'algorithme stoppait dès les premières itérations, le test (1.35) avec $\varepsilon = 10^{-5}$ étant vérifié, sans pour cela que le $\phi$ calculé soit solution du problème. La valeur optimale du paramètre de pénalisation se situe entre 0.6 et 0.8.

Si nous prenons maintenant une différence de potentiel de 0.82, (nous résolvons le problème toujours avec le même algorithme ALG2 ler schéma), nous obtenons les résultats présentés sur la figure VI.7.
La valeur optimale du paramètre de pénalisation se situe vers 0.5, la forme
de la courbe 4 semble plus aiguë, ce qui signifie que le problème "devient
plus difficile à résoudre". On constate le même phénomène de non convergence
si r est trop petit (mais peut-être est-ce simplement le test $\epsilon = 10^{-5}$ qui
est trop faible dans 1.35).

Dans ces deux premiers exemples on a pris $q_{max} = 1$, pour l'algorithme, mais comme nous l'avons déjà dit cette valeur n'est pas atteinte avec la solution du problème.

Le troisième exemple (fig. VI.8) correspond à une différence de potentiel de 0.85. On est à la limite de l'écoulement subsonique ; en deux endroits de la tuyère, on obtient $q = 1.01$. Nous avons traité ce problème avec $q_{max} = 1.25$, et donc les parties des courbes situées à gauche de r = 0,3 n'a pas grande signification (cela est lié au mode de résolution indiqué dans la remarque 6.1). La valeur optimale du paramètre de pénalisation rendant minimum le nombre d'itérations est voisine de 0.4 (donc très près de la limite de validité de l'algorithme : r = 0.3).

Sur la figure (VI.9), on donne les résultats obtenus pour une différence de potentiel de 0.865, et pour $q_{max} = 1$. On est donc grâce à cette contrainte assuré de l'unicité de la solution calculée. La valeur optimale du paramètre de pénalisation r se situe vers r = 0.7.

La figure suivante VI.10 donne le comportement de l'algorithme toujours pour une différence de potentiel de 0.865 mais pour $q_{max} = 1.35$. La "zone de validité" de la méthode est située à droite de r=0.5, mais même pour r > 0.5, la solution $\phi$ calculée n'est pas toujours la même et dépend apparemment de r ; les courbes représentées ne sont pas régulières (comparées à celles des figures précédentes), mais comme nous n'avons pris que des valeurs discrètes de r et interpolé ensuite, il est possible que ces courbes soient encore plus perturbées que celà !

Les figures VI.11, VI.12, VI.13 donnent le comportement de l'algorithme ALG2 ler schéma, pour une différence de potentiel de 0.9 et respectivement pour $q_{max} = 1., 1.25, 1.35$.
Pour $q_{max} = 1$, la valeur optimale de r semble se situer entre 1 et 2, alors que pour $q_{max} = 1.25$ ou $q_{max} = 1.35$, elle est plus petite que 1 et voisine de la limite de validité.

On a constaté le même phénomène de "fausse convergence" si le paramètre de pénalisation r est trop petit. Nous avons aussi fait quelques expériences numériques avec l'algorithmme ALG1 2ème schéma (1.25)-(1.30).

La figure VI.14 donne le comportement de cet algorithmme dans le cas d'écoulement subsonique (on a pris une différence de potentiel de 0.7) ; la valeur optimale du paramètre de pénalisation se situe vers 1, mais si on compare le nombre d'itérations avec celui de l'algorithmme ALG2 ler schéma, fig. VI.6, aucun doute, ce dernier algorithmme est plus performant.

On peut remarquer que le nombre d'itérations d'UZAWA diminue lorsque le paramètre de pénalisation r augmente, mais cela ne suffit pas pour juger de l'efficacité de l'algorithmme.

Les résultats présentés sur la figure VI.15, correspondent à une différence de potentiel de 0.85 et pour $q_{max} = 1.25$ (résultats à comparer à ceux de la figure VI.8).

Les résultats présentés sur la figure VI.16 correspondent à une différence de potentiel de 0.9 et pour $q_{max} = 1$ (résultats à comparer à ceux de la figure VI.11), et sur la figure VI.17, on a prix $q_{max} = 1.35$ (voir aussi fig. VI.13).

Sur ces trois dernières figures, le nombre d'itérations a considérablement augmenté, et les courbes ne sont pas toujours très régulières.

Sur les figures VI.18 à VI.21, nous présentons quelques résultats graphiques concernant l'écoulement. Sur la figure VI.18 a, on a représenté les équipoientielles de la solution calculée $\phi$. La différence de potentiel imposée est 0.7, le paramètre de pénalisation $r=1$, le paramètre $q_{max}$ (qui n'intervient pas ici) est pris à 1. Le pas entre deux équipoientielles consécutives est 0.035.
Sur la figure VI.18 b on a représenté les équipotentielles de $|∇φ|$, pour cela nous avons procédé de la façon suivante : pour chaque noeud de la triangulation $C_h$, on affecte une valeur $|∇φ|$ en faisant la moyenne des valeurs de $|∇φ|$ sur les triangles ayant ce noeud pour sommet ; on interpole ensuite linéairement la fonction ainsi obtenue et on trace les équipotentielles. Nous obtenons ainsi 0.209 pour la valeur minimum de la fonction, et 0.733 pour la valeur maximum ; l'écart entre 2 lignes consécutives est 0.0262.

Sur la figure VI.18 c, par le même procédé que précédemment nous avons tracé les équipotentielles de $|λ|$ (multiplicateur de Lagrange). Les valeurs extrêmes de la fonction sont 0.206 et 0.580 ; l'écart entre deux lignes consécutives est 0.0187.

La figure V.19 se rapporte aux résultats relatifs à une différence de potentiel de 0.85 (valeur limite) ; le paramètre de pénalisation est $r = 0.8$ et $q_{max}$ retenu pour le calcul est $q_{max} = 1.25$.

La figure VI.19 a, représente les équipotentielles de la fonction $φ_h$ calculée ; le pas entre 2 lignes consécutives est 0.01. La figure VI.19 b représente elle les équipotentielles de $|∇φ_h|$, interprétées comme précédemment. Les valeurs extrêmes de cette fonction sont 0.234 et 1.0. L'écart entre deux lignes consécutives est 0.02.

Les résultats des figures VI.20 et VI.21 correspondent à une différence de potentiel de 0.9. Pour la figure VI.20 on prend $q_{max}=1$ et pour la figure VI.21 on prend $q_{max} = 1.35$.

Sur les figures VI.20 a et VI.21 a on représente les équipotentielles de la fonction $φ_h$ avec un pas entre équipotentielles, respectivement de 0.02 et 0.045.

Les figures VI.20 b et VI.21 b représentent les équipotentielles de la fonction $|∇φ_h|$ ; les valeurs extrêmes de la fonction et le pas sont respectivement 0.248, 1.0, pas 0.02 et 0.222, 1.35, pas 0.0565.

Sur les figures VI.20 c et VI.21 c, on a représenté les zones saturées c'est-à-dire zones où $|∇φ_h| = q_{max}$.
<figure><img src="image_132.png" /><figcaption>Fig. VI.1</figcaption></figure>

Famille de courbes $g(q)$ pour le problème
de l'écoulement dans une tuyère

<figure><img src="image_133.png" /><figcaption>Fig. VI.3</figcaption></figure>

Dérivée seconde de $g(q)$

<figure><img src="image_134.png" /><figcaption>Fig. VI.2</figcaption></figure>

Famille des courbes dérivées

<figure><img src="image_135.png" /><figcaption>Fig. VI.4</figcaption></figure>

Courbe $\rho(q)$ donnée par (6.2)
<figure><img src="image_136.png" /><figcaption>Fig. VI.5
Maillage de la tuyère</figcaption></figure>

<figure><img src="image_137.png" /><figcaption>Fig. VI.6</figcaption></figure>

<figure><img src="image_138.png" /><figcaption>Fig. VI.7</figcaption></figure>

Comportement de l'algorithme en fonction du paramètre de pénalisation r
<figure><img src="image_139.png" /><figcaption>Fig. VI.8</figcaption></figure>

<figure><img src="image_140.png" /><figcaption>Fig. VI.9</figcaption></figure>

<figure><img src="image_141.png" /><figcaption>Fig. VI.10</figcaption></figure>

<figure><img src="image_142.png" /><figcaption>Fig. VI.11</figcaption></figure>

<figure><img src="image_143.png" /><figcaption>Fig. VI.12</figcaption></figure>

<figure><img src="image_144.png" /><figcaption>Fig. VI.13</figcaption></figure>
<figure><img src="image_145.png" /><figcaption>Fig. VI.14</figcaption></figure>

<figure><img src="image_146.png" /><figcaption>Fig. VI.15</figcaption></figure>

<figure><img src="image_147.png" /><figcaption>Fig. VI.16</figcaption></figure>

<figure><img src="image_148.png" /><figcaption>Fig. VI.17</figcaption></figure>

Comportement de l'algorithme en fonction du paramètre de pénalisation r
<figure><img src="image_149.png" /><figcaption>Fig. VI.18</figcaption></figure>

a) équipotentielles de la solution $\Phi$ - ddp = 0,7, pas entre 2 lignes 0.035

b) équipotentielles de $|\nabla\Phi|$, min = 0.209, max = 0.773, pas entre 2 lignes 0.0262

c) équipotentielles de $|\lambda|$, min = 0.206, max = 0.580, pas entre 2 lignes 0.0187

<figure><img src="image_150.png" /><figcaption>Fig. VI.19</figcaption></figure>

a) équipotentielles de la solution $\Phi$, ddp = 0.85, pas entre 2 lignes 0.01

b) équipotentielles de $|\nabla\Phi|$, min = 0.234, max = 1.0, pas entre 2 lignes 0.02.
<figure><img src="image_151.png" /><figcaption>Fig. VI.20</figcaption></figure>

a) équipotentielles de la solution $\Phi$, ddp = 0,9, $q_{max}$ = 1.
pas entre 2 lignes 0.02

b) équipotentielles de $|\nabla\Phi|$, min = 0.248, max = 1.0,
pas entre 2 lignes 0.02

c) visualisation des zones saturées $|\nabla\Phi| = q_{max}$

<figure><img src="image_152.png" /><figcaption>Fig. VI.21</figcaption></figure>

a) équipotentielles de la solution $\Phi$, ddp = 0.9, $q_{max}$ = 1.35
pas entre 2 lignes 0.045

b) équipotentielle de $|\nabla\Phi|$, min = 0.222, max = 1.35,
pas entre 2 lignes 0.0565

c) visualisation des zones saturées $|\nabla\Phi| = q_{max}$
VII - EXEMPLE 6 - PROBLEME DES SURFACES MINIMA

VII.1 Présentation du problème.

Soit Ω un ouvert borné de R<sup>N</sup> (N=2), de frontière Γ régulière. On veut minimiser l'aire d'une hypersurface s'appuyant sur un contour c R<sup>N+1</sup> et définie par

$$
(7.1) \quad \mathcal{C} = \{(x,z) \in \Gamma \times \mathbb{R} \mid z = f(x)\}
$$

où g est une certaine fonction donnée. Alors le problème peut se formuler

$$
(7.2) \quad \left\{
\begin{array}{ll}
\text{trouver } u \in V \text{ solution de} & \\
\ Inf\ J(v) & \\
v \in V & \\
v|_{\Gamma}=f &
\end{array}
\right.
$$

avec J donné par

$$
(7.3) \quad J(v) = 2 \int_{\Omega} \sqrt{1 + |\nabla v|^2} \, dx .
$$

L'espace V est du type $W_o^{1,1}(\Omega)$ (après avoir fait la modification habituelle
due à la présence des conditions aux limites de Dirichlet non homogènes),
cet espace est un Banach non réflexif et on rencontre les mêmes difficultés
théoriques que pour le cas du pseudo-laplacien. On pourra consulter [14] [15]
pour l'existence de solution généralisées pour (7.2). Bien entendu ces
problèmes théoriques disparaissent lorsque nous considérons le problème
discret en remplaçant V par $V_h$.

VII.2. Formulation discrète en Lagrangien augmenté.

Pour le problème des surfaces minima, le lagrangien augmenté s'écrit

$$
(7.4) \quad \mathcal{L}_r(v_h, q_h, \mu_h) = 2 \int_{\Omega} \sqrt{1 + |q_h|^2} \, dx + \frac{r}{2} \int_{\Omega} |\nabla v_h - q_h|^2 \, dx + \int_{\Omega} \mu_h \cdot (\nabla v_h - q_h) \, dx
$$

avec

$$
\mathcal{L}_r : V_h \times L_h \times L_h \to \mathbb{R}
$$

(L'espace $V_h$ tenant implicitement compte de la condition au bord).
Explicitons maintenant le système (1.11)-(1.12) que l'on doit résoudre

$$
(7.5) \quad \left\{
\begin{array}{l}
\displaystyle \int_{\Omega} \nabla \hat{u}_h \cdot \nabla v_h \, dx = \int_{\Omega} (\hat{p}_h - \frac{1}{r} \mu_h) \cdot \nabla v_h \, dx \\
\forall v_h \in V_h, \quad \hat{u}_h \in V_h
\end{array}
\right.
$$

$$
(7.6) \quad \left\{
\begin{array}{l}
\displaystyle \int_{\Omega} \left( \frac{\hat{p}_h}{\sqrt{1+|\hat{p}_h|^2}} + r \hat{p}_h \right) \cdot q_h \, dx = \int_{\Omega} (r \nabla \hat{u}_h + u_h) \cdot q_h \, dx \\
\forall q_h \in L_h, \quad \hat{p}_h \in L_h
\end{array}
\right.
$$

Comme dans les exemples précédents on peut ramener le système (7.6) à un
ensemble de NT equations non linéaires à 1 inconnue (NT est le nombre de
triangles de la triangulation $\mathcal{T}_h$).

L'équation qu'il faut alors résoudre est

$$
(7.7) \quad g(q) = \left( \frac{q}{\sqrt{1+q^2}} + rq \right) = |r \nabla \hat{u}_h + u_h|
$$

L'équation (7.7) est résolue par la méthode de Newton classique. On donne
sur la figure VII.1, la variation de la fonction g(q) (7.7) pour les valeurs
de r suivantes 10⁻², 5·10⁻², 10⁻¹, .5, 1., 10.

La fonction dérivée g'(q) est représentée sur la figure VII.2 pour les mêmes valeurs du paramètre r. La dérivée seconde est elle représentée sur la figure VII.3. La fonction g (7.7) est donc concave et strictement croissante.

VII.3. Résultats numériques.

Le domaine $\Omega$ est constitué par une couronne. Le petit cercle $\rho_1$ de rayon 1.
et le grand cercle $\rho_2$ de rayon 4.

Le domaine est découpé en 704 triangles et il y a 376 noeuds (voir figure VII.4).
Le bord Γ est aussi constitué de deux cercles concentriques. Pour les applications
numériques on a toujours pris f(x) = V1 pour le grand cercle ρ₂ et tout d'abord
f(x) = V2 sur le cercle ρ₁. Dans ce cas précis on connaît une solution analytique
au problème et la valeur critique pour la constante V2 est 2.07 (si la valeur
V1 = 0).
Donnons tout d'abord les résultats numériques concernant l'algorithme ALG2
1er schéma (1.31), (1.35).

Sur la figure VII.5 on donne le nombre d'itérations nécessaires à la conver-
gence de l'algorithme en fonction du paramètre de pénalisation r. La valeur
optimale du paramètre de pénalisation semble être légèrement au dessous de 1.
Il faut moins de 20 itérations pour obtenir la convergence avec $\epsilon = 10^{-5}$
dans 1.35.

Maintenant si on augmente la valeur de V2 (V2 = 2.07), on constate sur la
figure VII.6 que la valeur optimale du paramètre de pénalisation r a diminué,
que par contre le nombre d'itérations nécessaires à la convergence a augmenté
et que les courbes sont "moins plates" que sur la figure précédente : le pro-
blème est plus difficile à résoudre.

Si on augmente encore la valeur de V2 (V2=4.) ces divers phénomènes s'accen-
tuent. On pourra s'en rendre compte en observant les courbes de la figure VII.7.

Les résultats donnés par l'algorithme ALG1 2ème schéma (1.25)-(1.30), sont
nettement moins bons que ceux donnés par ALG2 1er schéma (en temps calcul).
Sur la figure VII.8 on représente le comportement de cet algorithme, la
valeur de la constante V2 est ici 1. Il faut presque 100 itérations (à l'op-
timum, r=1) pour résoudre le problème. Si on augmente la valeur de V2, le
nombre d'itérations croit fortement et par exemple pour V2=4, à l'optimum
il faut plus de 300 itérations.

Les figures qui suivent représentent les solutions calculées pour différents
exemples. Tout d'abord la figure VII.9 représente les lignes de niveau de la
solution u lorsque les paramètres V1 et V2 valent respectivement 0 et 1. ;
sur la figure suivante (VVI.10) nous avons représenté des coupes successives
(équidistantes) de la solution u obtenue.

La figure VII.11 corespond aux paramètres V1=0 et V2=2.07. Les gradients sont
plus forts au voisinage du petit cercle. La figure VII.12 donne les coupes
successives de la solution ainsi obtenue.
La figure VII.13 donne les lignes de niveau de la solution lorsque les paramètres valent $V1=0$ et $V2=4$. Les gradients sont vraiment très forts près du petit cercle. Théoriquement la paroi est verticale, mais cela est évidemment difficile à obtenir numériquement. On représente la solution par coupes équidistantes sur la figure VII.14.

Nous avons aussi fait d'autres expériences numériques, en prenant, toujours la valeur 0 sur le bord du grand cercle et en faisant passer un plan par le petit-cercle, l'équation étant donnée par

$$z = \left(\frac{V2-V1}{2}\right)x_1 + \left(\frac{V2+V1}{2}\right)$$

(V1 et V2 étant les valeurs imprimées sur les figures).

Si on prend $V1=3$. et $V2=1.5$ on obtient les résultats présentés sur les figures VII.15 et VII.16. Pour $V1=4$ et $V2=1$ les résultats sont sur les figures VII.17 et VII.18 et pour $V1=2$, $V2=-1$ les résultats sont sur les figures VII.19, VII.20.
<figure><img src="image_153.png" /><figcaption>Fig. VII.1
Famille de courbes g(q) pour le problème des surfaces minima.</figcaption></figure>

<figure><img src="image_154.png" /><figcaption>Fig. VII.2
Famille de fonctions dérivées</figcaption></figure>

<figure><img src="image_155.png" /><figcaption>Fig. VII.3
Dérivée seconde de g(q)</figcaption></figure>
<figure><img src="image_156.png" /><figcaption>Fig. VII.4
Maillage de Ω</figcaption></figure>

<figure><img src="image_157.png" /><figcaption>Fig. VII.5</figcaption></figure>

<figure><img src="image_158.png" /><figcaption>Fig. VII.6</figcaption></figure>

Comportement de l'algorithme en fonction du paramètre de pénalisation r.
<figure><img src="image_159.png" /><figcaption>Fig. VII.7</figcaption></figure>

<figure><img src="image_160.png" /><figcaption>Fig. VII.8</figcaption></figure>

Comportement de l'algorithme en fonction du paramètre de pénalisation r.

<figure><img src="image_161.png" /><figcaption>Fig. VII.9</figcaption></figure>

Equipotentielles de la solution

<figure><img src="image_162.png" /><figcaption>Fig. VII.10</figcaption></figure>

Coupes de la solution
<figure><img src="image_163.png" /><figcaption>Fig. VII.11
Equipotentielles de la solution</figcaption></figure>

<figure><img src="image_164.png" /><figcaption>Fig. VII.12
Coupes de la solution</figcaption></figure>

<figure><img src="image_165.png" /><figcaption>Fig. VII.13
Equipotentielles de la solution</figcaption></figure>

<figure><img src="image_166.png" /><figcaption>Fig. VII.14
Coupes de la solution</figcaption></figure>
<figure><img src="image_167.png" /><figcaption>Fig. VII.15
Equipotentielles de la solution</figcaption></figure>

<figure><img src="image_168.png" /><figcaption>Fig. VII.16
Coupes de la solution</figcaption></figure>

<figure><img src="image_169.png" /><figcaption>Fig. VII.17
Equipotentielles de la solution</figcaption></figure>

<figure><img src="image_170.png" /><figcaption>Fig. VII.18
Coupes de la solution</figcaption></figure>
<figure><img src="image_171.png" /></figure>

PAS ENTRE 2 LIGNES 0.150E+0

VALEURS EXTREMES DE LA FONCTION
-0.100E+1 0.200E+1

PARAMETRES V1,V2

2.00 -1.00

LABORIA- A.MARROCCO

Fig. VII.19

<figure><img src="image_172.png" /><figcaption>Equipotentielles de la solution</figcaption></figure>

Fig. VII.20

<caption>Coupes de la solution</caption>
VIII. CONCLUSION

La technique du lagrangien augmenté combinée à l'approximation par éléments finis permet de résoudre une assez grande variété de problèmes non-linéaires ; il suffit d'avoir à sa disposition un logiciel de résolution de système linéaire (symétrique). C'est le "gros du travail", mais il est commun à tous les problèmes. Quant à la non linéarité du problème, elle se ramène à la résolution d'une équation à l'inconnue.

On peut aussi envisager des problèmes évolutifs du type parabolique. Si le pas de temps est choisi constant, alors la matrice du problème est aussi constante et il suffit de la factoriser une seule fois. A chaque pas de temps nous avons un problème stationnaire du type de ceux étudiés dans les chapitres précédents.

Des essais numériques ont été faits sur le problème du pseudo laplacien (avec $p=1,1$) et sur le problème du fluide de Bingham.
BIBLIOGRAPHIE

[1] M. FORTIN, R. GLOWINSKI (Ed.), Résolution numérique de problèmes aux limites par des méthodes de Lagrangien Augmenté (à paraître).

[2] R. GLOWINSKI, A. MARROCCO, Sur l'approximation par éléments finis d'ordre 1 et la résolution par pénalisation-dualité d'une classe de problème de Dirichlet non linéaires R.A.I.R.O., R.2, p. 41-76, août 1975.

[3] J.L. LIONS, Contrôle optimal des sytèmes gouvernés par des équations aux dérivées partielles, DUNOD-GAUTHIER VILLARS, Paris, 1968.

[4] R. GLOWINSKI, A. MARROCCO, Sur l'approximation par éléments finis d'ordre 1 et la résolution par pénalisation-dualité d'une classe de problèmes de Dirichlet non linéaires, Rapport Laboria N° 115, avril 1975.

[5] R. GLOWINSKI, A. MARROCCO, Finite element approximation and iterative methods for 2.D non linear magnetostatic problems. Proceedings of COMPUMAG CONFERENCE, Oxford, 1976, pp. 112-125.

[6] R. GLOWINSKI, A. MARROCCO, Numerical solution of two-dimensional magnetostatic problems by augmented lagrangian methods. Computer methods in applied mechanics and engineering, 12 (1977), pp. 33-46.

[7] G. DUVAUT, J.L. LIONS, Les inéquations en mécanique et en physique, DUNOD, Paris (1972).

[8] D. GABAY, B. MERCIER, A dual algorithm for the solution of non linear variational problems via finite element approximation, Rapport Laboria N° 126, juin 1975.

[9] M. GOURSAT, Analyse numérique de problèmes d'élasto-plasticité et de visco-plasticité, Thèse de 3ème cycle, Paris, 1971.
[10] R. GLOWINSKI, J.L. LIONS, R. TREMOLIERES, Analyse numérique des inéquations variationnelles, DUNOD, Paris, 1976.

[11] R. GLOWINSKI, H. LANCHON, Torsion élasto-plastique d'une barre cylindrique de section multiconnexe, Journal de Mécanique, Vol. 12 N° 1, Mars 1973.

[12] O. PIRONNEAU, Sur des problèmes d'optimisation de structure en mécanique des fluides. Thèse Paris VI, 1976.

[13] M.O. BRISTEAU, Application of optimal control theory to transonic flow computations by finite elements methods. A paraître dans les proceedings du 3ème colloque international sur les méthodes de calcul scientifique et technique (Versailles, Dec. 1977).

[14] I. EKELAND, R. TEMAM, Analyse convexe et problèmes variationnels, Dunod, 1971.

[15] A. VANENDE-MORREEUW, Thèse 3ème cycle, Paris 1972.
<figure><img src="image_173.png" /></figure>
TABLE DES MATIERES

<table><thead><tr><th></th><th></th><th>Pages</th></tr></thead><tbody><tr><td>I.</td><td>INTRODUCTION</td><td>2</td></tr><tr><td></td><td>Alg 1</td><td>3</td></tr><tr><td></td><td>Alg 2</td><td>4</td></tr><tr><td>II.</td><td>EXEMPLE 1 - PROBLEME DU PSEUDO-LAPLACIEN</td><td>12</td></tr><tr><td></td><td>II.1 Présentation du problème</td><td>12</td></tr><tr><td></td><td>II.2 Formulation discrète en lagrangien augmenté</td><td>13</td></tr><tr><td></td><td>II.3 Résultats numériques</td><td>15</td></tr><tr><td></td><td>Figures II.1 à II.34</td><td>19</td></tr><tr><td>III.</td><td>EXEMPLE 2 - PROBLEME DE MAGNETOSTATIQUE</td><td>28</td></tr><tr><td></td><td>III.1 Présentation du problème</td><td>28</td></tr><tr><td></td><td>III.2 Approximation de v<sub>r</sub></td><td>29</td></tr><tr><td></td><td>III.3 Formulation du problème</td><td>30</td></tr><tr><td></td><td>III.4 Formulation du problème en lagrangien augmenté</td><td>30</td></tr><tr><td></td><td>III.5 Résultats numériques</td><td>32</td></tr><tr><td></td><td>Figures III.1 à III.39</td><td>40</td></tr><tr><td>IV.</td><td>EXEMPLE 3 - ECOULEMENT D'UN FLUIDE DE BINGHAM</td><td>49</td></tr><tr><td></td><td>IV.1 Présentation du problème</td><td>49</td></tr><tr><td></td><td>IV.2 Formulation discrète en Lagrangien augmenté</td><td>49</td></tr><tr><td></td><td>IV.3 Résultats numériques</td><td>52</td></tr><tr><td></td><td>Figures IV.2 à IV.19</td><td>55</td></tr><tr><td>V.</td><td>EXEMPLE 4 - PROBLEME DE LA TORSION ELASTO-PLASTIQUE</td><td>59</td></tr><tr><td></td><td>V.1. Présentation du problème</td><td>59</td></tr><tr><td></td><td>V.2. Formulation discrète en Lagrangien augmenté</td><td>61</td></tr><tr><td></td><td>V.3. Résultats numériques</td><td>62</td></tr><tr><td></td><td>Figures V.2 à V.38</td><td>66</td></tr><tr><td>VI.</td><td>EXEMPLE 5 - PROBLEME DE L'ECOULEMENT DANS UNE TUYERE</td><td>76</td></tr><tr><td></td><td>VI.1. Présentation du problème</td><td>76</td></tr><tr><td></td><td>VI.2. Résolution par Lagrangien augmenté</td><td>77</td></tr><tr><td></td><td>VI.3. Résultats numériques</td><td>79</td></tr><tr><td></td><td>Figures VI.1 à VI.21</td><td>83</td></tr></tbody></table>
<table><tr><td>VII. EXEMPLE 6 - PROBLEME DES SURFACES MINIMA</td><td>89</td></tr><tr><td>VII.1. Présentation du problème</td><td>89</td></tr><tr><td>VII.2. Formulation discrète en Lagrangien augmenté</td><td>89</td></tr><tr><td>VII.3. Résultats numériques</td><td>90</td></tr><tr><td>Figures VII.1 à VII.20</td><td>93</td></tr><tr><td>VIII. CONCLUSION</td><td>99</td></tr><tr><td>BIBLIOGRAPHIE</td><td>100</td></tr></table>