r7.01.33 Couplage élasto-plasticité-endommagement

Résumé:

Ce document présente une méthode de couplage d’un modèle élasto-plastique avec un modèle d’endommagement, laquelle est implantée dans le cadre du KIT_DDI. La méthode utilise une approche de partitionnement du modèle en trois modules: élasto-plasticité, endommagement et coupleur. C’est particulièrement cette modularité qui représente l’intérêt principal de la méthode. Elle permet d’utiliser le même module «coupleur» pour des couplages entre différents modèles d’élasto-plasticité et d’endommagement sans intervenir dans l’implantation de ces derniers. Exactement la même idée est utilisée pour les couplages mécanique/fluage (voir [R7.01.19]).

Pour l’instant, l’implantation permet de coupler que le modèle GLRC_DM pour l’endommagement avec les modèles VMIS_(ISOT/CINE)_LINE pour l’élasto-plasticité et ceci seulement pour des modélisations DKTG.

Modèles d’élasto-plasticité-endommagement en générale

Enjeux et difficultés

Dans la modélisation du comportement des solides, on constate que les phénomènes non-linéaires les plus importants à prendre en compte, lors d’une analyse de structure sous chargement extrême, concernent la présence et l’évolution de la déformation anélastique (plasticité) et l’affaiblissement de la raideur (endommagement). L’importance relative de ces deux phénomènes peut varier drastiquement d’un matériau à l’autre. Ainsi, le plus souvent on ne considère que la plasticité lors de la modélisation des métaux tout en négligeant l’endommagement, tandis que pour un matériau de type béton, l’endommagement est jugé plus important et la plasticité n’est pas prise en compte. Néanmoins, souvent cette simplification n’est pas justifiable et on est amené à modéliser les deux phénomènes, la plasticité et l’endommagement, à la fois.

Depuis déjà un certain temps le cadre thermodynamique pour la modélisation des deux phénomènes est bien établi (e.g. voir [bib1], [bib2]) et il n’existe aucun obstacle pour l’étendre aux couplages élasto-plastique-endommageables. Quant à la résolution numérique de ces modèles dans le cadre de la méthode des éléments finis (MEF), il existe des algorithmes très robustes notamment pour la plasticité (e.g. voir [bib13]) et aussi pour l’endommagement (e.g. voir [bib14]) à condition que le phénomène de localisation ne devienne pas prépondérant. Par contre, la résolution numérique d’un modèle couplé est beaucoup moins évidente.

Les difficultés numériques concernant l’algorithme de résolution d’un modèle couplé interviennent sur deux aspects essentiellement :

• satisfactions des conditions (de Kuhn-Tucker) des seuils de plasticité et d’endommagement,

• construction de la matrice tangente.

La première difficulté est due au fait que les variables internes plastiques et endommageables agissent très différemment sur le niveau des contraintes, ce qui rend le système des équations difficile à résoudre. Le deuxième point est une simple conséquence de la complexité augmentée du modèle couplé par rapport aux modèles à un seul phénomène.

Dernièrement, on a proposé dans [bib3] et [bib4] un algorithme de résolution d’un modèle élasto-plastique-endommageable basé sur le partitionnement strict en une partie plastique et à une partie endommageable. Il s’agit d’introduire un processus itératif supplémentaire afin de satisfaire l’admissibilité plastique et l’admissibilité endommageable simultanément. Ce processus itératif est piloté par ce qu’on appelle la déformation d’endommagement, une variable qui, contrairement à la déformation plastique, n’est pas une variable d’état, mais représente plutôt un moyen de couplage entre différents modèles. Ainsi, dans le modèle introduit dans [bib3] le choix du modèle élasto-plastique étant libre, on devait néanmoins modifier la formulation d’endommagement, qui ne prévoyait en l’occurrence que la modélisation de comportements non adoucissants. Dans [bib4] on a démontré l’intérêt d’intégrer ce type de lois dans le cadre d’éléments finis mixtes.

Dans ce document on généralise la démarche présentée dans [bib3] pour pouvoir coupler pratiquement n’importe quel modèle élasto-plastique avec n’importe quel modèle d’endommagement. L’intérêt d’une telle plate-forme va en réalité au-delà du couplage élasto-plastique-endommageable.

Cadre thermodynamique

Dans cette partie on rappelle la démarche habituelle pour construire des modèles d’endommagement et d’élasto-plasticité, basés sur la thermodynamique. La combinaison des deux types de modèle peut se faire d’une manière assez naturelle, pourvu que le couplage soit fait au niveau du calcul des contraintes. Toute la démarche présentée ici s’effectue sous l’hypothèse que les variables internes d’endommagement et de plasticité ne sont pas couplées directement entre elles. De plus, on se met dans le cadre des petites déformations et sans phénomène thermique.

Plasticité

L’énergie libre \(\psi\) peut s’écrire comme la somme de l’énergie élastique et de l’énergie dite stockée, due à la plasticité:

(3784)\[\psi (\varepsilon ,{\varepsilon}^{p},\alpha )={\psi}^{e}(\varepsilon -{\varepsilon}^{p})+{\psi}^{p}(\alpha )\]

où \(\varepsilon\) est la déformation totale, \({\varepsilon}^{p}\) la déformation plastique et \(\alpha\) un vecteur de variables internes. L’hypothèse essentielle qu’on fait dans () est que l’énergie élastique ne dépend que de la déformation élastique, \({\psi}^{e}={\psi}^{e}({\varepsilon}^{e})\) , où \({\varepsilon}^{e}=\varepsilon -{\varepsilon}^{p}\) .

Dans la deuxième étape on définit le taux de dissipation plastique et on s’assure qu’elle reste toujours positive,

(3784)\[{D}^{p}=\sigma :\dot{\varepsilon}-\rho \dot{\psi}=(\sigma -\rho \frac{\partial {\psi}^{e}}{\partial {\varepsilon}^{e}}):{\dot{\varepsilon}}^{e}+\sigma :{\dot{\varepsilon}}^{p}-A\mathrm{\ast }\dot{\alpha}\ge 0\]

