r5.03.28 Loi d’endommagement à gradient ENDO_FISS_EXP

Résumé:

Ce document décrit le modèle de comportement élastique fragile ENDO_FISS_EXP, disponible en modélisation non locale à gradient d’endommagement GRAD_VARI. L’endommagement y est modélisé de manière isotrope, ce qui implique préférentiellement de modéliser les fissures de manière individuelle et non pas par homogénéisation: l’échelle de discrétisation doit être adaptée à ce choix. Sur le plan phénoménologique, le critère d’endommagement est adapté à la modélisation des bétons, en distinguant notamment les états de traction et de compression. En outre, une restauration de rigidité en compression est introduite pour refléter la fermeture des fissures. En revanche, le modèle n’est a priori pas adapté à modéliser l’endommagement du béton en compression. Enfin, il est à noter que la réponse du modèle tend vers celle d’une loi cohésive lorsque la longueur caractéristique devient petite devant la taille de la structure. Il en résulte que la question de l’identification de cette longueur caractéristique ne se pose pas pour peu qu’elle soit suffisamment petite. Et le modèle peut donner accès aux grandeurs d’intérêt des lois cohésives, comme par exemple l’ouverture de fissure.

Modèle continu

L’ensemble des équations du modèle est réuni dans le tableau . Le lien entre les paramètres internes du modèle et les grandeurs plus accessibles à l’ingénieur y est également précisé. Les choix de modélisation appellent quelques remarques qu’on expose ci-dessous.

Tableau 2-1: Équations du modèle continu

Relation contraintes – déformations

La fonction \(A(a)\) qui intervient dans la relation contrainte – déformation mesure l’impact du niveau d’endommagement sur la rigidité résiduelle (en traction). Sa forme particulière est guidée par la cohérence avec une loi cohésive. Elle joue le même rôle que la forme plus classique en (1- a ) utilisée dans nombre de modèles de la littérature.

Par ailleurs, pour tenir compte de la restauration de rigidité en compression, la loi d’élasticité n’est plus linéaire, même si elle reste hyperélastique ( i.e. elle dérive d’une énergie) afin d’éviter toute production ou consommation parasites d’énergie. Elle est semblable à la formulation adoptée pour la loi ENDO_ISOT_BETON [R7.01.04]; elle tient compte non seulement du signe des déformations propres mais aussi de celui de la trace des déformations. Sur le plan pratique, elle nécessite donc de déterminer les valeurs propres du tenseur de déformation.

Enfin, pour améliorer la robustesse du modèle et la convergence du calcul, le saut de rigidité au passage traction – compression a été régularisé, au moyen de la fonction \(S(x)\) infiniment dérivable. Elle assure une transition progressive entre une valeur nulle et une fonction quadratique; la rapidité (ou brutalité) de la transition y est contrôlée par le paramètre \(\gamma\) . Une valeur par défaut est proposée dans la commande DEFI_MATER_GC qui stipule que 90% de la rigidité sans régularisation est retrouvée pour une déformation (en compression) de l’ordre de \(\mathit{fc}/E\) .

Evolution de l’endommagement

La présence d’un terme en Laplacien d’endommagement dans la fonction seuil ainsi qu’une condition aux limites sur le champ d’endommagement de type gradient normal nul au bord reflète le caractère non local du modèle. Plus précisément, cela correspond à la formulation non locale à gradient de variable interne telle qu’introduite dans le fascicule [R5.04.01]. Sur le plan numérique, elle autorise un traitement (presque) classique de la loi de comportement. Toutefois, le choix d’une surface d’endommagement phénoménologiquement acceptable conduit à renoncer à l’existence d’une unique fonctionnelle énergie pour décrire l’ensemble des équations du problème. Ici, l’évolution de l’endommagement et l’équilibre de la structure résultent de la minimisation de deux fonctionnelles énergies distinctes. Sur le plan pratique, l’opérateur tangent n’est alors plus symétrique.

Intégration numérique

Équations discrétisées

La discrétisation en temps des équations de comportement s’appuie sur un schéma d’Euler implicite, c’est-à-dire que les différentes variables du problème sont exprimées à l’instant final du pas de temps. L’expression de la contrainte en fin de pas de temps obéit ainsi à la même relation que celle dans le tableau, la déformation et l’endommagement y étant également exprimés à la fin du pas de temps. Reste la partie délicate qui consiste à calculer l’évolution de l’endommagement.

