# Un Algorithme de Gradient Conjugué Projeté Préconditionné pour la Résolution de Problèmes Unilatéraux

Nicolas Tardieu, Fabien Youbissi, Eric Chamberland

► To cite this version:

Nicolas Tardieu, Fabien Youbissi, Eric Chamberland. Un Algorithme de Gradient Conjugué Projeté Préconditionné pour la Résolution de Problèmes Unilatéraux. 8e Colloque national en calcul des structures, CSMA, May 2007, Giens, France. hal-01504125

HAL Id: hal-01504125

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

Submitted on 8 Apr 2017

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 CC0 1.0 - Universal - International License
# Un Algorithme de Gradient Conjugué Projeté Préconditionné pour la Résolution de Problèmes Unilatéraux

Nicolas Tardieu * — Fabien Youbissi ** — Éric Chamberland **

* EDF-R&D, 1 av du Gal de Gaulle, 91128 Clamart Cedex, France
nicolas.tardieu@edf.fr

** GIREF, Université Laval, Québec (Québec), Canada G1K7P4
{fabien.youbissi, eric.chamberland}@giref.ulaval.ca

**RÉSUMÉ.** L'algorithme du Gradient Conjugué Projeté (GCP) est issu du domaine de l'optimisation convexe où il est utilisé pour résoudre des problèmes de minimisation sous contraintes inégalités linéaires. Il se prête donc bien au cas du contact unilatéral. Appliqué à la forme duale du problème de contact et associé à un préconditionneur de Dirichlet, il exhibe un comportement très performant notamment quand on raffine le maillage.

**ABSTRACT.** The Projected Conjugate Gradient method comes from the field of convex optimization, where it is used to solve minimization problems with linear inequality constraints. Thus, it is well suited for unilateral contact. Utilized on the dual form of the contact problem and associated to a Dirichlet preconditionner, it exhibits a performing behaviour specially as the mesh size decreases.

**MOTS-CLÉS :** Contact, Gradient Conjugué, Préconditionneur.

**KEYWORDS:** Contact, Conjugate Gradient, Preconditionner.
# 1. Introduction

L'algorithme du Gradient Conjugué Projeté (GCP) est issu du domaine de l'optimisation convexe où il est utilisé pour résoudre des problèmes de minimisation sous contraintes inégalités linéaires. Des travaux précédents ont été réalisés portant sur l'utilisation de cet algorithme dans le cadre des problèmes de contact. La spécificité du présent travail réside dans l'application de cet algorithme à des cas de solides déformables (par opposition à (Renouf et al., 2005)), dans l'utilisation de la formulation duale du problème de contact (par opposition à (H.-O.May, 1986) et (G.Dilintas et al., 1988)) et surtout dans la définition d'un préconditionneur basé sur des problématiques mathématiques résolues par des considérations mécaniques.

On considère un problème de contact sur une partie $\Gamma_C$ de la frontière d'un solide élastique $\Omega$, dont la résolution s'exprime sous la minimisation suivante :

$$ \min_{u \cdot n \ge \text{Gap}} \mathcal{J}(u) \qquad [1] $$

Où :

- $\mathcal{J}$ est la fonctionnelle d'énergie potentielle

- $u$ est le champ de déplacement

- $n$ est la normale extérieure sur $\Gamma_C$

- Gap est la fonction écart

Il est classique d'introduire un champ de multiplicateurs de Lagrange $\lambda$ de manière à relaxer la contrainte sur le champ de déplacements. Après discrétisation par éléments finis, des manipulations algébriques permettent d'exprimer cette résolution sous sa forme duale :

$$ \min_{\lambda \ge 0} \mathcal{H}(\lambda) \qquad [2] $$

Où :

- $\mathcal{H}(\lambda) = \frac{1}{2}\lambda^T B K^{-1} B^T \lambda - \lambda^T B K^{-1} + \frac{1}{2} F^T K^{-1} F$

- $K$ est la matrice de rigidité du système

- $B$ est l'opérateur de jeu sur $\Gamma_C$

- $F$ est le second membre qui contient les chargements

L'intérêt d'utiliser cette expression est que l'opérateur obtenu présente un meilleur conditionnement que son pendant primal. Comme l'algorithme envisagé est de nature itérative, cet élément prend toute son importance. Le problème (2) peut être réécrit en supprimant la contrainte explicite de positivité sur $\lambda$ en ajoutant à la fonctionnelle l'indicatrice $\mathcal{I}_{\Delta^+}$ du convexe $\Delta^+ = \{\lambda \ge 0\}$. Pour la définition précise de la positivité du multiplicateur, nous renvoyons à (Hild, 2000).