Où

(3784)\[A=\rho \frac{\partial {\psi}^{p}}{\partial \alpha }\]

sont les forces thermodynamiques associées à des variables internes et le symbole \(\ast\) fait référence à l’opérateur de produit approprié entre \(A\) et \(\dot{\alpha}\) .

On suppose que la dissipation n’évolue pas, i.e. \(\dot{D}=0\) , lorsque les variables plastiques n’évoluent pas, i.e. \({\dot{\epsilon}}^{p}=\dot{\alpha}=0\) , ce qui nous permet de définir la contrainte comme:

(3784)\[\sigma =\rho \frac{\partial {\psi}^{e}}{\partial {\epsilon}^{e}}\]

Dans ce document, on se contente des matériaux qui ont un comportement linéaire et isotrope. Dans ce cas, la contribution élastique de l’énergie libre est donnée par:

(3784)\[\rho \psi ({\epsilon}^{e})=\frac{1}{2}{\epsilon}^{e}:{D}^{e}:{\epsilon}^{e}\]

Où \({D}^{e}\) est le tenseur d’élasticité isotrope standard. Des deux équations précédentes, on déduit l’expression du tenseur des contraintes:

(3784)\[\sigma ={D}^{e}:{\epsilon}^{e}\]

Un autre ingrédient essentiel d’un modèle élasto-plastique est la définition de la fonction du seuil d’écoulement qui dépend du tenseur des contraintes et des forces thermodynamiques associées aux variables internes. Cette fonction est négative quand les déformations sont purement élastiques et atteint zéro quand l’écoulement plastique est imminent:

(3784)\[{\Phi}^{p}={\Phi}^{p}\left(\sigma ,A\right)\le 0\]

La caractérisation complète du modèle général de plasticité nécessite la définition des lois d’évolution des variables internes, i.e. les varaibles associées aux phénomènes de dissipation. Dans notre cas, les variables internes sont le tenseur des déformations plastiques et l’ensemble de variables internes \(\alpha\) . Pour établir les lois d’évolution de ces variables internes ont utilises souvent le principe de la dissipation plastique maximale. Cependant ceci n’est pas essentiel dans notre démarche. On se servira d’une écriture plus générale:

(3784)\[\begin{split}\begin{array}{c}{\dot{\epsilon}}^{p}={\dot{\gamma}}^{p}N={\dot{\gamma}}^{p}\frac{\partial \psi }{\partial \sigma }\\ \dot{\alpha}={\dot{\gamma}}^{p}H=-{\dot{\gamma}}^{p}\frac{\partial \psi }{\partial A}\end{array}\end{split}\]

\(\dot{\gamma}\) est le multiplicateur plastique, \(N\) est appelé le vecteur d’écoulement et \(H\) est le module d’écrouissage. De plus toute solution doit satisfaire les conditions d” admissibilité plastiques , qui sont:

(3784)\[{\Phi}^{p}\le 0,\dot{\gamma}\ge 0,{\Phi}^{p}\dot{\gamma}=0\]

Endommagement

Pour les modèles d’endommagement il est plus difficile d’avoir une construction aussi générale que pour celle de la plasticité. Par conséquent on traitera un cas simple mais assez représentatif des caractéristiques souhaitées du modèle. Ainsi, on définit l’énergie libre comme:

(3784)\[\psi \left(\epsilon ,d\right)={\psi}^{e}\left(\epsilon \right)\zeta \left(d\right)\]

où \(\zeta\) est la fonction d’endommagement et \(d\) la variable d’endommagement, qui est un scalaire mais qui peut aussi bien être un tenseur. De même, la fonction \(\zeta\) peut en général dépendre de \(\epsilon\) , notamment du signe de sa trace pour distinguer des comportements typiquement très différents en traction et en compression.

En utilisant (), on définit la dissipation d’endommagement comme,

(3784)\[{D}^{d}=\sigma :\dot{\epsilon}-\rho \dot{\psi}=\left(\sigma -\rho \frac{\partial \psi }{\partial \epsilon }\right):\dot{\epsilon}-\rho \frac{\partial \psi }{\partial d}\dot{d}\ge 0\]

La dissipation n’évolue pas, i.e. \({\dot{D}}^{d}=0\) , lorsque la variable d’endommagement n’évolue pas, i.e. \(\dot{d}=0\) , ce qui nous permet de définir la contrainte comme:

(3784)\[\sigma =\rho \frac{\partial \psi }{\partial \epsilon }=\rho \frac{\partial {\psi}^{e}}{\partial \epsilon }\zeta \left(d\right)={D}^{e}\left(d\right):{\epsilon}^{e}\]

Où \({D}^{e}(d)\) représente le tenseur d’élasticité isotrope standard évoluant en fonction de l’endommagement. En faisant l’hypothèse d’un endommagement isotrope, en prenant \(\zeta =1-d\) la contrainte devient:

(3784)\[\sigma =(1-d){D}^{e}:{\varepsilon}^{e}\]

Comme pour la plasticité on introduit une fonction seuil d’endommagement:

(3784)\[{\Phi}^{d}={\Phi}^{d}(\varepsilon ,Y)\le 0\]

Où \(Y\) est la variable duale à \(d\) . Cette fois, le seuil est une fonction directe de la déformation totale, \(\varepsilon\) , et on cherche à satisfaire les conditions d’admissibilité d’endommagement:

(3784)\[{\Phi}^{d}\le 0,\dot{d}\ge 0,{\Phi}^{d}\dot{d}=0\]

On doit aussi formuler la loi d’évolution de \(d\) :

(3784)\[\dot{d}=\dot{d}\left(Y\right)\]

Couplage

Enfin, on combine les formulations d’endommagement et de plasticité pour arriver à un modèle couplé. L’énergie libre du modèle couplé s’écrit donc comme:

(3784)\[\psi (\varepsilon ,{\varepsilon}^{p},\alpha ,d)={\psi}^{e}(\varepsilon -{\varepsilon}^{p})\zeta (d)+{\psi}^{p}(\alpha )\]

La dissipation s’écrit:

(3784)\[D=\sigma :\dot{\varepsilon}-\rho \dot{\psi}=(\sigma -\rho \frac{\partial {\psi}^{e}}{\partial {\varepsilon}^{e}}):{\dot{\varepsilon}}^{e}+\sigma :{\dot{\varepsilon}}^{p}-A\mathrm{\ast }\dot{\alpha}-\frac{\partial {\psi}^{d}}{\partial d}\dot{d}\ge 0\]

Le tenseur des contraintes vaut:

(3784)\[\sigma =\rho \frac{\partial {\psi}^{d}}{\partial {\varepsilon}^{e}}=\rho \frac{\partial {\psi}^{e}}{\partial {\varepsilon}^{e}}\zeta (d)={D}^{e}(d):(\varepsilon -{\varepsilon}^{p})\]

Si l’élasticité est linéaire isotrope et l’endommagement est isotrope:

(3784)\[\sigma =(1-d){D}^{e}:(\varepsilon -{\varepsilon}^{p})\]

On a alors deux fonctions seuil: une pour l’écoulement plastique et l’autre pour l’endommagement:

(3784)\[\begin{split}\begin{array}{c}{\Phi}^{p}={\Phi}^{p}(\sigma ,A)\le 0\\ {\Phi}^{d}={\Phi}^{d}({\varepsilon}^{e},Y)\le 0\end{array}\end{split}\]

Les lois d’évolution des variables d’endommagement sont données par:

(3784)\[\begin{split}\begin{array}{c}{\dot{\varepsilon}}^{p}={\dot{\gamma}}^{p}N={\dot{\gamma}}^{p}\frac{\partial \psi }{\partial \sigma }\\ \dot{\alpha}={\dot{\gamma}}^{p}H=-{\dot{\gamma}}^{p}\frac{\partial \psi }{\partial A}\\ \dot{d}=\dot{d}(Y)\end{array}\end{split}\]

Les conditions d’admissibilité sont:

(3784)\[\begin{split}\begin{array}{c}{\Phi}^{p}\le 0,\dot{\gamma}\ge 0,{\Phi}^{p}\dot{\gamma}=0\\ {\Phi}^{d}\le 0,\dot{d}\ge 0,{\Phi}^{d}\dot{d}=0\end{array}\end{split}\]

Intégration numérique

Rappel de la résolution classique des problèmes non-linéaire par éléments finis

Dans la méthode proposée on rajoute un processus itératif supplémentaire par rapport à un calcul classique par la MEF, qu’on rappelle dans ce chapitre. Dans une résolution globale d’un problème de structure on cherche l’état d’équilibre, défini par:

(3784)\[{f}^{int}(u(t))={f}^{\mathit{ext}}(t)\]

Où \({f}^{int}\) est le vecteur des forces internes et \({f}^{\mathit{ext}}\) le vecteur des forces externes (le chargement). Le premier est une fonction du vecteur des déplacements, la variable principale du système, et le deuxième une fonction du paramètre t (pseudo-temps), définissant la notion de trajectoire aussi bien pour le chargement que pour l’évolution du champ des déplacements. Dans le cadre de la MEF la force interne est déterminée par:

(3784)\[{f}^{int}(u)={\int}_{\Omega}{B}^{T}\sigma (\varepsilon )d\Omega\]

où \(\Omega\) définit le domaine considéré, \(\varepsilon\) le tenseur de la déformation, \(\sigma\) le tenseur de la contrainte et \(B\) est la matrice classique correspondant à l’opérateur discrétisé du gradient symétrique qui permet de passer des déformations aux déplacements.

Comme nous avons considéré l’hypothèse des petites déformations, le tenseur \(\varepsilon\) est défini comme:

(3784)\[{\epsilon}_{ij}=\frac{1}{2}(\frac{\partial {u}_{i}}{\partial {x}_{j}}+\frac{\partial {u}_{j}}{\partial {x}_{i}})\]

En utilisant la définition de \(B\) , on a:

(3784)\[{\varepsilon}_{ij}={B}_{\mathit{ijk}}{u}_{k}\]

Après avoir discrétisé le champ des déplacements par la MEF. Lorsque le comportement du matériau est non-linéaire, on a le plus souvent recours à la méthode de Newton pour résoudre (). Pour ce faire, on est amené à linéariser (), ce qui nous donne:

(3784)\[{f}^{int,(k+1)}={f}^{int,(k)}+\frac{\partial {f}^{int,(k)}}{\partial {u}^{(k)}}\Delta {u}^{(k)}\]

Où l’indice \(k\) représente le numéro de l’itération en cours. À la convergence on doit avoir:

(3784)\[{f}^{int,(k)}\stackrel{k\to \infty }{\to }{f}^{\mathit{ext}}(t)\]

Afin de pouvoir calculer la correction de déplacements d’une itération à l’autre, \(\Delta {u}^{(k)}\) , on doit déterminer:

(3784)\[K=\frac{\partial {f}^{int,(k)}}{\partial {u}^{(k)}}=\underset{\Omega}{\int}{B}^{T}{D}^{\mathit{epd}}Bd\Omega\]

Où \({D}^{\mathit{epd}}\) est le module tangent cohérent pour la loi de comportement utilisée:

(3784)\[{D}^{\mathit{epd}}=\frac{\partial \sigma (\varepsilon )}{\partial \varepsilon }\]

La méthode détaillée ci-dessous se positionne au niveau d’un point matériel (voire au point de Gauss) et intervient au niveau du calcul de \(\sigma\) et \({D}^{\mathit{epd}}\) à partir de la déformation \(\epsilon\) , donnée pour chaque appel du processus global des itérations de Newton.

