# Solutions analytiques de problèmes de contact en grandes déformations

Eric Chamberland, André Fortin, Nicolas Tardieu

► To cite this version:

Eric Chamberland, André Fortin, Nicolas Tardieu. Solutions analytiques de problèmes de contact en grandes déformations. 9e Colloque national en calcul des structures, CSMA, May 2009, Giens, France. hal-01420721

HAL Id: hal-01420721

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

Submitted on 20 Dec 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
# Solutions analytiques de problèmes de contact en grandes déformations

E. Chamberland¹, A. Fortin¹, N. Tardieu²

¹ GIREF, Université Laval
Québec (Québec), G0V 0A6, CANADA
{Eric.Chamberland,Andre.Fortin}@giref.ulaval.ca

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

**Résumé** — Il est impératif pour l'ingénieur de s'assurer de la validité de ses outils de modélisation numérique et il lui revient impérativement de vérifier s'ils donnent la solution attendue pour un problème donné. Ceci peut être fait par l'utilisation de solutions analytiques lorsqu'elles existent ou par la méthode des solutions manufacturées (MSM). Nous proposons dans ce travail d'étendre la MSM aux problèmes de contact en grandes déformations pour lesquels, à notre connaissance, il n'existe pas de solutions analytiques non triviales. Nous montrerons comment l'utiliser comme un outil de vérification pour l'ingénieur.

**Mots clés** — Solutions manufacturées, contact, grandes déformations.

## 1 Sur la vérification de logiciels

La vérification de logiciel est un aspect très important du génie logiciel. On y inclut les tests de non-régression, où l'on compare les résultats d'un calcul dans l'état actuel du logiciel avec des résultats obtenus et vérifiés précédemment, ainsi que les compilations multi-architectures, multi-compilateurs effectuées sur une base régulière. Ce premier type de vérification assure un niveau de qualité minimale du logiciel : tant que les résultats des tests ne changent pas, le logiciel conserve la qualité qu'il avait lors de la production des résultats étalons. L'ensemble de ces vérifications est présenté comme faisant partie de l'assurance qualité logiciel par Roy [15].

Deux types de problèmes peuvent toutefois persister malgré ces tests quotidiens. Il se peut en effet que certaines parties du logiciel ne soient pas couvertes par ces tests. Il s'agit alors d'une partie du logiciel qui n'est pas vérifiée de façon régulière et qui peut arrêter de donner les bons résultats. Pire encore, et c'est le deuxième type de problèmes, il se peut qu'elle n'ait jamais donné de bons résultats. En modélisation numérique, et plus spécifiquement en éléments finis, les problèmes possèdent rarement des solutions analytiques connues. Comment peut-on alors valider des résultats pour lesquels il n'existe aucune solution analytique connue ?

La MSM, introduite par Steinberg et Roach [18], est utilisée depuis longtemps en mécanique des fluides (voir Pelletier et al. [12]). Salari et Knupp [16] présentent un rapport dans lequel plusieurs aspects de la vérification de code sont abordés. On y retrouve la MSM avec plusieurs détails intéressants, dont une série de 21 erreurs différentes pouvant ou non être détectées par cette méthode. Roache présente aussi la méthode en détails dans [14] et a publié un livre sur le sujet (voir [13]).
L'intérêt principal de la MSM est qu'elle constitue une méthodologie systématique pour vérifier les calculs obtenus par des modélisations pour lesquels aucune solution analytique n'est connue. Ceci vient combler un besoin énorme, notamment dans les problèmes de mécanique des solides en grandes déformations. Avec la MSM, il sera donc possible de vérifier que le logiciel donne le bon résultat. Cette technique a été jugée suffisamment fiable pour faire partie de la procédure standard de vérification de logiciels en mécanique des fluides (voir [1]) et en mécanique des solides (voir [2]). Il est intéressant de noter que l'« American Institute of Aeronautics and Astronautics » (AIAA) l'a incluse dans ses normes en 1998, alors que la « American Society of Mechanical Engineers » (ASME) l'a fait en 2006, soit beaucoup plus récemment. À notre connaissance, il n'existe aucune publication à ce jour qui fait usage de la MSM en mécanique des solides en 2D ou en 3D bien qu'il y ait un exemple qui paraît dans le résumé du guide de l'ASME pour un cas 1D ( [17]). Un article est en préparation sur le sujet (voir Chamberland et al. [7]), issu du travail de maîtrise de Chamberland [6].

## 2 Procédure de vérification de logiciels de contact en grandes déformations

