# Formulation variationnelle régularisée par équations intégrales en élastoplasticité tridimensionnelle

B. BURGARDT<fnref n="1" />, M. BONNET<fnref n="2" />, A. LE VAN<fnref n="1" />

## Résumé

Cette communication concerne la développement d'une formulation variationnelle symétrique pour la résolution de problèmes élastoplastiques HPP, reposant sur les éléments de frontière. Un maillage d'éléments volumiques de la région potentiellement élastoplastique permet l'interpolation des déformations inélastiques. L'ensemble de la formulation, y compris l'opérateur tangent résultant du traitement implicite de l'intégration du comportement, est symétrique. Des difficultés liées à l'intégration numérique subsistent en l'état actuel de ce travail.

## Présentation

La méthode des équations intégrales est applicable à la résolution de problèmes mécaniques présentant des non-linéarités matérielles ou géométriques. Les représentations intégrales font alors apparaître des opérateurs intégraux de volume dont les densités représentent les déformations plastiques $\varepsilon^p$ ou le champ gradient du déplacement $u$. La méthode des éléments de frontière s'avère ainsi très intéressante pour résoudre certains problèmes spécifiques (milieux infinis, zone plastique confinée, ...).

Dans ce travail nous traitons de la résolution numérique de problèmes élastoplastiques tridimensionnels en petites perturbations. L'intégration en temps est réalisé au moyen d'un schéma de type Newton implicite. L'approche variationnelle [5] proposée pour la résolution en espace diffère du cheminement classique et largement répandu, qui consiste à utiliser une méthode de collocation.

L'approche variationnelle appliquée à la résolution d'équations intégrales en élastoplastique a été proposée par Polizzotto [2] et permet, contrairement à la méthode de collocation, d'établir des résultats de convergence et de stabilité [4]. Son principal avantage pratique est de faire intervenir des opérateurs matriciels qui sont généralement symétriques et de signe défini.

L'étude proposée ici exploite la régularisation des formulations intégrales en élastostatique [3], qui clarifie et facilite la mise en œuvre numérique. Considérons un domaine

<footnote n="1">¹ Laboratoire Mécanique et Matériaux - Ecole Centrale de Nantes, BP 92101 - 44321 Nantes cedex 3</footnote>

<footnote n="2">² Laboratoire de Mécanique des Solides - Ecole Polytechnique, 91128 Palaiseau cedex</footnote>
tridimensionnel $\Omega$ de frontière extérieure $\partial\Omega = S_u \cup S_t$ et comportant une zone potentiellement plastique $\Omega_p$ à l'extérieur de laquelle la déformation plastique $\varepsilon^p$ est supposée nulle. On montre que le système variationnel portant sur la variable $\mathbf{u}$ (déplacement, sur $S_t$), $\mathbf{t}$ (vecteur-contrainte, sur $S_u$) et $\varepsilon^p$ (sur $\Omega_p$) associé à l'écriture des conditions aux limites peut être mis sous la forme suivante, symétrique :