Approche au problème couplé

Dans ce document on ne détaille pas l’intégration numérique dans le temps pour chacun des modèles. On suppose tout simplement qu’on sait le faire et ceci avec un schéma implicite d’Euler (arrière). Disons, que dans une résolution incrémentale par la méthode des éléments finis, pour chaque itération de Newton du processus itératif global, on veut calculer, pour chaque point de Gauss, les contraintes \({\sigma}_{n+1}\) et les variables internes \({\epsilon}_{n+1}^{p}\) , \({\alpha}_{n+1}\) et \({d}_{n+1}\) à partir des valeurs de la déformation totale, \({\epsilon}_{n+1}\) , tout en supposant qu’on a convergé à l’instant précédent et qu’on connaît \({\sigma}_{n}\) , \({\epsilon}_{n}^{p}\) , \({\alpha}_{n}\) et \({d}_{n}\) . Les indices \(n\) ou \(n+1\) indiquent que la valeur de la variable correspond à l’instant \({t}_{n}\) ou \({t}_{n+1}\) , respectivement.

Ainsi on suppose que les deux modèles peuvent être intégrés en mode découplé et qu’on est donc capable d’écrire:

(3784)\[{\sigma}_{n+1}^{p}={\sigma}^{p}\left({\stackrel{̃}{\epsilon}}_{n+1}\right)\]

et

(3784)\[{\sigma}_{n+1}^{d}={\sigma}^{d}\left({\stackrel{̃}{\epsilon}}_{n+1}\right)\]

Les exposants \(d\) ou \(p\) dans \({\sigma}^{d}\) et \({\sigma}^{p}\) signifient que la contrainte est calculée avec le modèle d’endommagement et avec le modèle de plasticité, respectivement. \(\stackrel{̃}{\epsilon}\) est une variable du type déformation.

Dans l’algorithme présenté dans la suite on introduit un processus itératif supplémentaire pour que les contraintes \({\sigma}^{d}\) et \({\sigma}^{p}\) convergent vers une seule valeur. Il est semblable à celui proposé dans [bib3] sauf qu’ici la construction du modèle d’endommagement est plus générale, ce qui exige de considérer une autre variable que la déformation d’endommagement (utilisée dans [bib3]) sur laquelle on itère.

Construction du système non-linéaire

Par rapport à une résolution classique de la loi de comportement, on doit faire converger ici deux valeurs de contraintes en introduisant une autre inconnue du type déformation. Ce système partitionné est construit plus facilement en introduisant la décomposition suivante de la déformation totale,

(3784)\[\epsilon ={\epsilon}^{e}+{\epsilon}^{p}+{\epsilon}^{d}\]

Où on a introduit la déformation d’endommagement, ce qui nous permet d’écrire la relation d’élasticité comme:

(3784)\[\sigma ={D}^{e}:\left(\epsilon -{\epsilon}^{p}-{\epsilon}^{d}\right)\]

ou encore:

(3784)\[{\epsilon}^{d}\left(\sigma ,\epsilon ,{\epsilon}^{p}\right)=\epsilon -{\epsilon}^{p}-{({D}^{e})}^{-1}:\sigma\]

On peut résumer les équations à résoudre en trois groupes :

Module plasticité

(3784)\[{\sigma}_{n+1}^{p}={\sigma}^{p}{\left({\epsilon}_{n+1}-{\epsilon}_{n+1}^{d}\right)}_{e} {\epsilon}_{n+1}^{p}={\epsilon}^{p}\left({\epsilon}_{n+1}-{\epsilon}_{n+1}^{d}\right)\]

L’indice \(e\) signifie que la loi élastoplastique est calculée avec l’opérateur d’élasticité \({D}^{e}\)

Module endommagement

(3784)\[{\sigma}_{n+1}^{d}={\sigma}^{d}{\left({\epsilon}_{n+1}-{\epsilon}_{n+1}^{p}\right)}_{e}\]

Module couplage

(3784)\[{\epsilon}_{n+1}^{d}={\epsilon}_{n+1}-{\epsilon}_{n+1}^{p}-{({D}^{e})}^{-1}:{\sigma}_{n+1}^{d} {\sigma}_{n+1}^{p}={\sigma}_{n+1}^{d}\]

Dans [bib3] on a choisi \({\epsilon}_{n+1}^{d}\) comme l’inconnue additionnelle, ce qui n’est pas possible dans notre cas, puisque \({\sigma}^{d}\) n’est pas une fonction de \({\epsilon}^{d}\) mais de \({\stackrel{̄}{\epsilon}}^{\mathit{ed}}=\epsilon -{\epsilon}^{p}\) , la variable qui est finalement utilisée pour piloter les processus itératifs décrits ci-dessous. Par ailleurs elle vaut aussi \({\stackrel{̄}{\epsilon}}^{\mathit{ed}}={\epsilon}^{e}+{\epsilon}^{d}=\epsilon -{\epsilon}^{p}\) .

Algorithme de résolution

Pour résoudre le système non-linéaire, défini dans les équations ( à ) en les attribuant aux différents modules (plasticité, endommagement et couplage), on applique la méthode classique de Newton. Afin de linéariser le système, on doit définir les dérivées suivantes:

(3784)\[{D}^{\mathit{ep}}\left({\stackrel{̄}{\epsilon}}^{\mathit{ed}}\right)=\frac{\partial {\sigma}_{n+1}^{p}}{\partial {\stackrel{̄}{\epsilon}}^{\mathit{ed}}} {D}^{d}\left({\stackrel{̄}{\epsilon}}^{\mathit{ed}}\right)=\frac{\partial {\sigma}_{n+1}^{d}}{\partial {\stackrel{̄}{\epsilon}}^{\mathit{ed}}} \frac{\partial {\epsilon}_{n+1}^{d}}{\partial {\stackrel{̄}{\epsilon}}^{\mathit{ed}}}=1-{\left({D}^{e}\right)}^{-1}{D}^{d}\left({\stackrel{̄}{\epsilon}}_{n+1}^{\mathit{ed}}\right)\]