Plusieurs développements informatiques sont nécessaires pour la prise en compte des problèmes de contact dans un code éléments finis. Cela inclut au minimum une intégrale de surface pour imposer la condition de non-pénétration. Un opérateur de projection sur les surfaces voissines admissibles au contact sera aussi nécessaire pour calculer le point projeté ($x_2$), la normale ($n_2$) et en obtenir la distance ($g = (x_1 - x_2) \cdot n_2$). Le dernier élément important à mettre en place est l'extraction de la pression de contact résultant du calcul. Une fois ces éléments informatiques programmés, vient le moment de vérifier leur exactitude.

Pour le contact déformable-déformable, Taylor et Papadopoulos [20] proposent un test permettant de vérifier la bonne transmission des pressions de contact entre deux domaines déformables. Crisfield [9] reprend ce test et propose une modification pour les problèmes de contact déformable-déformable avec une interface de contact courbe. Dans les deux cas, le problème proposé possède une solution analytique linéaire dont on peut extraire les composantes des contraintes de Cauchy. Toutefois, dans le cas des interfaces courbes, on ne montre pas comment calculer la contrainte normale ($\sigma_{nn}$) du bord en contact. C'est précisément cette contrainte qui est la plus intéressante. De plus, la solution linéaire proposée est insuffisante pour vérifier des éléments quadratiques. Il serait en effet intéressant de pouvoir construire une solution quadratique et en extraire les contraintes normales au bord de contact pour vérifier que les éléments quadratiques représentent de façon exacte cette pression de contact.

Dans plusieurs travaux, on utilise le maillage le plus fin que l'on puisse se permettre pour obtenir une solution de référence. Il est alors possible d'étudier le taux de convergence sur une une série de maillages plus grossiers. Il est aussi possible de comparer les résultats d'un logiciel de contact avec un logiciel reconnu. Par contre, comme le mentionne Roache [14], les éditeurs de logiciels ne publient pas leurs propres résultats de vérifications et peut-être avec raison !

Nous considèrerons donc ici les problèmes de contact en grandes déformations dans une formulation mixte en déplacements $u$ et en force de contact $\lambda$. Il serait donc intéressant de pouvoir comparer, sur une base neutre et libre d'erreur numérique, différents choix d'espaces de discrétisation $u_h - \lambda_h$. Ceci permetterait, par exemple, une comparaison juste entre les triangles et les quadrangles. Cette base neutre devrait aussi permettre de comparer les solutions obtenues avec des maillages adaptés. C'est dans une solution analytique que l'on retrouve tous les éléments intéressants pour faire ces comparaisons.

Il est important de mentionner ici qu'il n'a pas été possible de trouver une procédure de vérification de calcul de contact en grandes déformations avec solution analytique pour les éléments de
<figure><img src="image_2.png" /><figcaption>Figure 1 – (a) Domaine initial et déformé avec obstacle rigide, (b) Contrainte normale ($\sigma_{nn}J_s$) sur le bord $\Gamma_C^0$.</figcaption></figure>

tout ordre d'approximation. C'est pourquoi nous voudrions introduire la MSM qui nous procurera des solutions analytiques mais qui seront d'ordre plus élevé que celles présentées dans [20], [9] et [10] et pourront donc servir à vérifier l'exactitude des éléments de tout ordre.

Une procédure complète de vérification devrait inclure les items suivants :

1. Vérification de l'obtention de la solution exacte ($u_h, \sigma_h, u_h \in \Gamma_C, \lambda_h$) pour une solution analytique comprise dans l'espace d'approximation $u_h - \lambda_h$.

2. Vérification du taux de convergence sur $u_h, \sigma_h, u_h \in \Gamma_C$ et $\lambda_h$ pour une solution analytique qui n'est pas comprise dans l'espace d'approximation.

Pour les grandes déformations, une formulation Lagrangienne totale en déplacement $u$ sera utilisée. Le contact sera traité par l'introduction d'un multiplicateur de Lagrange $\lambda_h$ et résolu par la méthode du gradient conjugué projeté (voir Tardieu et al. [19]). On note $\sigma_{nn}^-$ et $\sigma_{nn}^+$ les parties respectivement négative et positive de $\sigma_{nn}$. À convergence, le multiplicateur approchera $\sigma_{nn}^- J_s$, ce qui est lié à notre choix d'intégrer la contribution du jacobien surfacique $J_s$ dans $\lambda_h$. Ceci évite de linéariser la variation de $J_s$, ce qui introduirait un terme non-symétrique dans la matrice tangente, comme dans le cas de la pression suiveuse (voir Bonet et Wood [4]). Sous forme variationnelle, on a :

