# Méthode de Newton généralisée en mécanique du contact

Pierre Alart

► To cite this version:

Pierre Alart. Méthode de Newton généralisée en mécanique du contact. Journal Africain de Mathématiques pures et Appliquées, 1997, 76, pp.83-108. 10.1016/S0021-7824(97)89946-1. hal-01352449

HAL Id: hal-01352449

https://hal.science/hal-01352449v1

Submitted on 8 Aug 2016

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 CC BY 4.0 - Attribution - International License
# MÉTHODE DE NEWTON GÉNÉRALISÉE EN MÉCANIQUE DU CONTACT

Par P. ALART

**RÉSUMÉ.** – Cet article présente une étude de la convergence globale de la méthode de Newton généralisée aux systèmes d'équations non différentiables issues de la mécanique du contact. La première partie fournit quelques résultats théoriques pour un problème d'obstacle type. Dans la seconde partie, seule une approche heuristique permet d'appréhender le comportement de l'algorithme de Newton pour des problèmes comportant, outre le contact unilatéral, des non linéarités comme le frottement et l'élasto-plasticité.

**ABSTRACT.** – This paper is devoted to studying the global convergence of the Newton method generalized to systems of non differentiable equations issued from the contact mechanics. In the first part, some theoretical results are presented for a typical obstacle problem. In the second part, the behaviour of the Newton algorithm is analyzed by an heuristic approach. In this way, others non linearities, like friction or elastoplasticity, together with unilateral contact, can be taken up.

## Introduction

Utiliser un algorithme de Newton est un gage de simplicité et, on l'espère, d'efficacité. C'est la méthode de prédilection de l'ingénieur dès qu'il se trouve confronté à un système d'équations non linéaires. Plusieurs variantes sont d'ailleurs souvent disponibles dans les codes standards de calcul des structures. La vitesse quadratique de la convergence au voisinage de la solution est son principal attrait mais, en contrepartie, les instabilités numériques sont fréquentes. La convergence globale n'est garantie que sous des hypothèses très restrictives sur l'opérateur (isotonie, ordre-convexité, faible non linéarité) [27, 28].

Sur les problèmes de contact (liaison unilatérale et frottement éventuel), l'usage de méthodes de type Newton peut sembler hérétique, le système à résoudre étant non différentiable. Ceci conduit d'ailleurs certains utilisateurs, pris de scrupules, à régulariser les équations, ce qui n'élimine pas pour autant les problèmes de stabilité. Mais une attitude plus répandue consiste à identifier différents sous problèmes « classiques », d'optimisation convexe sous contraintes ou de point fixe, disposant chacun de méthodes numériques de résolution dont la convergence est établie. Cette approche possède de nombreuses variantes selon les algorithmes adoptés et leur agencement pour résoudre le problème global [19, 20, 22, 30, 31]. Mais la convergence du schéma global n'est pas assurée, surtout si l'on adjoint au contact et au frottement des non linéarités issues, par exemples, des grandes déformations.
Le parti pris dans mes travaux depuis plusieurs années privilégie l'algorithm de Newton comme la méthode de résolution non linéaire standard. Associée dans un premier temps avec des méthodes de pénalisation pour traiter le contact et le frottement [10, 14], cette approche a donné, en pratique, toute satisfaction, quand elle fut couplée avec une formulation mixte de type Lagrangien augmenté [3], et ce sans recours à des techniques de stabilisation (amortissement, line-search...). Un des intérêts pratiques de la méthode de Newton est qu'une seule boucle d'itérations peut résoudre toutes les non linéarités, et ce, sans modifier la structure d'un logiciel standard d'éléments finis.

L'étude théorique est cependant délicate car peu de travaux concernent la méthode de Newton appliquée à des équations non différentiables (voir [29] dans un contexte différent); la différentiabilité est souvent essentielle dans les résultats de convergence. Pour cette raison, cet article est découpé en deux parties. La première rassemble les résultats théoriques et se limite à un problème d'obstacle type, sur lequel existent de nombreuses autres méthodes de résolution. La méthode de Newton généralisée n'est donc qu'une alternative, avec ses avantages et ses limites.

La seconde partie développe une approche heuristique en se limitant à des exemples élémentaires, par leur faible nombre d'inconnues, que l'on souhaite cependant instructifs sur le comportement de l'algorithm dans le cas général. En contre partie, on considère des problèmes plus complexes que dans la première partie comme ceux avec frottement, ou couplant plusieurs non linéarités (contact unilatéral et frottement, frottement et élasto-plasticité). Cette démarche, à défaut de théorèmes, permet de contrôler aisément la stabilité de la méthode de Newton généralisée par un choix peu contraignant du facteur de pénalisation spécifique des méthodes de pénalité-dualité, seul paramètre à « caler ».

Cet article rassemble un certain nombre de résultats obtenus au cours des dernières années; chacun d'entre eux ne constitue qu'une réponse partielle, mais pris dans leur ensemble, ils fournissent les raisons de la confiance que nous avons dans cette approche, entretenue par ailleurs par son efficacité sur des problèmes industriels [6, 17, 18, 24]. Des techniques multigrilles ou parallèles tout à fait spécifiques sont en cours de développement afin de réduire le temps calcul des algorithmes [7, 8].

# Partie I : Résultats théoriques

## 1. Formulation mixte d'un problème d'obstacle type

On considère un problème d'obstacle abstrait type exprimé en dimension finie, qui consiste à trouver un vecteur réalisant le minimum d'une fonction quadratique sur le cône positif,

$$ (1) \qquad u = \arg \min \left\{ \frac{1}{2} v \cdot A v - v \cdot f ; v \in \mathbb{R}_+^n \right\}. $$

La matrice A est supposée symétrique et définie positive. On se limite au cas discret car l'algorithm de Newton sera directement appliqué à ce type de problème. Cet exemple n'est
pas très représentatif d'un problème de contact, le nombre de contraintes (contacts nodaux) étant égal à la dimension de l'espace de l'inconnue (déplacements). On se reportera à la fin de la Section I.3 pour un commentaire plus détaillé sur ce sujet. En introduisant un multiplicateur, noté $\lambda$, on a conformément à [3] le système suivant à résoudre :

$$ (2) \qquad F(u, \lambda) = 0 \Leftrightarrow \begin{cases} Au + \lambda - f = 0, \\ \lambda - \operatorname{proj}(\lambda + ru; \mathbb{R}_{-}^{n}) = 0. \end{cases} $$

Le facteur $r$ est strictement positif et $\operatorname{proj}(\gamma; \mathbb{R}_{-}^{n})$ désigne la projection sur le cône négatif. Dans [3], le système à résoudre est obtenu en dérivant un lagrangien augmenté. Bien qu'équivalent à (2), il en est différent, les matrices jacobiennes étant symétriques. Mais l'avantage du système (2) réside, comme nous le verrons plus loin, dans la linéarité de la première ligne.

La démarche utilisée depuis plusieurs années pour résoudre ce système consiste à utiliser « brutalement » un algorithme de Newton. L'opérateur $F$ n'est de toute évidence pas différentiable et la méthode de Newton est dite alors généralisée aux systèmes lipschitziens non différentiables. L'intention sous-jacente est d'appliquer cette approche à des problèmes non linéaires comme ceux de contact avec frottement [25], pour lesquels on ne dispose pas de problème d'optimisation associé mais uniquement d'un système comparable à (2) [3].

Pour obtenir les résultats ci-après, on peut exploiter la linéarité de la première équation du système (2) et ainsi éliminer la variable $u$. Le système contenant $2n$ équations à $2n$ inconnues est remplacé par le système suivant à $n$ équations,

$$ (3) \qquad G(\gamma) = (I - r A^{-1})^{-1} \gamma + g - \operatorname{proj}(\gamma; \mathbb{R}_{-}^{n}) = 0, $$

où

$$ g = -r (A - r I)^{-1} f \quad \text{et} \quad \gamma = (I - r A^{-1}) (\lambda - g). $$

L'écriture du système (3) par changement de variables suppose implicitement que $I - r A^{-1}$ est inversible. En fait l'algorithme de Newton peut fonctionner sur (2) même avec une valeur rendant $I - r A^{-1}$ singulière. On peut exclure une telle situation en prenant $r$ différent des valeurs propres de $A^{-1}$ ou plus simplement plus petit que la plus petite valeur propre de $A^{-1}$. Le système (3) peut s'écrire synthétiquement, en séparant les parties linéaire et non différentiable,