Le système non-linéaire (à ) est tout d’abord ramené à la dernière des équations, c’est-à-dire à \({\sigma}_{n+1}^{p}={\sigma}_{n+1}^{d}\) , puisqu’on peut exprimer toutes les inconnues par rapport à la valeur de \({\stackrel{̄}{\epsilon}}_{n+1}^{\mathit{ed}}\) . On peut donc définir le résidu comme:

(3784)\[{R}_{\sigma}^{(k)}={\sigma}^{p}\left({\epsilon}_{n+1}-{\epsilon}_{n+1}^{d,(k)}\right)-{\sigma}^{d}\left({\stackrel{̄}{\epsilon}}_{n+1}^{\mathit{ed},(k)}\right)\]

où l’exposant \(k\) indique l’itération courante et où:

(3784)\[{\epsilon}_{n+1}^{d,(k)}={\epsilon}_{n+1}-{\epsilon}_{n+1}^{p,(k)}-{\left({D}^{e}\right)}^{-1}:{\sigma}^{d}\left({\stackrel{̄}{\epsilon}}_{n+1}^{\mathit{ed},(k)}\right)\]

En appliquant la méthode de Newton, on calcule la nouvelle approximation de \({\stackrel{̄}{\epsilon}}_{n+1}^{\mathit{ed},(k+1)}\) comme:

(3784)\[{\stackrel{̄}{\epsilon}}_{n+1}^{\mathit{ed},(k+1)}={\stackrel{̄}{\epsilon}}_{n+1}^{\mathit{ed},(k)}+\Delta {\stackrel{̄}{\epsilon}}_{n+1}^{\mathit{ed},(k)}\]

Où

(3784)\[\Delta {\stackrel{̄}{\epsilon}}_{n+1}^{\mathit{ed},(k)}=-{\left(\frac{\partial {R}_{\sigma}^{(k)}}{\partial {\stackrel{̄}{\epsilon}}_{n+1}^{\mathit{ed}}}\right)}^{-1}{R}_{\sigma}^{(k)}\]

Le module tangent correspondant peut être exprimé entièrement par les modules découplés ():

(3784)\[\left(\frac{\partial {R}_{\sigma}^{(k)}}{\partial {\stackrel{̄}{\epsilon}}_{n+1}^{\mathit{ed}}}\right)=-{D}^{\mathit{ep},(k)}\frac{\partial {\epsilon}_{n+1}^{d,(k)}}{\partial {\stackrel{̄}{\epsilon}}_{n+1}^{\mathit{ed},(k)}}-{D}^{d,(k)}=-\left({D}^{d,(k)}+{D}^{\mathit{ep},(k)}\left(1-{\left({D}^{e}\right)}^{-1}{D}^{d,(k)}\right)\right)\]

On rappelle ici que la déformation totale \({\epsilon}_{n+1}\) , issue du processus global des itérations de Newton, est constante pour le processus itératif local déterminant \({\stackrel{ˉ}{\varepsilon}}_{n+1}^{\mathit{ed}}\) . L’algorithme du calcul des variables à l’itération \(k+1\) s’effectue de la manière suivante:

0. Initialisation: \({\stackrel{ˉ}{\varepsilon}}_{n+1}^{\mathit{ed},(k=0)}={\varepsilon}_{n+1}-{\varepsilon}_{n}^{p}\)

1. Boucle sur les \(k=0,\mathrm{...},{N}_{max}\)

2. Calcul de l’incrément:

\({\stackrel{ˉ}{\varepsilon}}_{n+1}^{\mathit{ed},(k+1)}={\stackrel{ˉ}{\varepsilon}}_{n+1}^{\mathit{ed},(k)}+\Delta {\stackrel{ˉ}{\varepsilon}}_{n+1}^{\mathit{ed},(k)}\) où \(\Delta {\stackrel{ˉ}{\varepsilon}}_{n+1}^{\mathit{ed},(k)}={({D}^{d,(k)}+{D}^{\mathit{ep},(k)}(1-{({D}^{e})}^{-1}{D}^{d,(k)}))}^{-1}{R}_{\sigma}^{(k)}\)

3. \({\sigma}_{n+1}^{d,(k+1)}={\sigma}^{d}{({\stackrel{ˉ}{\varepsilon}}_{n+1}^{\mathit{ed},(k+1)})}_{e}\)

4. \({\varepsilon}_{n+1}^{d,(k+1)}={\stackrel{ˉ}{\varepsilon}}_{n+1}^{\mathit{ed},(k+1)}-{\varepsilon}_{n+1}^{\mathit{el},(k+1)}\) où \({\varepsilon}_{n+1}^{\mathit{el},(k+1)}={({D}^{e})}^{-1}:{\sigma}^{d}({\stackrel{ˉ}{\varepsilon}}_{n+1}^{\mathit{ed},(k+1)})\)

5. \({\sigma}_{n+1}^{p,(k+1)}={\sigma}^{p}{({\varepsilon}_{n+1}-{\varepsilon}_{n+1}^{d,(k+1)})}_{e}\)

6. \({R}_{\sigma}^{(k+1)}={\sigma}_{n+1}^{p,(k+1)}-{\sigma}_{n+1}^{d,(k+1)}\)

7. Si \(\mid \Delta {\stackrel{ˉ}{\varepsilon}}_{n+1}^{\mathit{ed},(k)}{R}_{\sigma}^{(k+1)}\mid >\mathit{tol}\) reprend la boucle au 1), sinon FIN

**Remarque 1: les seuls endroits où les spécificités physiques des modèles couplés sont prises en compte sont les points 3) et 5), où on fait appel aux modules de plasticité et d’endommagement, respectivement.*