$$ \left\{ \begin{aligned} a_{uu}(\tilde{\mathbf{u}}, \mathbf{u}) + a_{ut}(\tilde{\mathbf{u}}, \mathbf{t}) &= \ell_u(\tilde{\mathbf{u}}) + a_{u\varepsilon}(\tilde{\mathbf{u}}, \varepsilon^p) \\ a_{tu}(\tilde{\mathbf{t}}, \mathbf{u}) + a_{tt}(\tilde{\mathbf{t}}, \mathbf{t}) &= \ell_t(\tilde{\mathbf{t}}) + a_{t\varepsilon}(\tilde{\mathbf{t}}, \varepsilon^p) \end{aligned} \right. \qquad (1) $$

où les formes linéaires $\ell_u$, $\ell_t$ contiennent les données aux limites et les fonctions-test $\tilde{\mathbf{u}}, \tilde{\mathbf{t}}$ sont respectivement définies sur $S_t$ et $S_u$. Les divers opérateurs ci-dessus s'expriment sous forme d'intégrales doubles et en termes du déplacement fondamental de Kelvin $U(\mathbf{x}, \mathbf{y}) = U^T(\mathbf{y}, \mathbf{x})$ (déplacement créé en $\mathbf{y}$ par une force ponctuelle unitaire appliquée en $\mathbf{x}$) et de quantités dérivées. Par exemple :

$$ a_{tt}(\tilde{\mathbf{t}}, \mathbf{t}) = \int_{S_t} \int_{S_t} \tilde{\mathbf{t}}(\mathbf{x}).U(\mathbf{x}, \mathbf{y}).\mathbf{t}(\mathbf{y})dS_ydS_x \qquad (2) $$

La représentation intégrale du champ de déformation totale $\varepsilon$ sous forme faible donne alors :

$$ \langle \tilde{\varepsilon}, C : \varepsilon \rangle_{\Omega_P} = \ell_{\varepsilon}(\tilde{\varepsilon}) + a_{\varepsilon u}(\tilde{\varepsilon}, \mathbf{u}) + a_{\varepsilon t}(\tilde{\varepsilon}, \mathbf{t}) + a_{\varepsilon \varepsilon}(\tilde{\varepsilon}, \varepsilon^p) \qquad (3) $$

avec les propriétés de symétrie :

$$ a_{\varepsilon u}(\tilde{\varepsilon}, \mathbf{u}) = a_{u \varepsilon}(\tilde{\mathbf{u}}, \varepsilon^p) \qquad (4) $$

$$ a_{\varepsilon t}(\tilde{\varepsilon}, \mathbf{t}) = a_{t \varepsilon}(\tilde{\mathbf{t}}, \varepsilon^p) \qquad (5) $$

La structure de l'algorithmme de résolution dépend de la nature de l'interpolation (continue ou non) de la déformation plastique $\varepsilon^p$ et de la déformation totale $\varepsilon$. Pour une interpolation définie élément par élément (interpolation non-continue), l'algorithmme résultant peut être identique à celui rencontré pour la méthode de collocation et la matrice tangente est une matrice diagonale par blocs, chaque bloc représentant l'opérateur tangent en un point de Gauss. Mais dans ce cas l'obtention de la matrice associée à la forme bilinéaire $a_{\varepsilon\varepsilon}$ est complexe et fait intervenir les discontinuités des déformations plastiques aux interfaces des éléments.

Nous avons préféré mettre en œuvre une interpolation continue des déformations. Dans ce cas, l'opérateur tangent obtenu pour le schéma de Newton associé n'est symétrique qu'après écriture sous forme faible des conditions de normalité de la loi d'écoulement plastique [2]. La matrice tangente prend alors une expression plus complexe et il est nécessaire de faire appel à une procédure d'assemblage à chaque itération et une réactualisation de l'opérateur $C : \frac{\partial \varepsilon^p}{\partial \varepsilon_{n+1}}$ aux points de Gauss [1].

**Calcul des contraintes élastiques.** L'écriture sous forme faible du champ $C : \varepsilon$ conduit à l'expression suivante :

$$ \int_{\Omega_p} \delta\varepsilon : C : \varepsilon_{n+1} dV = \ell_{\varepsilon}(\delta\varepsilon) + a_{\varepsilon u}(\delta\varepsilon, \mathbf{u}_{n+1}) + a_{\varepsilon t}(\delta\varepsilon, \mathbf{t}_{n+1}) - a_{\varepsilon \varepsilon}(\delta\varepsilon, \varepsilon^p) \qquad (6) $$
La forme discrétisée de ces équations mène alors au système non-linéaire suivant :

$$ M \cdot E_{n+1} = G_{n+1} + Q^t \cdot Y_{n+1} - \Xi \cdot E_{n+1}^p (E_{n+1}) \quad (7) $$

où $Y_{n+1}$ désigne le vecteur des densités surfaciques inconnues solution de :

$$ K \cdot Y_{n+1} = F_{n+1} + Q \cdot E_{n+1}^p \quad (8) $$

Après avoir éliminé $Y_{n+1}$ on aboutit à :

$$ M \cdot E_{n+1} = G'_{n+1} - \Xi' \cdot E_{n+1}^p (E_{n+1}) \quad (9) $$

$$ G'_{n+1} = G_{n+1} + Q^t \cdot K^{-1} \cdot F_{n+1} \quad (10) $$

$$ \Xi' = \Xi - Q^t \cdot K^{-1} \cdot Q \quad (11) $$

En l'absence de déformations plastiques, le vecteur $E_{n+1}$ est donné par :

$$ E_{n+1} = M^{-1} \cdot G'_{n+1} \quad (12) $$

Le vecteur $M^{-1} \cdot G'_{n+1}$ représente ainsi le champ de déformation élastique du milieu. En fait, on a déterminé une solution EF collant au mieux au sens de la norme énergie, au champ de contrainte défini par la représentation intégrale. Le fait de trouver de bonnes valeurs pour ce vecteur, permet de valider le calcul des quantités $Q$ et $G$.

## Exemple numérique

La validation en cours de la méthode proposée repose notamment sur des comparaisons avec la solution exacte du problème de torsion élastoplastique. On considère un bloc cylindrique d'axe $e_z$, de section circulaire (rayon $R = 1$ mm, longueur $l = 1$ mm) sollicité en torsion. Le matériau constitutif est supposé élastoplastique parfait (module d'Young $E = 200000$ MPa, coefficient de poisson $\nu = 0,3$, limite élastique $\sigma_0 = 240$ MPa). La face $z = 0$ est bloquée, la face $z = l$ étant soumise à un déplacement imposé (rotation d'ensemble d'angle $\theta$) autour de $e_z$.

