r7.01.22 Loi de comportement viscoplastique VISC_DRUC_PRAG#
Résumé:
Ce document décrit la loi de comportement viscoplastique VISC_DRUC_PRAG basée sur le modèle élastoplastique de Drucker-Prager et prenant en compte la viscosité selon une loi puissance de type Perzyna. Son domaine d’application est l’argilite qui est la roche hôte du concept de stockage.
Le modèle proposé comporte un seul mécanisme viscoplastique. Le critère s’écrouit avec la déformation viscoplastique cumulée en passant par trois seuils: élastique, de pic et ultime. L’écoulement est non associé, le potentiel d’écoulement étant un potentiel de Drucker-Prager s’écrouissant selon trois niveaux : élastique, de pic et ultime. Entre les seuils, les écrouissages sont linéaires.
Cette loi peut être utilisée dans une modélisation mécanique pure comme elle peut être utilisée dans une modélisation THM. Elle est disponible en 3D, déformations planes et axisymétrique. Elle est intégrée par la résolution d’une seule équation scalaire non linéaire.
Table des matières
Introduction#
Ce document décrit l’intégration de la loi de comportement viscoplastique VISC_DRUC_PRAG dans Code_Aster. Cette loi comporte un seul mécanisme viscoplastique. Le critère viscoplastique s’écrouit avec la déformation viscoplastique déviatorique cumulée en passant par trois seuils : élastique pour une déformation viscoplastique nulle, un seuil dit de pic pour une déformation viscoplastique dite de pic (paramètre du modèle) et un seuil ultime pour une déformation viscoplastique dite ultime (paramètre du modèle). Entre les seuils, les fonctions d’écrouissage sont linéaires. Dans Code_Aster il existe une autre loi basée sur le modèle de Drücker_Prager et utilisée en élastoplasticité sous une forme associée sous le nom DRUCK_PRAGER ou non associée sous le nom DRUCK_PRAG_N_A (voir [R7.01.16]).
Formulation du modèle viscoplastique VISC_DRUC_PRAG#
Équations du modèle#
Ce modèle se base sur une formulation viscoplastique du type Drücker-Prager, où la surface de charge \(f(\sigma ,p)\) est définie par:
\(f=\sqrt{\frac{3}{2}}{s}_{\mathrm{II}}+\alpha (p){I}_{1}-R(p)\)
\(\alpha (p)\) et \(R(p)\) sont des fonctions de la déformation viscoplastique déviatorique cumulée \(p\) ,
On introduit un potentiel d’écoulement viscoplastique \(g(\sigma ,p)\) :
\(g=\sqrt{\frac{3}{2}}{s}_{\text{II}}+\beta (p){I}_{1}\)
Pour l’évolution du critère \(f\) et du potentiel \(g\) nous distinguons trois seuils distincts correspondants à trois valeurs de la variable d’écrouissage: un seuil élastique, un seuil de pic et un seuil ultime. Entre ces seuils, l’écrouissage est linéaire. Entre le seuil élastique et le seuil de pic, l’écrouissage est positif, après le pic l’écrouissage est négatif et devient constant après le seuil ultime.
Les fonctions liées à la cohésion s’écrivent sous la forme suivante:
\(\alpha (p)=(\frac{{\alpha}_{\mathit{pic}}-{\alpha}_{0}}{{p}_{\mathit{pic}}})p+{\alpha}_{0}\) pour \(0<p<{p}_{\mathit{pic}}\)
\(\alpha (p)=(\frac{{\alpha}_{\mathit{ult}}-{\alpha}_{\mathit{pic}}}{{p}_{\mathit{ult}}-{p}_{\mathit{pic}}})(p-{p}_{\mathit{pic}})+{\alpha}_{\mathit{pic}}\) pour \({p}_{\mathit{pic}}<p<{p}_{\mathit{ult}}\)
\(\alpha (p)={\alpha}_{\mathit{ult}}\) pour \(p>{p}_{\mathit{ult}}\)
Les fonctions liées à la dilatance s’écrivent sous la forme suivante:
\(\beta (p)=(\frac{{\beta}_{\mathit{pic}}-{\beta}_{0}}{{p}_{\mathit{pic}}})p+{\beta}_{0}\) pour \(0<p<{p}_{\mathit{pic}}\)
\(\beta (p)=(\frac{{\beta}_{\mathit{ult}}-{\beta}_{\mathit{pic}}}{{p}_{\mathit{ult}}-{p}_{\mathit{pic}}})(p-{p}_{\mathit{pic}})+{\beta}_{\mathit{pic}}\) pour \({p}_{\mathit{pic}}<p<{p}_{\mathit{ult}}\)
\(\beta (p)={\beta}_{\mathit{ult}}\) pour \(p>{p}_{\mathit{ult}}\)
Les fonctions d’écrouissage s’écrivent:
\(R(p)=(\frac{{R}_{\mathit{pic}}-{R}_{0}}{{p}_{\mathit{pic}}})p+{R}_{0}\) pour \(0<p<{p}_{\mathit{pic}}\)
\(R(p)=(\frac{{R}_{\mathit{ult}}-{R}_{\mathit{pic}}}{{p}_{\mathit{ult}}-{p}_{\mathit{pic}}})(p-{p}_{\mathit{pic}})+{R}_{\mathit{pic}}\) pour \({p}_{\mathit{pic}}<p<{p}_{\mathit{ult}}\)
\(R(p)={R}_{\mathit{ult}}\) pour \(p>{p}_{\mathit{ult}}\)
Les contraintes sont reliées aux déformations par la loi de Hooke:
\(\sigma ={D}^{e}(\varepsilon -{\varepsilon}^{\text{vp}})\)
Quand le seuil viscoplastique est atteint, des déformations irréversibles viscoplastiques sont générées et exprimées selon la théorie de Perzyna par:
\(d{\varepsilon}_{ij}^{\text{vp}}=A{\langle \frac{f}{{P}_{\text{ref}}}\rangle }^{n}\frac{\partial g}{\partial {\sigma}_{ij}}\text{dt}\)
\(f\) étant le critère de viscoplasticité; \(A\) et \(n\) sont des paramètres du modèle; \({P}_{\text{ref}}\) une pression de référence.
\(\frac{\partial g}{\partial {\sigma}_{ij}}=\sqrt{\frac{3}{2}}\frac{\partial {s}_{\text{II}}}{\partial {\sigma}_{ij}}+\beta (p)\frac{\partial {I}_{1}}{\partial {\sigma}_{ij}}\) et \(\dot{p}=\sqrt{\frac{2}{3}{\dot{\tilde{\varepsilon}}}_{ij}^{\text{vp}}{\dot{\tilde{\varepsilon}}}_{ij}^{\text{vp}}}\)
avec, \({\tilde{\varepsilon}}_{ij}^{\text{vp}}\) le déviateur du tenseur de déformation viscoplastique,
\(\frac{\partial {s}_{\text{II}}}{\partial {\sigma}_{ij}}=\frac{\partial {s}_{\text{II}}}{\partial {s}_{kl}}\frac{\partial {s}_{kl}}{\partial {\sigma}_{ij}}=\frac{{s}_{kl}}{{s}_{\text{II}}}({\delta}_{\text{ik}}{\delta}_{\text{jl}}-\frac{1}{3}{\delta}_{ij}{\delta}_{kl})=\frac{{s}_{ij}}{{s}_{\text{II}}}\)
et \(\frac{\partial {I}_{1}}{\partial {\sigma}_{ij}}=\frac{\partial \text{tr}({\sigma}_{ij})}{\partial {\sigma}_{ij}}={\delta}_{ij}\)
d’où \(\frac{\partial g}{\partial {\sigma}_{ij}}=\sqrt{\frac{3}{2}}\frac{{s}_{ij}}{{s}_{\text{II}}}+\beta (p){\delta}_{ij}\)
Résumé des équations: Le critère: \(f=\sqrt{\frac{3}{2}}{s}_{\text{II}}+\left[(\frac{{\alpha}_{\text{pic}}-{\alpha}_{0}}{{p}_{\text{pic}}})p+{\alpha}_{0}\right]{I}_{1}-\left[(\frac{{R}_{\text{pic}}-{R}_{0}}{{p}_{\text{pic}}})p+{R}_{0}\right]\) pour \(0<p<{p}_{\text{pic}}\) \(f=\sqrt{\frac{3}{2}}{s}_{\text{II}}+\left[(\frac{{\alpha}_{\text{ult}}-{\alpha}_{\text{pic}}}{{p}_{\text{ult}}-{p}_{\text{pic}}})(p-{p}_{\text{pic}})+{\alpha}_{\text{pic}}\right]{I}_{1}-\left[(\frac{{R}_{\text{ult}}-{R}_{\text{pic}}}{{p}_{\text{ult}}-{p}_{\text{pic}}})(p-{p}_{\text{pic}})+{R}_{\text{pic}}\right]\) pour \({p}_{\text{pic}}<p<{p}_{\text{ult}}\) \(f=\sqrt{\frac{3}{2}}{s}_{\text{II}}+{\alpha}_{\text{ult}}{I}_{1}-{R}_{\text{ult}}\) pour \(p\ge {p}_{\text{ult}}\) : Le potentield’écoulement: \(g=\sqrt{\frac{3}{2}}{s}_{\text{II}}+\left[(\frac{{\beta}_{\text{pic}}-{\beta}_{0}}{{p}_{\text{pic}}})p+{\beta}_{0}\right]{I}_{1}\) pour \(0<p<{p}_{\text{pic}}\) \(g=\sqrt{\frac{3}{2}}{s}_{\text{II}}+\left[(\frac{{\beta}_{\text{ult}}-{\beta}_{\text{pic}}}{{p}_{\text{ult}}-{p}_{\text{pic}}})(p-{p}_{\text{pic}})+{\beta}_{\text{pic}}\right]{I}_{1}\) pour \({p}_{\text{pic}}<p<{p}_{\text{ult}}\) \(g=\sqrt{\frac{3}{2}}{s}_{\text{II}}+{\beta}_{\text{ult}}{I}_{1}\) pour \(p\ge {p}_{\text{ult}}\) \({\alpha}_{0}\) , \({R}_{0}\) et \({\beta}_{0}\) : paramètres d’écrouissage correspondant au seuil d’élasticité (\(p=0\) ) \({\alpha}_{\text{pic}}\) , \({R}_{\text{pic}}\) et \({\beta}_{\text{pic}}\) : paramètres d’écrouissage correspondant au paramètre \({p}_{\text{pic}}\) \({\alpha}_{\text{ult}}\) , \({R}_{\text{ult}}\) et \({\beta}_{\text{ult}}\) : paramètres d’écrouissage correspondant au paramètre \({p}_{\text{ult}}\) La loi de Hooke: \(\sigma ={D}^{e}(\varepsilon -{\varepsilon}^{\text{vp}})\) \(f(\sigma ,p)\le 0\) domaine d’élasticité; \({\dot{\varepsilon}}_{ij}^{\text{vp}}=0\) \(f(\sigma ,p)>0\) viscoplasticité ; \({\dot{\varepsilon}}_{ij}^{\text{vp}}=A{\langle \frac{f}{{P}_{\text{ref}}}\rangle }^{n}\frac{\partial g}{\partial {\sigma}_{ij}}\) ; \(\dot{p}=\sqrt{\frac{2}{3}{\dot{\tilde{\varepsilon}}}_{ij}^{\text{vp}}{\dot{\tilde{\varepsilon}}}_{ij}^{\text{vp}}}\) |
Intégration dans Code_Aster#
Décomposition du tenseur de déformation#
La décomposition de l’incrément de déformation totale s’écrit:
\(\Delta \varepsilon =\Delta {\varepsilon}^{e}+\Delta {\varepsilon}^{\text{vp}}\)
où \(\Delta {\varepsilon}^{e}\) et \(\Delta {\varepsilon}^{\text{vp}}\) sont les incréments des tenseurs élastiques et viscoplastiques.
Mise à jour des contraintes#
Les notations suivantes sont adoptées: \({A}^{-}\) , \(A\) et \(\Delta A\) désignant respectivement les quantités au début, à la fin du pas de temps et son incrément durant le pas.
On exprime les contraintes actualisées à l’instant \({t}^{+}\) par rapport à celles calculées à l’instant \({t}^{-}\) :
\(\sigma ={\sigma}^{-}+{D}^{e}\Delta {\varepsilon}^{e}\) ; \(s={s}^{-}+2\mu \Delta {\tilde{\varepsilon}}^{e}\) ; \({I}_{1}={I}_{1}^{-}+3K\Delta {\varepsilon}_{v}^{e}\)
\({\sigma}_{ij}={s}_{ij}+\frac{{I}_{1}}{3}{\delta}_{ij}\)
\(\Delta {\varepsilon}_{ij}=\Delta \tilde{\varepsilon}+\text{tr}\frac{(\Delta \varepsilon )}{3}{\delta}_{ij}=\Delta \tilde{{\varepsilon}_{ij}}+\frac{\Delta {\varepsilon}_{v}}{3}{\delta}_{ij}\)
\({I}_{1}=\text{tr}(\sigma )\) ; \({\varepsilon}_{v}=\text{tr}(\Delta \varepsilon )\)
Prédiction élastique:
\({\sigma}^{\mathrm{el}}={\sigma}^{-}+{D}^{e}\Delta \varepsilon\) ; \({s}^{\mathrm{el}}={s}^{-}+2\mu \Delta \tilde{\varepsilon}\) ; \({I}_{1}^{\mathrm{el}}={I}_{1}^{-}+3K\Delta {\varepsilon}_{v}\)
Solution élastique
Calcul de l’incrément des contraintes en régime élastique:
\(\Delta {\sigma}_{ij}=\Delta {s}_{ij}+\frac{\Delta {I}_{1}}{3}{\delta}_{ij}\) ; \(\Delta {\varepsilon}_{ij}=\Delta {\tilde{\varepsilon}}_{ij}+\frac{\Delta {\varepsilon}_{v}}{3}{\delta}_{ij}\)
\(\Delta {\sigma}_{ij}=2\mu \Delta \tilde{{\varepsilon}_{ij}}+3K\frac{\Delta {\varepsilon}_{v}}{3}{\delta}_{ij}=2\mu \Delta \tilde{{\varepsilon}_{ij}}+K\Delta {\varepsilon}_{v}{\delta}_{ij}=2\mu (\Delta {\epsilon}_{ij}-\frac{\text{tr}(\Delta \varepsilon )}{3}{\delta}_{ij})+K\text{tr}(\Delta \varepsilon ){\delta}_{ij}\)
\(\Delta {\sigma}_{ij}=2\mu \Delta {\varepsilon}_{ij}+(K-\frac{2G}{3})\text{tr}(\Delta \varepsilon ){\delta}_{ij}\)
\(\left\lbrace \begin{array}{c}\Delta {\sigma}_{11}\\ \Delta {\sigma}_{22}\\ \Delta {\sigma}_{33}\\ \sqrt{2}\Delta {\sigma}_{12}\\ \sqrt{2}\Delta {\sigma}_{13}\\ \sqrt{2}\Delta {\sigma}_{23}\end{array}\right\rbrace =\underset{{D}^{e}}{\underset{\underbrace{}}{\left[\begin{array}{cccccc}\frac{4\mu }{3}+K& K-\frac{2\mu }{3}& K-\frac{2\mu }{3}& 0& 0& 0\\ K-\frac{2\mu }{3}& \frac{4\mu }{3}+K& K-\frac{2\mu }{3}& 0& 0& 0\\ K-\frac{2\mu }{3}& K-\frac{2\mu }{3}& \frac{4\mu }{3}+K& 0& 0& 0\\ 0& 0& 0& 2\mu & 0& 0\\ 0& 0& 0& 0& 2\mu & 0\\ 0& 0& 0& 0& 0& 2\mu \end{array}\right]}}.\left\lbrace \begin{array}{c}\Delta {\varepsilon}_{11}\\ \Delta {\varepsilon}_{22}\\ \Delta {\varepsilon}_{33}\\ \sqrt{2}\Delta {\varepsilon}_{12}\\ \sqrt{2}\Delta {\varepsilon}_{13}\\ \sqrt{2}\Delta {\varepsilon}_{23}\end{array}\right\rbrace\)
Solution viscoplastique
On exprime le champ de contraintes à l’instant \({t}^{+}\) :
\({\sigma}_{ij}={\sigma}_{ij}^{-}+{D}_{ijkl}^{e}\Delta {\varepsilon}_{kl}^{e}={\sigma}_{ij}^{-}+{D}_{ijkl}^{e}(\Delta {\varepsilon}_{kl}-\Delta {\varepsilon}_{kl}^{\text{vp}})={\sigma}_{ij}^{\mathit{el}}-{D}_{ijkl}^{e}\Delta {\varepsilon}_{kl}^{\text{vp}}\)
\({s}_{ij}={s}_{ij}^{\mathit{el}}-2\mu \Delta {\tilde{\varepsilon}}_{ij}^{\text{vp}}\) et \({I}_{1}={I}_{1}^{\mathrm{el}}-3K\Delta {\varepsilon}_{v}^{\text{vp}}\)
\({\sigma}_{ij}={s}_{ij}+\frac{{I}_{1}}{3}{\delta}_{ij}\)
qui s’écrit en remplaçant l’accroissement des déformations visqueuses par leurs expressions sous la forme:
\({\sigma}_{ij}={\sigma}_{ij}^{\mathrm{el}}-{D}_{ijkl}^{\mathrm{el}}\langle \Phi \rangle \frac{\partial G}{\partial {\sigma}_{ij}}(\sigma ,p)\Delta t\) avec \(\Phi =A{(\frac{f(\sigma ,p)}{{P}_{\text{ref}}})}^{n}\) où
\(\Phi\) et \(\frac{\partial g}{\partial \sigma }\) caractérisent l’amplitude et la direction de la vitesse des déformations irréversibles.
\(f(\sigma ,p)\) étant le critère de viscoplasticité, \(A\) et \(n\) sont des paramètres du modèle.
Le critère viscoplastique à l’instant \({t}^{+}\) s’écrit:
\(f(\sigma ,p)=f({\sigma}_{ij}^{\mathrm{el}}-{D}_{ijkl}^{e}\langle \Phi \rangle \frac{\partial g}{\partial {\sigma}_{ij}}(\sigma ,p)\Delta t,p)\)
L’incrément de la déformation viscoplastique étant:
\(\Delta {\varepsilon}_{ij}^{\text{vp}}=\langle \Phi \rangle \frac{\partial g}{\partial {\sigma}_{ij}}\Delta t=\langle \Phi \rangle (\sqrt{\frac{3}{2}}\frac{{s}_{ij}}{{s}_{\text{II}}}+\beta (p){\delta}_{ij})\Delta t\)
La déformation volumique viscoplastique étant:
\(\Delta {\varepsilon}_{v}^{\text{vp}}=3\langle \Phi \rangle \beta (p)\Delta t\)
La composante déviatorique de la déformation viscoplastique s’écrit sous la forme:
\(\Delta {\tilde{\varepsilon}}_{ij}^{\text{vp}}=\langle \Phi \rangle \sqrt{\frac{3}{2}}\frac{{s}_{ij}}{{s}_{\text{II}}}\Delta t\) ou \(\Delta {\tilde{\varepsilon}}_{ij}^{\text{vp}}=\langle \Phi \rangle \frac{3}{2}\frac{{s}_{ij}}{{\sigma}_{\text{eq}}}\Delta t\)
comme \({\sigma}_{\text{eq}}=\sqrt{\frac{3}{2}}{s}_{\text{II}}\) , \({s}_{\text{II}}=\sqrt{{s}_{ij}{s}_{ij}}\) et \({\sigma}_{\text{eq}}^{\mathit{el}}=\sqrt{\frac{3}{2}{s}_{ij}^{\mathit{el}}{s}_{ij}^{\mathit{el}}}\)
On écrit aussi les égalités suivantes:
\({s}_{ij}\frac{{\sigma}_{\text{eq}}^{\mathit{el}}}{{\sigma}_{\text{eq}}}={s}_{ij}^{\mathit{el}}\)
\(\Delta p=\sqrt{(\frac{2}{3}\Delta {\tilde{\epsilon}}_{ij}^{\mathit{vp}}\Delta {\tilde{\epsilon}}_{ij}^{\mathit{vp}})}\)
\(\frac{\Delta p}{\Delta t}=\langle \Phi \rangle =A{\langle \frac{f(\sigma ,p)}{{P}_{\mathrm{ref}}}\rangle }^{n}\) éq 1
d’où : \(\Delta p=\langle \phi \rangle \Delta t\)
En utilisant ces égalités on peut trouver une expression pour \({s}_{ij}\) , \({\sigma}_{\mathit{eq}}\) et \({I}_{1}\) en fonction de \({s}_{ij}^{\mathit{el}}\) , \({\sigma}_{\mathit{eq}}^{\mathit{el}}\) , \({I}_{1}^{\mathit{el}}\) et \(\Delta p\) :
\({s}_{ij}={s}_{ij}^{\mathit{el}}-2\mu \Delta {\tilde{\varepsilon}}_{ij}^{\text{vp}}={s}_{ij}^{\mathit{el}}-3\mu \langle \Phi \rangle \frac{{s}_{ij}}{{\sigma}_{\text{eq}}}\Delta t={s}_{ij}^{\mathit{el}}-3\mu \langle \Phi \rangle \frac{{s}_{ij}^{\mathit{el}}}{{\sigma}_{\text{eq}}^{\mathit{el}}}\Delta t\)
\({s}_{ij}={s}_{ij}^{\mathit{el}}(1-\frac{3\mu \langle \Phi \rangle }{{\sigma}_{\text{eq}}^{\mathit{el}}}\Delta t)={s}_{ij}^{\mathit{el}}(1-\frac{3\mu }{{\sigma}_{\text{eq}}^{\mathit{el}}}\Delta p)\)
\({\sigma}_{\text{eq}}={\sigma}_{\text{eq}}^{\mathit{el}}-3\mu \langle \Phi \rangle \Delta t={\sigma}_{\text{eq}}^{\mathit{el}}-3\mu \Delta p\) éq 2
\({I}_{1}={I}_{1}^{\mathrm{el}}-3K\Delta {\varepsilon}_{v}^{\text{vp}}={I}_{1}^{\mathrm{el}}-9K\beta \langle \Phi \rangle \Delta t={I}_{1}^{\mathrm{el}}-9K\beta \Delta p\) éq 3
Calcul de l’inconnue
L’incrément de déformation viscoplastique cumulée \(\Delta p\) est la seule inconnue du problème. Pour le déterminer, on écrit la loi d’écoulement viscoplastique (éq 1):
\(\frac{\Delta p}{\Delta t}=A{\langle \frac{{\sigma}^{\text{eq}}+\alpha (p){I}_{1}-R(p)}{{P}_{\text{ref}}}\rangle }^{n}\)
\(R(p)=R({p}^{-}+\Delta p)={R}^{-}+{R}_{\text{const}}\Delta p\) ; \({R}_{\text{const}}=\frac{\partial R}{\partial p}\)
\(\alpha (p)=\alpha ({p}^{-}+\Delta p)={\alpha}^{-}+{\alpha}_{\text{const}}\Delta p\) ; \({\alpha}_{\text{const}}=\frac{\partial \alpha }{\partial p}\)
\(\beta (p)=\beta ({p}^{-}+\Delta p)={\beta}^{-}+{\beta}_{\text{const}}\Delta p\) ; \({\beta}_{\text{const}}=\frac{\partial \beta }{\partial p}\)
Par souci de simplification de l’écriture de l’équation en \(\Delta p\) , on pose:
\(C=\frac{(A\Delta t)}{{{P}_{\mathit{ref}}}^{n}}\)
Soit, en remplaçant \({\sigma}_{\mathit{eq}}\) et \({I}_{1}\) par leurs expressions (éq 2 et éq 3), on obtient:
\(F(\Delta p)=C{\langle \begin{array}{c}({\sigma}_{\mathrm{eq}}^{\mathrm{el}}+\alpha {I}_{1}^{\mathrm{el}}-{R}^{-})-(3\mu +{R}_{\mathrm{const}}-{\alpha}_{\mathrm{const}}{I}_{1}^{\mathrm{el}}+9K{\alpha}^{-}{\beta}^{-})\Delta p\\ \\ (9K{\alpha}^{-}{\beta}_{\mathrm{const}}+9K{\alpha}_{\mathrm{const}}{\beta}^{-})\Delta {p}^{2}-(9K{\alpha}_{\mathrm{const}}{\beta}_{\mathrm{const}})\Delta {p}^{3}\end{array}\rangle }^{n}-\Delta p=0\)
On cherche \(\Delta p\) tel que \(F(\Delta p)=0\) .
\(F(\Delta p)=0\) est une équation scalaire non linéaire. La borne inférieure étant \({x}_{\inf}=0\) et la borne supérieure peut être fixée à:
\({x}_{\sup}=A{\langle \frac{{\sigma}_{\text{eq}}^{\mathit{el}}+\alpha {I}_{1}^{\mathit{el}}-{R}^{-}}{{P}_{\text{ref}}}\rangle }^{n}\Delta t\)
On utilise la méthode des cordes avec un contrôle de l’intervalle de recherche en s’inspirant du document [R5.03.04].
\(\Delta p\in \left[{x}_{\inf},{x}_{\sup}\right]\) ;
\(x=\Delta p\)
Si \(\mid F({x}_{\inf})\mid <\eta\) alors \(\Delta p={x}_{\inf}\)
Si \(\mid F({x}_{\sup})\mid <\eta\) alors \(\Delta p={x}_{\sup}\)
Si \(F({x}_{\inf})>0\) alors \({x}_{2}={x}_{\inf}\) et \({y}_{2}=F({x}_{\inf})\)
Si \(F({x}_{\sup})<0\) alors on fait une boucle en découpant \({x}_{\sup}\) par 10 jusqu’à obtenir une valeur de \({x}_{\sup}\) pour laquelle \(F({x}_{\sup})>0\) dans ce cas on multiplie la dernière solution par 10 et on fixe \({x}_{1}={x}_{\sup}\) et \({y}_{1}=F({x}_{\sup})\)
Si \(F({x}_{\sup})>0\) alors on fait une boucle en multipliant \({x}_{\sup}\) par 10 jusqu’à obtenir une valeur de \({x}_{\sup}\) pour laquelle \(F({x}_{\sup})<0\) et on fixe \({x}_{1}={x}_{\sup}\) et \({y}_{1}=F({x}_{\sup})\)
Si \(F({x}_{\inf})<0\) alors \({x}_{1}={x}_{\inf}\) et \({y}_{1}=F({x}_{\inf})\)
Si \(F({x}_{\sup})>0\) alors on fait une boucle en découpant \({x}_{\sup}\) par 10 jusqu’à obtenir une valeur de \({x}_{\sup}\) pour laquelle \(F({x}_{\sup})<0\) dans ce cas on multiplie la dernière solution par 10 et on fixe \({x}_{2}={x}_{\sup}\) et \({y}_{2}=F({x}_{\sup})\)
Si \(F({x}_{\sup})<0\) alors on fait une boucle en multipliant \({x}_{\sup}\) par 10 jusqu’à obtenir une valeur de \({x}_{\sup}\) pour laquelle \(F({x}_{\sup})>0\) et on fixe \({x}_{2}={x}_{\sup}\) et \({y}_{2}=F({x}_{\sup})\)
Des vérifications sont faites sur les valeurs que peuvent prendre les bornes et notamment si elles sont plus faibles qu’une tolérance fixée à 1.E-12, elles seront considérées égales à 0. et donc la solution \(\Delta p\) également. Si les bornes sont égales, on fait un re-découpage du pas de temps.
Les valeurs \({x}_{1}\) , \({x}_{2}\) , \({y}_{1}\) et \({y}_{2}\) seront les valeurs à donner en entrée à la routine zeroco qui se base sur la méthode des cordes. La solution est calculée par la formule suivante:
\({x}^{n+1}={x}^{n-1}-F({x}^{n-1})\frac{{x}^{n}-{x}^{n-1}}{F({x}^{n})-F({x}^{n-1})}\)
Avec les valeurs suivantes, on représente la fonction scalaire à résoudre.
\({\sigma}_{\mathit{eq}}^{\mathit{el}}\) |
\(6,315\mathit{MPa}\) |
\({\alpha}^{-}\) |
\(6,86{10}^{-2}\) |
\({I}_{1}^{\mathit{el}}\) |
\(-21,061\mathit{MPa}\) |
\({\beta}^{-}\) |
\(-0,147\) |
\(N\) |
\(4,5\) |
\({R}^{-}\) |
\(1,394\mathit{MPa}\) |
\(\Delta t\) |
\(10s\) |
\({\alpha}_{\mathit{const}}\) |
\(13.\) |
\(A\) |
\(1,5{10}^{-12}\) |
\({\beta}_{\mathit{const}}\) |
\(10.\) |
\({P}_{\mathit{ref}}\) |
\(0,1\mathit{MPa}\) |
\({R}_{\mathit{const}}\) |
\(329,732\mathit{MPa}\) |
L’inconnue \(x\) pour laquelle \(F(x)\) s’annule se situe entre \(6.{10}^{-5}\) et \(7.{10}^{-5}\) qui se situe bien entre la borne inférieure \({x}_{inf}\) et la borne supérieure \({x}_{\sup}\) qui vaut dans ce cas précis \(1,2913{10}^{-4}\) .
Figure 4-1: Allure de la fonction scalaire
Opérateur tangent cohérent#
On cherche à calculer: \(\frac{\partial \sigma }{\partial \varepsilon }=\frac{\partial s}{\partial \varepsilon }+\frac{1}{3}\mathrm{Id}\otimes \frac{\partial {I}_{1}}{\partial \varepsilon }\)
Avec:
\(\frac{\partial s}{\partial \varepsilon }=\frac{\partial {s}^{\mathrm{el}}}{\partial \varepsilon }(1-\frac{3\mu }{{\sigma}_{\text{eq}}^{\mathrm{el}}}.\Delta p)+\frac{3\mu }{{({\sigma}_{\text{eq}}^{\mathrm{el}})}^{2}}.\Delta p({s}^{\mathrm{el}}\otimes \frac{\partial {\sigma}_{\text{eq}}^{\mathrm{el}}}{\partial \varepsilon })-\frac{3\mu }{{\sigma}_{\text{eq}}^{\mathrm{el}}}.({s}^{\mathrm{el}}\otimes \frac{\partial \Delta p}{\partial \varepsilon })\)
\(\frac{\partial {I}_{1}}{\partial \varepsilon }=\frac{\partial {I}_{1}^{\mathrm{el}}}{\partial \varepsilon }-9K\beta (p)\frac{\partial \Delta p}{\partial \varepsilon }\)
Calcul de \(\frac{\partial {s}^{\mathrm{el}}}{\partial \varepsilon }\) :
\(\frac{\partial {s}^{\mathrm{el}}}{\partial \varepsilon }=2\mu ({\mathrm{Id}}_{4}-\frac{1}{3}\mathrm{Id}\otimes \mathrm{Id})\)
\(\frac{\partial {s}_{ij}}{\partial {\varepsilon}_{\text{pq}}}=2\mu ({\delta}_{\text{ip}}{\delta}_{\text{jq}}-\frac{1}{3}{\delta}_{ij}{\delta}_{\text{pq}})\)
Calcul de \(\frac{\partial {I}_{1}^{\mathrm{el}}}{\partial \varepsilon }\) :
\(\frac{\partial {I}_{1}^{\mathrm{el}}}{\partial \varepsilon }=3K\mathrm{Id}\) soit:
\(\frac{\partial {I}_{1}^{\mathrm{el}}}{\partial {\varepsilon}_{\text{pq}}}=3K{\delta}_{\text{pq}}\)
Calcul de \(\frac{\partial {\sigma}_{\text{eq}}^{\mathrm{el}}}{\partial \varepsilon }\) :
\(\frac{\partial {\sigma}_{\text{eq}}^{\mathrm{el}}}{\partial \varepsilon }=\frac{\partial {\sigma}_{\text{eq}}^{\mathrm{el}}}{\partial {\sigma}^{\mathrm{el}}}\frac{\partial {\sigma}^{\mathrm{el}}}{\partial \varepsilon }=\frac{\partial {\sigma}_{\text{eq}}^{\mathrm{el}}}{\partial {s}^{\mathrm{el}}}\frac{\partial {s}^{\mathrm{el}}}{\partial {\sigma}^{\mathrm{el}}}\frac{\partial {\sigma}^{\mathrm{el}}}{\partial \varepsilon }=\frac{3}{2}\frac{{s}^{\mathrm{el}}}{{\sigma}_{\text{eq}}^{\mathrm{el}}}({\mathrm{Id}}_{4}-\frac{1}{3}\mathrm{Id}\otimes \mathrm{Id}){D}^{e}=\frac{3}{2}\frac{{s}^{\mathrm{el}}}{{\sigma}_{\text{eq}}^{\mathrm{el}}}{D}^{e}\)
Calcul de \(\frac{\partial \Delta p}{\partial \varepsilon }\) :
\(\frac{\Delta p}{\Delta t}=A{\langle \frac{f(\sigma ,p)}{{P}_{\text{ref}}}\rangle }^{n}\)
soit \(F(\Delta p)=\frac{A\Delta t}{{P}_{\mathrm{ref}}^{n}}{\langle f(\sigma ,p)\rangle }^{n}-\Delta p\)
\(\frac{\partial \Delta p}{\partial \varepsilon }=\frac{\partial \Delta p}{\partial {\sigma}^{\mathrm{el}}}\frac{\partial {\sigma}^{\mathrm{el}}}{\partial \varepsilon }\)
pour calculer \(\frac{\partial \Delta p}{\partial {\sigma}^{\mathrm{el}}}\) , on utilise \(F({\sigma}^{\mathrm{el}},p)=0\)
\(\frac{\partial F({\sigma}^{\mathrm{el}},p)}{\partial {\sigma}^{\mathrm{el}}}\delta {\sigma}^{\mathrm{el}}+\frac{\partial F({\sigma}^{\mathrm{el}},p)}{\partial \Delta p}\delta \Delta p=0\) d’où:
\(\frac{\delta \Delta p}{\delta {\sigma}^{\mathrm{el}}}=-\frac{\frac{\partial F({\sigma}^{\mathrm{el}},p)}{\partial {\sigma}^{\mathrm{el}}}}{\frac{\partial F({\sigma}^{\mathrm{el}},p)}{\partial \Delta p}}\)
\(F(\Delta p)=C{\langle \begin{array}{c}({\sigma}_{\mathrm{eq}}^{\mathrm{el}}+\alpha {I}_{1}^{\mathrm{el}}-{R}^{-})-(3\mu +{R}_{\mathrm{const}}-{\alpha}_{\mathrm{const}}{I}_{1}^{\mathrm{el}}+9K{\alpha}^{-}{\beta}^{-})\Delta p\\ \\ (9K{\alpha}^{-}{\beta}_{\mathrm{const}}+9K{\alpha}_{\mathrm{const}}{\beta}^{-})\Delta {p}^{2}-(9K{\alpha}_{\mathrm{const}}{\beta}_{\mathrm{const}})\Delta {p}^{3}\end{array}\rangle }^{n}-\Delta p=0\)
\(\frac{\partial F({\sigma}^{\mathrm{el}},p)}{\partial {\sigma}^{\mathrm{el}}}=\mathrm{C.}n.{\langle f({\sigma}^{\mathrm{el}},p)\rangle }^{n-1}\frac{\partial f({\sigma}^{\mathrm{el}},p)}{\partial {\sigma}^{\mathrm{el}}}\)
où \(\frac{\partial f({\sigma}^{\mathrm{el}},p)}{\partial {\sigma}^{\mathrm{el}}}=(\frac{\partial {\sigma}_{\text{eq}}^{\mathrm{el}}}{\partial {\sigma}^{\mathrm{el}}}+\alpha \frac{\partial {I}_{1}^{\mathrm{el}}}{\partial {\sigma}^{\mathrm{el}}})+{\alpha}_{\text{const}}(\frac{\partial {I}_{1}^{\mathrm{el}}}{\partial {\sigma}^{e}})\Delta p\)
\(\frac{\partial F({\sigma}^{\mathrm{el}},p)}{\partial \Delta p}=\mathrm{C.}n.{\langle f({\sigma}^{\mathrm{el}},p)\rangle }^{n-1}\frac{\partial f({\sigma}^{\mathrm{el}},p)}{\partial \Delta p}-1\)
\(\begin{array}{}\frac{\partial f({\sigma}^{\mathrm{el}},p)}{\partial \Delta p}=-(3\mu +{R}_{\mathrm{const}}-{\alpha}_{\mathrm{const}}{I}_{1}^{\mathrm{el}}+9K{\alpha}^{-}{\beta}^{-})\\ -2\Delta p9K({\alpha}^{-}{\beta}_{\mathrm{const}}+{\alpha}_{\mathrm{const}}{\beta}^{-})-3\Delta {p}^{2}9K({\alpha}_{\mathrm{const}}{\beta}_{\mathrm{const}})\end{array}\)
Calcul de \(\frac{\partial {\sigma}_{\text{eq}}^{\mathrm{el}}}{\partial {\sigma}^{\mathrm{el}}}\) :
\(\frac{\partial {\sigma}_{\text{eq}}^{\mathrm{el}}}{\partial {\sigma}^{\mathrm{el}}}=\frac{\partial {\sigma}_{\text{eq}}^{\mathrm{el}}}{\partial {s}^{\mathrm{el}}}\frac{\partial {s}^{\mathrm{el}}}{\partial {\sigma}^{\mathrm{el}}}=\frac{3}{2}\frac{{s}^{\mathrm{el}}}{{\sigma}_{\text{eq}}^{\mathrm{el}}}.({I}_{4}-\frac{1}{3}\mathrm{Id}\otimes \mathrm{Id})=\frac{3}{2}\frac{{s}^{\mathrm{el}}}{{\sigma}_{\text{eq}}^{\mathrm{el}}}\)
Calcul de \(\frac{\partial {I}_{1}^{\mathrm{el}}}{\partial {\sigma}^{\mathrm{el}}}\) :
\(\frac{\partial {I}_{1}^{\mathrm{el}}}{\partial {\sigma}^{\mathrm{el}}}=\frac{\partial \text{tr}({\sigma}^{\mathrm{el}})}{\partial {\sigma}^{\mathrm{el}}}=\mathrm{Id}\)
Données matériaux#
Les 16 paramètres du modèle sont:
sous ELAS
\(E\) : module de Young (\(\mathit{Pa}\) ou \(\mathit{MPa}\) )
\(\nu\) : coefficient de Poisson
sous VISC_DRUC_PRAG
\({P}_{\text{ref}}\) : pression de référence (\(\mathit{Pa}\) ou \(\mathit{MPa}\) )
\(A\) : paramètre viscoplastique (en \({s}^{-1}\) )
\(n\) : puissance de la loi fluage
\({p}_{\text{pic}}\) : taux d’écrouissage au niveau du seuil de pic
\({p}_{\text{ult}}\) : taux d’écrouissage au niveau du seuil ultime
\({\alpha}_{0}\) , \({\alpha}_{\text{pic}}\) et \({\alpha}_{\text{ult}}\) : paramètres de la fonction de cohésion \(\alpha (p)\)
\({R}_{0}\) , \({R}_{\text{pic}}\) et \({R}_{\text{ult}}\) : paramètres de la fonction d’écrouissage \(R(p)\)
\({\beta}_{0}\) , \({\beta}_{\text{pic}}\) et \({\beta}_{\text{ult}}\) : paramètres de la fonction de dilatance \(\beta (p)\)
Les variables internes#
\({v}_{1}=p\) : déformation viscoplastique déviatorique cumulée;
\({v}_{2}=(\text{0 ou 1})\) : indicateur de plasticité;
\({v}_{3}=\mathit{pos}\) : position du point de charge par rapport au seuil:
(\(\text{pos}=\text{1 si 0}<p<{p}_{\text{pic}}\) ; \(\text{pos}={\text{2 si p}}_{\text{pic}}<p<{p}_{\text{ult}}\) ; \(\text{pos}=\text{3 si p}>{p}_{\text{ult}}\) )
\({v}_{4}\) : nombre d’itérations locales.
Résumé de l’algorithme de résolution#
L’algorithme de résolution tel qu’il est implanté dans Code_Aster:
\({\sigma}^{\mathrm{el}}={\sigma}^{-}+{D}^{e}\Delta \varepsilon\)
Le critère: \(f({\sigma}^{\mathrm{el}},{p}^{-})={\sigma}_{\text{eq}}^{\mathrm{el}}+\alpha ({p}^{-}){I}_{1}^{\mathrm{el}}-R({p}^{-})\)
Élasticité: si \(f({\sigma}^{\mathrm{el}},{p}^{-})\le 0\) alors \(\Delta p=0\) ;
Viscoplasticité: si \(f({\sigma}^{\mathrm{el}},{p}^{-})=0\) alors \(\Delta p\ge 0\) avec \(\Delta p\) solution de l’équation \(F(\Delta p)=0\)
où:
\(\frac{\Delta p}{\Delta t}=A{\langle \frac{{\sigma}^{\text{eq}}+\alpha (p){I}_{1}-R(p)}{{P}_{\text{ref}}}\rangle }^{n}=\frac{A}{{P}_{\text{ref}}^{n}}{\langle f(\sigma ,p)\rangle }^{n}\)
et
\(F=\frac{A\Delta t}{{P}_{\mathrm{ref}}^{n}}{\langle f(\sigma ,p)\rangle }^{n}-\Delta p\)
Mise à jour du tenseur des contraintes:
\(\sigma ={\sigma}^{\mathrm{el}}-{D}^{e}\Delta {\varepsilon}^{\mathrm{vp}}\)
\(s={s}^{\mathrm{el}}(1-3\frac{\mu}{{\sigma}_{\mathrm{eq}}^{\mathrm{el}}}\Delta p)\)
\({\sigma}_{\mathrm{eq}}={\sigma}_{\mathrm{eq}}^{\mathrm{el}}-3G\Delta p\)
\({I}_{1}={{I}_{1}}^{\mathrm{el}}-9K\beta \Delta p\)
\(\sigma =s+\frac{1}{3}{I}_{1}\mathrm{Id}\)
Une fois que \(\Delta p\) est calculée, les contraintes et les variables internes mises à jour, on vérifie la position de \(p\) par rapport à \({p}^{-}\) et le signe de \(f(\sigma ,p)\) :
Si \(0<{p}^{-}<{p}_{\text{pic}}\) ; tester 1) sinon 2) sinon 3)
Si \({p}_{\text{pic}}<{p}^{-}<{p}_{\text{ult}}\) ; tester 2) sinon 3)
Si \({p}^{-}>{p}_{\text{ult}}\) ; tester 3)
Si \({p}^{-}+\Delta p<{p}_{\text{pic}}\) ;
on vérifie \(f(\sigma ,p)>0\) avec \(R\) , \(\alpha ` et :math:\)mathrm{beta }` correspondantes à \(0<p<{p}_{\text{pic}}\) ,
si \(f(\sigma ,p)>0\) alors on met à jour les champs de contraintes
et de variables internes,
sinon, on considère que \(\Delta p\) n’est pas valable et on re-découpe le pas de temps
on vérifie \(f(\sigma ,p)>0\) avec \(R\) , \(\alpha ` et :math:\)beta` correspondantes à \({p}_{\text{pic}}<p<{p}_{\text{ult}}\)
si \(f(\sigma ,p)>0\) alors on met à jour les champs de contraintes
et de variables internes,
sinon, on considère que \(\Delta p\) n’est pas valable et on redécoupe le pas de temps
on vérifie \(f(\sigma ,p)>0\) avec \(R\) , \(\alpha\) et \(\beta\) correspondantes à \(p\ge {p}_{\text{ult}}\)
si \(f(\sigma ,p)>0\) alors on met à jour les champs de contraintes
et de variables internes,
sinon, on considère que \(\Delta p\) n’est pas valable et on redécoupe le pas de temps
Résultats d’un essai triaxial#
Il s’agit de simuler un essai triaxial (voir le cas test ssnv211) en imposant comme une contrainte de confinement de \(5\mathit{MPa}\) . Une déformation uniaxiale est imposée en compression et qui évolue dans le temps.
La vitesse du chargement est fixée à \({10}^{-5}m/s\) . Le déviateur des contraintes et la déformation volumique en fonction de la déformation axiale imposée sont représentées ci contre.
Figure 5-1: Déviateur des contraintes en fonction de la déformation uniaxiale
Figure 5-2: Déformation volumique en fonction de la déformation uniaxiale
Fonctionnalités et vérification#
La loi de comportement peut être définie par le mot-clé VISC_DRUC_PRAG (commande STAT_NON_LINE, mot clé facteur COMPORTEMENT). Elle est associée au matériau VISC_DRUC_PRAG (commande DEFI_MATERIAU).
La loi VISC_DRUC_PRAG est vérifiée par les cas tests suivants:
Références#
EL GHARIB et C. CHAVANT, “ Calage sur des essais triaxiaux d’une loi de comportement viscoplastique pour l’argilite basée sur le modèle Drucker_Prager », H-T64-2008-04194-FR,
EL GHARIB et C. CHAVANT, «Mise en oeuvre dans Code_Aster d’un modèle viscoplastique simplifié», H-T64-2007-01800-FR,
Descriptif des versions du document#
Version Aster |
Auteur(s) Organisme(s) |
Description des modifications |
10.0 |
|
Texte initial |
14.4 |
F.Voldoire EDF R&D / ERMES |
Petites corrections de forme des équations |