r5.03.29 Lois de comportement d’endommagement ductile plastique etviscoplastique GTN et VISC_GTN

Résumé:

Ce document décrit le modèle d’endommagement ductilede Gurson – Tvergaard – Needleman (GTN). Il est disponible en version locale et à gradient de variables internes (GRAD_VARI ou GRAD_INCO). Lamodélisation des grandes déformations est accessible via la cinématique GDEF_LOG. Le comportement peut être indépendant du temps (GTN) ou viscoplastique (VISC_GTN). Dans ce dernier cas, la plasticité évolue selon un modèle de Norton dans lequel intervient la fonction seuil.

Modèle continu

Équations de comportement : modèle plastique local

On se place en grandes déformations logarithmiques, en notant \(E\) la déformation mécanique (logarithmique) et \(T\) la contrainte associée (par dualité). Les variables internes du modèle consistent en la déformation plastique \({E}^{p}\) , la variable d’écrouissage \(\kappa\) , l’endommagement \(d\) (qui dépend de la porosité \(f\) ) et la déformation plastique cumulée (qui gouverne une partie de la porosité de germination).

Tout d’abord, la relation contrainte – déformation s’appuie sur une décomposition additive en une partie plastique et une partie élastique:

(2064)\[ \begin{align}\begin{aligned}T=K\mathit{tr}(E-{E}^{p})\mathit{Id}+2\mu \mathit{dev}(E-{E}^{p})\\\textrm{avec}\\3K=\frac{E}{1-2\nu }` et :math:`2\mu =\frac{E}{1+\nu }\end{aligned}\end{align} \]

où :math:`` et :math:`` désignent respectivement la trace et le déviateur d’un tenseur d’ordre deux.

La fonction seuil de plasticité notée \(F\) dépend du niveau d’écrouissage \(\kappa\) à traversla fonction d’écrouissage :math:`` ,d’une mesure (isotrope) de l’intensité des contraintes notée \({T}_{\text{*}}\) qui dépend aussi de l’endommagement et de la déformation à travers le Jacobien de la transformation \(J=\exp(\mathrm{tr}E)\) . Elle s’écrit comme suit:

(2065)\[F(T,\kappa ;d,E)=\frac{{T}_{\text{*}}(T,d)}{J(E)}-R(\kappa )\]

La fonction d’écrouissage est identique à celle déployée dans la loi VMIS_ISOT_NL [R5.03.33], soit:

(2066)\[R(\kappa )={R}_{0}+{R}_{H}\kappa +{R}_{1}(1-{e}^{-{\gamma}_{1}\kappa })+{R}_{2}(1-{e}^{-{\gamma}_{2}\kappa })+{R}_{K}{({p}_{0}+\kappa )}^{{\gamma}_{K}}\]

Elle dépend au plus de neuf paramètres dont certains peuvent rester nuls si on ne souhaite pas activer tous les termes de l’expression (choix par défaut).

Remarque: on rappelle qu’il n’est pas possible d’introduire un plateau de Lüders contrairement à la loi VMIS_ISOT_NL.

La définition de la contrainte de Gurson \({T}_{\text{*}}\) est implicite. À contrainte \(T\) et endommagement \(d\) donnés, \({T}_{\text{*}}\) est la solution de l’équation suivante:

(2067)\[G(T,{T}_{\text{*}},d)=0;G(T,{T}_{\text{*}},f)\equiv {\left(\frac{{T}_{\mathit{eq}}}{{T}_{\text{*}}}\right)}^{2}+2d\mathrm{chm}\left(\frac{3}{2}{q}_{2}\frac{{T}_{H}}{{T}_{\text{*}}}\right)-{(1-d)}^{2};\mathrm{chm}(x)\equiv \mathrm{ch}(x)-1\]

où \({q}_{1}\) et \({q}_{2}\) sont deux paramètres matériau et où on a introduit les invariants de la contrainte:

(2068)\[ \begin{align}\begin{aligned}{T}_{H}=\frac{1}{3}\mathrm{tr}(T)\\\textrm{et}\\{T}_{\mathit{eq}}=\sqrt{\frac{3}{2}\mathrm{dev}(T):\mathrm{dev}(T)}\end{aligned}\end{align} \]

L’évolution de la déformation plastique est gouvernée par l’équation d’écoulement. Lorsque la fonction seuil est dérivable par rapport à \(T\) , c’est-à-dire quand \(T\ne 0\) , la direction d’écoulement est normale à la surface seuil:

(2069)\[\text{si}T\ne 0\text{}{\dot{E}}^{p}=\dot{\kappa}N\text{};\text{}N=\frac{\partial F}{\partial T}=-\frac{1}{J}{\left(\frac{\partial G}{\partial {T}_{\text{*}}}\right)}^{-1}\left(\frac{\partial G}{\partial T}\right)\]

Dans le cas singulier où \(T=0\) , la dérivée est généralisée via des notions d’analyse convexe (sous-gradient). La direction d’écoulement n’est alors assujettie qu’à la condition suivante:

(2070)\[\text{si}T=0\text{}{\Pi}_{N}(\dot{{E}^{p}},d)\le \frac{\dot{\kappa}}{J}\text{}\]

La fonction d’appui \({\Pi}_{N}\) a l’expression suivante, avec des notations similaires à () :