$$ \min_{\lambda} \mathcal{L}(\lambda) = \min_{\lambda} (\mathcal{H}(\lambda) + \mathcal{I}_{\Delta^{+}}(\lambda)) \qquad [3] $$
## 2. Algorithme de résolution

Comme (3) est un problème de minimisation sans contraintes, on peut le résoudre par l'algorithm du GCP en prenant garde à la non-différentiabilité de la fonctionnelle qui nécessite la manipulation de sous-gradients. En outre, de manière à obtenir une direction de descente efficace, on cherche à déterminer l'élément de norme minimale du sous-gradient $\partial \mathcal{L}$ (M.S.Bazaraa et al., 1993). On obtient alors :

$$ \underset{j \in \partial \mathcal{L}}{\operatorname{argmin}} \|j\| = \begin{cases} \max(Bu - \text{Gap}, 0) & \text{si } \lambda = 0 \\ Bu - \text{Gap} & \text{sinon} \end{cases} \quad [4] $$

On peut utiliser ce résultat pour écrire un algorithme de descente très classique où l'introduction de la conjugaison des directions de recherche accélère sensiblement la convergence. On propose en Figure 1 l'algorithm en question dont l'écriture est basée sur la formulation primale du problème de contact. On insiste sur le fait que cet algorithme réalise bien la minimisation de la fonctionnelle $\mathcal{L}(\lambda)$ mais il manipule les grandeurs primales qui sont communément accessibles dans un code de calcul par éléments finis. A la lecture de cet algorithme, on peut remarquer que :

- la phase de recherche linéaire est non-admissible au sens où la mise à jour du multiplicateur de Lagrange $\lambda$ peut le faire sortir du domaine admissible ce qui est corrigé par une phase de projection. Il est très facile de définir une recherche linéaire admissible en limitant la valeur du pas d'avancement $\rho$.

- il faut prendre garde à appliquer la conjugaison uniquement si l'ensemble des contraintes actives ne change pas d'une itération à l'autre sans quoi la conjugaison perd sa pertinence mathématique.

- l'algorithm trouve sa performance dans l'utilisation d'un solveur direct pour la résolution des systèmes linéaires. La matrice $K$ est alors déjà factorisée et chaque résolution se limite donc à une descente-remontée.

## 3. Préconditionneur

Comme il est usuel dans le cas de l'utilisation de méthodes itératives, il est souhaitable de préconditionner le problème afin d'en accélérer la résolution (Saad, 2005). Le préconditionneur proposé est basé sur des remarques sur les espaces fonctionnels auxquels appartiennent les différents objets manipulés dans l'algorithm. Durant la phase de recherche linéaire, le multiplicateur de Lagrange $\lambda$ est mis à jour par la formule $\lambda^k + \alpha W^k$. Or, ces deux vecteurs n'appartiennent pas aux mêmes espaces fonctionnels : $\lambda^k \in H^{-1/2}(\Gamma_C)$ tandis que $W^k \in H^{1/2}(\Gamma_C)$. Il manque ici un opérateur de Steklov-Poincaré pour transférer la direction de recherche $W^k$ d'un espace à l'autre.

En s'inspirant de la technique de relèvement développée dans (Fortin et al., 2005), on propose de résoudre le problème annexe suivant où $\tilde{W}^k$ est le préconditionné de la direction de descente $W^k$ du GCP :

