r5.03.25 Loi d’endommagement régularisée ENDO_SCALAIRE

Résumé:

Ce document décrit le modèle de comportement élastique fragile ENDO_SCALAIRE disponible seulement pour la modélisation non locale à gradient d’endommagement GRAD_VARI. L’endommagement est modélisé de manière scalaire; les chargements en compression et en traction ne sont pas distingués. A la différence des autres lois d’endommagement introduites précédemment, cette dernière se comporte d’une manière régulière (pas de snap-back, allongement fini à la rupture) au moins dans les cas unidimensionnels.

Formulation variationnelle du problème d’endommagement

Cas d’une loi générique

Deux approches équivalentes peuvent être utilisées pour décrire le processus d’endommagement d’un matériau isotrope fragile. D’un côté il est possible de dériver la loi d’endommagement dans le cadre de description standard généralisée. Dans ce cas il est nécessaire de définir une énergie libre du système, ainsi que le potentiel de dissipation. La règle d’écoulement établi alors l’évolution des variables internes.

Comme pour la description d’endommagement on a seulement besoin d’une variable scalaire, la description précédente se simplifie et peut être ramenée vers un problème variationnel sous contrainte d’accroissement d’endommagement [bib2].

Pour définir une loi de comportement à gradient d’endommagement [R5.04.01] il suffit donc d’exprimer la densité d’énergie libre totale (élastique+dissipation) en fonction de tenseur de déformation \(\varepsilon\) et de variable d’endommagement \(0\le a\le 1\) . La répartition spatiale de l’endommagement est donnée alors par un champ \(a(x)\) . La densité d’énergie libre se présente en générale sous la forme suivante :

\(\Phi (\varepsilon ,a)=A(a)w(\varepsilon )+\omega (a)+c/2{(\nabla a)}^{2}\) éq 2.1-1

Ici \(c\) est le paramètre de non localité (C_GRAD_VARI) \(w(\varepsilon )\) l’énergie de déformation élastique, \(\omega (a)\) l’énergie de dissipation et \(A(a)\) la fonction de rigidité. \(a=0\) correspond au matériau sain et \(a=1\) correspond au matériau complètement endommagé : \(A(1)=0,A(0)=1\) . Le problème d’évolution est désormais une simple problème de minimisation d’énergie libre de Helmholtz \(F\equiv \int\Phi (\varepsilon ,a)d\Omega\) sous contrainte \(\dot{a}\ge 0\) [1] .

\({\min}_{(\varepsilon ,a)}F(\varepsilon ,a),\text{où}F(\varepsilon ,a)=\int\left[A(a)\varepsilon :E:\varepsilon +\omega (a)+c/2{(\nabla a)}^{2}\right]d\Omega\)

où on a remplacé \(w(\varepsilon )=\varepsilon :E:\varepsilon /2\) en utilisant la définition du tenseur de Hooke. Deux équations dérivent du problème variationnel de minimisation : \(\delta F(\varepsilon ,a)/\delta \varepsilon =0\) [2] et \(\delta F(\varepsilon ,a)/\delta a\ge 0\) . L’inégalité dans la deuxième équation est liée à la présence de contrainte imposée. Ces deux équations doivent être satisfaites partout dans le domaine d’intégration \(\Omega\) . Elles sont complétées par une équation de cohérence de Kuhn-Tucker \(\dot{a}\cdot \delta F(\varepsilon ,a)/\delta a=0\) . Sur les bords \(\partial \Omega\) nous obtenons une condition de normalité supplémentaire \(\nabla a\cdot n=0\) , où \(n\) est un vecteur-normal. Enfin la variable d’endommagement et son gradient doivent être continus à intérieur du domaine d’intégration pour réaliser le minimum de la fonctionnelle en question (voir [bib2,4] pour plus de détails).

Relations de comportement

Le lien entre la formulation variationnelle et les lois d’évolution habituelles est direct. L’état du matériau est caractérisé par la déformation \(\varepsilon\) et l’endommagement \(a\) , compris entre 0 et 1. On définit la relation contrainte-déformation, qui reste élastique, et la rigidité est affectée par l’endommagement:

\(\sigma =\delta F(\varepsilon ,a)/\delta \varepsilon =A(a)E:\varepsilon\) éq 2.2-1

