# Méthodes des éléments de frontière appliquées à un problème élastoplastique

Bertrand Burgardt, Anh Le Van

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

## To cite this version:

Bertrand Burgardt, Anh Le Van. Méthodes des éléments de frontière appliquées à un problème élastoplastique. 14e Congrès Français de Mécanique, Aug 1999, Toulouse, France. hal-04670823

HAL Id: hal-04670823

https://minesparis-psl.hal.science/hal-04670823v1

Submitted on 13 Aug 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
# MÉTHODES DES ÉLÉMENTS DE FRONTIÈRE APPLIQUÉES À UN PROBLÈME ÉLASTOPLASTIQUE

B. Burgardt & A. Le van

Laboratoire Mécanique et Matériaux,
Ecole Centrale de Nantes, BP 92101,
44321 Nantes Cedex 3, France

## RESUME

Ce papier propose une formulation implicite par éléments de frontière d'un problème élastoplastique enптites perturbations, utilisant le concept d'opérateur tangent cohérent. Une nouvelle représentation intégrale du champ de déformation est incorporée dans un schéma par collocation. Parallèlement une méthode de Bubnov-Galerkin entièrement régularisée et aboutissant à un opérateur tangent symétrique, est présentée.

## ABSTRACT

This paper deals with an implicit boundary element method formulation for elastoplastic problems in small strains. A new strain boundary integral representation is incorporated in an collocation scheme. A fully regularized Galerkin formulation, which leads to a symmetric Consistent Tangent Operator, is also presented.

## MOTS-CLEFS

Elastoplasticité, Opérateur Tangent Cohérent, éléments de frontière.

## I. Introduction

La méthode des éléments de frontière représente une alternative intéressante aux éléments finis de domaine standard pour des problèmes élastiques, dans la mesure où seul le maillage de la frontière du milieu étudié est nécessaire. Pour certains problèmes spécifiques (milieux infinis, mécanique de la rupture), cette méthode permet d'éviter un maillage 3D souvent laborieux et d'aboutir à un système matriciel de taille modeste. Les formulations par équations intégrales ont été étendues à la résolution de problèmes pouvant présenter des non-linéarités matérielles (élastoplasticité, endommagement) et géométriques, mais leur place occupée dans la littérature demeure assez faible. Ceci s'explique principalement par la relative complexité des formulations intégrales du champ gradient des déplacements (ou des contraintes) et par l'introduction d'un maillage volumique. Il convient cependant de nuancer ce dernier point, en remarquant que la méthode est d'autant plus adaptée que la zone des non-linéarités est confinée comme c'est la cas pour des problèmes de rupture, d'endommagement ou de contact.

Ce papier traite d'un problème élastoplastique 3D avec écrouissage isotrope et cinématique linéaire sous l'hypothèse des petites perturbations. Actuellement la méthode de collocation est principalement utilisée pour résoudre les équations intégrales du problème. La présente étude propose une généralisation des processus de régularisation de la représentation intégrale du champ de déformation. Parallèlement on y présente une extension à l'élastoplasticité de la formulation variationnelle des équations intégrales menant à un opérateur tangent cohérent (CTO) symétrique.
## II. Modèle élastoplastique et méthodologie

Considérons un solide borné $\Omega$, de frontière Lipschitzienne $S$ et comportant une zone susceptible de plastifier $\Omega_p$ à l'intérieur de laquelle la déformation plastique $\varepsilon^p$ est supposée non nulle. On effectue une partition de $S$ en $S_u$ et $S_t$, supportant respectivement des conditions aux limites de type Dirichlet et Neumann.

Les champs de contrainte et de déformation sont notés $\boldsymbol{\sigma}$ et $\boldsymbol{\varepsilon}$ et nous appelons $\boldsymbol{\beta}$, la variable d'écrouissage cinématique et $p$ la déformation plastique cumulée. Un schéma d'intégration en temps de type implicite est choisi. Les forces de volume sont négligées. En suivant une démarche analogue à [Simo et al.89], les équations du problème élastoplastique peuvent alors être mises sous la forme ci-dessous :