$$ (4) \qquad G(\gamma) = B\gamma + \phi(\gamma) = 0. $$

Le changement de variable effectué en (3) rend la partie non différentiable diagonale, propriété utile par la suite.
## 2. Méthode de Newton généralisée

Pour définir rigoureusement l'algorithm de Newton généralisé, il convient d'introduire le jacobien généralisé [13] qui constitue une extension aux opérateurs lipschitziens de la notion de différentielle.

**DÉFINITION.** – Si $F$ est une application lipschitzienne de $\mathbb{R}^n$ dans $\mathbb{R}^n$, le jacobien généralisé, noté $\partial F(x)$ est l'enveloppe convexe de toutes les matrices obtenues comme limite d'une suite de matrices jacobiennes $\nabla F(x^i)$, où $x^i$ tend vers $x$ et $F$ est différentiable en $x^i$ :

$$ (5) \qquad \partial F(x) = \co \delta F(x) = \co \{ \lim \nabla F(x^i); x^i \to x, F \text{ différentiable en } x^i \} $$

L'ensemble $\delta F(x)$ est appelé base du jacobien généralisé et joue un rôle important dans l'étude de l'existence et de l'unicité du problème de contact avec frottement [1, 9]. Une suite d'itérés de la méthode de Newton généralisée pour résoudre $F(x) = 0$ est alors définie de la manière suivante,

$$ (6) \qquad x^{k+1} = x^k - J_k^{-1} F(x^k) \quad \text{avec} \quad J_k \in \partial F(x^k). $$

Cette extension de la méthode de Newton est en définitive assez formelle car, en pratique, les itérés n'ont que peu de chance de se situer en un point de non différentiabilité puisque l'ensemble de tels points est de mesure nulle. Il n'y a donc pas lieu, en général, de préciser le choix de la matrice dans l'ensemble $\partial F(x^k)$. On peut cependant se limiter de façon naturelle à $\delta F(x^k)$ qui est facilement accessible comme nous le verrons plus loin. Cette notion d'algorithm de Newton généralisée diffère de celle introduite par Ciarlet [11], qui regroupe sous ce vocable de nombreuses variantes de la méthode de Newton (Quasi-Newton, Newton modifié...).

A ce stade, il convient de discuter les différents résultats classiques sur la convergence globale de la méthode de Newton. Il est possible d'obtenir la convergence globale pour de faibles non linéarités [27]. Ceci est inutilisable ici où la non linéarité de type projection est sévère. Nous ferons cependant référence au résultat de Necas [27] dans la partie réservée à l'approche heuristique.

Des hypothèses d'ordre partiel conduisent également à la convergence globale monotone [28]. L'exigence de différentiabilité supposée dans [28] est aisément surmontable. La section suivante exploite cette démarche sur des problèmes du type (4).

Les méthodes d'amortissement classiques [28] utilisent une évaluation d'une fonctionnelle à minimiser. Par méthode d'amortissement il faut entendre une technique de réduction de l'incrément dans la direction de Newton, qui assure la convergence en garantissant la décroissance de la fonctionnelle à minimiser. Si le système d'équations dérive d'un potentiel, ce dernier doit être convexe pour que les méthodes de line-search convergent. Sinon on a recours aux algorithmes de Gauss-Newton basés sur un potentiel aux moindres carrés. Dans les deux cas la différentiabilité est essentielle.
Il existe cependant une technique d'amortissement spécifique aux systèmes linéaires par morceaux initialement développée pour des circuits électriques non linéaires [21]. Cette procédure consiste, dans la direction de Newton, à s'arrêter sur l'hyperplan séparant deux zones de linéarité. L'itération suivante utilise, pour calculer la nouvelle direction, la matrice jacobienne de la région adjacente. Si dans l'espace objet le trajet des itérés est une ligne brisée, son image est un segment de droite dans l'espace image, remarque qui constitue l'esprit de la preuve de convergence. La Figure 1 illustre le comportement de l'algorithme sur un opérateur de $\mathbb{R}^2$ dans $\mathbb{R}^2$ linéaire par cônes au nombre de quatre. Contrairement à ce qui est écrit plus haut le choix de la matrice dans $\delta F(x^k)$ est ici essentiel. Il est évident si l'on se trouve à l'interface de deux zones de linéarité car $\delta F(x^k)$ contient deux éléments, celui associé à la région d'où l'on vient, et celui de la zone dans laquelle on entre. Pour un itéré se situant à l'interface de plus de deux régions, on ne peut conclure et la convergence n'est plus nécessairement assurée. Notons enfin que, comme toute méthode d'amortissement, la convergence est garantie au prix d'une stratégie globale assez onéreuse en temps calcul [4].

<figure><img src="image_2.png" /><figcaption>Fig. 1. – Méthode d'amortissement spécifique.</figcaption></figure>

### 3. Un résultat de convergence monotone

La convergence globale peut être également obtenue sous des hypothèses d'ordre partiel dans $\mathbb{R}^n$. Rappelons qu'un vecteur est positif si chacune de ses composantes est positive. Une telle notion présente l'inconvénient de ne pas être invariante par changement de base. Cependant, seule une telle approche fournit des résultats de convergence, en outre monotone, faisant abstraction de la différentiabilité. Le premier théorème est une généralisation immédiate d'un résultat classique [28].

**THÉORÈME 1.** – On considère une application lipschitzienne $F$ de $\mathbb{R}^n$ dans $\mathbb{R}^n$, convexe, telle que

(i) l'équation $F(s) = 0$ admette une unique solution,

(ii) $\forall x \in \mathbb{R}^n, \forall J \in \partial F(x), J^{-1}$ existe et $J^{-1} \ge 0$.
Alors la suite engendrée par (6) converge pour toute valeur initiale $x^0$ et la convergence est monotone ($s \le \dots \le x^k \le \dots \le x^1$).

Preuve. – a) On démontre d'abord que $F(x^1) \ge 0$. En effet la convexité fournit l'inégalité,

$$ (7) \qquad F(x^1) \ge F(x^0) + J(x^1 - x^0), \quad \forall J \in \partial F(x^0). $$

En remplaçant dans (7) $x^1$ par sa valeur ($x^1 = x^0 - J_0^{-1}F(x^0)$) et en prenant $J = J_0$, il vient,

$$ (8) \qquad F(x^1) \ge F(x^0) + J_0 (x^0 - J_0^{-1} F(x^0) - x^0) = 0. $$

b) On vérifie que $x^1 \ge s$. En effet, la convexité entraîne,

$$ (9) \qquad 0 = F(s) \ge F(x^1) + J_1(s - x^1), $$

ou encore, en pré-multipliant par $J_1^{-1}$,

$$ (10) \qquad x^1 \ge s + J_1^{-1} F(x^1). $$

Le second terme du deuxième membre est positif d'après l'hypothèse (ii) et d'après l'inégalité (8).

c) Par récurrence et en procédant comme dans la deuxième étape, la suite $(x^k)$ s'avère décroissante monotone et minorée par $s$,

$$ (11) \qquad x^k \ge x^k - J_k^{-1} F(x^k) = x^{k+1} \ge s. $$

Elle converge donc vers une valeur $x^*$ satisfaisant $F(x^*) = 0$, car la multiapplication $\partial F(\cdot)$ est semi-continue inférieurement [13],

$$ (12) \qquad x^* = x^* - J_*^{-1} F(x^*), \quad J_*^{-1} \in \partial F(x^*). $$

L'hypothèse (i) montre alors que $x^*$ et $s$ sont confondus.

**COROLLAIRE.** – Si $F(x) = Bx + \varphi(x)$ où $B$ est une $M$-matrice et $\varphi$ est diagonal, convexe et isotone ($x \ge y \Rightarrow \varphi(x) \ge \varphi(y)$), la suite $(x^k)$ converge.

Ce dernier résultat n'est malheureusement pas applicable au système (3), car $\varphi$ n'est pas isotone mais antitone ($x \ge y \Rightarrow \varphi(x) \le \varphi(y)$). Pour obtenir la convergence de l'algorithme (6) pour résoudre le système (3), il est donc nécessaire de vérifier directement les hypothèses du théorème 1. Mais il convient au préalable d'exprimer le jacobien
généralisé en tout point $\gamma$ de l'application $G$, définie en (3). Comme $\varphi$ est linéaire par cônes de sommet l'origine, il suffit d'écrire $\partial G(0)$. En effet [1],