**Remarque 2: la valeur de* \({\varepsilon}_{n+1}\) est considérée constante dans le processus d’équilibrage des contraintes et lors de la variation de \({\stackrel{ˉ}{\varepsilon}}_{n+1}^{\mathit{ed}}\) on varie en réalité la valeur de \({\varepsilon}_{n+1}^{p}\) . Cette déformation plastique est en principe déterminée par le module d’élasto-plasticité comme une fonction de \({\varepsilon}_{n+1}\) (voir ()). En réalité, \({\stackrel{ˉ}{\varepsilon}}_{n+1}^{\mathit{ed}}\) est considérée ici comme une variable et son lien avec \({\varepsilon}_{n+1}^{p}\) n’est rétabli qu’à la convergence du processus itératif. Ainsi:

(3784)\[{\stackrel{ˉ}{\varepsilon}}_{n+1}^{\mathit{ed},(k)}\mathrm{\ne }{\varepsilon}_{n+1}-{\varepsilon}^{p}({\varepsilon}_{n+1}-{\varepsilon}_{n+1}^{d,(k)})\]

Construction du module tangent cohérent

Enfin, afin d’obtenir une performance optimale dans le processus itératif global on doit calculer le module tangent cohérent, une fois le processus décrit ci-dessus convergé, i.e. \({\sigma}_{n+1}={\sigma}_{n+1}^{p}={\sigma}_{n+1}^{d}\) :

(3784)\[{D}_{n+1}^{\mathit{epd}}=\frac{\partial {\sigma}_{n+1}}{\partial {\varepsilon}_{n+1}}\]

D’abord on observe que si on dérive () par rapport à \({\sigma}_{n+1}^{p}\) on peut écrire:

(3784)\[I={D}_{n+1}^{\mathit{ep}}\left({D}_{n+1}^{\mathit{epd}}-\frac{\partial {\epsilon}_{n+1}^{d}}{\partial {\sigma}_{n+1}^{p}}\right)\]

et que de la même manière, en dérivant () par rapport à \({\sigma}_{n+1}^{d}\) on obtient:

(3784)\[I={D}_{n+1}^{d}\left({D}_{n+1}^{\mathit{epd}}-\frac{\partial {\epsilon}_{n+1}^{p}}{\partial {\sigma}_{n+1}^{d}}\right)\]

avec, en fin du processus itératif \({\sigma}_{n+1}={\sigma}_{n+1}^{p}={\sigma}_{n+1}^{d}\) , on obtient finalement:

(3784)\[\frac{\partial {\epsilon}_{n+1}^{d}}{\partial {\sigma}_{n+1}}={({D}_{n+1}^{\mathit{epd}})}^{-1}-{({D}_{n+1}^{\mathit{ep}})}^{-1}\]

En utilisant ensuite \({\epsilon}_{n+1}={\epsilon}_{n+1}^{e}+{\epsilon}_{n+1}^{p}+{\epsilon}_{n+1}^{d}\) et donc que:

(3784)\[\frac{\partial {\epsilon}_{n+1}}{\partial {\sigma}_{n+1}}=\frac{\partial {\epsilon}_{n+1}^{e}}{\partial {\sigma}_{n+1}}+\frac{\partial {\epsilon}_{n+1}^{p}}{\partial {\sigma}_{n+1}}+\frac{\partial {\epsilon}_{n+1}^{d}}{\partial {\sigma}_{n+1}}\]

on obtient le module tangent cohérent directement comme:

(3784)\[{D}_{n+1}^{\mathit{epd}}={\left({\left({D}_{n+1}^{\mathit{epd}}\right)}^{-1}+{\left({D}_{n+1}^{\mathit{epd}}\right)}^{-1}+{\left({D}_{n+1}^{\mathit{epd}}\right)}^{-1}\right)}^{-1}\]

Limitations du coupleur

Il existe quelques conditions à satisfaire pour l’utilisation de cette méthode de couplage concernant les aspects à la fois physiques et numériques:

• Variables internes distinctes:

Comme précisé dans l’introduction du cadre thermodynamique, les deux modèles couplés ne doivent pas partager les mêmes variables internes. La méthode est conçue pour pouvoir combiner des modèles représentant des phénomènes physiques découplés, ce qui est jugé être le cas dans la plupart des applications cibles.

• Existence des modules tangents cohérents:

Pour que l’approche fonctionne, il faut absolument qu’on dispose des modules tangents cohérents pour les deux modules, élasto-plastique et endommageable. Sous «module tangent cohérent» on sous-entend qu’il est la dérivée cohérente avec la discrétisation temporelle employée du tenseur de contraintes à l’instant donné par rapport au tenseur de déformation au même instant.

• Non-singularité des modules tangents:

On est obligé à plusieurs reprise d’inverser les matrices tangentes de modules individuels. Si le comportement d’un module a une pente non-linéaire plate (élasto-plasticité parfaite par exemple) il est nécessaire de lui donner une valeur non-nulle. Cette valeur non-nulle doit être choisie telle qu’elle soit suffisamment grande pour que l’algorithme fonctionne et suffisamment petite pour qu’elle soit cohérente avec le modèle physique étudié.

Modèles spécifiques

Modèles globaux

Rappel des modèles globaux

Dans ce contexte, un modèle de comportement de plaque dit global ou d’élément de structure en général, signifie que la loi de comportement s’écrit directement en termes de relation entre les contraintes généralisées et les déformations généralisées. L’approche globale de modélisation du comportement des structures s’applique notamment aux structures composites, par exemple le béton armé (voir Figure 1), et représente une alternative aux approches dites locales ou semi-globales, qui sont des modélisations plus fines mais plus coûteuses (voir [bib6] et [bib7]). Dans l’approche locale on utilise une modélisation fine pour chacune des phases (acier, béton) et leurs interactions (adhérence) et dans l’approche semi-globale on exploite l’élancement de la structure pour simplifier la description de la cinématique, ce qui aboutit aux modèles PMF (Poutre Multi-Fibre) ou coques multi-couches.