$$ 
\begin{aligned}
\int_{\Omega^0} \Pi : \nabla_X w \, dV - \int_{\Gamma_C^0} \lambda_h w \cdot n_2 \, dA &= \int_{\Omega^0} r_0 \cdot w \, dV + \int_{\Gamma_N^0} h_0 \cdot w \, dA + \int_{\Gamma_C^0} k_0 \cdot w \, dA \quad \forall w \\
\int_{\Gamma_C^0} (\mu_h - \lambda_h) g(u) \, dA &\le 0 \qquad \forall \mu_h \le 0
\end{aligned}
\quad (1) $$

où $g(u)$ est la distance entre les corps (le « gap »). De manière équivalente, on a :

$$ 
\left\{ 
\begin{array}{rcl@{\quad}l}
-\nabla_X \Pi & = & r_0 & \text{dans } \Omega^0 \\
u_h & = & u_0 & \text{sur } \Gamma_D^0 \\
\Pi \cdot N & = & h_0 & \text{sur } \Gamma_N^0
\end{array}
\right. 
\qquad \left. 
\begin{array}{rcl}
\Pi \cdot N - \lambda_h n_2 & = & k_0 \\
\lambda_h & \le & 0 \\
g & \ge & 0 \\
\lambda_h g & = & 0
\end{array}
\right\} \text{ sur } \Gamma_C^0 $$

Le choix de l'espace de discrétisation pour $u_h$ et $\lambda_h$ est conditionné par une condition de type $inf - sup$ (voir Bathe et Brezzi [3], Coorevits et al. [8]). Plusieurs choix d'espaces de discrétisation satisfaisant ou non la condition $inf - sup$ seront présentés dans cette étude.

### 3 Solution manufacturée pour le contact

Pour utiliser la méthode des solutions manufacturées pour le contact en grandes déformations, il faut suivre les étapes suivantes, tout en observant certains critères :
1. Choisir un domaine de modélisation $\Omega^0$ ainsi qu'une solution manufacturée $u$;