♦ Equations de champ,

$$ \int_{\Omega} \nabla \delta \mathbf{u}(x) : \boldsymbol{\sigma}_{n+1}(x) \, dV(x) = \int_{S} \delta \mathbf{u}(x) \cdot \boldsymbol{\sigma}_{n+1}(x) \cdot \mathbf{n}(x) \, dS(x) \quad (1) $$

$$ \int_{\Omega} \delta \varepsilon(x) : \mathbf{C} : \varepsilon_{n+1}(x) \, dV(x) = \int_{\Omega} \delta \varepsilon(x) : \mathbf{C} : \nabla^s \mathbf{u}_{n+1}(x) \, dV(x) \quad (2) $$

$$ \int_{\Omega} \delta \varepsilon^p(x) : \boldsymbol{\sigma}_{n+1}(x) \, dV(x) = \int_{\Omega} \delta \varepsilon^p(x) : \mathbf{C} : [\varepsilon_{n+1}(x) - \varepsilon_{n+1}^p(x)] \, dV(x) \quad (3) $$

♦ Conditions aux limites,

$$ \int_{S_u} \delta \mathbf{t}(x) \cdot \mathbf{u}_{n+1}(x) \, dS(x) = \int_{S_u} \delta \mathbf{t}(x) \cdot \mathbf{u}_{n+1}^0(x) \, dS(x) \quad (4) $$

$$ \int_{S_u} \delta \mathbf{u}(x) \cdot \boldsymbol{\sigma}_{n+1}(x) \cdot \mathbf{n}(x) \, dS(x) = \int_{S_u} \delta \mathbf{u}(x) \cdot \mathbf{t}_{n+1}^0(x) \, dS(x) \quad (5) $$

♦ Lois d'évolution.

$$ \int_{\Omega_p} \delta \varepsilon(x) : \mathbf{C} : [ \varepsilon_{n+1}^p(x) - \varepsilon_n^p(x) ] \, dV(x) = \int_{\Omega_p} \delta \varepsilon(x) : \mathbf{C} : \lambda_{n+1}(x) \hat{\mathbf{n}}_{n+1}(x) \, dV(x) \quad (6) $$

$$ \lambda_{n+1}(x) = \frac{1}{2\theta\left(G + \frac{H+h}{3}\right)} \langle \phi(\mathbf{s}_{n+1}^E(x), \boldsymbol{\beta}_n(x), p_n(x)) \rangle \quad (7) $$

$\mathbf{s}_{n+1}^E$, $h$ et $H$ représentant respectivement le déviateur d'essai des contraintes et les modules d'écrouissage isotrope et cinématique, le critère de plasticité étant pris sous la forme :

$$ \phi_{n+1} = \| \mathbf{s}_{n+1} - \boldsymbol{\beta}_{n+1} \| - \sqrt{\frac{2}{3}(\sigma_0 + hp_{n+1})} \quad (8) $$

Alors que le traitement classique du caractère non-linéaire des équations est actuellement bien établi [Simo et al.89] et engendre l'algorithme dit de Retour Radial associé à la notion de CTO, l'introduction de ces concepts aux éléments de frontière, n'a été proposée que récemment dans [Bonnet et al.96]. L'écriture des équations sous forme faible permet de regrouper dans le même formalisme les formulations obtenues par une méthode de collocation et par une méthode de Bubnov-Galerkin.

L'idée maîtresse des éléments de frontière, consiste à résoudre de manière implicite les équations de domaine en ramenant le problème sur la frontière, les équations du problème se résumant alors aux conditions aux limites et portant sur les traces du champ de déplacement sur $S$. Introduisons l'ensemble $\mathcal{V}$ défini par :

