r7.01.42 Loi d’endommagement homogène ENDO_LOCA_EXP#

Résumé:

Ce document décrit le modèle de comportement élastique fragile ENDO_LOCA_EXP, disponible en modélisation locale uniquement (par opposition à une formulation non locale type GRAD_VARI). Cela implique de restreindre l’usage du modèle à des situations dans lesquelles l’endommagement est diffus et non pas localisé. La version non locale correspondante est fournie par le modèle ENDO_FISS_EXP [R5.03.28]. Si l’endommagement y est modélisé de manière isotrope, le critère d’endommagement est quant à lui bien adapté à la modélisation des bétons, en distinguant notamment les états de traction et de compression, à la différence du modèle ENDO_ISOT_BETON [R7.01.04] par exemple. 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.

Modèle continu#

Équations de comportement#

La relation contrainte -déformation prend la forme suivante, dans laquelle \({⟨x⟩}_{-}\) et \({⟨x⟩}_{+}\) désignent respectivement les parties négatives et positives d’un scalaire ou d’un tenseur symétrique:

(4065)#\[\sigma (\epsilon ,b)\text{}=\text{}B(b)\left[\lambda {⟨\mathit{tr}(\epsilon )⟩}_{+}+2\mu {⟨\epsilon ⟩}_{+}\right]\text{}+\text{}\left[\lambda {⟨\mathit{tr}(\epsilon )⟩}_{-}+2\mu {⟨\epsilon ⟩}_{-}\right]\]

Quant à la fonction de rigidité \(B(b)\) , elle vaut:

(4066)#\[B(b)\text{}=\text{}\frac{1-b}{\left(1-b\right)+\frac{\kappa}{2}b\left({m}_{0}+({D}_{1}-{m}_{0}){b}^{r-1}\right)}\]

L’évolution de l’endommagement est gouvernée par les conditions de Kuhn et Tucker suivantes:

(4067)#\[\dot{b}⩾0;Q(\epsilon )-\varphi (b)⩽0;\dot{b}\left[Q(\epsilon )-\varphi (b)\right]=0\]

La fonction seuil \(\varphi (b)\) a pour expression:

(4068)#\[\varphi (b)={\left[(1-b)+\frac{\kappa}{2}b\left({m}_{0}+({D}_{1}-{m}_{0}){b}^{r-1}\right)\right]}^{2}\]

Quant à la fonction de la déformation \(Q(\epsilon )\) , elle s’écrit:

(4069)#\[Q(\epsilon )={\epsilon}_{\text{*}}^{2}\]

Dans cette expression, la pseudo-norme \({\epsilon}_{\text{*}}\) est la solution de l’équation scalaire suivante:

(4070)#\[{f}_{\sigma}\left(\frac{\lambda \mathit{tr}(\epsilon )\text{Id}+2\mu \epsilon }{{\epsilon}_{\text{*}}}\right)=0\]

Enfin, le seuil d’initiation a pour expression:

(4071)#\[{f}_{\sigma}(\sigma )=\Vert \frac{\sigma}{{\sigma}_{0}}+\left({\beta}_{0}-\frac{1}{3}\right)\mathrm{tr}\left(\frac{\sigma}{{\sigma}_{0}}\right)\text{Id}\Vert +\Vert \exp\left(\frac{\sigma}{{\sigma}_{0}}\right)\Vert -{\gamma}_{0}\]

La fonction de rigidité \(B(b)\) 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 ainsi que celle du seuil \(\varphi (b)\) sont guidées par la cohérence avec la réponse de la loi ENDO_FISS_EXP dans un essai de traction uniaxiale confinée. Elle joue le même rôle que la forme plus classique \((1-b)\) utilisée dans nombre de modèles de la littérature.

Énergie élastique et transition traction – compression#

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 déformations propres. Sa valeur est fournie parmiles variables internes (ENERELAS).

Ses propriétés sont étudiées de manière détaillée dans [Badel et al., 2007]. On retiendra qu’elle est continûment dérivable, ce qui assure la continuité de la relation contrainte – déformation. En revanche, elle n’est pas deux fois continûment dérivable, ce qui signifie que la rigidité présente un saut au franchissement de la transition traction – compression. Cela peut nuire à la convergence du calcul; c’est pourquoi on propose de régulariser cette transition.

L’idée, explicitée dans [Lorentz, 2017] et [Lorentz, 2020] consiste à régulariser l’énergie élastique. Cela a comme conséquence sur la relation contrainte – déformation () que les parties positive et négative y sont approchées par des fonctions infiniment dérivables:

(4072)#\[{⟨x⟩}_{-}\simeq \text{}{N}_{\gamma}(x);{⟨x⟩}_{+}\simeq \text{}x-{N}_{\gamma}(x)\]

L’approximation \({N}_{\gamma}\) dépend d’un paramètre \(\gamma >0\) qui tend vers l’infini en l’absence de régularisation. Elle a pour expression:

(4073)#\[\begin{split}{N}_{\gamma}(x)=\lbrace \begin{array}{cc}\left(x+\frac{1}{2\gamma }\right)\exp\left(-\text{}\frac{1}{\gamma x}\right)& \text{si}x<0\\ 0& \text{si}x⩾0\end{array}\end{split}\]

Paramètres du modèle#

On peut recenser les paramètres utilisés dans les équations de comportement. Il s’agit de \(\lambda ` , :math:\)mu ` , \({\sigma}_{0}\) , \({\gamma}_{0}\) , \({\beta}_{0}\) , \(\kappa ` , :math:`{m}_{0}\) , \({D}_{1}\) et \(r\) . On se propose de donner leur signification et leur expression, pour lister finalement ceux devant être renseignés par l’utilisateur.

Elasticité – Les caractéristiques élastiques, qui interviennent dans l’expression de l’énergie de déformation, se réduisent classiquement aux coefficients de Lamé \(\lambda\) et \(\mu\) ou, de manière équivalente, au module de Young \(E\) et au coefficient de Poisson \(\nu\) . On en déduit en particulier la rigidité confinée \({E}_{c}=\lambda +2\mu\) .

Seuil d’initiation de l’endommagement – Le seuil d’initiation est défini via la fonction \({f}_{\sigma}\) qui dépend de trois paramètres: \({\sigma}_{0}\) , \({\gamma}_{0}\) et \({\beta}_{0}\) . En pratique, on fixe \({\beta}_{0}\) à 0.1 et on détermine \({\sigma}_{0}\) et \({\tau}_{0}\) en fonction des résistances du béton en traction \({f}_{t}\) en en compression \({f}_{c}\) , selon la méthodologie présentée dans [Lorentz, 2017]. Ces données permettent alors de calculer \(Q\left(\epsilon \right)\) pour toute déformation. En particulier, on peut déterminer la contrainte d’initiation en traction uniaxiale confinée \({\sigma}_{c}\) par la résolution de l’équation suivante, dans laquelle \(n\) désigne la direction (arbitraire) de traction:

(4074)#\[{f}_{\sigma}\left(\left(\frac{\lambda}{{E}_{c}}\mathit{Id}+\frac{2\mu }{{E}_{c}}n\otimes n\right){\sigma}_{c}\right)=0\]

On accède ensuite à \({w}_{c}\) l’énergie de déformation lorsque le seuil d’initiation est atteint en traction uniaxiale confinée via \({w}_{c}={\sigma}_{c}^{2}/2{E}_{c}\) .

  • Energie d’endommagement – L’énergie de fissuration \({G}_{F}\) mesure l’énergie consommée par unité de surface fissurée. Elle est à rapprocher de la distance moyenne \(L\) qui sépare deux fissures pour construire l’énergie volumique \(k={G}_{F}/L\) consommée par le modèle d’endommagement; on choisit de normaliser cette dernière via \(\kappa ={G}_{F}/\left(L{w}_{c}\right)\) . Les deux paramètres \({G}_{F}\) et \(L\) n’interviennent donc dans le modèle d’endommagement qu’à travers leur ratio même si tous deux ont une signification physique propre.

  • Facteur de forme – Un dernier paramètre, \(p⩾1\) , influe sur la forme de la réponse en phase d’endommagement (réponse d’autant plus courbée dans un diagramme \(\sigma -\epsilon\) que \(p\) est grand). De la valeur de \(p\) , on déduit les paramètres \({m}_{0}\) , \({D}_{1}\) et \(r\) selon les expressions suivantes:

(4075)#\[{m}_{0}=\frac{3\pi }{2}{\left(p+2\right)}^{-\text{}\frac{3}{2}}\]
(4076)#\[{D}_{1}=\frac{3\pi }{4}\sqrt{p+1}\]
(4077)#\[r=\frac{2({D}_{1}-1)-{m}_{0}}{2-{m}_{0}}\]

Et de là, connaissant également l’énergie volumique à rupture normalisée \(\kappa\) , on peut exprimer pour toute valeur d’endommagement la fonction de rigidité \(B(b)\) et le seuil \(\varphi \left(b\right)\) . La donnée de \(p\) limite également les valeurs possibles pour \(\kappa ` pour des questions de stabilité du point matériel ce qui, indirectement, fixe une borne supérieure :math:`{L}_{c}\) à la distance inter-fissures \(L\) via l’inégalité suivante:

(4078)#\[L\text{}⩽\text{}{L}_{c}=\frac{{G}_{F}}{{w}_{c}}\frac{3\pi }{4}{(p+2)}^{-\text{}\frac{3}{2}}\]

  • Régularisation de la transition traction – compression – Le dernier paramètre, \(\gamma ` , conditionne la régularisation des fonctions parties positive et négative intervenant dans la relation contrainte – déformation, comme expliqué au § :ref:`2.2 <RefHeading___Toc 2097_2561604576>\) . Alternativement à la donnée d’une valeur pour \(\gamma ` , on peut préférer fixer la proportion de la rigidité (sans régularisation) qui est retrouvée pour une déformation (en compression) de l'ordre de :math:`{f}_{c}/E\) , proportion qu’on note \({\rho}_{\gamma}\) .

Au final, en s’appuyant sur la commande DEFI_MATERIAU, on est amené à renseigner les neuf paramètres suivants: \(E\) , \(\nu ` , :math:`{\sigma}_{\mathrm{c}}\) , \({\sigma}_{0}\) , \({\gamma}_{0}\) , \({\beta}_{0}\) , \(\kappa ` , :math:`p\) et \(\gamma ` . Des valeurs par défaut sont proposées: :math:`{\beta}_{0}=0.1\) et \(\gamma =0\) (pas de régularisation). Alternativement, la commande DEFI_MATER_GC permet de renseigner ces paramètres de manière plus commune. On y fournit en effet les dix paramètres suivants: \(E\) , \(\nu ` , :math:`{f}_{t}\) , \({f}_{c}\) , \({\gamma}_{0}\) , \({\beta}_{0}\) , \({G}_{F}\) , \(L\) , \(p\) et \({\rho}_{\gamma}\) . Là aussi, deux valeurs par défaut sont proposées: \({\beta}_{0}=0.1\) et \({\rho}_{\gamma}=0.95\) (une régularisation jugée raisonnable).

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 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 scalaire non linéaire. Il résulte de la discrétisation temporelle des conditions de Kuhn et Tucker (), en indiçant par \(n\) les quantités au début du pas de temps et en omettant tout indice pour celle en fin de pas de temps:

(4079)#\[b-{b}_{n}⩾0;Q(\epsilon )-\varphi (b)⩽0;(b-{b}_{n})\left[Q(\epsilon )-\varphi (b)\right]=0\]

Il s’agit alors de déterminer la valeur de \(b\) , :math:`epsilon ` étant connu; c’estce qu’on appelle l’intégration de la loi de comportement.

Il faut aussi bien noter que la déformation :math:`epsilon ` correspond à la partie mécanique des déformations, c’est-à-dire la déformation totale de laquelle on a retiré les déformations de retrait, notamment thermique (dilatation). Dans la terminologie propre au code Aster (Code_Aster), il s’agit d’une loi exprimée en déformation mécanique.

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

Pour traiter l’alternative (), on reprend la méthode décrite dans [Lorentz, 2020] où on a également démontré l’existence et l’unicité de la solution \(b\) . On commence par évaluer la fonction \(Q(\epsilon )\) ; cela passe par la résolution de l’équation scalaire (). On renvoie le lecteur au document [R5.03.28] pour la méthode de la résolution, la fonction \(Q\) étant égale à un facteur près à la fonction :math:`Gamma ` définie dans ce fascicule.

On procède ensuite à un tir élastique en supposant que \(b={b}_{n}\) . Si \(Q(\epsilon )-\varphi ({b}_{n})⩽0\) , alors il s’agit de la solution; il s’agit du régime élastique. Dans le cas contraire, on examine la valeur ultime \(b=1\) . Si \(Q(\epsilon )-\varphi (1)⩾0\) ,alors il s’agit de la solution; il s’agit du régime saturé. Enfin, dernier cas de figure, le régime endommageant, la solution \(b\) est la racine de l’équation suivante, comprise entre \({b}_{n}\) et 1:

(4080)#\[\varphi (b)-Q(\epsilon )=0\]

Le membre de gauche est croissant avec \(b\) et de signe opposé aux deux bornes. La résolution est réalisée par une méthode de Newton à bornes contrôlées.

La précision avec laquelle les deux équations scalaires sont résolues, celle pour évaluer \(Q(\epsilon )\) et () est indirectement fixée par l’utilisateur, de sorte que l’erreur commise sur la contrainte soit inférieure à la limite en traction \({\sigma}_{c}\) multipliée par le paramètre RESI_INTE_RELA de STAT_NON_LINE. On pourra se reporter à [Lorentz,2020] pour le calcul des précisions attendues pour ces deux équations scalaires.

Matrice tangente#

La matrice tangente est obtenue par variation de la contrainte dans () par rapport à la déformation. Pour cela, on commence par noter \({E}_{+}(\epsilon )\) et \({E}_{-}(\epsilon )\) les parties positive et négative de l’opérateur d’élasticité:

(4081)#\[{E}_{+}(\epsilon )=\lambda {⟨\mathit{tr}(\epsilon )⟩}_{+}+2\mu {⟨\epsilon ⟩}_{+};{E}_{-}(\epsilon )=\lambda {⟨\mathit{tr}(\epsilon )⟩}_{-}+2\mu {⟨\epsilon ⟩}_{-}\]

En différentiant (), il vient alors:

(4082)#\[\delta \sigma =\left[B(b)\frac{{\mathit{dE}}_{+}}{d\epsilon }(\epsilon )+\frac{{\mathit{dE}}_{-}}{d\epsilon }(\epsilon )\right]:\delta \epsilon +{E}_{+}(\epsilon )\frac{\mathit{dB}}{\mathit{db}}\delta b\]

En régime élastique ou saturé, le dernier terme est nul (aucune variation de \(b\) ). En régime endommageant, la variation de \(b\) est obtenue en différentiant l’équation (). Cela conduit à:

(4083)#\[\delta b={\left(\frac{d\varphi }{\mathit{db}}(b)\right)}^{-1}\frac{\mathit{dQ}}{d\epsilon }:d\epsilon\]

La partie délicate consiste alors à dériver \(\Gamma ` , :math:`{E}_{+}\) et \({E}_{-}\) par rapport à \(\epsilon ` . On renvoie pour cela au document [:ref:`R5.03.28 <R5.03.28>\)] dans lequel ces opérations sont détaillées. Bien que ces quantités dépendent des déformations propres, on peut noter que la méthodologie de dérivation qui est mise en œuvre s’affranchit 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.