Intégrer l’équation d’évolution de l’endommagement est un problème certes scalaire mais non linéaire. Tout d’abord, en suivant la démarche de résolution proposée dans [Lorentz and Godard, 2011] et reprise dans le fascicule [R5.04.01], un terme affine en endommagement vient se substituer au terme en Laplacien du champ d’endommagement dans la fonction seuil \(g\) :

(2053)\[\stackrel{̃}{g}(a)=-A'(a)\Gamma (\epsilon )+\lambda +r(\alpha -a)-k\]

A l’étape d’intégration de la loi de comportement, les grandeurs \(k\) et \(\alpha\) sont considérées connues (ainsi que le coefficient de pénalisation \(r\) positif), à l’instar des déformations à la fin du pas de temps car elles s’expriment en fonction des inconnues nodales du problème. De même, la valeur de la fonction \(\Gamma (\epsilon )\) est également connue, même si son calcul nécessite lui aussi la résolution d’une équation non linéaire sur laquelle on reviendra plus loin. C’est pourquoi la fonction seuil algorithmique ci-dessus ne dépend que de \(a\) . Au final, le problème à résoudre consiste à trouver \(a\le 1\) tel que:

(2054)\[\Delta a\ge 0\phantom{\rule{3em}{0ex}};\phantom{\rule{3em}{0ex}}\stackrel{~}{g}(a)\le 0\phantom{\rule{3em}{0ex}};\phantom{\rule{3em}{0ex}}\Delta a\stackrel{~}{g}(a)=0\]

où \(\Delta a=a–{a}_{n}\) désigne l’incrément d’endommagement pendant le pas de temps. Et s’il n’existe pas de solution (inférieure à \(1\) ), alors \(a=1\) (endommagement saturé).

Résolution de l’équation d’évolution algorithmique

Comme \(A\) est convexe, \(\Gamma (\epsilon )\) positif et \(r\) positif, la fonction seuil algorithmique est décroissante. Il y a donc existence et unicité de la solution (important pour la robustesse de l’intégration numérique). En pratique, on résout le problème en trois étapes:

• Si \(\stackrel{~}{g}({a}_{n})\le 0\) alors \(a={a}_{n}\)

• Si \(\stackrel{~}{g}(1)\ge 0\) alors \(a=1\)

• Sinon résoudre \(\stackrel{~}{g}(a)=0\) avec \({a}_{n}\le a\le 1\) par une méthode de Newton à bornes contrôlées

Calcul effectif de la fonction gamma

D’après le tableau , le calcul de \(\Gamma (\epsilon )\) renvoie à celui de \(\chi (\epsilon )\) qui n’est défini que de manière implicite: il faut donc résoudre une équation non linéaire dont on sait qu’elle admet une solution unique, cf. [Lorentz, 2016]. Pour cela, on commence par introduire quelques notations et un changement de variable dans lequel \(x\) est la nouvelle inconnue:

(2055)\[{\sigma}^{e}=\lambda \mathit{tr}(\epsilon )\mathit{Id}+2\mu \epsilon\]

L’espace du tenseur \(s\) est généré par:

(2056)\[\mathit{span}(s)=\lbrace {s}_{1},{s}_{2},{s}_{3}\rbrace\]

Avec:

(2057)\[{s}_{\beta}=\parallel s+({\beta}_{0}-\frac{1}{3})\mathit{tr}(s)\mathit{Id}\parallel\]

L’équation à résoudre pour déterminer \(x\) et donc \({\rm X}\) s’écrit ainsi:

(2058)\[{s}_{\beta}x+\sqrt{\sum_{i=1}^{3}\exp(2{s}_{i}x)}-{\gamma}_{0}=0\]

Le premier membre est convexe. On estime un majorant \({x}_{0}\) de la solution puis on résout l’équation par une méthode de Newton à partir de \({x}_{0}\) , ce qui conduit à une suite de valeurs décroissantes et qui convergent vers la solution de l’équation du fait de la convexité annoncée. A noter que le majorant est obtenu par l’une deux relations suivantes, où \({s}_{M}\) désigne la plus grande des trois valeurs propres de \(s\) :

(2059)\[{s}_{\beta}x+\sqrt{\sum_{i=1}^{3}\exp(2{s}_{i}x)}\ge {s}_{\beta}x\Rightarrow x\le \frac{{\gamma}_{0}}{{s}_{\beta}}\]

Ou:

(2060)\[{s}_{\beta}x+\sqrt{\sum_{i=1}^{3}\exp(2{s}_{i}x)}\ge \exp({s}_{M}x)\phantom{\rule{2em}{0ex}}\Rightarrow \phantom{\rule{2em}{0ex}}x\le \frac{\ln{\gamma}_{0}}{{s}_{M}}\]

On prendra le plus petit des deux majorants par souci de performance (ou le seul si \({s}_{M}\le 0\) ).

Matrice tangente

Déterminer la matrice tangente cohérente au schéma d’intégration choisi ne soulève pas de difficulté particulière (dérivation de fonctions composées et dérivation de fonctions implicites) à ceci prèsque les valeurs propres de la déformation interviennent dans le calcul de \(\Gamma (\epsilon )\) et dans celui de la loi élastique au moment de distinguer la traction de la compression. Nous examinons ici chacun de ces points plus précisément, dans l’optique de s’affranchir de la dérivée des vecteurs propres associés qui, on le rappelle, ne sont pas définis lorsque les valeurs propres ne sont pas distinctes.

Pour la dérivée de \(\Gamma (\epsilon )\) , seul le terme qui à \(x\) associe la norme de \(\exp(xs)\) nécessite un traitement particulier. Or:

(2061)\[\delta \parallel \exp(xs)\parallel =\frac{\delta \exp(xs):\exp(xs)}{\parallel \exp(xs)\parallel }=\parallel \exp(xs)\parallel \delta x+\frac{\exp(xs)}{\parallel \exp(xs)\parallel }:\exp'(xs):(x\delta s)\]

A ce stade, même si la dérivée de l’exponentielle (tenseur d’ordre quatre) apparaît formellement dans l’expression, il n’est pas utile de procéder effectivement à la dérivée. En effet, la dérivée de l’exponentielle est symétrique car c’est la dérivée de la fonction scalaire isotrope suivante (cf. [Silhavy, 2000] pour les propriétés de dérivabilité):

(2062)\[f(S)=\sum_{i=1}^{3}{e}^{{A}_{i}}\]

Et donc:

(2063)\[\exp(A):\exp'(A)=\exp'(A):\exp(A)=\exp(2A)\]

Ce dernier résultat se montre aisément en remarquant que \(A\) et \(\exp(A)\) commutent. Au final, le calcul de la matrice tangente passera par celui de l’exponentielle tensorielle, d’après la relation ci-dessus, sans avoir besoin de dériver les vecteurs propres. Et pour le calcul de l’exponentielle tensorielle, on s’appuiera sur la méthode dite de «scaling and squaring», voir par exemple [Higham, 2005], qui s’affranchit du calcul des vecteurs propres.

Pour ce qui concerne la relation contraintes – déformations, elle exhibe un terme de trace qui ne pose pas de problème de dérivation particulier et une relation fonctionnelle infiniment dérivable entre la contrainte propre \({\sigma}_{i}\) et la déformation propre \({\epsilon}_{i}\) . Pour celle-ci, on entre dans le cadre de la dérivée de fonctions tensorielles basées sur une fonction scalaire des valeurs propres: l’article [Carlson and Hoger, 1986] fournit alors les expressions de la dérivée (tensorielle) que les valeurs propres soient distinctes ou non.

Description des versions du document

Version Aster	Auteur(s) Organisme(s)	Description des modifications
13.2	E.LORENTZ, EDF-R&D/AMA	Texte initial

Bibliographie

• Carlson D. E., Hoger A. (1986) The derivative of a tensor-valued function of a tensor. Quart. Appl. Math. 44, 409-423.

• Higham N. J. (2005) The scaling and squaring method for the matrix exponential revisited. SIAM J. Matrix Anal. Appl. 26, 1179-1193.

• Lorentz E. (2016) A nonlocal damage model for plain concrete consistent with cohesive fracture. Submitted.

• Lorentz E., Godard V. (2011). Gradient damage models: toward full-scale computations. Comput. Methods Appl. Mech. Engrg. 200, 1927-1944.

• Silhavy M. (2000) Differentiability properties of isotropic functions. Duke Mathematical Journal 104, 367-373.