(2071)\[{\Pi}_{N}(\dot{{E}^{p}},d)=\frac{2}{{q}_{2}}\mathrm{arcosh}(1+{C}^{\text{*}})\left|\dot{{E}_{H}^{p}}\right|+\frac{2}{3}\left(1-d\right)\sqrt{1-\gamma {C}^{\text{*}}}\dot{{E}_{\mathit{eq}}^{p}}\]

(2072)\[{C}^{\text{*}}=-\mathrm{\chi }+\sqrt{{\mathrm{\chi }}^{2}+\frac{2}{\gamma}(\mathrm{\chi }-1)};\mathrm{\chi }=\frac{{({q}_{2}\dot{{E}_{\mathit{eq}}^{p}})}^{2}+\frac{{(3\dot{{E}_{H}^{p}})}^{2}}{d}}{{({q}_{2}\dot{{E}_{\mathit{eq}}^{p}})}^{2}};\gamma =\frac{2d}{{\left(1-d\right)}^{2}}\]

L’évolution de la variable d’écrouissage est fixée par la condition de cohérence:

(2073)\[\dot{\kappa}\ge 0;F(T,\kappa )\le 0;\dot{\kappa}F(T,\kappa )=0\]

Il reste enfin à définir l’évolution de l’endommagement. Le modèle est spécifique à EDF; il coïncide avec le modèle de GTN classique en charge mais il précise également ce qui se passe en décharge. Plus exactement, l’endommagement est défini comme suit:

(2074)\[d={q}_{1}f+{d}^{\mathit{co}};f={f}^{nu}+{f}^{\mathit{gr}}\]

Où \(f\) est la porosité, somme de la porosité de germination (nucleation en anglais) \({f}^{nu}\) et de la porosité de croissance (growth) \({f}^{\mathit{gr}}\) , et \({d}^{\mathit{co}}\) est l’endommagement de coalescence. A l’état initial (naturel), on a:

(2075)\[{f}^{nu}(0)=0;{f}^{\mathit{gr}}(0)={f}_{0};{d}^{\mathit{co}}(0)=0\text{}\Rightarrow \text{}f(0)={f}_{0}\]

Avec \({f}_{0}\) la porosité initiale.

La porosité decroissance découle de la conservation du volume plastique:

(2076)\[\dot{{f}^{\mathit{gr}}}=(1-f)\mathrm{tr}\dot{{E}^{p}};{f}^{\mathit{gr}}\ge {f}_{0}\]

Elle ne peut descendre en-dessous de \({f}_{0}\) mais elle est réversible (on peut refermer les cavités).

Pour ce qui concerne la porosité de germination, deux familles de modèles sont implantées:

(2077)\[{f}^{nu}={f}_{(1)}^{nu}+{f}_{(2)}^{nu}\]

La premièrefamille adopte une évolution de la porosité liée à celle de la variable d’écrouissage:

(2078)\[{\dot{f}}_{(1)}^{nu}={B}_{n}(\kappa )\dot{\kappa}\text{}\Rightarrow \text{}{f}_{(1)}^{nu}(\kappa )=\underset{0}{\overset{\kappa}{\int}}{B}_{n}(s)\mathit{ds}\]

Pour recouvrir deux modèles utilisés dans la littérature, on choisit la fonction \({B}_{n}\) suivante:

(2079)\[{B}_{n}(\kappa )=\frac{{f}_{N}}{{s}_{N}\sqrt{2\pi }}\exp\left[-\frac{1}{2}{\left(\frac{\kappa -{\kappa}_{N}}{{s}_{N}}\right)}^{2}\right]+\frac{{c}_{0}}{{\kappa}_{f}-{\kappa}_{i}}{\mathrm{\chi }}_{[{\kappa}_{i};{\kappa}_{f}]}\]

où \({\chi }_{I}\) vaut 1 dans l’intervalle \(I\) et 0 à l’extérieur. Compte tenu des primitives de ces deux termes, on en déduit:

(2080)\[{f}_{(1)}^{nu}=\frac{{f}_{N}}{2}\left[\mathrm{erf}\left(\frac{\kappa -{\kappa}_{N}}{\sqrt{2}{s}_{N}}\right)+\mathrm{erf}\left(\frac{{\kappa}_{N}}{\sqrt{2}{s}_{N}}\right)\right]+{c}_{0}\phantom{\rule{0.5em}{0ex}}\min\left(\frac{⟨\kappa -{\kappa}_{i}⟩}{{\kappa}_{f}-{\kappa}_{i}},1\right)\]

où les crochets de McCauley \(⟨x⟩\) désignent la partie positive de \(x\) .

La seconde famille de modèles de germination dépend de la déformation plastique cumulée qui est définie comme suit:

(2081)\[{\dot{E}}_{\mathit{eq}}^{p}\equiv \sqrt{\frac{2}{3}{\dot{E}}^{p}:{\dot{E}}^{p}}\]

La fonction de germination implantée est proportionnelle à \({E}_{\mathit{eq}}^{p}\) à partir d’un seuil:

(2082)\[{f}_{(2)}^{nu}({E}_{\mathit{eq}}^{p})={b}_{0}⟨{E}_{\mathit{eq}}^{p}-{E}_{c}^{p}⟩\]

Comme la variable d’écrouissage \(\kappa\) et la déformation plastique cumulée \({E}_{\mathit{eq}}^{p}\) sont croissantes par définition, il en va de même de la porosité de germination; son évolution est irréversible.