Dans l'exemple traité, nous avons choisi une valeur de $T = 94,25$ N mm correspondant à une contrainte de cisaillement de peau (valeur maximale) théorique de $\tau(r = R) = 60$ MPa. Le maillage considéré est composé de 626 nœuds surfaciques et 656 nœuds volumiques correspondant à 208 éléments de frontière de type $Q8$ et 96 éléments de domaine de type $H20$.

Les résultats numériques obtenus en régime élastique sont concluants. En effet la “relation de comportement” macroscopique (entre angle de torsion appliqué et moment de torsion) est bien retrouvée numériquement. Le champ de contrainte calculé admet bien une composante principale de cisaillement radial.

L'évaluation de la composante circonférentielle à l'aide de la procédure décrite plus haut est également correcte, ce qui valide donc les opérateurs matriciels $Q$ et $G$ qui sont également utilisés pour le calcul élastoplastique.
<figure><img src="image_1.png" /><figcaption>FIGURE 1 – Maillage par éléments de frontière (les côtés des éléments sont en fait courbes puisque l’interpolation unidimensionnelle est quadratique). Les cellules volumiques occupent un cylindre correspondant aux deux couronnes externes d’éléments de frontière.</figcaption></figure>

En revanche, les résultats actuellement obtenus pour le régime élastoplastique ne sont pas satisfaisants : les courbes moment-angle de torsion numérique et analytique divergent dès l’apparition de plasticité. Cette première tentative de validation souffre manifestement (pour la phase élastoplastique) de défauts de jeunesse. En particulier, l’intégration numérique, très délicate pour certains termes, est probablement perfectible. Il n’est pas non plus exclu que des erreurs subsistent au niveau du calcul de la matrice Ξ ou de la mise en œuvre de l’algorithme de retour radial. Des résultats plus étoffés seront disponibles lors du colloque.

<table><thead><tr><th></th><th>Val. théorique</th><th>Val. num. max.</th><th>Val. num. min.</th><th>Érreur</th></tr></thead><tbody><tr><td>τ(r = R) en MPa</td><td>60,0</td><td>63,1</td><td>59,5</td><td>5,2 %</td></tr><tr><td>θ en 10<sup>-3</sup> rad</td><td>7,80</td><td>7,82</td><td>7,78</td><td>2,6 %</td></tr></tbody></table>

Références

[1] SIMO J.C., KENNEDY J.G., TAYLOR R.L., Complementary mixed finite element formulations for elastoplasticity, *Comput. Meth. in Appl. Mech. Engng.*, **74**, pp. 177-206, 1989.

[2] POLIZZOTTO C., PANZECA T., ZITO M., A consistent Boundary/Interior Element Method for evolutive elastic plastic structural analysis, In *Avances in BE Techniques*, eds. J.H. Kane, G. Maier, N. Tosaka, S.N. Atluri, Springer, 1993.

[3] BONNET M., Regularized direct and indirect symmetric variational BIE formulations for three-dimensional elasticity, *Eng. Anal. with Bound. Elem.*, **15**, pp. 93-102, 1995.
<figure><img src="image_2.png" /><figcaption>FIGURE 2 – Comparaison des relations entre couple et angle de torsion analytique et numérique</figcaption></figure>

[4] MAIER G., MICCOLI S., NOVATI G., PEREGO U., Symmetric Galerkin boundary element method in plasticity and gradient plasticity, *Comput. Mech.*, **17**, pp. 115-129, 1995.

[5] BONNET M., MAIER G., POLIZZOTTO C., Symmetric Galerkin boundary element method, *Appl. Mech. Review*, **51**, pp. 669-704, 1998.