L’intérêt du modèle global réside dans le fait que l’élément fini correspondant ne nécessite qu’un point d’intégration dans l’épaisseur et surtout dans l’obtention d’un comportement homogénéisé. Cet avantage est encore plus important dans l’analyse du béton armé, puisqu’on contourne les problèmes de localisation rencontrés lors de la modélisation du béton sans armatures. Évidemment, un modèle global représente les phénomènes locaux d’une manière grossière.

Figure 3.1-a. Dalle en béton armé.

Implantation du couplage GLRC_DM/VMIS_

L’hypothèse la plus importante qui est faite dans ce couplage particulier est que l’élasto-plasticité ne peut être activée qu’en sollicitant la plaque en membrane. Autrement dit, en flexion pure la structure ne plastifiera jamais. D’une part, cette simplification se justifie par les applications cibles, où les sollicitations sont plus importantes en membrane qu’en flexion et où on ne s’attend donc pas à des influences dominantes de la flexion. D’autre part, l’hypothèse s’impose pour des raisons techniques, car actuellement on ne dispose pas d’un modèle «global» élasto-plastique convenable pour le couplage avec un modèle «global» d’endommagement.

Les modèles globaux sont formulés en terme de contraintes et déformations généralisées, ce qui nous impose, dans le coupleur (voir §2), de remplacer les contraintes, \(\sigma\) , et déformations, \(\epsilon\) , par:

(3792)\[\begin{split}\sigma \to \Sigma =\left(\begin{array}{c}{N}_{xx}\\ {N}_{yy}\\ {N}_{xy}\\ {M}_{xx}\\ {M}_{yy}\\ {M}_{xy}\end{array}\right)\end{split}\]

où les déformations généralisées, \({\rm E}\) , se décomposent en extensions membranaires \(ϵ\) et les courbures \(\kappa\) tandis que les contraintes généralisées sont composées des efforts membranaires \(N\) et des moments fléchissants \(M\) .

Il est important de noter que dans ce cadre on fait l’hypothèse de plaques minces où la distorsion transverse et l’effort tranchant sont négligeables. Ci-dessous on rappelle la cinématique de Hencky-Mindlin (voir [bib9] pour les détails) pour les coques et les plaques ainsi que la définition des contraintes généralisées:

(3792)\[\begin{split}(\begin{array}{}{U}_{1}({x}_{1,}{x}_{2,}z)\\ {U}_{2}({x}_{1,}{x}_{2,}z)\\ {U}_{z}({x}_{1,}{x}_{2,}z)\end{array})=\underset{\begin{array}{}\mathrm{cinématique}\mathrm{de}\mathrm{plaque}\\ {u}^{s}\in {V}_{s}\end{array}}{\underset{\underbrace{}}{(\begin{array}{}{u}_{1}({x}_{1,}{x}_{2})\\ {u}_{2}({x}_{1,}{x}_{2})\\ {u}_{z}({x}_{1,}{x}_{2})\end{array})+z(\begin{array}{}{\theta}_{2}({x}_{1,}{x}_{2})\\ -{\theta}_{1}({x}_{1,}{x}_{2})\\ 0\end{array})}}+\underset{\begin{array}{}\mathrm{déplacement}\mathrm{complémentaire}\\ {u}^{c}\in {V}_{c}\end{array}}{\underset{\underbrace{}}{(\begin{array}{}{u}_{1}^{c}({x}_{1,}{x}_{2,}z)\\ {u}_{2}^{c}({x}_{1,}{x}_{2,}z)\\ {u}_{z}^{c}({x}_{1,}{x}_{2,}z)\end{array})}}\end{split}\]

où \(U={({U}_{1}{U}_{2}{U}_{z})}^{T}\) est le champ de déplacement en 3D, \(u={({u}_{1}{u}_{2}{u}_{z})}^{T}\) le déplacement du feuillet moyen et \({\theta}_{x}\) , \({\theta}_{y}\) ses rotations. Ainsi, le tenseur de déformation, défini comme:

(3792)\[{\varepsilon}_{ij}=\frac{1}{2}(\frac{\partial {U}_{i}}{\partial {x}_{j}}+\frac{\partial {U}_{j}}{\partial {x}_{i}}),i,j=\mathrm{1..3}\]

s’écrit aussi comme:

(3792)\[\begin{split}\begin{array}{}{\varepsilon}_{11}=\underset{{\varepsilon}_{11}^{s}}{\underset{\underbrace{}}{{\epsilon}_{11}+z{\kappa}_{11}}}+{u}_{1,1}^{c}\\ {\varepsilon}_{22}=\underset{{\varepsilon}_{22}^{s}}{\underset{\underbrace{}}{{\epsilon}_{22}+z{\kappa}_{22}}}+{u}_{2,2}^{c}\\ {\varepsilon}_{12}=\underset{{\varepsilon}_{12}^{s}}{\underset{\underbrace{}}{{\epsilon}_{12}+z{\kappa}_{12}}}+\frac{1}{2}({u}_{2,1}^{c}+{u}_{1,2}^{c})\\ {\varepsilon}_{\mathrm{1z}}={\varepsilon}_{\mathrm{1z}}^{c}=\frac{1}{2}{u}_{3,1}^{c}\\ {\varepsilon}_{\mathrm{2z}}={\varepsilon}_{\mathrm{2z}}^{c}=\frac{1}{2}{u}_{3,2}^{c}\\ {\varepsilon}_{zz}={\varepsilon}_{zz}^{c}={u}_{3,z}^{c}\end{array}\end{split}\]

où \(\epsilon\) est le tenseur de l’extension membranaire, défini dans le plan:

(3792)\[{\epsilon}_{ij}=\frac{1}{2}(\frac{\partial {u}_{i}}{\partial {x}_{j}}+\frac{\partial {u}_{j}}{\partial {x}_{i}}),i,j=\mathrm{1..2}\]