$$ (13) \qquad \partial G(0) = \co \delta G(0) = \co \cup \{\delta G(\gamma); \gamma \in \mathbb{R}^n\}. $$

Le jacobien généralisé et sa base ont la forme suivante,

$$ (14) \qquad \partial G(0) = \left\{ J \in \mathbb{R}^{n \times n}; J = (I - r A^{-1})^{-1} - \sum_{i=1}^{n} \alpha_i e_i \otimes e_i, \alpha_i \in [0, 1] \right\}, $$

$$ (15) \qquad \delta G(0) = \left\{ J \in \mathbb{R}^{n \times n}; J = (I - r A^{-1})^{-1} - \sum_{i=1}^{n} \alpha_i e_i \otimes e_i, \alpha_i \in \{0, 1\} \right\}, $$

où $(e_i)_i$ dénote la base canonique de $\mathbb{R}^n$ et $\otimes$ le produit tensoriel ou dyadique. La positivité de $J^{-1}$ est la condition la plus délicate à vérifier. Le lemme suivant fournit un outil technique dans cette voie. Pour une matrice carrée $C$, on considère la famille $(J_p)_{p=0,n}$, définie par,

$$ (16) \qquad J_p = (I - rC)^{-1} - \sum_{i=1}^{p} e_i \otimes e_i. $$

$|C_p|$ désigne le mineur principal d'ordre $p$ défini par :

$$ (17) \qquad |C_p| = \det C_p = \det E_p^T C E_p \quad \text{où } E_p = [e_1, \dots, e_p] \in \mathbb{R}^{n \times p}. $$

LEMME (conjecture). – *Pour les valeurs de $p$ de 1 à $n$,*

$$ (18) \qquad 1 - e_p \cdot J_{p-1}^{-1} e_p = r \frac{|C_p|}{|C_{p-1}|}. $$

Le produit scalaire est noté avec un point. On peut vérifier la relation (18) pour $p$ petit ($p \le 3$ par exemple) en utilisant la formule suivante donnant l'inverse de la somme d'une matrice inversible et d'une matrice de rang 1 :

$$ (19) \qquad (A + u \otimes v)^{-1} = A^{-1} - \frac{A^{-1} u \otimes v A^{-1}}{1 + v \cdot A^{-1} u}. $$

Mais un raisonnement par récurrence s'est avéré, jusqu'à présent, infructueux. Dans le théorème suivant, le choix de $J$ dans (6) est restreint à $\delta F(x^k)$.
### THÉORÈME 2. – La méthode de Newton généralisée (6) appliquée au système d'équations (3) ou (4) converge si :

(i) $a_{ij}^{-1} \le 0$ pour tout $i$ et $j$ avec $i \neq j$,

(ii) $r \le \min \frac{1}{a_{ii}^{-1}}$.

*Preuve.* – Il convient de vérifier les hypothèses du théorème 1. G est convexe car chaque composante $\phi_i$ est une fonction réelle convexe de la seule variable $\gamma_i$. De plus, la matrice A étant symétrique et définie positive, la solution du problème (1), (2) ou (3) est unique. D'après (13) et (14), il suffit de vérifier la condition (ii) du théorème 1 pour toute matrice J de $\delta G(0)$. Quitte à réarranger les vecteurs de base $(e_i)_i$, l'étude peut se limiter aux éléments de la famille $(J_p)_p$ définie en (16), avec C égal à $A^{-1}$. On procède alors par récurrence. Pour p égal à 0,

$$J_0^{-1} = I - r A^{-1} \ge 0 \text{ dès que (i) et (ii) sont satisfaits.}$$

On suppose qu'à l'ordre $p-1$, $J_{p-1}^{-1} \ge 0$. A l'ordre $p$ on a, par (19),

$$(20) \qquad J_p^{-1} = (J_{p-1} - e_p \otimes e_p)^{-1} = J_{p-1}^{-1} + \frac{J_{p-1}^{-1} e_p \otimes e_p J_{p-1}^{-1}}{1 - e_p \cdot J_{p-1}^{-1} e_p}.$$

Pour conclure, il suffit de montrer que $(1 - e_p \cdot J_{p-1}^{-1} e_p)$ est positif. Or en appliquant le lemme,

$$(21) \qquad 1 - e_p \cdot J_{p-1}^{-1} e_p = r \frac{|A_p^{-1}|}{|A_{p-1}^{-1}|}.$$

Ce coefficient est positif car les mineurs principaux d'une matrice définie positive sont strictement positifs.

La convergence de l'algorithme sur (3) se traduit par la convergence sur le système (2), mais le comportement monotone n'est pas transféré sur (2). Ce résultat intéressant est cependant d'un usage limité. En effet, la condition sur A est extrêmement exigeante. En particulier, elle n'est pas vérifiée dans les problèmes de contact unilatéral à l'origine de ce travail. Pour illustrer ce constat, on considère un problème d'équilibre d'un corps élastique, discrétisé dans une modélisation bidimensionnelle, en contact éventuel avec un obstacle rigide (Fig. 2) que l'on peut écrire sous la forme :

$$(22) \qquad \begin{cases} K u + N \lambda \\ \lambda - \operatorname{proj}(\lambda + r N^T u, \mathbb{R}^n) \end{cases} = \begin{cases} f, \\ 0, \end{cases}$$

avec $u \in \mathbb{R}^p$, $\lambda \in \mathbb{R}^n$, $N \in \mathbb{R}^{p \times n}$, $n < \frac{p}{2}$.

Remarquons que nous nous plaçons dans le cadre de l'élasticité linéarisée, non seulement pour des raisons de simplicité mais afin de ne pas mélanger les non linéarités issues des grandes déformations et celles provenant du contact. Pour une approche théorique du contact unilatéral en élasticité tridimensionnelle se référer à [12]. Le nombre de degrés
<figure><img src="image_3.png" /><figcaption>Fig. 2. – Problème de contact.</figcaption></figure>

de liberté est égal à *p* soit deux fois le nombre de nœuds libres du maillage. La lettre *n* désigne le nombre de nœuds de la surface en contact éventuel. *N* est une matrice donnant tous les vecteurs normaux de contact. Le système réduit (3) s'écrit alors :

$$ (23) \quad G(\gamma) = (I - r N^T K^{-1} N)^{-1} (\gamma - r N^T K^{-1} f) - \text{proj}(\gamma, \mathbb{R}^n) = 0, $$

$$ \text{avec } \gamma = \lambda + r N^T K^{-1} (f - N \lambda). $$

On est donc amené à examiner l'hypothèse *i* du théorème 2 sur $N^T K^{-1} N$. Or cette condition est violée pour des exemples élémentaires tel que le troisième de la Figure 3. Il y a, en effet deux coefficients non diagonaux positifs (gras). On peut cependant remarquer que le module de ces termes est petit au regard des autres. Le coefficient de Poisson est 0.3 pour les trois exemples, les éléments sont linéaires à quatre nœuds.

En conséquence, si la convergence monotone n'est pas établie rigoureusement pour les problèmes de contact unilatéral, en pratique, l'algorithme suit un tel comportement (2 à 3 itérations même pour un grand nombre de nœuds).

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

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

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

$$ N^T K^{-1} N = \begin{pmatrix} 2.582 & -0.878 \\ - & 2.582 \end{pmatrix} $$

$$ N^T K^{-1} N = \begin{pmatrix} 2.214 & -0.1424 & -0.0714 \\ - & 0.8916 & -0.1424 \\ - & - & 2.214 \end{pmatrix} $$

$$ N^T K^{-1} N = \begin{pmatrix} 2.2 & -0.145 & \mathbf{0.0246} & -0.041 \\ - & 0.89 & -0.075 & \mathbf{0.0246} \\ - & - & 0.89 & -0.145 \\ - & - & - & 2.2 \end{pmatrix} $$

<caption>Fig. 3. – Problèmes de contact et matrices associées.</caption>
## 4. Complexité et stabilité