Il reste à définir l’évolution de l’endommagement de coalescence. Pour cela, on introduit la définition (usuelle pour GTN) de la porosité en régime de coalescence \({f}^{\text{*}}\) ,où \({f}_{c}\) désigne le seuil de porosité à coalescence et \({f}_{F}\) la porosité à rupture:

(2083)\[\begin{split}{f}^{\text{*}}=\lbrace \begin{array}{cc}f& \text{si}f<{f}_{c}\\ {f}_{c}+\delta (f-{f}_{c})& \text{si}f\ge {f}_{c}\end{array}\end{split}\]

L’endommagement de coalescence obéit alors à la loi d’évolution suivante:

(2084)\[{d}^{\mathit{co}}(t)=\underset{\tau \le t}{\max}[{q}_{1}{f}^{\text{*}}(\tau )-{q}_{1}f(\tau )]\]

Il est lui aussi croissant par définition, donc irréversible.

Pour des raisons numériques, on introduit également un seuil d’endommagement à rupture \({d}_{f}\) au-delà duquel on considère que le point est cassé, c’est-à-dire \(T=0\) (il y a donc un saut de contrainte qui peut perturber la résolution des équations d’équilibre mais facilite l’intégration du comportement). Une fois qu’un point est cassé, il le reste définitivement. Ce mécanisme n’est pas activé par défaut (\({d}_{f}=1\) ).

Prise en compte de la viscosité

En présence de viscosité, l’écoulement plastique est différé. Plus exactement, la vitesse d’évolution de l’écrouissage n’est plus fixée par la condition de cohérence () mais elle est dorénavant fonction de l’intensité de la fonction seuil suivant la loi de Norton:

(2085)\[\dot{\kappa}={\left(\frac{⟨F(T,\kappa )⟩}{\stackrel{~}{K}}\right)}^{n}\]

Les autres équations du modèlerestent inchangées. Lorsque \(\stackrel{~}{K}\to 0\) , on retrouve le modèle élastoplastique.

En général, le paramètre \(\stackrel{~}{K}\) coïncide avec le paramètre \(K\) du modèle de Norton. Toutefois, lorsque l’endommagement \(d\) s’approche de 1, l’écart entre \({T}_{\text{*}}\) et \(T\) s’accroît (puisque \(T\) tend vers zéro) et la viscoplasticité peut devenir trop brutale pour conserver un rôle efficace, notamment en présence d’erreurs numériques d’intégration. C’est pourquoi on propose une modification de la loi de Norton usuelle en couplant la viscoplasticité avec l’endommagement, au-delà d’un seuil \({d}_{v}\) :

(2086)\[\begin{split}\begin{array}{cc}\text{si}d\le {d}_{v}& \stackrel{~}{K}=K\\ \text{si}d\ge {d}_{v}& \text{}\stackrel{~}{K}=K\frac{1-{d}_{v}}{1-d}\end{array}\end{split}\]

Enfin, la viscoplasticité permet d’accommoder des sauts de contraintes potentiels en les lissant dans le temps. Il est important de préserver ce mécanisme y compris pour les points cassés. Dans ce cas, l’annulation des contraintes (point cassé) est différée dans le temps, de manière similaire à un mécanisme de relaxation, à condition d’avoir choisi \({d}_{v}<1\) . Plus précisément, en régime viscoplastique de point cassé, la plasticité et les contraintes obéissent à la loi d’évolution suivante:

(2087)\[\dot{\kappa}=0;\dot{{E}^{p}}=0\]

(2088)\[\begin{split}\begin{array}{cc}\text{si}\Vert T\Vert \ge K(1-{d}_{v})& \text{}\frac{\dot{T}}{E}=-{\left[\frac{\Vert T\Vert }{K(1-{d}_{v})}\right]}^{n}\frac{T}{\Vert T\Vert }\\ \text{si}\Vert T\Vert \le K(1-{d}_{v})& \frac{\dot{T}}{E}=-\frac{T}{K(1-{d}_{v})}\end{array}\end{split}\]

Le passage à un exposant égal à 1 lorsque \(\Vert T\Vert \le K(1-{d}_{v})\) évite de trop différer l’effacement de la contrainte.

Prise en compte du gradient d’écrouissage

Une formulation non locale à gradient de variable interne est adaptée lorsque le gradient d’écrouissage \(\nabla \kappa\) devient important. Elle consiste à introduire un terme supplémentaire \({\Phi}^{\mathit{grad}}\) dans l’énergie libre qui reflète les interactions entre points matériels voisins, voir [R5.04.01]:

(2089)\[{\Phi}^{\mathit{grad}}(\nabla \kappa )=\frac{1}{2}c{\left(\nabla \kappa \right)}^{2}\]

En utilisant la formulation relaxée augmentée présentée dans [R5.04.01], l’impact sur la relation de comportement se traduit par une modification de la force thermodynamique associée à la variable d’écrouissage donc, en pratique, par une modification de la fonction d’écrouissage :

(2090)\[\stackrel{~}{R}(\kappa )=R(\kappa )+r\kappa -\varphi ;\varphi =\lambda +\mathit{ra}\]