$$ \mathcal{V} = \{ \mathbf{u} \in H^1(\mathcal{L}, \Omega) / \mathcal{L} \mathbf{u}_{n+1} = 0 \} \quad (9) $$
avec $\mathcal{L} . = \text{Div}[C : (\nabla^s . - \varepsilon_{n+1}^p)]$ où $\varepsilon_{n+1}^p$ est considéré comme un champ de déformation initiale. L'ensemble $\mathcal{V}$ est alors exactement connu et chaque élément est donné sous la forme d'une représentation intégrale.

$$ u_k^{(n+1)}(x) = \int_S \tilde{t}_i^{(n+1)}(y) U_i^k(y-x) \, dS(y) - \int_S \sum_{ij}^k (y-x)n_j(y) \, \tilde{u}_i^{(n+1)}(y) \, dS(y) + \int_{\Omega_P} \sum_{ij}^k (y-x) \varepsilon_{ij}^{p,(n+1)}(y) \, dV(y) \quad (10) $$

où $\tilde{u}$ et $\tilde{t}$ désignent respectivement la trace et le vecteur contrainte sur la frontière $S$. Les termes $U_i^k(r)$ représentent les composantes du tenseur de Kelvin et $\sum_{ij}^k(r)$, le tenseur des contraintes correspondant à $U_i^k(r)$, avec $U(r) = O(r^{-1})$ et $\Sigma(r) = O(r^{-2})$. Le champ gradient de $u^{(n+1)}$ peut s'exprimer sous forme régularisée, après dérivation sous le signe somme de (10).

$$ \begin{aligned} u_{k,l}^{(n+1)}(x) = & \int_S \sum_{ij,l}^k (y-x)n_j(y) \tilde{u}_i^{(n+1)}(y) \, dS(y) - \int_S \tilde{t}_i^{(n+1)}(y) U_{i,l}^k(y-x) \, dS(y) \\ & - \int_{\Omega_P} \sum_{ij,l}^k (y-x) \left[ \varepsilon_{ij}^{p,(n+1)}(y) - \varepsilon_{ij}^{p,(n+1)}(x) \right] \, dV(y) \\ & - \varepsilon_{ij}^{p,(n+1)}(x) \int_S \sum_{ij}^k (y-x)n_l(y) \, dS(y) \end{aligned} \quad (11) $$

L'algorithme élastoplastique par éléments de frontière s'apparente à une formulation à trois champs à savoir, densités surfaciques inconnues ($\tilde{u}, \tilde{t}$) notées $y$, déformation totale $\varepsilon$ et déformation plastique $\varepsilon^p$. Introduisons une discrétisation de type éléments finis :

$$ \varepsilon_{n+1}(x) = N(x) \cdot E_{n+1} \qquad \varepsilon_{n+1}^p(x) = N(x) \cdot E_{n+1}^p \qquad y_{n+1}(x) = N(x) \cdot Y_{n+1} \quad (12) $$

L'écriture des conditions aux limites, utilisant la représentation intégrale (10), permet d'obtenir le système classique d'équations intégrales, dans lequel apparaît un terme de volume faisant intervenir les déformations plastiques, sous la forme d'une matrice d'influence $Q$ :

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

L'évaluation des déformations dans le domaine $\Omega_p$ à partir de la représentation intégrale (11), permet d'obtenir une relation directe entre le champ de déformation totale et le champ de déformation plastique.

$$ M \cdot E_{n+1} = G_{n+1} + R \cdot Y_{n+1} - H \cdot E_{n+1}^p \quad (14) $$

En éliminant les densités surfaciques inconnues $Y_{n+1}$ au moyen de (13), on aboutit finalement au système d'équations non-linéaires à résoudre, à savoir :

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

A ce stade de l'étude, nous disposons d'une équation linéaire reliant la déformation plastique $E_{n+1}^p$ à la déformation totale $E_{n+1}$. Il reste à discrétiser les lois d'évolution (6) pour obtenir l'équation non-linéaire caractérisant le comportement du matériau.