avec \(E\) le tenseur de Hooke. L’évolution de l’endommagement, toujours croissante, est gouvernée par la fonction seuil suivante:

\(f(\varepsilon ,a)=-\delta \Phi (\varepsilon ,a)/\delta a=-\frac{1}{2}A'(a)\varepsilon :E:\varepsilon -\omega '(a)+c\Delta a\) éq 2.2-2

La condition de cohérence prend alors sa forme habituelle:

\(f(\varepsilon ,a)\le 0\dot{a}\ge 0\dot{a}f(\varepsilon ,a)=0\) éq 2.2-3

On note deux particularités de cette formulation. Premièrement, la fonction seuil est non-locale à cause de la présence du laplacien d’endommagement. Ensuite, l’absence de condition d’écoulement se justifie par le double rôle de l’endommagement \(a\) , d’un côté elle se présente comme une variable d’évolution interne, de l’autre côté elle remplit la mission du paramètre de Lagrange \(\lambda \equiv a\) .

On voit aussi l’avantage de présentation des lois d’endommagement sous leur forme variationnelle. Il suffit de décrire la densité d’énergie libre totale (éq.2.1-1), qui inclut la dissipation, pour définir complètement la loi d’évolution.

Identification des paramètres pour la loi ENDO_SCALAIRE

Dans la loi ENDO_SCALAIRE les fonctions de rigidité et dissipation sont choisies comme suit :

\(\omega (a)=\mathrm{ka},A(a)={(\frac{1-a}{1+\gamma a})}^{2}\)

Les paramètres de cette loi de comportement sont alors au nombre de cinq. D’une part, le module de Young \(E\) et le coefficient de Poisson \(\nu\) qui déterminent le tenseur de Hooke par:

\({E}^{-1}\cdot \sigma =\frac{1+\nu }{E}\sigma -\frac{\nu}{E}(\text{tr}\sigma )\text{Id}\) éq 2.2-1

D’autre part, \(k,\gamma ,c\) qui définissent le comportement adoucissant, ainsi que la largeur caractéristique de la bande d’endommagement. Ces derniers peuvent être recalés aux paramètres macroscopiques à partir du modèle unidimensionnel, qui admet une solution semi-analytique (bib 6,7). En notant la contrainte au pic par \({\sigma}_{y}\) , l’énergie de la rupture de Griffith par \({G}_{f}\) et la taille de la zone endommagée à la rupture par \(D\) on obtient :

\(k=\frac{3{G}_{f}}{4D},c=\frac{3}{8}D{G}_{f},\gamma =\frac{3E{G}_{f}}{4{\sigma}_{y}^{2}D}-1\)

Les essais numériques ont montré que pour éviter la présence de snap-back dans la réponse force-déplacement en 1D, il faudrait avoir \(\gamma \ge 2.8\) . Pour cette raison le choix a été fait de simplifier l’entrée des données du modèle, on renseigne non pas le jeu complet de paramètres macroscopiques \({\sigma}_{y},{G}_{f},D\) , mais directement les paramètres du modèle \(\gamma ,c\) et la contrainte au pic \({\sigma}_{y}\) , donnés sous les mots-clé facteurs ENDO_SCALAIRE (GAMMA, SY) et NON_LOCAL (C_GRAD_VARI) de l’opérateur DEFI_MATERIAU. Quant à \(E\) et \(\nu\) , ils sont donnés classiquement sous le mot-clé facteur ELAS ou ELAS_FO. Le raisonnement qui suit, n’est valable au sens strict que pour la modélisation 1D, mais peut se trouver utile pour les utilisateurs non avertis. Si les paramètres \(E,\nu ,{G}_{f},{\sigma}_{y}\) sont définis a priori, l’utilisateur peut faire varier le paramètre \(D\) afin de satisfaire la condition d’absence de snap-backs locaux \(\gamma \ge 2.8\) . Il doit par la suite s’assurer que la taille du système considéré est supérieure à la largeur de la bande d’endommagement \(D\) .

Exemple du béton en traction

\(\begin{array}{}E=30\text{GPa,}\nu =0.2\\ {G}_{f}=100\text{N/m}\\ {\sigma}_{y}=3\text{MPa}\end{array}\)