où \(\lambda\) est le multiplicateur de Lagrange associé à la relaxation, \(r>0\) le coefficient d’augmentation et \(a\) l’interprétation de la variable d’écrouissage à l’échelle de la structure. Ces trois quantités sont des données pour ce qui concerne la relation de comportement. Au final, la fonction d’écrouissage est corrigée d’un terme affine, ce qui ne rajoute donc pas de complexité par rapport au traitement numérique du modèle. En revanche, le signe et l’amplitude de \(\varphi\) n’étant pas fixés, cela explique qu’il faille prendre en compte le cas d’un écoulement singulier pour lequel \({T}_{\text{*}}=0\) .

Intégration numérique

Pour ce qui concerne l’intégration de la loi de comportement, c’est-à-dire le calcul des variables internes et des contraintes à histoire de déformation (mécanique) donnée, il n’y a plus lieu de distinguer s’il s’agit d’une formulation locale ou à gradient. En effet, l’impact de cette dernière se limite à modifier la fonction d’écrouissage selon (). On continuera ici à la noter \(R(\kappa )\) sans distinction.

En revanche, on distinguera les cas plastique et viscoplastique, même si on peut montrer que ce dernier se limite également à une correction de la fonction seuil dans les équations discrétisées.

Équations discrétisées dans le cas sans viscosité

La discrétisation en temps des équations de comportement s’appuie sur un schéma étagé: à endommagementfixé, on applique un schéma d’Euler implicite (même en viscoplasticité) pour calculer les variables plastiques (déformation plastique, variable d’écrouissage) ainsi que la contrainte, c’est-à-dire que les différentes variables du problème sont exprimées à l’instant final du pas de temps considéré. Puis on évalue l’évolution des différentes contributions à l’endommagementen connaissant l’évolution des variables plastiques. On renvoie à [Zhang et al., 2018] pour plus d’informations sur la discrétisation temporelle et la résolution du système algébrique correspondant.

On note \({Q}_{n}\) la valeur d’une quantité \(Q\) au début du pas de temps, \(\Delta Q\) son incrément pendant le pas de temps et (simplement) \(Q\) sa valeur en fin de pas de temps. L’état mécanique au début du pas de temps \(({E}_{n},{{E}^{p}}_{n},{\kappa}_{n},{f}_{n}^{nu},{f}_{n}^{\mathit{gr}},{d}_{n}^{\mathit{co}},{d}_{n})\) est supposé connu ainsi que l’incrément de déformation \(\Delta E\) (et donc également la déformation \(E\) ). Il s’agit alors de calculer les incréments de variables internes ainsi que la contrainte à la fin du pas de temps \(T\) .

Concernant la valeur de l’endommagement choisiepour intégrer la plasticité, [Zhang et al., 2018] proposaient d’adopter celle en début de pas de temps, soit \({d}_{\mathit{ex}}={d}_{n}\) . Pour des questions de précision d’intégration (et donc de taille de pas de temps), on enrichit ce schéma en extrapolant l’endommagementen cours de pas de temps en s’appuyant sur sa vitesse d’évolution au pas précédent. Plus précisément, on définit:

(2091)\[\dot{{d}_{n}}=\frac{{d}_{n}-{d}_{n-1}}{{t}_{n}-{t}_{n-1}};{d}_{\mathit{ex}}={d}_{n}+\theta (t-{t}_{n})\dot{{d}_{n}}\]

où la valeur de \(\theta\) est fixée par l’utilisateur sous le mot-clé COMPORTEMENT/ PARM_THETA; on recommande la valeur \(\theta =1/2\) .

Connaissant maintenant la valeur d’endommagement extrapolée, on peut procéder à l’intégration implicite de la plasticité. Si \({d}_{\mathit{ex}}\ge {d}_{f}\) , on bascule simplement en régime de point cassé et \(T=0\) . Dans le cas contraire, la discrétisation de la relation contrainte – déformation () s’écrit comme suit:

(2092)\[T=K\mathrm{tr}(E-{{E}^{p}}_{n}-\Delta {E}^{p})\mathit{Id}+2\mu \mathrm{dev}(E-{{E}^{p}}_{n}-\Delta {E}^{p})\]

On choisit de faireapparaître une contrainte élastique (ou contrainte d’essai élastique), notée \({T}^{e}\) et qui est connue:

(2093)\[T={T}^{e}-K\mathrm{tr}(\Delta {E}^{p})-2\mu \mathrm{dev}(\Delta {E}^{p});{T}^{e}=K\mathrm{tr}(E-{{E}^{p}}_{n})\mathit{Id}+2\mu \mathrm{dev}(E-{{E}^{p}}_{n})\]

La discrétisation de la fonction seuil s’écrit simplement:

(2094)\[F(T,\kappa )=\frac{{T}_{\text{*}}}{J}-R(\kappa )\]

Et la condition de cohérence vaut:

(2095)\[\Delta \kappa \ge 0;F(T,\kappa )\le 0;\Delta \kappa \phantom{\rule{0.5em}{0ex}}F(T,\kappa )=0\]

En particulier, en charge plastique, \(F=0\) , si bien qu’on peut en déduire la contrainte équivalente de GTN en fonction du niveau d’écrouissage:

(2096)\[{T}_{\text{*}}=\widehat{T}(\kappa )\equiv JR(\kappa )\]

Dans le cas d’un écoulement régulier, la discrétisation temporelle de () s’écrit:

(2097)\[\mathrm{tr}(\Delta {E}^{p})=3\frac{\lambda}{J}\Theta \left(\frac{{T}_{H}}{{T}_{\text{*}}},\frac{{T}_{\mathit{eq}}}{{T}_{\text{*}}}\right)\Lambda \left(\frac{{T}_{H}}{{T}_{\text{*}}}\right)\]

(2098)\[\mathrm{dev}(\Delta {E}^{p})=\frac{3}{2}\frac{\lambda}{J}\Theta \left(\frac{{T}_{H}}{{T}_{\text{*}}},\frac{{T}_{\mathit{eq}}}{{T}_{\text{*}}}\right)\frac{\mathrm{dev}(T)}{{T}_{\text{*}}}\]

On y a introduit les fonctions suivantes:

(2099)\[\Lambda (x)=\frac{1}{2}{q}_{2}{d}_{\mathit{ex}}\phantom{\rule{0.5em}{0ex}}\sinh\left(\frac{3}{2}{q}_{2}x\right);\Theta (x,y)={\left[{y}^{2}+3x\Lambda (x)\right]}^{-1}\]

Par substitution de () et () dans (), on en déduit:

(2100)\[\frac{{T}_{H}^{e}-{T}_{H}}{3K}=\frac{\lambda}{J}\Theta \left(\frac{{T}_{H}}{{T}_{\text{*}}},\frac{{T}_{\mathit{eq}}}{{T}_{\text{*}}}\right)\Lambda \left(\frac{{T}_{H}}{{T}_{\text{*}}}\right)\]

(2101)\[\mathrm{dev}(T)={\left[1+3\mu \frac{\lambda}{J}\Theta \left(\frac{{T}_{H}}{{T}_{\text{*}}},\frac{{T}_{\mathit{eq}}}{{T}_{\text{*}}}\right)\frac{1}{{T}_{\text{*}}}\right]}^{-1}\mathrm{dev}({T}^{e})\]

Cette dernière expression indique que \(T\) a la même direction et le même sens que \({T}^{e}\) . En outre, en substituant \(\Theta ` par son expression issue de (), on obtient la valeur de :math:`{T}_{\mathit{eq}}\) en fonction de \({T}_{H}\) et \({T}_{\text{*}}\) :

(2102)\[{T}_{\mathit{eq}}={T}_{\mathit{eq}}^{e}{\left[1+\frac{\mu}{K}\frac{{T}_{H}^{e}-{T}_{H}}{{T}_{\text{*}}\Lambda (\frac{{T}_{H}}{{T}_{\text{*}}})}\right]}^{-1}\equiv {\widehat{T}}_{\mathit{eq}}({T}_{H},{T}_{\text{*}})\]

Compte tenu de la définition de la contrainte équivalente de GTN (), on peut exprimer unlien implicite entre \({T}_{H}\) et \({T}_{\text{*}}\) en y substituant ():

(2103)\[\widehat{G}({T}_{H},{T}_{\text{*}})\equiv {\left[\frac{{\widehat{T}}_{\mathit{eq}}({T}_{H},{T}_{\text{*}})}{{T}_{\text{*}}}\right]}^{2}+2{d}_{\mathit{ex}}\cosh\left[\frac{3}{2}{q}_{2}\frac{{T}_{H}}{{T}_{\text{*}}}\right]-1-{d}_{\mathit{ex}}^{2}=0\]

Le multiplicateur plastique, égal à \(\Delta \kappa\) , se déduit de ():

(2104)\[\Delta \kappa =\widehat{\lambda}({T}_{H},{T}_{\text{*}})\equiv J\frac{{T}_{H}^{e}-{T}_{H}}{3K\Lambda \left(\frac{{T}_{H}}{{T}_{\text{*}}}\right)}\frac{1}{\Theta \left(\frac{{T}_{H}}{{T}_{\text{*}}},\frac{{\widehat{T}}_{\mathit{eq}}({T}_{H},{T}_{\text{*}})}{{T}_{\text{*}}}\right)}\]

Lorsque l’écoulement est singulier, \({T}_{\text{*}}=0\) donc \(T=0\) , ce qui permet de calculer \(\Delta {E}^{p}\) via (). Il reste alors à vérifier que l’écoulement appartient bien au cône prescrit par () soit, après discrétisation:

(2105)\[{\Pi}_{N}(\Delta {E}^{p},{d}_{\mathit{ex}})\le \frac{\Delta \kappa }{J}\]

Connaissant l’évolution de la plasticité \(\Delta {E}^{p}\) et \(\Delta \kappa\) , on peut maintenant procéder à l’intégration de la porosité. La porosité de germination \({f}^{nu}\) est simplement calculéevia les fonctions () et (); la déformation plastique cumulée est évaluée comme suit:

(2106)\[\Delta {E}_{\mathit{eq}}^{p}=\sqrt{\frac{2}{3}\Delta {E}^{p}:\Delta {E}^{p}}\]

Concernant la porosité de croissance, on applique un thêta-schéma, avec la même valeur de \(\theta\) que pendant la phase de prédiction. La discrétisation correspondante s’écrit:

(2107)\[\Delta f-\Delta {f}^{nu}=(1-\theta f-(1-\theta ){f}_{n})\text{}\mathrm{tr}(\Delta {E}_{p})\]

On en déduit:

(2108)\[\left[1+\mathrm{tr}\Delta {E}^{p}\right]{f}^{\mathit{gr}}={f}_{n}^{\mathit{gr}}+(1-\theta {f}^{nu}-(1-\theta ){f}_{n})\mathrm{tr}\Delta {E}^{p}\]