Il nous faut à présent expliciter les matrices intervenant dans les systèmes ci-dessus. Leurs expressions dependent de la nature des champs variés choisis lors de l'écriture sous forme faible des équations et font l'objet des prochaines sections.
### III. Application de la méthode de collocation

La méthode de collocation consiste à prendre des mesures de Dirac en un nombre suffisant de points, comme champ varié dans l'écriture des équations sous forme faible. Généralement pour des raisons de facilité de mise en oeuvre numérique, il est préférable de choisir les noeuds du maillage. L'écriture des conditions aux limites en déplacement aux noeuds du maillage de la frontière S donne lieu à la DBIE (Displacement Boundary Integral Equation) et est donnée sous forme régularisée ou non, par exemple dans [Bonnet94] et [Sladek et al.98].

La discrétisation de (11) permet d'exprimer le tenseur de déformation totale au point de collocation considéré. Certains problèmes surviennent néanmoins, lorsque le point de collocation est situé sur la frontière S puisque $\Sigma(r) = O(r^{-2})$. Certains articles récents donnent une valeur limite de (11) pour une point régulier de la surface, moyennant un processus de régularisation lourd à mettre en oeuvre numériquement ([Bonnet94],[Sladek et al.98],[Huber et al.96]). La démarche que nous proposons allège la formulation et permet de choisir un point de collocation sur une arête ou sur un coin de la frontière. Introduisons alors $\tilde{h}$ et $\tilde{\sigma}$, représentant respectivement le gradient du déplacement et la contrainte, reconstitués sur S.

$$ 
\begin{aligned}
\tilde{h}_{ij}(y) &= D_j \tilde{u}_i(y) + \frac{1}{G} \left[ \delta_{ik} - \frac{n_i(y)n_k(y)}{2(1-\nu)} \right] \left[ \tilde{t}_k(y) + C_{klab} \varepsilon_{ab}^p(y) n_l(y) \right] n_j(y) \\
&\quad - \frac{\nu}{1-\nu} D_k \tilde{u}_k(y) n_i(y) n_j(y) - n_k(y) D_i \tilde{u}_k(y) n_j(y) \\
\tilde{\sigma}_{ij}(y) &= C_{ijab} \left[ \tilde{h}_{ab}(y) - \varepsilon_{ab}^p(y) \right]
\end{aligned}
$$

La représentation intégrale (11) peut alors être mise sous la forme suivante pour tout point intérieur au domaine :

$$ 
\begin{aligned}
\forall x \in \Omega, u_{k,l}(x) = & \int_S U_i^k(x,y) D_{bl} \tilde{\sigma}_{ab}(y) \, dS(y) - \int_S T_i^k(x,y) \tilde{h}_{il}(y) \, dS(y) \\
& + \int_{\Omega} \Sigma_{ij}^k(x,y) \varepsilon_{ij,l}^p(y) \, dV(y)
\end{aligned}
\quad (16) $$

Il est à noter, que les noyaux des opérateurs intégraux, intervenant dans (16), sont identiques à ceux, apparaissant dans (10). La trace de (16) sur S est alors déterminée par la relation intégrale ci-dessous, moyennant la condition de régularité sur $\tilde{h}$ au point de collocation $\tilde{h} \in C^{0,\alpha}(S)$.

$$ 
\begin{aligned}
\forall z \in S, u_{k,l}(z) = & \int_S U_i^k(z,y) D_{bl} \tilde{\sigma}_{ab}(y) \, dS(y) - \int_S T_i^k(z,y) [\tilde{h}_{il}(y) - \tilde{h}_{il}(z)] \, dS(y) \\
& + \tilde{h}_{kl}(z) + \int_{\Omega} \Sigma_{ij}^k(z,y) \varepsilon_{ij,l}^p(y) \, dV(y)
\end{aligned}
\quad (17) $$