\(\begin{array}{}\text{ELAS}(\text{E}=\mathrm{3e10,}\text{NU}=0.2)\\ \text{ENDO\_SCALAIRE}(\text{GAMMA}=1/(\mathrm{4D})-1,\text{SY}=3e6)\\ \text{NON\_LOCAL}(\text{C\_GRAD\_VARI}=\mathrm{37.5D})\end{array}\)

La largeur de la bande d’endommagement est à choisir en respectant \(\gamma \ge 2.8\) <=> \(D\le 66\text{mm}\)

Intégration de la loi de comportement en locale

Nous présentons ici la méthode d’intégration de la loi ENDO_SCALAIRE dans sa version locale (\(c=0\) ), afin que l’utilisateur puisse faire une généralisation pour le cas non-local, qui elle est générique et se repose entièrement sur l’algorithme présenté dans la doc. [R5.04.01].

La discrétisation temporelle des équations [éq 2.2-1] à [éq 2.1-3] sur un pas de temps

est réalisée par un schéma d’Euler implicite. Intégrer en temps la loi de comportement consiste à déterminer l’état de contrainte et d’endommagement de la solution du système non linéaire suivant:

\(\sigma =A(a)E:\varepsilon\) éq 2.4-1

\({f}_{\text{loc}}(\epsilon ,a)\le 0a-{a}^{-}\ge 0(a-{a}^{-})\cdot {f}_{\mathrm{loc}}(\epsilon ,a)=0\) éq 2.4-2

où les variables sans indices correspondent au pas de temps final \(t\) , comme par exemple la déformation \(\varepsilon\) ; l’état du matériau au début du pas de temps \(({\varepsilon}^{-},{a}^{-})\) est indiqué par l’indice «-». La fonction seuil locale est donnée par (éq. 2.2-2) : \({f}_{\text{loc}}(\varepsilon ,a)=(1+\gamma )\frac{(1-a)}{{(1+\gamma a)}^{3}}\varepsilon :E:\varepsilon -k\)

Une méthode de résolution a été proposée par [bib3]. Elle commence par examiner la solution sans évolution de l’endommagement (aussi appelée essai élastique) puis, si nécessaire, procède à une correction pour vérifier la condition de cohérence. Dans le cas présent, l’existence et l’unicité de la solution garantissent le bon fonctionnement de la méthode. Considérons l’essai élastique:

\(a={a}^{-}\) solution si \({f}^{\text{el}}(\varepsilon )\equiv {f}_{\text{loc}}(\varepsilon ,{a}^{-})\le 0\) éq 2.4-3

Dans le cas contraire, l’endommagement est obtenu en résolvant \({f}_{\text{loc}}(\varepsilon ,a)=0\) (polynôme d’ordre 3).

\((1-a)(1+\gamma )\varepsilon :E:\varepsilon /k={(1+\gamma a)}^{3}\) éq 2.4-4

C’est la plus grande racine qui est choisie parmi les trois existant.

Il reste encore à s’assurer que l’endommagement ne dépasse pas la valeur 1. En fait, lorsque \(a=1\) , la rigidité du point matériel considéré s’annule \(A(1)=0\) . Dans la mesure où aucune technique de suppression des éléments finis «cassés» n’est mise en œuvre (technique éventuellement délicate lorsque les éléments finis possèdent plusieurs points de Gauss), des pivots nuls peuvent apparaître dans la matrice de rigidité. C’est pourquoi on introduit un seuil numérique de rigidité résiduelle élastique pour la matrice tangente, qui peut être renseigné sous le mot-clé facteur COEF_RIGI_MINI de l’opérateur DEFI_MATERIAU. Cette valeur sans dimension est un coefficient multiplicateur du module élastique d’un modèle linéaire isotrope. Pour préserver un conditionnement raisonnable de la matrice de rigidité, on choisit la valeur par défaut \(\minA(a)={10}^{-5}\) .

Un indicateur \(\chi\) , rangé dans la deuxième variable interne, précise alors le comportement pendant le pas de temps courant:

• 

comportement élastique (énergie de déformation inférieure au seuil)

• 

évolution de l’endommagement

• 

endommagement saturé (\(a=1\) ).

Intégration de la loi de comportement en non locale