Enfin, l’endommagement de coalescence est simplement défini via:

(2109)\[{d}^{\mathit{co}}=\max\left({d}_{n}^{\mathit{co}},{q}_{1}({f}^{\text{*}}-f)\right)\]

Résolution du problème discrétisé sans viscosité

Remarque préliminaire : les équations sont résolues à une certaine précision près. Si on note \(ϵ\) la précision fournie par l’utilisateur via le mot-clé RESI_INTE_RELA, on peut introduire les précisions suivantes, respectivement en termes de contrainte, de déformation et sans unité pour l’équation ():

(2110)\[{ϵ}_{\sigma}=R(0)ϵ;{ϵ}_{\kappa}=\frac{1}{10}\frac{R(0)}{E}ϵ;{ϵ}_{G}=(1-{d}_{\mathit{ex}})ϵ\]

Si \(1-{d}_{\mathit{ex}}<{ϵ}^{1/2}\) , on considère que le point est cassé: \(T=0\) , \(\Delta \kappa =0\) , \(f={f}_{F}\) (même si \({d}_{\mathit{ex}}<{d}_{f}\) , car à la précision numérique demandée, la solution «point cassé» est effectivement acceptable).Dans le cas contraire, il faut résoudre le problème discrétisé.

Ce dernier admet une solution unique. On commence par déterminer le régime duquel elle relève: réponse élastique, réponse plastique à écoulement singulier ou réponse plastique à écoulement régulier. La solution est élastique si:

(2111)\[{T}_{\text{*}}^{e}<\widehat{T}({\kappa}_{n})\]

Dans ce cas, \(\Delta \kappa =0\) , \(\Delta {E}^{p}=0\) et \(T={T}^{e}\) .

Si la solution n’est pas élastique, on examine le cas d’un écoulement plastique singulier. Dans ce cas, \({T}_{\text{*}}=0\) , \(T=0\) , \(\Delta {E}^{p}\) est obtenu via () et on note \({\Pi}_{N}^{s}\) la valeur correspondante de \({\Pi}_{N}\) selon (). On calcule également \({\kappa}^{z}\) tel que (résolution par une méthode de Newton):

(2112)\[\widehat{T}({\kappa}^{z})=0\]

L’écoulement est singulier si \({\kappa}^{z}\ge {\kappa}_{n}\) et si \(J\phantom{\rule{0.5em}{0ex}}{\Pi}_{N}^{s}\le \Delta {\kappa}^{z}\) d’après (); on a alors \(\kappa ={\kappa}^{z}\) .

Enfin, dans le cas contraire, l’écoulement plastique est régulier et on a en particulier \(0<{T}_{\text{*}}<{T}_{\text{*}}^{e}\) . On peut remarquer qu’à \({T}_{\text{*}}\) donnée, l’équation () permet de remonter à la valeur de \({T}_{H}\) . En notant \({T}_{H}=p\phantom{\rule{0.5em}{0ex}}{T}_{H}^{e}\) , la fonction \(\widehat{g}:p\to \widehat{G}({T}_{H},{T}_{\text{*}})\) est croissante et la solution vérifie l’encadrement \(0\le p\le 1\) . Pour faciliter la résolution, on peut calculer une borne minimale \({p}_{min}>0\) en majorant la fonction \(\widehat{g}\) . En notant \(r\equiv {T}_{H}/{T}_{\text{*}}\) et \(\alpha \equiv \mu /K\) , on a:

(2113)\[\widehat{g}(p)={\left(\frac{{T}_{\mathit{eq}}^{e}}{{T}_{\text{*}}}\right)}^{2}{\left(\frac{\Lambda (r)}{\alpha \frac{{T}_{H}^{e}}{{T}_{\text{*}}}+\Lambda (r)-\alpha \phantom{\rule{0.5em}{0ex}}r}\right)}^{2}+{d}_{\mathit{ex}}\phantom{\rule{0.5em}{0ex}}2\mathrm{chm}(\frac{3}{2}{q}_{2}r)-{(1-{d}_{\mathit{ex}})}^{2}\]

On peut minorer analytiquement le dénominateur du second facteur et on le note \({D}_{min}\) . Comme en outre :math:`` , il vient:

(2114)\[\widehat{g}(p)\le {\left(\frac{{T}_{\mathit{eq}}^{e}}{{T}_{\text{*}}}\right)}^{2}{\left(\frac{\Lambda (r)}{{D}_{min}}\right)}^{2}+{d}_{\mathit{ex}}\phantom{\rule{0.5em}{0ex}}2\mathrm{sh}(\frac{3}{2}{q}_{2}r)-{(1-{d}_{\mathit{ex}})}^{2}\]

La racine du membre de droite peut être obtenue analytiquement et permet de remonter à la borne \({p}_{min}\) . Une résolution de l’équation \(\widehat{g}(p)=0\) par une méthode de Newton dans l’intervalle \({p}_{min}\le p\le 1\) permet de trouver la solution \(p({T}_{\text{*}})\) . On en contrôle la précision via la condition de convergence \(\left|\widehat{g}\right|\le {ϵ}_{G}\) .