A défaut d'un résultat de convergence applicable, nous sommes en mesure de fournir des résultats sur le comportement de la méthode de Newton généralisée, en particulier sur la complexité maximale et sur les formes d'instabilités envisageables. L'opérateur $F$ introduit en (2) comporte $2^n$ cônes de linéarité et la convergence du schéma (6), si elle a lieu, s'effectue alors au plus en $2^n$ itérations. En effet, le nombre d'itérés possibles est égal au nombre de régions de linéarité. D'autre part, si pour une valeur initiale donnée, l'algorithme passe deux fois non consécutives par la même zone, il y aura nécessairement divergence sous forme d'un cycle sur un nombre fini d'itérés. Deux itérés consécutifs identiques traduisent la convergence. En fait, le résultat suivant permet de faire mieux.

A partir de la formulation (2), on montre aisément que tout itéré $x^k = (u^k, \lambda^k)$ vérifie,

$$ (24) \qquad u_i^k = 0 \quad \text{ou} \quad \lambda_i^k = 0, \quad i = 1, n. $$

Cette relation traduit le fait qu'à chaque itération l'algorithme résout un problème avec un certain nombre de contraintes d'égalité ($u_i^k = 0, i \in Ic$), les autres composantes de $u^k$ n'étant point contraintes ($\lambda_i^k = 0, i \in \{1, n\}/Ic$). La condition de complémentarité (produit scalaire $u^k \cdot \lambda^k = 0$) est donc bien satisfaite à chaque itération, mais les conditions d'inégalité ($u^k \ge 0, \lambda^k \le 0$), ne sont remplies qu'à la convergence par la solution. On peut remarquer que le facteur $r$ ne joue aucun rôle sur ce comportement.

**THÉORÈME 3.** – La convergence de la suite $(x^k)$ engendrée par (6) sur le système (2), si elle a lieu, s'effectue en au plus $2^{n-1}$ itérations.

**Preuve.** – On montre que seule la moitié des régions de linéarité contient un itéré de Newton éventuel. Les cônes de linéarité sont des quadrants de $\mathbb{R}^n$, pour la variable $\theta = \lambda + ru$. Une propriété caractéristique concerne le cône polaire $X^0$ d'un quadrant $X$ [26].

$$ (25) \qquad X^0 = \{z \in \mathbb{R}^n; z \cdot \theta \le 0, \theta \in X\} = -X = \{z \in \mathbb{R}^n; z = -\theta, \theta \in X\}. $$

Supposons donc que deux itérés appartiennent à deux quadrants mutuellement polaires,

$$ (26) \qquad \theta^k = \lambda^k + ru^k \in X \quad \text{et} \quad \theta^1 = \lambda^1 + ru^1 \in X^0 = -X. $$

Les relations (24) permettent de préciser les choses,

$$ (27) \qquad (\lambda^k \in X \text{ et } u^k \in X) \quad \text{et} \quad (\lambda^1 \in -X \text{ et } u^1 \in -X). $$

La définition des cônes mutuellement polaires fournit les inégalités suivantes,

$$ (28) \qquad \lambda^1 \cdot \lambda^k \le 0, \quad u^k \cdot u^1 \le 0, \quad \lambda^1 \cdot u^k \le 0, \quad \lambda^k \cdot u^1 \le 0. $$
La première équation de (2) exprimée pour les deux itérations 1 et $k$ et le caractère défini positif de la matrice $A$ conduisent à la relation,

$$ (29) \qquad (u^k - u^1) \cdot (\lambda^k - \lambda^1) = -(u^k - u^1) \cdot A(u^k - u^1) < 0. $$

D'autre part, en développant le premier membre de l'équation précédente, on obtient,

$$ (30) \qquad (u^k - u^1) \cdot (\lambda^k - \lambda^1) = u^k \cdot \lambda^k + u^1 \cdot \lambda^1 - u^k \cdot \lambda^1 - u^1 \cdot \lambda^k \geq 0. $$

En effet, les deux premiers termes sont nuls car la condition de complémentarité est satisfaite à chaque itération et les deux derniers sont négatifs de par les inégalités (28). Les relations (29) et (30) sont contradictoires, ce qui termine la démonstration.

Ce niveau maximal de complexité pour l’algorithme de Newton généralisée est deux fois moins élevé que celui obtenu pour la méthode de Lemke [23] égal au maximum $2^n$ [15]. Cette comparaison est licite en terme de nombre d’itérations mais pas nécessairement en temps de calcul, car, à chaque itération la méthode de Lemke ne résout pas tout un système linéaire.

Le théorème précédent n’élimine pas cependant des instabilités sous forme de cycles sur plusieurs itérés. C’est, rappelons-le, le seul mode de divergence de la méthode de Newton appliquée à des systèmes linéaires par morceaux. Il est pourtant possible d’exclure les cycles d’ordre deux, c’est-à-dire portant sur deux itérés successifs.

**THÉORÈME 4.** – La suite $(x^k)$ ne peut diverger via un cycle d’ordre deux.

*Preuve.* – Considérons $k$ tel que $x^{k+2} = x^k$, avec $x^{k+1} \neq x^k$, et écrivons (24) sous la forme plus détaillée suivante,

$$ (31) \qquad \theta_1^{k+1} = \begin{cases} \lambda_i^{k+1} & \text{si } \theta_i^k < 0 \quad (u_i^{k+1} = 0) & (a) \\ r u_i^{k+1} & \text{si } \theta_i^k \ge 0 \quad (\lambda_i^{k+1} = 0) & (b) \end{cases} $$

On pose $E = \{i \in \{1, \dots, n\}; \theta_i^{k+1} \ge 0\}$. Si $i$ n’appartient pas à $E$, alors il existe deux possibilités :

$$ (32) \qquad * \theta_i^k \ge 0 \Rightarrow \theta_i^{k+1} = r u_i^{k+1} \ge 0 \quad (\text{par 31 b}) \Rightarrow \theta_i^{k+2} = r u_i^{k+2} = \theta_i^k, $$

$$ (33) \qquad * \theta_i^k < 0 \Rightarrow \theta_i^{k+1} = \lambda_i^{k+1} < 0 \quad (\text{par 31 a}) \Rightarrow \theta_i^{k+2} = \lambda_i^{k+2} = \theta_i^k. $$

On peut synthétiser les relations (32) et (33) sous la forme d’une inclusion sous-différentielle utilisant le sous-différentiel de la fonction indicatrice des réels positifs :

$$ (34) \qquad \lambda_i^1 \in \partial \Psi_+(u_i^1), \quad i \in E, \quad 1 = k, k+1. $$
Pour i n'appartenant pas à E, par un raisonnement analogue, on obtient :

$$
(35) \qquad \lambda_i^1 \in \partial \Psi_-(u_i^1), \quad i \notin E, \quad 1 = k, k+1.
$$

Les inclusions (34) et (35) se résument en l'unique inclusion :

$$
(36) \quad \lambda_i^1 \in \partial \Psi_C(u_i^1), \quad 1 = k, k+1, C = \{z \in \mathbb{R}^n; z_i \ge 0 \text{ si } i \in E, z_i \le 0 \text{ si } i \notin E\}.
$$

Pour conclure, la monotonie du sous-différentiel [26] conduit à l'inégalité suivante,
contradictoire avec la relation (29) en prenant 1 égal à k + 1,

$$
(37) \quad (\lambda^{k+1} - \lambda^k) \cdot (u^{k+1} - u^k) \ge 0
$$

Ce résultat, d'une portée limitée, garantit cependant la convergence pour *n* égal à 1 ou 2 (ou *p* quelconque et *n* égal à 1 ou 2 si l'on se réfère à l'exemple développé en Section I.3). D'autre part, ce type de cycle est celui que l'on rencontre pour les problèmes de frottement ou de contact avec frottement que l'on examine dans la seconde partie.

**Partie II : Approche heuristique**

Par approche heuristique, nous entendons une étude portant sur des exemples élémentaires
par leur faible nombre d'inconnues, doublée d'une analyse du comportement numérique de
l'algorithme sur des problèmes de plus grande taille. Une telle démarche n'apporte pas de
renseignements supplémentaires pour le problème d'obstacle (1), car, en pratique, aucune
instabilité n'a jamais été constatée. L'étude menée dans la première partie a apporté une
réponse partielle à ce constat. D'autre part le théorème 4 assure la convergence pour les
exemples élémentaires (n = 1 ou 2) bases de cette deuxième partie.