Nous présentons ici seulement la méthode d’intégration de la loi ENDO_SCALAIRE dans sa version locale (\(c=0\) ), car la généralisation pour le cas non-local est générique et se repose entièrement sur l’algorithme présenté dans la doc. [R5.04.01]. On note que pour la version non locale la fonction seuil est décalée, nous obtenons donc un polynôme d’ordre 4 à résoudre. Quant à la contrainte, elle est donnée par [éq 2.4-1] dans tous les cas.

Description des variables internes

Les variables internes sont au nombre de trois:

• \(\mathrm{VI}(1)\) endommagement \(a\)

• \(\mathrm{VI}(2)\) indicateur \(\chi\)

• \(\mathrm{VI}(3)\) rigidité résiduelle \(1-A(a)\)

[1]

On note par \(\nabla a`la dérivée spatiale du champ d'endommagement et par :math:\)dot{a}`celle liée à l’évolution temporelle

[2]

\(\delta F(\varepsilon ,a)/\delta \varepsilon`est la dérivée variationnelle partielle selon la direction du champ spatial :math:\)varepsilon (x)`, champ :math:`a(x)`restant fixé.

Pilotage par prédiction élastique

Le pilotage de type PRED_ELAS standard contrôle l’intensité du chargement pour satisfaire une certaine équation liée à la valeur de la fonction seuil \({f}^{\mathrm{el}}\) lors de l’essai élastique [bib5]. Par conséquent, seuls les points où l’endommagement n’est pas saturé sont pris en compte. L’algorithme qui prend en charge ce mode de pilotage, cf. [R5.03.80], requiert la résolution en chacun de ces points de Gauss de l’équation scalaire suivante dans laquelle \(\Delta \tau\) est une donnée et \(\eta\) l’inconnue:

\({f}^{\mathrm{el}}({\varepsilon}_{\text{impo}}+\eta {\varepsilon}_{\text{pilo}},{a}^{\text{\_}})=\Delta \tau\) éq 3-1

Notons que cette équation est modifiée pour le pilotage PRED_ELAS en ENDO_SCALAIRE afin d’avoir le paramètre \(\Delta \tau\) qui correspond à l’incrément d’endommagement que l’on cherche à obtenir pour au moins un point de la structure. On ne cherche alors plus un paramètre de pilotage \(\eta\) qui fasse sortir le critère d’une valeur \(\Delta \tau\) avec l’endommagement issu du pas de temps précédent (cf. Eq 3-1), mais un paramètre \(\eta\) qui nous ramène sur le critère avec un endommagement augmenté de \(\Delta \tau\) :

\({f}^{\mathrm{el}}({\varepsilon}_{\text{impo}}+\eta {\varepsilon}_{\text{pilo}},{a}^{\text{\_}})=\Delta \tau \text{}\Rightarrow \text{}{f}^{\mathrm{el}}({\varepsilon}_{\text{impo}}+\eta {\varepsilon}_{\text{pilo}},{a}^{\text{\_}}+\Delta \tau )=0\) éq 3-2

Description des versions du document

Version Aster	Auteur(s) Organisme(s)	Description des modifications
10.0	K.KAZYMYRENKO, E.LORENTZ, S.CUVILLIEZ EDF-R&D/AMA	Texte initial
10.2	K.KAZYMYRENKO, EDF-R&D/AMA	Corrections mineures des notations

,

Bibliographie

1. LEMAITRE J., CHABOCHE J.L.: Mécanique des matériaux solides. Dunod: Paris, 1988.

2. MARIGO J.J.: Formulation d’une loi d’endommagement d’un matériau élastique. Compte‑rendu de l’Académie des Sciences, Paris 1981; série II, 292(19): 1309-1312.

3. SIMO J.C., TAYLOR R.L.: Consistent tangent operators for rate-independent elastoplasticity. Computer Methods in Applied Mechanics and Engineering 1985 ; 48 : 101-118.

4. LORENTZ E., ANDRIEUX S.: A variational formulation for nonlocal damage models. International Journal of Plasticity 1999; 15: 119-138.

5. LORENTZ E., BADEL P.: A load control method for damage finite element simulations. International Journal for Numerical Methods in Engineering 2004, 60:499–526

6. LORENTZ E., GODARD V. : Gradient damage models: toward full-scale computations. Comput. Methods Appl. Mech. Eng. (2010) in press

7. LORENTZ E., CUVILLIEZ S., KAZYMYRENKO K. Convergence of a gradient damage model toward a cohesive zone model. CRAS 339 (2011) 20–26