$$ \left\{ \begin{array}{l} Ku + B_{act}^T \tilde{W}^k = 0 \\ B_{act}u = W^k \end{array} \right. \qquad [5] $$
◇ $k = 0, \lambda = 0, \|r^0\| = 2 * \text{Critère}$

◇ Tant que $\|r^k\| > \text{Critère}$

—— DIRECTION DE RECHERCHE ET PROJECTION ——

$$- \text{Résoudre } K u^k = F - B^T \lambda^k$$

$$- \text{Soit } r^k = B u^k - \text{Gap, si } \lambda_i^k > 0 \text{ alors } r_i^k = \max(r_i^k, 0)$$

—— CONJUGAISON ——

$$- \text{Si } k > 0 \text{ et si l'état de contact est identique au précédent, alors } \gamma = \frac{r^{k+1} - r^k}{r^k}$$

$$- \text{Si } k > 0, W^k = r^k + \gamma W^{k-1}$$

—— RECHERCHE LINÉAIRE ——

$$- \text{Résoudre } K w^k = B^T W^k$$

$$- \alpha = \frac{W^k - W^{k-1}}{B^T w^k}$$

$$- \lambda_i^k = \max(\lambda_i^k + \alpha W_i^k, 0)$$

$$- k = k + 1$$

◇ Fin tant que

**Figure 1. Algorithme de GCP**

Où $B_{act}$ est l'opérateur de jeu limité aux liaisons actives de l'itération courante.

Ce problème linéaire est résolu en suivant la même démarche présentée précédemment en exploitant sa forme duale. On manipule alors des gradients et non des sous-gradients pour alimenter un algorithme de Gradient Conjugué. On remarque que les résolutions de systèmes linéaires impliquent la matrice $K$ qui est déjà factorisée.

On notera que ce préconditionneur est de nature Dirichlet tandis que le solveur développé est de nature Neumann, comme c'est le cas des techniques FETI en décomposition de domaine (Farhat et al., 1991).

**4. Illustration**

On propose d'illustrer le comportement de l'algorithme proposé sur le cas d'un problème de Hertz (écrasement d'une sphère sur un plan rigide) modélisé en axisymétrie. On compare le nombre d'itérations ainsi que le rapport des temps CPU requis pour obtenir la solution avec et sans préconditionneur quand on raffine le maillage. On utilise pour cela un maillage à base d'éléments triangulaires linéaires que l'on subdivise en découpant leurs arrêtes en 2. On passe ainsi de $10^3$ à $2 \cdot 10^6$ degrés de liberté. Dans les figures 2(a) et 2(b), H représente la taille initiale des éléments tandis que h est leur taille courante. On remarque sur cette figure que :
<figure><img src="image_2.png" /><figcaption>Figure 2. Comparaison des algorithmes avec et sans préconditionneur</figcaption></figure>

– l’utilisation du préconditionneur réduit toujours le nombre d’itérations requis pour atteindre la convergence.

– l’augmentation du nombre d’itérations avec la taille du problème quand on utilise le préconditionneur est très limitée : on passe de 5 à 12 itérations tandis que sans préconditionneur, on passe de 79 à 1823 itérations, ce qui autorise à suppposer son optimalité.

– en terme de temps CPU, l’utilisation du préconditionneur est toujours intéressante et cet intérêt croit avec la taille du modèle.
**5. Conclusions**

Dans ce papier a été présenté un algorithme de Gradient Conjugué Projeté Préconditionné basé sur la formulation duale d'un problème de contact unilatéral. Associé à un préconditionneur de nature Dirichlet, il exhibe des performances intéressantes quand la taille du maillage augmente. Une autre qualité majeure de cet algorithme est le très faible surcoût mémoire et disque du traitement du contact (pas de matrice tangente ou de factorisation supplémentaire). Si le problème sans contact "passe" sur une machine donné, il est acquis qu'il en sera de même pour le problème avec contact. Cette constatation a été vérifiée lors de l'implantation de cet algorithme et de son préconditionneur dans le logiciel de simulation *Code_Aster* (Code_Aster, 2007), où des simulations industrielles impliquant plusieurs millions de degrés de liberté ont pu être menées à bien.

**6. Bibliographie**

Code_Aster, « http://www.code-aster.org », 2007.

Farhat C., Roux F.-X., « A method of finite element tearing and interconnecting and its parallel solution algorithm », *Int. J. Numer. Meth. Engng.*, vol. 32, n° 6, p. 1205-1227, 1991.

Fortin M., Tardieu N., Chamberland E., « Un algorithme de sous-gradient pour le traitement du contact frottant », *Actes du Colloque National en Calcul de Structures de Giens*, Éditions Hermès, Paris, 2005.

G.Dilintas, P.Laurent-Gengoux, D.Trystram, « A conjugate projected gradient method with preconditioning for unilateral contact problems », *Comput. Struct.*, vol. 29, n° 4, p. 675-680, 1988.

H.-O.May, « The conjugate gradient method for unilateral problems », *Comput. Struct.*, vol. 22, n° 4, p. 595-598, 1986.

Hild P., « Numerical implementation of two nonconforming finite element methods for unilateral contact », *Comput. Methods Appl. Mech. Engrg.*, vol. 184, n° 1, p. 99-123, 2000.

M.S.Bazaraa, H.D.Sherali, M.Shetty, *Nonlinear Programming - Theory and Algorithms*, John Wiley and Sons, Philadelphia, 1993.

Renouf M., Alart P., « Conjugate gradient type algorithms for frictional multi-contact problems : applications to granular materials », *Comput. Methods Appl. Mech. Engrg.*, vol. 194, n° 18-20, p. 2019-2041, 2005.

Saad Y., *Iterative Methods for Sparse Linear Systems*, SIAM, Philadelphia, 2005.