Pilotage par prédiction élastique#

S’agissant d’une loi d’endommagement, des instabilités matérielles peuvent apparaître. Pour permettre de les suivre, on recommande l’usage d’une méthode de pilotage par prédiction élastique (PRED_ELAS) en l’absence de variables de commande (telles qu’un transitoire thermique) et de chargements dépendant du temps physique, voir document [R5.03.80].

L’incrément de pilotage, noté \(\Delta \tau\) , s’interprète ici comme un incrément d’endommagement maximal cible. Au niveau d’un point d’intégration, il s’agit alors de déterminer la ou les valeurs du paramètre de pilotage \(\eta ` qui placent la déformation sur le seuil d’endommagement. A ce stade, la déformation s’exprime comme une fonction affine de :math:\)eta ` , c’est-à-dire que \(\epsilon ={\epsilon}_{0}+\eta {\epsilon}_{1}\) . Compte tenu de l’expression de la condition de cohérence (), il s’en suit que :math:`eta ` est solution de l’équation scalaire suivante:

(4084)#\[Q({\epsilon}_{0}+\eta {\epsilon}_{1})=\varphi ({b}_{n}+\tau )\]

Compte-tenu de la convexité de \(Q\) , cette équation admet au plus deux solutions. On détaille dans [Lorentz, 2020] comment établir un intervalle borné dans lequel chercher ces solutions et l’application d’une méthode de Newton pour les déterminer.

Description des versions du document#

Version Aster

Auteur(s) Organisme(s)

Description des modifications

15.1

E.LORENTZ, EDF-R&D/ERMES

Texte initial

Bibliographie#

  • Badel P., Godard V., Leblond J.-B. (2007) Application of some anisotropic damage model to the prediction of the failure of some complex industrial concrete structure. Int. J. Solids Struct. 44, 5848-5874.

  • Lorentz E. (2017) A nonlocal damage model for plain concrete consistent with cohesive fracture. Int. J. Fract. 207, 123-159.

  • Lorentz E. (2020) CIWAP3 – Fissuration du béton: proposition d’un modèle local isotrope représentatif de l’endommagement en traction et cisaillement pour simuler des zones fissurées homogènes. EDF R&D 6125-1724-2020-01604.