2. Choisir la loi de comportement et les propriétés matérielles (le module d'Young $E_Y$ et le coefficient de Poisson $\nu$ pour un matériau hookéen);

3. Calculer $F$, le tenseur de Green-Lagrange $E = \frac{1}{2}(F^t \cdot F - I)$ et les tenseurs de Piola-Kirchoff $S$ et $\Pi = F \cdot S$;

4. Calculer la divergence du premier tenseur de Piola-Kirchoff : $r_0 = -\nabla_X \cdot \Pi$;

5. Avec la normale $N$ à la géométrie non déformée, calculer $h_0 = \Pi \cdot N$ pour l'imposition des conditions aux limites de type Neumann sur $\Gamma_N^0$;

6. Calculer de tenseur des contraintes de Cauchy $\sigma = \frac{1}{J} F \cdot S \cdot F^\top$;

7. Calculer la normale à l'état déformé $n = \frac{F^{-\top} \cdot N}{\|F^{-\top} \cdot N\|}$ et $J_s = \det(F) \frac{\|F^{-\top} \cdot N\|}{\|F^{-\top} \cdot N\|}$

8. Calculer les contraintes normale $\sigma_{nn} = (\sigma \cdot n) \cdot n$ et tangentielle $\sigma_t = \sigma \cdot n - \sigma_{nn} n$;

9. Calculer $k_0 = (\sigma_{nn}^+ n + \sigma_t) J_s$ pour l'imposition des conditions aux limites sur $\Gamma_C^0$

10. Choisir $\Gamma_C^0$ et vérifier que $\sigma_{nn} < 0$ sur une partie de celui-ci;

11. Créer l'obstacle rigide ($\Omega_R$) correspondant à la déformée analytique de $\Gamma_C^0$ où $\sigma_{nn} < 0$;

La figure 1(a) présente un exemple de problème assez complexe pour lequel une solution manufacturée a été construite. Le domaine intial est le carré ($\Omega^0 = [0, 1] \times [0, 1]$). La solution manufacturée choisie est celle de l'équation (2) et une feuille de calcul (voir [5]) a été utilisée pour générer l'ensemble des variables nécessaires au calcul ($r_0, h_0, k_0$) et aux comparaisons ($\sigma_{nn}, J_s$).

$$ 
\begin{aligned}
u_x &= \frac{1}{50} \left( X - \frac{101}{4} + 7/2 \sqrt{51} \right)^2 + \frac{1}{50} Y^5 X^3 & E_Y &= X + 1 \\
u_y &= \frac{1}{50} Y^4 \left( X - 1/2 \right)^3 - \frac{1}{50} X & v &= \frac{3}{10}
\end{aligned} 
\quad (2) $$

Cette solution déformera le carré en une parabole sur la partie inférieure du carré ($\Gamma_C^0$). Sur cette partie parabolique, la contrainte normale ($\sigma_{nn}$) est inférieure à 0 pour $X \in [0, a]$ (zone où il y aura contact) et plus grande que zéro pour $X \in [a, 1]$ avec $a = 1/4$. La valeur analytique de ($\sigma_{nn}, J_s$) sur la frontière $\Gamma_C^0$ est présentée à la figure 1(b). L'obstacle rigide ($\Omega_R$) a été représenté par un seul élément $P_2$ sur l'intervalle $X = [-1, a]$ et déformé par les composantes $u$ de la solution manufacturée (2). Une série de maillages composés de quadrangles a été générée en divisant le domaine initial en $n$ subdivisions pour $n = [4, 8, 16, 32, 64, 128, 256, 512]$. Les sommets des maillages résultants ont été aléatoirement perturbés d'un facteur $h/10$ pour ne pas avoir d'effets de maillages trop réguliers. Tous les maillages de quadrangles ont été ensuite divisés en triangles. Tout est maintenant en place pour lancer les calculs et faire les comparaisons.

# 4 Résultats numériques

Il est intéressant de regarder la convergence de plusieurs variables. Nous nous contenterons de commenter ici la convergence de $\lambda_h$ sur $\Gamma_C^0$, car c'est la quantité d'intérêt pour un problème de contact. Bien entendu, grâce à la solution manufacturée, l'erreur sur $u_h$ peut être calculée en norme $H^1$ (voir figure 2) et $L^2$ dans $\Omega^0$ et aussi en norme $L^2$ sur le bord $\Gamma_C^0$ (voir figure 3) afin de compléter les vérifications du logiciel.

Sur la figure 4, on peut voir l'erreur relative en norme $L^2$ sur la pression de contact $\lambda_h$. On peut observer que $\lambda_h$ converge bien vers $\sigma_{nn}^- J_s$ (et non vers $\sigma_{nn}^-$). On distingue aussi trois paires de courbes :

1. Les éléments les moins précis sont les $P_2 - P_0$ et $Q_2 - P_0$ pour lesquels la convergence est limitée à $O(h)$.
<figure><img src="image_3.png" /><figcaption>Figure 2 – Errur relative en norme H¹ (||u − u_h||₁,Ω_h / ||u||₁,Ω_h) en fonction de la taille des éléments (h).</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3 – Errur relative en norme L² (||u − u_h||₀,Γ_c^0 / ||u||₀,Γ_c^0) en fonction de la taille des éléments
(h).</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 4 – Erreur relative en norme $L^2 (||\sigma_{nn}^- J_s - \lambda_h||_{0,\Gamma_C^0} / ||\sigma_{nn}^- J_s||_{0,\Gamma_C^0})$ en fonction de la taille des éléments (h).</figcaption></figure>

2. Les éléments $P_1 - P_1$ et $Q_1 - P_1$ ont une convergence limitée à $O(h^{1.24})$ et sont d'un ordre de grandeur plus précis que les précédents.

3. Les éléments $P_2 - P_1$ et $Q_2 - P_1$ sont limités dans leur convergence à $O(h^{1.42})$, mais restent les plus précis de tous.

Dans presque tous les calculs, on ne voit aucune différence marquée entre les triangles et les quadrangles sauf pour le $Q_1 - P_1$ en comparaison avec le $P_1 - P_1$ si on regarde l'erreur sur le déplacement.

Tel que mentionné dans [10], une difficulté majeure des problèmes de contact est l’intégration exacte du gap. On voit dans les taux de convergence présentés ici l’impact d’une telle difficulté. Si on s’assure d’avoir un sommet placé exactement là où le changement de signe de $\sigma_{nn}$ se produit (soit à $X = 1/4$), on retrouve un taux de convergence $O(h^2)$ (voir figure 5) sur $\lambda_h$ pour les éléments $P_2 - P_1$ et $Q_2 - P_1$.

Concernant les taux de convergence, Hild et Laborde [11] présentent des résultats théoriques et numériques pour quelques espaces d'approximations. Cependant, le taux de convergence nécessite le calcul d'erreur sur $\lambda_h$ en norme $H^{-1/2}$ et il est difficile de donner l'équivalent en norme $L^2$.

# 5 Conclusion

L'originalité de ce travail réside dans l'utilisation de la MSM pour les problèmes de contact en grandes déformations. La vérification de logiciel étant une partie incontournable de tout développement informatique, il est maintenant possible de le faire d'une manière plus complète qu'avec les méthodes actuellement utilisées, grâce à la MSM.

L'utilisation de la MSM est aussi possible pour les problèmes de contact déformable-déformable. Il serait donc possible de compléter le patch test de contact proposé dans [20] et [9] avec des solutions analytiques $P_2$. Un autre intérêt de la MSM est de comparer différentes méthodes de résolution du contact (pénalisation, Lagrangien augmenté, GCP) entre elles. Il est aussi possible de faire une comparaison pour les formulations incompressibles, où l'on verrait apparaître la pression (p) comme inconnue et donc un triplet de variables à vérifier : $u_h - p_h - \lambda_h$. La suite de ce travail sera de présenter comment choisir une solution permettant de vérifier les algorithmes de contact frottant.
<figure><img src="image_6.png" /></figure>

Figure 5 – Erreur relative en norme $L^2 (||\sigma_{nn}^- J_s - \lambda_h||_{0,\Gamma_C^0} / ||\sigma_{nn}^- J_s||_{0,\Gamma_C^0})$ en fonction de la taille des éléments (h) pour les maillages dont un sommet est à $X = 1/4$.

Références

[1] AIAA. *Guide for the Verification and Validation of Computational Fluid Dynamics Simulations*. American Institute for Aeronautics and Astronautics, 1998.

[2] ASME. *Guide for Verification and Validation in Computational Solid Mechanics*. American Society of Mechanics Engineers, 2006.

[3] K.J. Bathe et F. Brezzi. Stability of finite element mixed interpolation for contact problems. *Rend. Mat. Acc. Lincei s.*, 9(12):159-166, 2001.

[4] J. Bonet et D. R. Wood. *Nonlinear continuum mechanics for finite element analysis*. Cambridge, 1997.

[5] É. Chamberland. http://www.giref.ulaval.ca/~ericc/mms.

[6] É. Chamberland. Solutions manufacturées en grandes déformations hyper-élastiques. Master's thesis, Université Laval, 2009.

[7] É. Chamberland, A. Fortin, et M. Fortin. A performance study of finite element discretizations using manufactured solutions for hyperelastic materials in large deformations. 2008. En préparation.

[8] P. Coorevits, P. Hild, K. Lhalouani, et T. Sassi. Mixed finite element methods for unilateral problems : convergence analysis and numerical studies. *Math. Comput.*, 71(237) :1-25, 2002.

[9] M. A. Crisfield. Re-visiting the contact patch test. *International Journal for Numerical Methods in Engineering*, 48(3) :435-449, 2000.

[10] N. El-Abbasi et K.-J. Bathe. Stability and patch test performance of contact discretizations and a new solution algorithm. *Computers & Structures*, 79(16) :1473-1486, 2001.

[11] P. Hild et P. Laborde. Quadratic finite element methods for unilateral contact problems. *Appl. Numer. Math.*, 41(3) :401-421, 2002.

[12] D. Pelletier, É. Turgeon, D. Lacasse, et J. Borggaard. Adaptivity, Sensivity and Uncertainty : Towards Standards in CFD. In *39th AIAA Aerospace Sciences Meeting and Exibit*, 2001.

[13] P. J. Roache. *Verification and validation in computational science and engineering*. Hermosa, Albuquerque, 1998.
[14] P. J. Roache. Code verification by the method of manufactured solutions. ASME, 124:4-10, 2002.

[15] C. J. Roy. Review of code and solution verification procedures for computational simulation. *J. Comput. Phys.*, 205(1):131-156, 2005.

[16] K. Salari et P. Knupp. Code verification by the method of manufactured solutions. *Sandia National Laboratories*, 2000.

[17] L. E. Schwer. An overview of the PTC 60/V&V 10: guide for verification and validation in computational solid mechanics: Transmitted by l. e. schwer, chair ptc 60/v&v 10. *Eng. with Comput.*, 23(4):245-252, 2007.

[18] S. Steinberg et P.J. Roache. Symbolic manipulation and computational fluid dynamics. *Journal of Computational Physics*, 57(2):251-284, 1985.

[19] N. Tardieu, F. Youbissi, et É. Chamberland. Un algorithme de gradient conjugué projeté préconditionné pour la résolution de problèmes unilatéraux. *Comptes Rendus Mécanique*, 336(11-12):840-845, 2008.

[20] R.L. Taylor et P. Papadopoulos. In: *Wriggers P and Wagner W, editors, Computational methods in nonlinear mechanics*, chapter On a patch test for contact problems in two dimensions, pages 690-702. Springer, 1991.