On peut ainsi ramenerle problème incrémental à la seule inconnue \(\Delta \kappa\) . En effet, à \(\kappa\) donné, on peut calculer \({T}_{\text{*}}=\widehat{T}(\kappa )\) . De là, on en déduit \(p({T}_{\text{*}})\) par résolution de l’équation précédente. Connaissant \(p\) donc \({T}_{H}\) , on peut calculer \(\lambda =\widehat{\lambda}({T}_{H},{T}_{\text{*}})\) via (), ce qui ramène le problème à la résolution de l’équation scalaire \(\Delta \kappa =\lambda\) .

Grâce aux propriétés de convexité de la contrainte équivalente \({T}_{\text{*}}\) par rapport à \(T\) , l’analyse du schéma de «return mapping» utilisé montre qu’à \(\lambda\) croissant de 0 à \(J\phantom{\rule{0.5em}{0ex}}{\Pi}_{N}^{s}\) , \({T}_{\text{*}}\) est décroissante de \({T}_{\text{*}}^{e}\) à 0. La résolution du problème revient donc à trouver l’intersection de la courbe croissante \(\Delta \kappa \to {T}_{\text{*}}\) avec la courbe décroissante \({T}_{\text{*}}\to \lambda\) . Une borne minimale sur \(\Delta \kappa\) est fournie par le maximum de0 et de \(\Delta {\kappa}_{z}\) , car \({T}_{\text{*}}\) solution doit êtrepositive et \(\Delta \kappa\) aussi. Une borne maximale est fournie par le minimumde \(J\phantom{\rule{0.5em}{0ex}}{\Pi}_{N}^{s}\) et de \(\Delta {\kappa}^{e}\) tel que \(\widehat{T}({\kappa}^{e})={T}_{\text{*}}^{e}\) (s’il existe), car \({T}_{\text{*}}\) solution doit êtrepositive et inférieure à \({T}_{\text{*}}^{e}\) .

L’équation est résolue par une méthode de Newton entre ces deux bornes. On contrôle la précision de la solution en considérant que l’algorithme a convergé dès que l’une des conditions \(\left|\Delta \kappa -\lambda \right|\le {ϵ}_{\kappa}\) ou \(\left|\widehat{T}(\Delta \kappa )-\widehat{T}(\lambda )\right|\le {ϵ}_{\sigma}\) est atteinte. On choisit d’initialiser la méthode de Newton par la solution \(\Delta {\kappa}_{\mathit{VM}}\) du problème équivalent de von Mises (i.e. en fixant \({T}_{H}=0\) ) lorsqu’elle respecte les bornes ou sinon par une méthode de corde vis-à-vis des valeurs aux bornes.

On remonte enfin à la contrainte via la relation suivante:

(2115)\[T={T}_{H}\mathrm{Id}+{T}_{\mathit{eq}}\frac{\mathrm{dev}\phantom{\rule{0.5em}{0ex}}{T}^{e}}{{T}_{\mathit{eq}}^{e}}\]

Résolution du problème discrétisé dans le cas viscoplastique

En présence de viscosité, la condition de cohérence () est remplacée par la loi d’évolution de Norton (). Une fois discrétisée, cette dernière s’écrit:

(2116)\[\frac{\Delta \kappa }{\Delta t}={\left(\frac{⟨F(T,\kappa )⟩}{\stackrel{~}{K}}\right)}^{n}\]

Et en inversant cette relation:

(2117)\[{C}_{\Delta t}{\Delta \kappa }^{q}=⟨F(T,\kappa )⟩;{C}_{\Delta t}=\frac{\stackrel{~}{K}}{{\Delta t}^{1/n}};q=\frac{1}{n}\]

Comme \(n>1\) , \(0<q<1\) .

Elle peut encore s’exprimer de manière équivalente sous la forme d’une relation de Kuhn et Tucker:

(2118)\[\Delta \kappa \ge 0;{F}_{v}(T,\kappa )\le 0;{F}_{v}(T,\kappa )\Delta \kappa =0\]

(2119)\[\text{où}{F}_{v}(T,\kappa )=F(T,\kappa )-{C}_{\Delta t}{\Delta \kappa }^{q}=\frac{{T}_{\text{*}}}{J}-{R}_{v}(\kappa );{R}_{v}(\kappa )=R(\kappa )+{C}_{\Delta t}{(\kappa -{\kappa}_{n})}^{q}\]

Après discrétisation, la prise en compte de la viscosité se réduit à l’ajout d’un terme (positif et croissant) dans la fonction d’écrouissage.