Les lois d'évolution sont écrites en chaque noeud du maillage de $\Omega_P$ pour se ramener à une équation non-linéaire en $E_{n+1}$ (pour de plus amples précisions, voir [Bonnet et al.96]).

### IV. Application de la méthode de Bubnov-Galerkin

Les équations du problème sont prises sous forme faible en choisissant pour chaque inconnue du système, le même espace pour le champ réel et le champs varié. La méthode proposée ci-dessous est principalement inspirée de [Polizzotto et al.94] et étend la formulation de [Bonnet95] à des problèmes élastoplastiques. En suivant une démarche analogue à [Bonnet95], on montre que
le système d'équations relatif aux conditions aux limites du problème, peut se mettre sous la forme suivante :

$$ 
\begin{aligned}
a_{uu}(\delta\tilde{u}, \tilde{u}) + a_{ut}(\delta\tilde{u}, t) &= \ell_u(\delta\tilde{u}) + a_{u\epsilon}(\delta\tilde{u}, \epsilon^p) \\
a_{tu}(\delta\tilde{t}, \tilde{u}) + a_{tt}(\delta\tilde{t}, t) &= \ell_t(\delta\tilde{t}) + a_{t\epsilon}(\delta\tilde{t}, \epsilon^p)
\end{aligned} 
\qquad (18) $$

où les formes bilinéaires introduites vérifient les propriétés de symétrie,