et \(\kappa\) le tenseur de variation de courbure, défini dans le plan:

(3792)\[{\kappa}_{11}=\frac{\partial {\theta}_{2}}{\partial {x}_{1}}\]

relations auxquelles se rajoute l’hypothèse de contraintes planes \({\sigma}_{zz}=0\) , \({\sigma}_{\mathrm{1z}}=0\) , \({\sigma}_{\mathrm{2z}}=0\) qui déterminera le champ de déplacement complémentaire \({\mathrm{u}}^{c}\in {V}_{c}\) . Dans la théorie utilisée ici, on n’introduit que deux composantes de rotation \({\theta}_{x}\) et \({\theta}_{y}\) , ce qui implique que le tenseur de variation de courbure est 2D et n’a que 3 composantes indépendantes.

Parmi les deux modèles couplés, GLRC_DM est un modèle global et formulé donc directement en terme de \(\Sigma =\Sigma (E)\) , mais les modèles VMIS_ISOT_LINE et VMIS_CINE_LINE sont les lois de comportement 3D classiques. Vu qu’on souhaite surtout que la partie élasto-plastique soit représentée en membrane, on n’applique la loi de comportement VMIS_ que sur la relation \(N=N(\epsilon )\) , puis on suppose que le comportement en flexion reste élastique linéaire. On a donc:

Endommagement :

(3792)\[{\Sigma}^{d}(\epsilon ,\kappa )={\Sigma}_{\text{GLRC\_DM}}(\epsilon ,\kappa )\]

Elasto-plasticité :

(3792)\[\begin{split}{\Sigma}^{p}(\epsilon ,\kappa )=(\begin{array}{c}{N}_{\mathit{VMIS}}(\epsilon )\\ {H}_{\mathit{ELAS}}:\kappa \end{array})\end{split}\]

où \({\Sigma}_{\text{GLRC\_DM}}(\epsilon ,\kappa )\) représente la loi de comportement globale GLRC_DM, \({N}_{\mathit{VMIS}}(\epsilon )\) la loi de Von Mises en contraintes planes et \({H}_{\mathit{ELAS}}\) le tenseur élastique agissant sur la courbure (voir [bib9] pour sa forme exacte). Ces deux modules, d’endommagement et d’élasto-plasticité, sont ensuite utilisés dans le coupleur (voir §2) en remplaçant \({\sigma}^{p}\to {\Sigma}^{p}\) et \({\sigma}^{d}\to {\Sigma}^{d}\) .

Variables internes

Les variables internes du modèle couplé sont stockées en série: les premières correspondent au modèle d’endommagement, qui sont suivies de celle du modèle élasto-plastique. Le coupleur a besoin de six (6) variables internes pour stocker le tenseur du type déformation, \({\stackrel{ˉ}{\varepsilon}}^{\mathit{ed}}\) , qui pilote la boucle interne (voir §2.3). Dans le cas du couplage GLRC_DM/VMIS_ on en rajoute quatre (4) à cause de l’algorithme permettant de satisfaire la condition des contraintes planes. Ces quatre variables internes ont la même signification que celles utilisées pour la méthode DEBORST (voir [bib12]).

Pour GLRC_DM/VMIS_ on a donc 24 variables internes: V1-V7: GLRC_DM; V8-V14: VMIS_ (on ne se sert que de V5 et V6 pour VMIS_ISOT_LINE); V15-V20: tenseur de déformation \({\stackrel{ˉ}{\epsilon}}^{\text{ed}}\) ;V21-V24: méthode DEBORST.

Validation

Ce modèle est validé par les tests SSNS106F,G (voir [bib8]).

Bibliographie

1. 10. LUBLINER, «Plasticity theory», MacMillan, New York, 1990.

2. J.LEMAITRE, J.L.CHABOCHE, «Mécanique d’endommagement»,

3. 1. IBRAHIMBEGOVIC, D. MARKOVIC, F. GATUINGT, «Constitutive model of coupled damage-plasticity and its finite element implementation», Eur. J. Finite Elem., 12 (4), 2003, p. 381-405.

4. 4. MARKOVIC, A. IBRAHIMBEGOVIC, «Complementary energy based FE modelling of coupled elasto-plastic and damage behavior for continuum microstructure computations», Comp. Methods Appl. Engrg., 195, 2006, p. 5077-5093.

5. D. MARKOVIC , «Modélisation multi-échelles de structures hétérogènes aux comportements anélastiques non-linéaires», thèse ENS de Cachan, 2004, disponible sur `http://tel.archives-ouvertes.fr/tel-00133643 <http://tel.archives-ouvertes.fr/tel-00133643>`_(version pdf) et dans le bureau V015 (version papier).

6. S.Moulin. Modélisation des structures en béton armé sous chargement sismique. Note HT-62/04/025/A, 12/2004.

7. S.Moulin. F. VOLDOIRE Étude d’une poutre en béton armé sous chargement de flexion. Note HT-62/05/013/A, 9/2006.

8. [V6.05.106] SSNS106 – Endommagement d’une plaque plane sous sollicitations variées avec la loi de comportement GLRC_DM.

9. [R3.07.03]– Éléments de plaque DKT, DST, DKQ, DSQ et Q4G.

10. [R7.01.31]– Loi de comportement de plaque en béton armé GLRC_DAMAGE.

11. [R7.01.32]– Loi de comportement de plaque en béton armé GLRC_DM.

12. [R5.03.03]– Prise en compte de l’hypothèse de contraintes planes dans les comportements non-linéaires.

13. [R5.03.02]– Intégration des relations de comportement élasto-plastique de von Mises.

14. [R7.01.09]– Loi de comportement ENDO_ORTH_BETON.

15. [R7.01.19]– Modélisation du couplage fluage/plasticité pour le béton.

Historique des versions du document

Version Aster	Auteur(s) ou contributeur(s), organisme	Description des modifications
9.3	D.MARKOVIC EDF-R&D/AMA	Texte initial