1. Frottement

Il en va différemment pour un autre problème d'optimisation non différentiable (de
type frottement) :

$$
(38) \quad u = \arg \min \left\{ \frac{1}{2} v \cdot Av - v \cdot f + \alpha \|v\|_{\infty} ; v \in \mathbb{R}^n \right\} \quad \text{où} \quad \|v\|_{\infty} = \max_i |v_i|.
$$

En introduisant un multiplicateur, le système à résoudre s’écrit formellement comme (2)
avec un opérateur projection quelque peu différent :

$$
(39) \quad F(u, \lambda) = 0 \iff \left\{
\begin{array}{ll}
Au + \lambda - f = 0 & \mathcal{C} = \prod_{i=1}^{n} [-\alpha_i, \alpha_i] \\
\lambda - \operatorname{proj}(\lambda + ru; \mathcal{C}) = 0, &
\end{array}
\right.
$$
Il en va de même pour la forme simplifiée analogue à (3),

$$ (40) \qquad G(\gamma) = (I - r A^{-1})^{-1} \gamma + g - \operatorname{proj}(\gamma; C) = 0. $$

L'application $G$ n'est plus linéaire par cônes mais par morceaux. La méthode d'amortissement présentée en Section I.1 est applicable. Mais il est instructif d'examiner le comportement de la méthode de Newton généralisée standard (non amortie) sur le problème le plus simple ($n = 1$). Il s'agit d'un oscillateur élémentaire disposant d'un seul degré de liberté présenté en médaillon sur la Figure 4 et pour lequel $\alpha = \mu q_n$ ($\mu$ est le coefficient de frottement). On met alors en évidence sur les Figures 4 et 5 deux types de graphes de $G$ (40) selon la valeur du facteur $r$. Il est clair qu'un cycle d'ordre 2 peut intervenir avec la seconde catégorie.

<figure><img src="image_7.png" /><figcaption>Fig. 4. – Problème de frottement élémentaire et graphe de $G$ (1re espèce).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 5. – Problème de frottement élémentaire et graphe de $G$ (2e espèce).</figcaption></figure>
On peut classer de telles fonctions de la manière suivante. Une fonction réelle continue présente un point d'inflexion en un point *x* s'il y a changement de convexité. Cette définition ne nécessite pas l'existence d'une dérivée seconde en tout point, hypothèse non satisfaite ici. Nous ne considérons dorénavant que des fonctions réelles lipschitziennes monotones et n'admettant qu'un changement de convexité. Autrement dit l'ensemble des points d'inflexion est un intervalle ouvert. Une notion faible de différentiabilité est accessible : le sous-différentiel de Clarke [13], dont la définition est comparable à celle du jacobien généralisé (5). On distingue alors deux types de fonctions selon que la valeur absolue de la dérivée en un point d'inflexion réalise le minimum (1<sup>re</sup> espèce, Fig. 6) ou le maximum (2<sup>e</sup> espèce, Fig. 7) sur R. Un tel énoncé n'est pas tout à fait correct, car la dérivée est remplacée par le sous-différentiel de Clarke qui est éventuellement multivoque, mais plus simple à exprimer. Des théorèmes de convergence sont alors disponibles pour chaque espèce de fonctions.

<figure><img src="image_9.png" /><figcaption>Fig. 6. – Fonction de 1<sup>re</sup> espèce.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 7. – Fonction de 2<sup>e</sup> espèce.</figcaption></figure>
### THÉORÈME 5.
Si $f$ est une fonction réelle de 1re espèce, la suite $(x^k)$ définie comme suit converge,

$$ (41) \qquad x^{k+1} = x^k - \frac{f(x^k)}{y^k} \quad \text{avec} \quad y^k \in \partial f(x^k). $$

*Preuve.* – Sans restreindre la généralité, on peut supposer $f$ croissante et 0 point d'inflexion. Il s'en suit que $f$ est concave sur $(-\infty, 0[$ et convexe sur $]0, +\infty)$. La solution est notée $s(f(s) = 0)$. Il est clair qu'alors,

$$ (42) \qquad x^k > s \Rightarrow x^{k+1} \le x^k \quad \text{et} \quad x^k < s \Rightarrow x^{k+1} \ge x^k. $$

C'est-à-dire que $(x^k)$ décroît tant que $x^k$ est supérieur à $s$ et croît dès que $x^k$ est inférieur à $s$. Supposons encore que la solution $s$ soit positive. Deux cas se présentent selon la valeur initiale $x^0$.

– Si $x^0$ est supérieur à $s$, $(x^k)$ décroît vers $s$. En effet, la convexité de $f$ donne :

$$ (43) \qquad f(x^1) \ge f(x^0) + y^0 (x^1 - x^0) = f(x^0) - y^0 \frac{f(x^0)}{y^0} = 0 $$

Donc $f(x^1) \ge 0$ c'est-à-dire $x^1 \ge s$. La suite $(x^k)$ est alors décroissante et minorée par $s$. Elle converge vers une valeur $q$ supérieure ou égale à $s$. La multifonction $\partial f$ étant semi-continue supérieurement, $\partial f(x^k)$ tend vers $\partial f(q)$ et $q$ vérifie,

$$ (44) \qquad q = q - \frac{f(q)}{y}, \quad y \in \partial f(q). $$

La valeur $q$ satisfait alors l'équation $(f(q) = 0)$ et est donc égale à $s$ par unicité.

– Si $x^0$ est inférieur à $s$, soit la suite est majorée par $s$ et converge donc vers $s$ ; soit il existe un rang $m$ pour lequel $x^k$ est supérieur à $s$ et l'on retrouve le premier cas.

Pour les fonctions de 2e espèce, la convergence n'est assurée que sous une hypothèse restrictive comme l'énonce le théorème suivant.

### THÉORÈME 6.
Si $f$ est une fonction réelle de 2e espèce et de plus vérifie la condition,

$$ (45) \qquad C_M < 2C_m $$

où

$$ C_M = \max \{|y|; y \in \partial f(x), x \in \mathbb{R}\} $$

$$ (46) \qquad \text{et} $$

$$ C_m = \min \{|y|; y \in \partial f(x), x \in \mathbb{R}\}, $$

alors la suite $(x^k)$ définie en (41) converge.
Preuve. – L’expression (46) sous-entend que les extrema sont atteints, et (45) que le minimum est strictement positif. Sans restreindre la généralité, on peut supposer $f$ croissante, 0 point d’inflexion et la solution $s$ positive. Il s’en suit que $f$ est convexe sur $(-\infty, 0[$ et concave sur $]0, +\infty)$ et que, si $x^0$ (ou $x^m$ à partir d’un certain rang) appartient à $[0, s]$, alors la suite $(x^k)$ est croissante et tend vers $s$.

Montrons au préalable que si deux itérés successifs n’appartiennent pas à l’intervalle $[0, s]$, ils sont alors de signe opposé,

$$ (47) \qquad (x^k, x^{k+1}) \notin [0, s] \times [0, s] \Rightarrow x^k x^{k+1} < 0. $$

En effet, considérons dans un premier temps que $x^k$ est supérieur à $s$, et supposons que $x^k$ est positif. Alors la concavité de $f$ sur $[x^{k+1}, x^k]$ conduit à,

$$ (48) \qquad f(x^{k+1}) \le f(x^k) + y(x^{k+1} - x^k) \quad \text{pour tout } y \in \partial f(x^k), $$

$$ (49) \qquad f(x^{k+1}) \le f(x^k) + y^k (x^{k+1} - x^k) \quad \text{pour } y^k \text{ définissant } x^{k+1}, $$

$$ (50) \qquad f(x^{k+1}) \le f(x^k) - f(x^k) = 0, \text{ ce qui est contraire à l'hypothèse.} $$

Pour le cas $x^k$ négatif, on procède comme précédemment en utilisant la convexité de $f$ sur $[x^k, x^{k+1}]$. Démontrons maintenant l’implication suivante,

$$ (51) \qquad (x^k, x^{k+1}) \notin [0, s] \times [0, s] \Rightarrow |x^{k+1} - s| \le c|x^k - s|, \quad 0 < c < 1. $$

En considérant, comme premier cas, $x^k$ négatif, la convexité sur $[x^k, 0]$ et la concavité sur $[0, s]$ donnent les inégalités suivantes,

$$ (52) \qquad f(0) - f(x^k) \le y(0 - x^k), \quad \text{pour tout } y \in \partial f(0), $$

$$ (53) \qquad -f(0) = f(s) - f(0) \le y(s-0), \quad \text{pour tout } y \in \partial f(0). $$

En additionnant (52) et (53) et en remplaçant dans (41), il vient,

$$ (54) \qquad x^{k+1} - s = x^k - s - \frac{f(x^k)}{y^k} \le \left( \frac{y}{y^k} - 1 \right) (s - x^k), \quad \text{pour tout } y \in \partial f(0). $$

Selon l'hypothèse (45), il existe un $\varepsilon$ positif tel que $C_M = (2 - \varepsilon) C_m$ et l'inégalité (51) est satisfaite pour $c = 1 - \varepsilon$. Remarquons que le supremum $C_M$ sera bien atteint sur le sous-différentiel de $f$ en 0.
Dans le cas où $x^k$ est supérieur à $s$, la concavité sur $[s, x^k]$ implique,

$$
(55) \quad f(x^k) = f(x^k) - f(s) \le y^s (x^k - s), \quad \text{pour tout } y^s \in \partial f(s),
$$

$$
(56) \quad s - x^{k+1} \le \left( \frac{y^s}{y^k} - 1 \right) (x^k - s) \le (1 - \varepsilon) (x^k - s).
$$

On en déduit alors la convergence de la suite (x^k) vers s.

La non différentiabilité alourdit les notations sans introduire pour autant de difficultés supplémentaires. L'hypothèse (45) signifie que le module de la dérivée ne peut varier que dans un rapport de 1 à 2. Ceci suppose donc une faible non linéarité. On retrouve ici un résultat plus général (différentiabilité mise à part) de Necas [27], qui exige que la plus grande valeur propre de toutes les matrices jacobiennes n'excède deux fois la plus petite. Dans notre cas, cette condition devrait être satisfaite au préalable par la matrice A, ce qui est extrêmement restrictif. Il est clair qu'une telle hypothèse, comme extrapolation de (45), ne tient pas compte du caractère directionnel de (45). Ceci explique que nous n'avons pas essayé de généraliser le théorème de Necas aux systèmes non différentiables. D'autant qu'un tel résultat ne peut rendre compte du théorème 5 pour les fonctions de 1re espèce. On comprend mieux aussi le recours que nous avons ici à une approche heuristique.

L'application des deux théorèmes précédents fournit une valeur maximale pour le facteur $r$ en deça de laquelle la convergence est assurée. En effet, pour le graphe de 2e espèce (Fig. 5), la condition (45) donne,

$$
(57) \quad \left| (1 - r a^{-1})^{-1} - 1 \right| < 2 \left| (1 - r a^{-1})^{-1} \right| \Rightarrow r < 2a.
$$

Si cette hypothèse n'est pas satisfaite, un cycle d'ordre 2 intervient entre les deux régions de glissement si la solution est adhérente (zone centrale).

A ce stade, il est intéressant de comparer cette estimation à celle obtenue pour l’algorithme d’Uzawa que l’on identifie à une méthode d’approximations successives sous la forme suivante appliquée au système (39) ou (40),

$$
(58) \quad \lambda^{k+1} = L(\lambda^k) \equiv \operatorname{proj}[\lambda^k + r A^{-1}(f - \lambda^k); \mathcal{C}].
$$

La condition de contraction de $L$ assurant la convergence s'exprime alors,

$$
(59) \quad \rho (I - r A^{-1}) < 1 \Leftrightarrow r < \frac{2}{\lambda_{\max} (A^{-1})} = r_{uz},
$$

où $\rho$ désigne le rayon spectral de la matrice et $\lambda_{\max}$ la plus grande valeur propre.

Pour l'équation élémentaire (n = 1), on obtient la condition (57) ! Cette coïncidence est assez troublante même si la méthode d'Uzawa peut apparaître comme un algorithme
<figure><img src="image_11.png" /><figcaption>Fig. 8. - Oscillateur à 2 degrés de liberté avec frottement.</figcaption></figure>

de Newton modifié [3], c'est-à-dire où seule est prise en compte la dérivée de la partie différentiable (ici linéaire).

La condition (59) peut donc être vue comme une extrapolation de l'inégalité (57) pour des systèmes de taille quelconque. L'exemple qui suit confirme le bien fondé de cette estimation. La Figure 8 représente un oscillateur à deux degrés de liberté pour lequel sont également précisés la matrice A et le second membre f. Sur les Figures 9, 10, 11 et 12 sont dessinés, dans le plan $(\gamma_1, \gamma_2)$ les 9 zones de linéarité, les 9 itérés potentiels de la méthode de Newton et l'ensemble des chemins possibles selon la région de linéarité à laquelle appartient la valeur initiale (trait pointillé). Un chemin particulier est représenté en trait continu. Le choix du second membre f permet de mettre en évidence un cycle d'ordre 3 donc plus complexe que le cycle élémentaire d'ordre 2 précédent. Cependant la convergence est établie pour toute valeur initiale si

$$ (60) \qquad r < 2 \max \left( \frac{1}{a_{11}^{-1}}, \frac{1}{a_{22}^{-1}} \right) = r_{\max}. $$

On constate que $r_{\max}$ est supérieur à $r_{uz} = 2\lambda_{min}(A) = 3 - \sqrt{5} \approx 0.76$ et que l'estimation $r_{uz}$ peut être retenue en pratique tant pour le problème de frottement que pour les problèmes des sections suivantes.

<figure><img src="image_12.png" /><figcaption>Fig. 9. - Oscillateur double, comportement de l'algorithme ($r = 0.25 < 0.5$).</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Fig. 10. – Oscillateur double, comportement de l’algorithme ($r = 0.75$).</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 11. – Oscillateur double, comportement de l’algorithme ($r = 1.5$).</figcaption></figure>

## 2. Contact avec frottement

Il s'agit donc d'un problème couplant contact unilatéral et frottement sec. Seul sera étudié un exemple élémentaire représentant un treillis de ressorts et un seul nœud en contact affleurant avec un plan rigide (Fig. 13). Le système à résoudre a alors la forme suivante,

$$ (61) \qquad \begin{cases} A u + \lambda - f = 0, \\ \lambda - \gamma_n^{-} n - \operatorname{proj}(\gamma_t; C(\gamma_n^{-})) = 0, \end{cases} \quad t \ge 0, $$

$$ (62) \quad \text{avec } \gamma = \lambda + ru = \gamma_t t + \gamma_n n, \ u \in \mathbb{R}^2, \ \lambda \in \mathbb{R}^2, \ (t, n) \text{ repère orthonormé,} $$

$$ (63) \quad \gamma_n^{-} = \operatorname{proj}(\gamma_n; \mathbb{R}_{-}) = \min(0, \gamma_n), \ C(\gamma_n^{-}) = [\mu \gamma_n^{-}, -\mu \gamma_n^{-}]. $$
<figure><img src="image_15.png" /><figcaption>Fig. 12. – Oscillateur double, comportement de l’algorithme ($r = 2.1 > 2$).</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>Fig. 13. – Problème élémentaire de contact avec frottement.</figcaption></figure>

Le système réduit s’exprime alors,

$$ (64) \qquad G(\gamma) = (I - r A^{-1})^{-1} \gamma + g - \gamma_n^{-} n - \operatorname{proj}(\gamma_t; C(\gamma_n^{-})) t = 0. $$

Il est plus commode pour la suite d’utiliser l’expression suivante en pré-multipliant le système précédent par la matrice $B = I - r A^{-1}$,

$$ (65) \qquad H(\gamma) = \gamma - \gamma_n^{-} (I - r A^{-1}) n - \operatorname{proj}(\gamma_t; C(\gamma_n^{-})) (I - r A^{-1}) t = -(I - r A^{-1}) g = h. $$

L’application $H$ est alors linéaire par cônes et égale à l’identité sur un demi-espace. Les matrices jacobiennes ($H_{\text{écart}}$, $H_{\text{adh}}$, $H_{s\epsilon}$) associées à chaque région sont données sur la Figure 14. La solution du système (65) est unique si [1],

$$ (66) \quad \det(H_{s\epsilon}) = \det(I - B v_{\epsilon} \otimes n) = 1 - n \cdot B v_{\epsilon} > 0 (\epsilon = \pm 1) \Leftrightarrow \mu < \frac{n \cdot A^{-1} n}{|n \cdot A^{-1} t|}. $$
<figure><img src="image_17.png" /><figcaption>Fig. 14. – Cônes de linéarité de H (65) et itérés de Newton.</figcaption></figure>

Nous notons $\gamma_g$, $\gamma_\varepsilon$, et $\gamma_a$ les itérés de Newton potentiels associés aux cônes Ecart (non contact), $S\varepsilon$ (glissement + et -) et Adh (contact adhérant). Nous discutons la convergence selon que la solution, notée $\gamma_s$, appartient ou non au cône Ecart.

Premier cas : $\gamma_s = \gamma_g$. C'est-à-dire $H(\gamma_s) = \gamma_s = h$ tel que $(\gamma_s)_n = h_n \ge 0$.

En exprimant les itérés $\gamma_\varepsilon$, et en tenant compte de la condition d'unicité (66) ils s'avèrent appartenir au demi-espace Ecart,

$$ (67) \qquad \gamma_{\varepsilon} = \left( I + \frac{1}{1-n.Bv_{\varepsilon}} B v_{\varepsilon} \otimes n \right) h \quad \text{avec} \quad (\gamma_{\varepsilon})_n = \frac{h_n}{1-n.Bv_{\varepsilon}} \ge 0. $$

Il y a donc trois itérés éventuels sur quatre dans le cône de la solution ; il ne peut y avoir cycle !

Deuxième cas : $\gamma_s \neq \gamma_g$. Aucun itéré ne pouvant appartenir au cône Ecart, on peut donc supprimer la projection sur les réels négatifs en (63) et résoudre le système :

$$ (68) \qquad H(\gamma) = \gamma - \gamma_n (I - r A^{-1}) n - \text{proj}(\gamma_t; C(\gamma_n)) (I - r A^{-1}) t = h. $$

A ce stade, il est nécessaire de décomposer $H(\gamma)$ dans la base $(t, n)$ en posant pour simplifier l'écriture,

$$ (69) \qquad c_{nn} = n \cdot A^{-1} n, \quad c_{tt} = t \cdot A^{-1} t, \quad c_{nt} = n \cdot A^{-1} t = t \cdot A^{-1} n, $$

$$ (70) \qquad \begin{cases} \gamma_t + r c_{nt} \gamma_n - \text{proj}(\gamma_t; C(\gamma_n)) (1 - r c_{tt}) = h_t, \\ r c_{nn} \gamma_n - \text{proj}(\gamma_t; C(\gamma_n)) r c_{nt} = h_n. \end{cases} $$

En éliminant $\gamma_n$, le système (70) se réduit à l'équation,

$$ (71) \qquad G(\gamma_t) = \frac{c_{nn}}{c_{nn} - r \det A^{-1}} \gamma_t - \text{proj}(\gamma_t; C(\gamma_n)) = \frac{c_{nn}}{c_{nn} - r \det A^{-1}} \left( h_t - \frac{c_{nt}}{c_{nn}} h_n \right), $$
qui a donc la forme du problème de frottement (40). On peut donc distinguer trois régions de linéarité sur $\mathbb{R}$ et deux valeurs possibles des pentes,

$$ (72) \qquad G'(\gamma_t) = \frac{c_{nn}}{c_{nn} - r \det A^{-1}} - 1 = \frac{r \det A^{-1}}{c_{nn} - r \det A^{-1}} \quad \text{sur } ]\mu\gamma_n, -\mu\gamma_n[. $$

$$ (73) \qquad G'(\gamma_t) = \frac{c_{nn}}{c_{nn} - r \det A^{-1}} \quad \text{sur } (-\infty, \mu\gamma_n[\cup] - \mu\gamma_n, +\infty). $$

En exprimant que ces dérivées doivent rester dans un rapport de 1 à 2, on obtient la condition de convergence suivante,

$$ (74) \qquad r < 2 \frac{c_{nn}}{\det A^{-1}} = \frac{2}{c_{tt} - \frac{c_{nx}^2}{c_{nn}}} = r_{\max}. $$

On constate que,

$$ (75) \qquad r_{\max} \ge \frac{2}{c_{tt}} \ge \frac{2}{\lambda_{\max}(A^{-1})}, $$

et donc que la condition (74) est relaxée par rapport à la relation (59) par le couplage du contact unilatéral et du frottement. Cette démonstration est bien entendu difficilement transférable à des problèmes de plus grande taille, c'est-à-dire avec plus d'un seul nœud en contact éventuel. D'ailleurs, avec seulement deux nœuds, on peut trouver des exemples de cycles complexes d'ordre 3 voire 4 ($r$ assez grand) [4]. Cependant la relation (75) est en pratique très largement suffisante pour assurer la convergence.

Le problème d'un seul contact nodal en discrétisation tridimensionnelle est plus difficile car les équations ne sont plus linéaires par cônes mais seulement linéaires par rayons (raywise linear) selon la terminologie développée dans [1].

### 3. Frottement et élasto-plasticité

La encore il faut se contenter d'une étude sur un modèle élémentaire. Une telle démarche a été motivée par certaines instabilités constatées sur des logiciels de simulation de l'emboutissage utilisant la méthode de Newton comme solveur non linéaire standard. Le modèle analogique d'oscillateur non linéaire est présenté sur la Figure 15. La loi de comportement n'est donnée qu'en traction étant entendu que le graphe est symétrique par rapport à l'origine pour la compression. On peut donc distinguer trois régions de linéarité pour l'application $A(.)$: élasticité (E) autour de l'origine, plasticité en traction (P+) et plasticité en compression (P−). Remarquons qu'est envisagé un coefficient d'écrouissage $k$ non nul. Le système à résoudre en présence de frottement (coefficient $\mu$) et d'une force normale $f_n$, est le suivant,

$$ (76) \qquad F(u, \lambda) = 0 \iff \left\{ \begin{array}{ll} A(u) + \lambda - f_t = 0, & C = f_n[-\mu, \mu] \\ \lambda - \operatorname{proj}(\lambda + ru; C) = 0, & \end{array} \right. $$
<figure><img src="image_18.png" /><figcaption>Fig. 15. – Oscillateur non linéaire modèle.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 16. – Régions de linéarité et itérés de Newton éventuels.</figcaption></figure>

A(.) étant non linéaire, on ne peut analyser le comportement de la méthode de Newton sur un système réduit à une seule équation comme en (40). En effet, même s'il est aisé de définir ici la fonction réciproque de A(.) et par conséquent d'effectuer formellement l'élimination de la première équation du système (76), le comportement de la méthode de Newton généralisée sur l'équation réduite ne traduit pas celui sur le système entier. On est donc amené à faire une étude directe sur l'application F qui présente 9 zones de linéarité détaillées sur la Figure 16.

L'expérience acquise sur les problèmes précédents nous conduit à considérer le cycle entre les régions (P-, S-) et (P+, S+) comme l'instabilité fondamentale à éviter. Trouver une parade à une telle situation ne garantit pas totalement la convergence de l'algorithme. Une discussion graphique permet de montrer que pour un facteur *r* assez petit, le cycle précédent peut être éliminé. Comme on ne peut représenter le graphe d'une application de $\mathbb{R}^2$ dans $\mathbb{R}^2$, sont représentés à la fois sur la Figure 16 le plan $(u, \lambda)$, sur lequel on distingue les 9 régions de linéarité (4 cônes, 4 demi-bandes et 1 quadrilatère en gris), et le plan $(u, f_t)$ sur lequel on peut suivre le cheminement de l'algorithme de Newton. Pour $f_t = 0$ et un couple initial $(u^0, \lambda^0)$ situé dans la zone (P+, S+), on peut déterminer graphiquement l'itéré suivant $(u^1, \lambda^1)$. Avec la valeur de *r* choisie sur la Figure, cet itéré n'est pas dans la région (P-, S-), mais pourrait s'y trouver pour une valeur plus grande, la bande caractérisant les régions associées à l'adhérence devenant plus oblique. En écrivant que $(u^1, \lambda^1)$ appartient à la zone (P-, A) $(\lambda^1 + ru^1 > \mu f_n)$, on obtient la condition suivante,

$$ (77) \qquad r < 2k \frac{1}{1 + \frac{SoK}{-\mu f_n (K+k)}} = r_{\max}. $$
On remarque que $r_{max}$ tend vers zéro quand $f_n$ tend vers zéro. En extrapolant à un problème couplant contact unilatéral et frottement ($f_n$ devient $\lambda_n$, nouvelle inconnue du problème), cette dernière estimation du facteur $r$ devient exigeante quand l'effort normal $\lambda_n$ devient faible. Or une telle situation se produit au cours de la simulation quasi-statique d'un processus d'emboutissage dans certaines zones de contact sur les flancs du poinçon et de la matrice signalée sur la Figure 17. Il s'agit d'une phase délicate du processus où l'ensemble de la tôle située sous le serre-flan commence à glisser caractérisant le début de l'avalement. En considérant que la géométrie représentée sur la Figure est axisymétrique, les nœuds situés dans les régions mentionnées plus haut ont pour le moment glissé vers l'extérieur. Avec l'avalement, le sens de glissement va s'inverser en passant momentanément par une phase d'adhérence. Fort de l'analyse effectuée ici, en diminuant de manière drastique le coefficient $r$, nous avons réussi à franchir cette étape [24].

<figure><img src="image_20.png" /><figcaption>Fig. 17. – Processus d'emboutissage.</figcaption></figure>

Conclusion

Cet article constitue un faisceau de résultats sur la convergence de la méthode de Newton généralisée dans le contexte de la mécanique du contact. Une telle étude était nécessaire car l'usage d'une méthode de Newton sur un système non différentiable n'est pas standard. Le mathématicien peut être déçu par le manque de généralité des résultats obtenus, mais l'ingénieur y trouvera des justifications théoriques au bon comportement constaté en pratique et les moyens de remédier à certaines instabilités. Des domaines autres que la mécanique du contact pourront peut-être exploiter le travail développé dans cet article. Mais l'extension des théorèmes présentés ici exigera sans doute des outils mathématiques nouveaux.
[1] P. ALART, Critères d'injectivité et de surjectivité pour certaines applications de $R^n$ dans lui-même : application à la mécanique du contact, *RAIRO, Modélisation Mathématique et Analyse Numérique*, 27, n° 2, 1993, p. 203-222.

[2] P. ALART, Multiplicateurs « augmentés » et méthode de Newton généralisée pour contact avec frottement, *Document LMA-DME-EPFL*, Lausanne, 1988.

[3] P. ALART et A. CURNIER, A Mixed Formulation for Frictional Contact Problems Prone to Newton Like Solution Methods, *Computer Methods in Applied Mechanics and Engineering*, 92, n° 3, 1991, p. 353-375.

[4] P. ALART et A. CURNIER, Contact discret avec frottement : unicité de la solution, convergence de l’algorithme, *Document LMA-DME-EPFL*, Lausanne, 1987.

[5] P. ALART et A. HEEGE, Consistent Tangent Matrices of Hybrid Curved Contact Operators Involving Anisotropic Friction, *European Journal of Finite Elements*, 4, n° 2, 1995, p. 183-207.

[6] P. ALART, A. HEEGE et G. KAPELSKI, Développement et validation d'un logiciel de simulation de pression entre plaquettes et disque de frein, Rapport de fin de contrat REGIENOV, Mars 1992.

[7] P. ALART et F. LEBON, Solution of Frictional Contact Problems Using ILU and Coarse/Fine Preconditioner, *Computational Mechanics*, 16, 1995, p. 98-105.

[8] P. ALART, M. BARBOTEU et F. LEBON, Solution of Frictional Contact Problems Using EBE preconditionners, en préparation.

[9] P. ALART, F. LEBON, F. QUITTAU et K. REY, Frictional Contact Problem in Elastostatics: Revisiting the Uniqueness Condition, Proceeding of 2nd Contact Mechanics International Symposium, 94, September 19-23 1994, Carry Le Rouet.

[10] R. CHARLIER et S. CESCOTTO, Modélisation du phénomène de contact unilatéral avec frottement dans un contexte de grandes déformations, *Journal de Mécanique Théorique et Appliquée*, Special issue: Numerical Method in Mechanics of Contact Involving Friction, 1988, p. 177-192.

[11] P. G. CIARLET, *Introduction à l'analyse numérique matricielle et à l'optimisation*, Masson, Paris, 1985.

[12] P. G. CIARLET et J. NECAS, Unilateral Problems in Nonlinear Three-Dimensional Elasticity, *Arch. for Rat. Mech. and Anal.*, 87, n° 4, 1985, p. 319-338.

[13] F. H. CLARKE, *Optimization and Nonsmooth Analysis*, Wiley, 1983.

[14] A. CURNIER et P. ALART, A Generalized Newton Method for Contact Problems with Friction, *Journal de Mécanique Théorique et Appliquée*, Special Issue: Numerical Method in Mechanics of Contact Involving Friction, 1988, p. 67-82.

[15] Y. FATHI, Computational complexity of LCP's associated with positive definite symmetric matrices, *Math. Prog. St. 17*, 1979, p. 335-344.

[16] M. FORTIN et G. GLOWINSKI, *Méthodes de Lagrangien augmenté*, Dunod, Collection Méthodes Mathématiques de l'Informatique, 1982.

[17] A. HEEGE, Simulation numérique 3D du contact avec frottement et application à la mise en forme, *Thèse de Doctorat*, INPG, Grenoble, 1992.

[18] A. HEEGE, P. ALART et E. ONATE, Numerical modelling and simulation of frictional contact using a generalized Coulomb law, *Engineering Computations*, 12, n° 7, 1995, p. 641-656.

[19] M. JEAN, F. JOURDAN et B. TATHI, Numerical dynamics for the simulation of deep drawing, *Proc. of IDDRG 94-18th*, Biennial Congress, May 16-17, 1994, Lisboa.

[20] M. JEAN et G. TOUZOT, Implementation of Unilateral Contact and Dry Friction in Computer Codes Dealing with Large Strain Stress Elastoplastic Problems, *Journal de Mécanique Théorique et Appliquée*, Special Issue: Numerical Method in Mechanics of Contact Involving Friction, 1988, p. 145-160.

[21] J. KATZENELSON, An Algorithm for Solving Nonlinear Resistive Networks, *Bell System Tech. J.*, 44, 1965.

[22] T. LADREYT, Modélisation bidimensionnelle et simulation numérique des processus de mise en forme des tôles dans l'industrie automobile, *Thèse de Doctorat*, INPG, Grenoble, 1992.

[23] C. E. LEMKE, Bimatrix Equilibrium Points and Mathematical Programming, *Management Science*, 4, 1978, p. 61-73.

[24] L. MENESES L., Modélisation tridimensionnelle et simulation numériques des processus de mise en forme ; application à l'emboutissage des tôles métalliques, *Thèse*, Université de Coïmbra, 1995.

[25] J. J. MOREAU, Application of convex analysis to some problems of dry friction, *Trends of Pure Mathematics Applied to Mechanics*, Zorski ed., 1979.
[26] J. J. MOREAU, Fonctionnelles convexes, Séminaire sur les Équations aux Dérivées Partielles, Collège de France, Paris, 1967.

[27] J. NECAS, Introduction to the Theory of Nonlinear Elliptic Equations, *Teubner-Texte zur Mathematik*, Band 52, 1982.

[28] J. ORTEGA et W. RHEINBOLT, Iterative Solutions of Nonlinear Equations in Several Variables, *Academic Press*, New York, 1970.

[29] J. S. PANG, Newton's Methods for B-Differentiable Equations, *Math. Op. Res.*, 15, n° 2, 1990, p. 311-341.

[30] M. RAOUS, P. CHABRAND et F. LEBON, Numerical Methods for Frictional Contact Problems and Applications, *Journal de Mécanique Théorique et Appliquée*, Special Issue: Numerical Method in Mechanics of Contact Involving Friction, 1988, p. 111-128.

[31] C. TEODOSIU, P. ALART et T. LADREYT, Développement d'une version implicite d'un logiciel Éléments Finis dédié à la simulation bidimensionnelle du retour élastique en emboutissage, Rapport de fin de contrat Renault, Février 1992.