$$ 
\left\{ 
\begin{aligned}
a_{uu}(\delta\tilde{u}, \tilde{u}) &= a_{uu}(\tilde{u}, \delta\tilde{u}) \\
a_{tt}(\delta\tilde{t}, \tilde{t}) &= a_{tt}(\tilde{t}, \delta\tilde{t}) \\
a_{ut}(\tilde{u}, \tilde{t}) &= a_{tu}(\tilde{t}, \tilde{u})
\end{aligned}
\right. 
$$

entrainant une matrice K symétrique dans (13). La représentation intégrale sous forme faible du champ de déformation totale et des lois d'évolution donnent :

$$ 
\begin{aligned}
\int_{\Omega_p} \delta\epsilon^p : C : \epsilon_{n+1} \, dV &= \ell_{\epsilon}(\delta\epsilon^p) + a_{\epsilon u}(\delta\epsilon^p, \mathbf{u}_{n+1}) + a_{\epsilon t}(\delta\epsilon^p, \mathbf{t}_{n+1}) - a_{\epsilon \epsilon}(\delta\epsilon^p, \epsilon_{n+1}^p) \\
\int_{\Omega_p} \delta\epsilon : C : (\epsilon_{n+1}^p - \epsilon_n^p) \, dV &= \int_{\Omega_p} \delta\epsilon : C : \lambda_{n+1} \hat{\mathbf{n}}_{n+\theta} \, dV
\end{aligned} 
$$

avec les propriétés suivantes :

$$ 
\left\{ 
\begin{aligned}
a_{u\epsilon}(\tilde{u}, \epsilon) &= a_{\epsilon u}(\epsilon, \tilde{u}) \\
a_{t\epsilon}(\tilde{t}, \epsilon) &= a_{\epsilon t}(\epsilon, \tilde{t}) \\
a_{\epsilon\epsilon}(\delta\epsilon^p, \epsilon^p) &= a_{\epsilon\epsilon}(\epsilon^p, \delta\epsilon^p)
\end{aligned}
\right. 
$$

On montre que ces équations résultent de l'écriture de la stationnarité d'une fonctionnelle. Celle-ci proposée dans [Polizzotto et al.94] pour un matériau élastoplastique parfait a été étendue à un matériau écrouissable, dans la mesure où il reste standard généralisé, en suivant [Simo et al.89]. Il en résulte certaines propriétés liées à l'existence, l'unicité et la stabilité du schéma numérique proposé [Maier et al.95].

La forme discrétisée de ces équations mène alors au système non-linéaire suivant :

$$ M.E_{n+1} = G_{n+1} - \Xi.E_{n+1}^p (E_{n+1}) \qquad (19) $$

$$ M.E_{n+1}^p = M.E_n^p + \Psi(E_{n+1}) \qquad (20) $$

où Ξ est une matrice symétrique. En éliminant le vecteur $E_{n+1}^p$, on obtient un système non-linéaire à résoudre, prenant la forme ci-dessous :

$$ R(E_{n+1}) = E_{n+1} + M^{-1} \Xi M^{-1} (\Psi(E_{n+1}) + M.E_n^p) - M^{-1} G_{n+1} = 0 \qquad (21) $$

On montre que l'opérateur cohérent tangent associé au schéma de Newton-Raphson du système est symétrique. Une forme réduite $\Lambda$, est donnée par l'expression suivante :

$$ \Lambda = I + M^{-1} \Xi M^{-1} . 2G \int_{\Omega_p} N^t(x) \frac{\partial(\lambda_{n+1} \hat{n}_{n+1})}{\partial \epsilon_{n+1}} N(x) \, dV \qquad (22) $$

$$ \frac{\partial}{\partial \epsilon_{n+1}} (\lambda_{n+1} \hat{n}_{n+1}) = 2G \frac{\lambda_{n+1}}{\| \xi_{n+1}^E \|} (J - \hat{n}_{n+1} \otimes \hat{n}_{n+1}) + \frac{\hat{n}_{n+\theta} \otimes \hat{n}_{n+\theta}}{1 + \frac{H+h}{3G}} \qquad (23) $$

où J est l'opérateur déviatorique.
V. Conclusion

L'utilisation d'algorithmes élastoplastiques implicites dans les formulations par éléments finis de frontière est très récente. Une approche en collocation permettant la prise en compte de points irréguliers de la géométrie pour l'évolution élastoplastique du système a été proposée. Parallèlement une méthode de Bubnov-Galerkin, aboutissant à un opérateur tangent symétrique a été obtenue. Des exemples numériques permettant de valider les formulations sont en cours.

REFERENCES

Bonnet M. et Mukherjee S. ; 1996, Implicit bem formulations for usual and sensitivity problems
in elastoplasticity using the consistent tangent operator concept, *Int. J. Solids and Structures*
vol. 33 pp. 4461-4480.

Bonnet M.;, *Équations intégrales et éléments de frontière*, CNRS Éditions - Eyrolles 1994.

Bonnet M.; 1995, Regularized direct and indirect symmetric variational bie formulations for three-dimensional elasticity, *Eng. Anal. with Bound. Elem.* vol. 15 pp. 93–102.

Huber O., Dallner R., Partheymueller P. et Kuhn G. ; 1996, Evaluation of the stress tensor in 3-d elastoplasticity by direct solving of hypersingular integrals, *Int. J. Num. Meth. Engng*. vol. 39 pp. 2555–2573.

Maier G., Miccoli S., Novati G. et Perego U. ; 1995, Symmetric galerkin boundary element method in plasticity and gradient plasticity, *Comput. Mech.* vol. 17 pp. 115–129.

Polizzotto C. et Zito M. ; 1994, A step-wise variational approach to elastic-plastic analysis by boundary-interior elements, *Int. J. of Plasticity* vol. 10 pp. 81–101.

Simo J. C., Kennedy J. G. et Taylor R. L. ; 1989, Complementary mixed finite element formula-
tions for elastoplasticity, *Comput. Meth. in Appl. Mech. Engng*. vol. 74 pp. 177–206.

Sladek V., Sladek J. et Levan A. ; 1998, Completely regularized integral representations and integral equations for anisotropic bodies with initial strains, Z.A.M.M. vol. 78 pp. 771–780.