On pourrait être tenté d’appliquer la même méthode de résolution que précédemment mais, malheureusement, ce nouveau terme n’est pas dérivable en \(\Delta \kappa =0\) (pente infinie). C’est pourquoi on propose de distinguer deux régimes, l’un à forte pente, l’autre à pente faible. La pente «frontière» est fournie forfaitairement par \(R'(0)\) .Plus précisément, on note \(V(\Delta \kappa )={C}_{\Delta t}{\Delta \kappa }^{q}\) le terme visqueux. La valeur \(\Delta {\kappa}_{v}\) frontière correspond à:

(2120)\[V'(\Delta {\kappa}_{v})=R'(0)\]

Si \(\Delta \kappa <\Delta {\kappa}_{v}\) , on se situe dans le régime à pente forte, i.e. dominé par la viscosité. A l’inverse, si \(\Delta \kappa >\Delta {\kappa}_{v}\) , on est dans le régime à pente faible, dominé par la plasticité. Compte tenue de la monotonie de l’équation scalaire \(\Delta \kappa -\widehat{\lambda}=0\) à résoudre par rapport à \(\Delta \kappa\) , on peut déterminer dès le départ le régime de la solution et réduire le cas échéant l’intervalle de recherche.

Dans le régime dominé par la plasticité, on résout l’équation en \(\Delta \kappa\) comme au paragraphe précédent; le terme visqueux ne présente pas de difficultés numériques particulières.

Au contraire, dans le régime dominé par la viscosité, une méthode de Newton portant sur \(\Delta \kappa\) n’est pas efficace (pentes trop fortes donc convergence très lente). On procède alors au changement de variable \(v=V(\Delta \kappa )\) (qui élimine les fortes pentes au voisinage de 0) et on applique la méthode de Newton par rapport à \(v\) . Les critères de convergence restent inchangés.

Matrice tangente

Dans le cas local, il s’agit de déterminer le tenseur d’ordre 4 \({H}_{\mathit{TE}}=\delta T/\delta E\) . Dans le cas non local, il faut également construire les tenseurs d’ordre 2 \({H}_{T\varphi }=\delta T/\delta \varphi\) et \({H}_{\kappa E}=\delta \kappa /\delta E\) ainsi que le scalaire \({H}_{\kappa \varphi }=\delta \kappa /\delta \varphi\) .

L’opérateur tangent dépend du régime d’écoulement. Dans le cas élastique, il se réduit à:

(2121)\[\delta T=K\mathrm{tr}(\delta E)\mathit{Id}+2\mu \mathit{dev}(\delta E);\delta \kappa =0\]

Dans le cas du régime singulier, l’expression est là aussi très simple:

(2122)\[\delta T=0;\delta \kappa =-{\left(\frac{\partial \widehat{T}}{\partial \kappa }\right)}^{-1}\left[\frac{\partial \widehat{T}}{\partial \varphi }\delta \varphi +\frac{\partial \widehat{T}}{\partial J}\frac{\partial J}{\partial E}\delta E\right]\]

où la fonction \(\widehat{T}\) contient le cas échéant les termes liés au non local et à la viscosité.

Le régime d’écoulement régulier conduit à un opérateur tangent plus complexe. On se contente ici de décrire les grandes lignes de la méthode pour le calculer. Ce régime est caractérisé par la résolution du système d’équation \(\widehat{G}({T}_{H},{T}_{\text{*}};{T}_{H}^{e},{T}_{\mathit{eq}}^{e},J)=0\) et \(\Delta \kappa -\widehat{\lambda}({T}_{H},{T}_{\text{*}};{T}_{H}^{e},{T}_{\mathit{eq}}^{e},J)=0\) avec la relation \({T}_{\text{*}}=\widehat{T}(\Delta \kappa ;J,\varphi )\) . La dérivée des fonctions implicites permet d’exprimer les variations des inconnues \(\delta {T}_{H}\) et \(\delta \kappa\) en fonction des variations des données d’entrée du système (système linéaire 2x2 à résoudre). De là, on peut remonter à la variation de la contrainte en différentiant ().

Enfin, il reste à lier les variations des données d’entrée à celles de la déformation \(E\) et de \(\varphi\) . La seule difficulté (minime) vient de la contrainte de von Mises. Pour cela, on rappelle la relation suivante:

(2123)\[ \begin{align}\begin{aligned}\delta {T}_{\mathit{eq}}^{e}=\frac{3}{2}\frac{\mathrm{dev}{T}^{e}}{{T}_{\mathit{eq}}^{e}}\\ \delta {T}_{\mathit{eq}}^{e}=\frac{3}{2}\frac{\mathrm{dev}{T}^{e}}{{T}_{\mathit{eq}}^{e}}\end{aligned}\end{align} \]

Variables calculées en post-traitement

Quatre des variables internes calculées en post-traitement nécessitent encore d’être définies à ce stade:

• La part des contraintes liée à l’écrouissage:

(2124)\[\text{SIEQ\_ECR =}J(1-{d}_{\mathit{ex}})R(\kappa )\]

• La part des contraintes liée à la viscosité:

(2125)\[\text{SIEQ\_VSC =}J(1-{d}_{\mathit{ex}}){C}_{\Delta t}\Delta {\kappa}^{q}\]

• La part des contraintes liée au non local:

(2126)\[\text{SIEQ\_NLC =}J(1-{d}_{\mathit{ex}})(r\kappa -\varphi )\]

• L’erreur d’extrapolation exprimée en contraintes:

(2127)\[\text{SIEQ\_ERX =}\left(1-\frac{{T}^{\text{*}}(T,{d}_{\mathit{ex}})}{{T}^{\text{*}}(T,d)}\right)\Vert T\Vert\]

Bibliographie

• Gurson A. L. (1977) Continuum theory of ductile rupture by void nucleation and growth: part I – Yield criteria and flow rules for porous ductile media. J. eng. Mater. Technol. 99, 2-15.

• Lorentz E., Andrieux S. (1999) A variational formulation for nonlocal damage models. Int. J. Plas. 15, 119-138.

• Tvergaard V., Needleman A. (1984) Analysis of the cup-cone fracture in a round tensile bar. Acta Metall. 32, 157-169.

• Zhang Y., Lorentz E., Besson J. (2018) Ductile damage modelling with locking-free regularised GTN model. Int. J. Numer. Methods Eng. 113, 1871-1903.