r7.01.23 Loi de comportement cyclique de HUJEUX pour les sols#
Résumé:
Le modèle de comportement dit de «Hujeux», conçu au laboratoire MSSMat de l’ECP [bib5], est un des modèles élastoplastiques cycliques de mécanique des sols (géomatériaux granulaires: argiles sableuses, normalement consolidées ou sur-consolidées, graves…) le plus adapté pour des simulations d’ouvrages géotechniques en séisme. Il est de plus exploité depuis de nombreuses années, son paramétrage étant donc bien maîtrisé.
Ce modèle multi-mécanismes (sphériques – pour un trajet de consolidation – et déviatoires) à variables mémoratrices est caractérisé par huit surfaces de charge avec écrouissage, définies pour des trajets monotones et pour des trajets cycliques. Les mécanismes sont définis par plans fixes, ce qui induit une orthotropie de comportement du sol. A l’intérieur de ces surfaces de réversibilité, le matériau est élastique non linéaire. L’écrouissage est régi par plusieurs variables et la règle d’écoulement normal est adoptée pour les mécanismes de consolidation, tandis que la règle d’écoulement pour les mécanismes déviatoires est non associée, suivant la règle de dilatance de Roscoe. Comme d’autres modèles de comportement de sols, l’écrouissage est positif en phase pré-pic et négatif en phase post-pic, qui correspond à l’effet de la dilatance; ces effets induisent le comportement de «liquéfaction» du sol. Le tenseur de déformation plastique résulte du cumul des contributions de divers mécanismes actifs. La déformation plastique volumique couple les mécanismes.
On décrit les équations du modèle, son paramétrage, puis son intégration numérique selon un schéma général de Newton implicite. Les utilisateurs ont la possibilité d’accéder à trois schémas d’intégration implicite pour le modèle de Hujeux : “SPECIFIQUE”, “SEMI_EXPLICITE” et “BASCULE_EXPLICITE”, les deux derniers étant des variantes du schéma “SPECIFIQUE” introduites pour fournir plus de souplesse à l’intégration locale du comportement lors des études en dynamique.
Intégration numérique de la relation de comportement#
Rappel du problème#
On emploie les notations suivantes : \({A}^{-},A,\Delta A\) pour une quantité évaluée à l’instant connu \({t}^{-}{\text{=t}}_{i-1}\) , à l’instant \({t}_{i}{\text{=t}}^{-}+\Delta t\) et son incrément \(\Delta t\) respectivement. Les équations sont discrétisées de manière implicite, c’est-à-dire exprimées en fonction des variables inconnues à l’instant \({t}_{i}{\text{=t}}^{-}+\Delta t\) .
Pour un incrément de chargement donné et un ensemble de variables données (champ initial de déplacement, contrainte et variable interne), on résout le système global discrétisé ([éq 2.2.2.2-1] de [bib1]) qui cherche à satisfaire les équations d’équilibre.
code_aster utilise une méthode de Newton [bib1], initiée par une étape de prédiction (tir d’Euler) qui fournit une estimation \(\Delta {\boldsymbol{U}}_{i}^{0}\) de l’incrément de déplacement, d’où la prédiction \('\boldsymbol{U}_{i}^{0}={\boldsymbol{U}}_{i-1}+\Delta {\boldsymbol{U}}_{i}^{0}\) , suivie d’itérations \(n=1,\mathrm{...}\) de correction, qui donnent des corrections \(\delta {\boldsymbol{U}}_{i}^{n+1}\) de l’incrément de déplacement, d’où la mise à jour \({\boldsymbol{U}}_{i}^{\text{n+}1}={\boldsymbol{U}}_{i}^{n}+\delta {\boldsymbol{U}}_{i}^{\text{n+}1}\) , qui doit converger vers la solution \({\boldsymbol{U}}_{i}={\boldsymbol{U}}^{+}\) .
L’option RIGI_MECA_TANG calcule l’opérateur tangent \({\boldsymbol{K}}_{i-1}\) , pour l’étape de prédiction (linéarisation des équations d’équilibre autour de l’équilibre à l’instant \({t}_{i-1}\) ), à partir de la relation vitesses de contrainte rapportées aux vitesses de déformation, à l’aide de l’opérateur déformation assemblé \({}^{T}\boldsymbol{Q}\) :
Remarque:
On peut décider, pour économiser du temps calcul, de ne pas réactualiser cette matrice, et de prendre la matrice de rigidité élastique, voir [U4.51.03] , mot-clé Newton, opérande PREDICTION, valeur ‘ELASTIQUE’, plutôt que ‘ TANGENTE’. Mais cela peut augmenter le nombre d’itérations de correction, voire rendre la convergence délicate, notamment à cause de la loi élastique non linéaire constituant la loi de Hujeux [éq 1.1.1-10].
Les options RAPH_MECA et FULL_MECA calculent les contraintes \({\boldsymbol{\sigma}}_{i}^{n}\) , les variables internes et les forces nodales \(\boldsymbol{R}({\boldsymbol{U}}_{i}^{n})={}^{T}\boldsymbol{Q}.{\boldsymbol{\sigma}}_{i}^{n}\) , à l’itération globale \(n\).
L’option FULL_MECA calcule de plus l’opérateur tangent \({\boldsymbol{K}}_{i}^{n}\) , à l’itération de correction globale \(n\) (sur demande via la commande STAT_NON_LINE, mot-clé facteur NEWTON, mot-clé MATRICE=”TANGENTE” et REAC_ITER), à partir des contraintes \({\boldsymbol{\sigma}}_{i}^{n}\) et de la matrice tangente du problème local discrétisé en temps :
La résolution de ces systèmes globaux d’équations d’équilibre nous donne des incréments \(\Delta {\boldsymbol{U}}_{i}^{n}\) , donc des incréments de déformation \(\Delta {\boldsymbol{\varepsilon}}_{i}^{n}\) . On cherche donc localement (en chaque point de Gauss) l’incrément des contraintes et de variables internes correspondant à \(\Delta \boldsymbol{\varepsilon}\) et qui satisfont la loi de comportement. L’évolution des contraintes et des variables internes est obtenue par résolution d’un système d’équations différentielles, par intégration locale, sous les conditions initiales au début du pas de temps, décrite au [§ 2.2].
Les équations de comportement élastoplastique multi-mécanisme présentées au [§ 1] ne permettent pas de démontrer par une analyse simple l’unicité de la solution du problème d’équilibre de structure, discrétisé en temps. Le schéma d’intégration doit avoir une influence certaine sur les solutions calculées.
Schéma général d’intégration locale#
Les schémas d’intégration locale choisis pour l’implantation de la loi de Hujeux dans code_aster sont entièrement implicites formulés sur le problème incrémental, pour des raisons de précision de calcul.
Ils utilisent une prédiction élastique puis des itérations de correction. Ils ont pour but de produire, un incrément de déformation \(\Delta \boldsymbol{\varepsilon}\) étant fourni, la valeur des contraintes et variables internes à l’instant \({t}_{i}{\text{=t}}^{-}+\Delta t\) . La [fig. 2.2-a] présente l’organisation générale commune des algorithmes d’intégration locale.
Remarque:
Les évolutions des grandeurs calculées du modèle de Hujeux à travers une intégration implicite sont réputées être sensibles à des perturbations des données de l’algorithme d’intégration (par exemple selon le choix de la plate-forme de calcul ou du choix de discrétisation spatiale…). C’est pourquoi diverses dispositions ont été prises afin de limiter cette sensibilité : voir ci-après la description des tolérances introduites sur les incréments. On doit noter qu’à cause du caractère multi-mécanisme avec variables mémoratrices du modèle de Hujeux, l’activation effective des mécanismes n’étant pas un processus différentiable, la formulation en vitesse ne constitue pas la limite de la formulation incrémentale en temps.
code_aster ne propose pas d‘intégration explicite pour ce modèle, cf. [R5.03.14], bien que cette méthode soit utilisée dans la littérature, cf. [bib3]. Les algorithmes d’intégration implicite disponibles sont au nombre de 3 et se distinguent selon les termes suivants :
ALGO_INTE = “SPECIFIQUE”, il s’agit du 1er schéma d’intégration implicite développé dans code_aster. Il est le schéma d’intégration le plus robuste, car il réalise pour l’utilisateur les subdivisions locales nécessaires pour assurer au maximum une convergence du schéma d’intégration.
ALGO_INTE=”SEMI_EXPLICITE”, il s’agit de s’appuyer sur la stratégie d’intégration du schéma “SPECIFIQUE” mais d’accepter, pendant les itérations du schéma de Newton global, des points de Gauss n’ayant pas ateint le critère de convergence local imposé par l’utilisateur. Néanmoins, la convergence globale n’est accepté que lorsque l’intégration locale est assurée pour l’ensemble des points de Gauss.
ALGO_INTE= »BASCULE_EXPLICIT », il s’agit dans ce cas d’accepter la non convergence locale du comportement si convergence du Newton global. Cette stratégie a été implémenté à des fins d’accélération de calcul en dynamique mais ne s’avère pas payante pour des modélisation fortement non linéaires.
Le processus de subdivision locale du pas de temps est disponible afin de réduire significativement les situations possibles de non-convergence, et d’économiser des itérations de correction. Il consiste à résoudre le comportement local sur des subdivisions linéaires de l’incrément de déformation \(\Delta \boldsymbol{\varepsilon}\) proposé par la résolution globale. La subdivision locale est activable dans l’opérateur de mécanique non linéaire de code_aster (mot-clé facteur COMPORTEMENT et mot-clé ITER_INTE_PAS, cf. [U4.51.11].
Phase de prédiction#
On commence par établir la prédiction élastique \({{\sigma}_{ij}^{+}}_{é}\) (tir élastique tangent approché – Euler):
Le tenseur \({C}_{\text{ijrs}}({\stackrel{ˆ}{p}}_{é})\) , défini par [éq 1.1.1-11], est obtenu à partir d’une expression discrétisée par un schéma approché de la pression de confinement:
Ce choix permet de limiter le risque d’avoir un prédicteur de pression de confinement produisant des contraintes trop éloignées des valeurs espérées, du fait de la forme non linéaire exponentielle de la loi élastique [éq 1.1.1-11].
On commence par une première étape de «reprédiction». Un re-découpage de l’incrément de déformation \(\Delta \boldsymbol{\varepsilon}\) – supplémentaire de celui éventuellement demandé par l’utilisateur via les mots-clés de l’opérateur de mécanique non linéaire de code_aster – est introduit de manière à borner les incréments suivants, associés à la prédiction \({\boldsymbol{\sigma}}_{é}\) , via une simple homothétie:
La finalité est d’éviter de possibles incursions au cours des itérations de correction vers des trajets de chargements trop éloignés, pouvant conduire à diverger (en sortant du domaine de convergence de l’algorithme de Newton).
Mécanismes plastiques potentiellement activés \({M}^{\text{pot}}\)#
Les mécanismes élastoplastiques \({M}_{k}^{K}\) , pour \(\text{k=}1,...,4\) et \(K=m,c\) , sont activés ou inactivés selon le type de chargement suivi. Une famille de mécanismes potentiellement actifs est ainsi déterminée au préalable au début du pas courant à partir de l’état atteint précédemment, et sert à calculer l’évolution réellement suivie. Le choix initial des mécanismes potentiellement activés au démarrage du pas de temps considéré résulte de l’état mécanique convergé atteint à l’instant \({t}^{-}\) : l’état d’activation des différents mécanismes (monotones ou cycliques) à l’instant \({t}^{-}\) et valeurs des variables \({\boldsymbol{\sigma}}^{-}\) et \({\boldsymbol{\alpha}}^{-}\) .
On note par \({M}^{\text{pot}}\) l’ensemble des mécanismes plastiques potentiellement activés cf. [fig. 2.2-b], et par \({M}^{\mathrm{act}}\) l’ensemble des mécanismes plastiques réellement activés lors de l’intégration, cf. [fig. 2.2-c]. Cette définition permet de réduire au strict minimum la taille du système d’équations non linéaires locales à résoudre à chaque pas. On utilise la prédiction élastique [éq 2.2.1-1] pour établir une ré-estimation de \({M}^{\text{pot}}\) à l’aide de l’évolution \(\dot{{\mathrm{f}}_{k}^{K}}\) pour \(\text{k=}1,...,4\) et \(K=m,c\) , des seuils des mécanismes et définir ainsi les situations de chargement en régime plastique ou de déchargement.
On initialise donc \({M}^{\text{pot}}:={M}^{\text{act}}({t}^{-})\) , et on enregistre les variables d’histoire des mécanismes cycliques «pères» \({\alpha}_{H}^{-}\) , à partir de la situation à l’instant précédent \({t}^{-}\) . On calcule pour tous les mécanismes \({M}_{k}^{K}\) deux grandeurs aidant à la décision :
S’il s’agit d’un mécanisme cyclique \({M}_{k}^{c}\) alors, on calcule aussi \({c}_{k}^{m}\) pour le mécanisme monotone \({M}_{k}^{m}\) associé. En cas de difficultés, on procède à un re-découpage du pas de temps \(\Delta t\) . Si la valeur de \({c}_{k}^{m}\ge 0\) , et si les deux surfaces de charge de \({M}_{k}^{c}\) et \({M}_{k}^{m}\) sont «proches», alors \({M}^{\text{pot}}:={M}^{\text{pot}}\cup \left\lbrace {M}_{k}^{m}\right\rbrace \left\lbrace {M}_{k}^{c}\right\rbrace\) et les propriétés d’écrouissage de \({M}_{k}^{c}\) sont réinitialisées à \({r}_{\text{éla}}^{c}\) . Cette proximité est établie en calculant la différence entre le facteur de mobilisation atteint pour le seuil monotone \({r}_{k}^{m}\) et la valeur du facteur de mobilisation à l’état convergé précédent \({r}_{k}^{c}\) ,
Puis on passe en revue tous les mécanismes possibles \({M}_{k}^{K}\) en testant \({c}_{k}^{K}\) et \({g}_{k}^{K}\) . On identifie ainsi les situations:
\({M}_{k}^{m}\in {M}_{\text{pot}}\) . Alors, si \({c}_{k}^{m}\ge 0\) , \({M}_{k}^{m}\) reste dans \({M}^{\text{pot}}\) ; sinon, on mémorise les variables d’histoire correspondantes \({\alpha}_{H}^{+}\) et si \({g}_{k}^{c}<0\) pour le mécanisme cyclique \({M}_{k}^{c}\) associé à \({M}_{k}^{m}\) , on désactive \({M}_{k}^{c}\) et \({M}_{k}^{m}\) ; mais si \({c}_{k}^{m}\ge 0\) on désactive \({M}_{k}^{m}\) au profit de \({M}_{k}^{c}\) mis dans \({M}^{\text{pot}}\) .
\({M}_{k}^{m}\notin {M}_{\text{pot}}\) . Alors, si \({c}_{k}^{m}\ge 0\) et \({g}_{k}^{m}\ge 0\) , alors \({M}^{\text{pot}}:{\text{=M}}^{\text{pot}}\cup \left\lbrace {M}_{k}^{m}\right\rbrace\) , mais si \({g}_{k}^{m}<0\) \({M}_{k}^{m}\) reste hors de \({M}^{\text{pot}}\) .
\({M}_{k}^{c}\in {M}_{\text{pot}}\) . Alors, si \({c}_{k}^{c}\ge 0\) , \({M}_{k}^{c}\) reste dans \({M}^{\text{pot}}\) ; sinon, on mémorise les variables d’histoire correspondantes \({\alpha}_{H}^{+}\) et si \({g}_{k}^{c}\ge 0\) , on maintient \({M}_{k}^{c}\) dans \({M}^{\text{pot}}\) , mais si \({g}_{k}^{c}<0\) , \({M}_{k}^{c}\) est retiré de \({M}^{\text{pot}}\) .
\({M}_{k}^{c}\notin {M}_{\text{pot}}\) . Alors, si \({c}_{k}^{c}\ge 0\) et \({g}_{k}^{c}\ge 0\) , alors \({M}^{\text{pot}}:{\text{=M}}^{\text{pot}}\cup \left\lbrace {M}_{k}^{c}\right\rbrace\) et on gère la potentialité de micro-décharges pour les seuls mécanismes déviatoires : on teste le seuil du mécanisme «père» de \({M}_{k}^{c}\) : si \({g}_{k}^{{c}^{'}}\ge 0\) alors le mécanisme «fils» \({M}_{k}^{c}\) est remplacé par ce mécanisme «père» \({M}_{k}^{{c}^{'}}\) avec mise à jour des variables d’histoire \({\alpha}_{H}^{+}\) . Si \({c}_{k}^{c}<0\) , on mémorise les variables d’histoire correspondantes et si \({r}_{k}^{c}\ge {r}_{\text{éla}}^{\text{dc}}\) on crée un mécanisme «fils» \({M}_{k}^{{c}^{'}}\) , dont la surface de charge est tangente à celle du mécanisme «père», qui est mis dans \({M}^{\text{pot}}\) si son seuil vérifie \({g}_{k}^{{c}^{'}}\ge 0\) , sinon, \({M}_{k}^{{c}^{'}}\) est désactivé.
Bien entendu, si \({M}_{k}^{K}\in {M}^{\text{pot}}\) pour \(K=m\) ou \(c\) , alors le mécanisme «frère» pour \(K=m\) ou \(c\) n’est pas dans \({M}^{\text{pot}}\) . Si au terme de l’analyse des diverses situations ci-dessus, aucun mécanisme ne subsiste dans \({M}^{\text{pot}}\) alors le régime sera élastique, voir cf. [fig. 2.2-b].
Le régime élastique advient lorsque aucun mécanisme n’est activé.
On impose à la prédiction \({\tilde{p}}_{{k}_{é}}\) de chaque mécanisme d’être strictement négatif pour pouvoir respecter les domaines de définition des diverses dérivées qui interviennent dans les équations ci-dessous: ainsi on obtient la prédiction en contrainte \({\tilde{\sigma}}_{\mathrm{ijé}}^{+}={\sigma}_{ij}^{-}+{C}_{\text{ijrs}}({\stackrel{ˆ}{p}}_{é}).\eta .\Delta {\varepsilon}_{\text{rs}}\) où \(\stackrel{ˆ}{{p}_{é}}\) est défini en [éq 2.2.1-2] et avec \(\eta\) tel que \({\tilde{p}}_{{k}_{é}}\le {10}^{-6}.{P}_{\text{réf}}\) . On note par la suite \({\tilde{Y}}_{é}=({\tilde{\sigma}}_{\mathrm{ijé}}^{+},{\varepsilon}_{v}^{{p}^{-}},{r}_{k}^{{K}^{-}},0)\) .
On cherche ensuite une solution d’essai explicite en calculant un candidat \({(\Delta {\lambda}_{k}^{K})}_{0}\) à l’aide des mécanismes potentiellement actifs \({M}_{k}^{K}\) , en résolvant les équations linéarisées autour de l’état sans évolution plastique \({\tilde{\sigma}}_{é}\) :
avec :
les dérivées étant données par [éq 7-13] à [éq 7-24]. D’où ensuite la solution d’essai \(\Delta {\boldsymbol{Y}}_{0}\) :
De plus, on impose la condition \({(\Delta {\lambda}_{k}^{K})}_{0}\ge 0\) ; les facteurs de mobilisation sont aussi bornés sur la solution d’essai: \({({r}_{k}^{m})}_{0}\le 1\) et \({({r}_{k}^{c})}_{0}\le {({r}_{k}^{m})}_{0}\) , de telle manière que la solution d’essai respecte la vraisemblance de l’évolution des variables internes. Sur la base de ces critères, on introduit également une restriction sur l’évolution de la déformation volumique plastique également, \(\lvert{(\Delta {\varepsilon}_{v}^{p})}_{0}\rvert\le 10\text{%}\) .
Remarque:
L’expression de la pression critique fonction de la déformation plastique volumique faisant intervenir une exponentielle [éq 1.1.1-13], on prévient le risque de valeurs trop élevées débordant du domaine de convergence en enclenchant un découpage de l’incrément de déformation \(\Delta \boldsymbol{\varepsilon}\) donné si nécessaire.
Phase de correction : équations non linéaires à résoudre#
Cette étape consiste à résoudre le système d’équations locales non linéaires établi sur la base des mécanismes potentiellement activés \({M}^{\text{pot}}\) . Après convergence, on réévalue l’ensemble des mécanismes réellement activés \({M}^{\text{act}}\) , cf. [§ 2.2.3.1], et s’il y a une différence avec \({M}^{\text{pot}}\) , on reprend la résolution locale non linéaire avec \({M}^{\text{act}}\) .
Itérations de correction de Newton#
Les itérations de correction de Newton consistent à résoudre les équations suivantes, pour \(\Delta \boldsymbol{\varepsilon}\) donné, voir [éq 6-1]:
l’équation d’état incrémentale [éq 1.1.1-9], notée \({\text{LE}}_{ij}\) (6 équations scalaires),
l’évolution incrémentale de la déformation plastique volumique \(\Delta {\varepsilon}_{v}^{p}\) [éq 1.1.5-1]: \(\text{LEVP}\) ,
l’évolution incrémentale des facteurs de mobilisation \(\Delta {r}_{k}^{K}\) [éq 1.1.2-5], [éq 1.1.2-15], [éq 1.1.3-3], et [éq 1.1.3-11], notés \(\text{LR.1}\) et \(\text{LR.2}\) (1 pour les mécanismes déviatoires, 2 pour les mécanismes de consolidation; autant d’équations que de mécanismes actifs),
les critères des divers mécanismes potentiellement actifs \({\mathrm{f}}_{k}^{K}=0\) [éq 1.1.2-1], [éq 1.1.2-9], [éq 1.1.3-1], et [éq 1.1.3-7], notés \(\text{LF.1}\) et \(\text{LF.2}\) (1 pour les mécanismes déviatoires, 2 pour les mécanismes de consolidation; autant d’équations que de mécanismes actifs).
Elles constituent un système carré \(\mathrm{\boldsymbol{R}}\left(\Delta \mathrm{\boldsymbol{Y}}\right)=0\) , où les inconnues sont \(\Delta \boldsymbol{Y}=(\Delta {\sigma}_{ij},\Delta {\varepsilon}_{v}^{p},\Delta {r}_{k}^{K},\Delta (\Delta {\lambda}_{k}^{K}))\) , qui couplent les équations (au nombre de 15 au maximum et 9 au minimum). On résout de manière implicite le système \(\boldsymbol{R}(\Delta \boldsymbol{Y})=0\) par une méthode de Newton, pour \(K=m,c\) et \(k=1,\mathrm{...},4\) .
A l’itération \(j\) de la boucle de correction locale de Newton, on résout l’équation matricielle:
où la matrice tangente \(\frac{d\mathrm{\boldsymbol{R}}}{d\left(\Delta \mathrm{\boldsymbol{Y}}\right)}{\rvert}_{\Delta {\mathrm{\boldsymbol{Y}}}_{j}}\) , non symétrique, est calculée de la manière présentée au [§ 6], selon [éq 6-2], avec une mise à jour du tenseur d’élasticité avec la valeur actualisée de la pression de confinement \(\frac{1}{3}\text{tr}({\boldsymbol{\sigma}}^{-}+\Delta {\boldsymbol{\sigma}}_{j})\) à l’itération précédente, selon [éq 1.1.1-11]:
On a mis au préalable les diverses lignes de ce système «à l’échelle», en divisant les équations d’état \({\text{LE}}_{ij}\) et les seuils \(\mathrm{LF}\) par le module de Young initial \({E}_{0}\) , les lois d’évolution \(\mathrm{LR}\) par \({P}_{\text{réf}}/{E}_{0}\) et les contraintes ainsi que les facteurs de mobilisation \({r}_{k}^{K}\) dans le vecteur inconnu \(\Delta {\boldsymbol{Y}}_{j}\) sont mises à l’échelle des déformations via le module de Young initial \({E}_{0}\) et la pression \({P}_{\mathrm{réf}}\) . Ce choix permet d’avoir des équations du même ordre de grandeur donc d’assurer une convergence plus «uniforme» sur l’ensemble du système [éq 2.2.2-3].
La matrice jacobienne présentée ci-dessus est évaluée analytiquement pour tous les schémas d’intégration du modèle.
Remarque:
Le déterminant de la matrice tangente \(\frac{d\boldsymbol{R}}{d(\Delta \boldsymbol{Y})}{\rvert}_{\Delta {\boldsymbol{Y}}_{j}}\) est a priori positif en phase d’écrouissage positif. Cependant, il pourrait devenir négatif; aucun traitement particulier n’est prévu.
La convergence est réputée acquise dès lors que (\(\mathrm{tol}\) est donnée par le mot-clé RESI_INTE_RELA du mot-clé facteur COMPORTEMENT de la commande STAT_NON_LINE, cf. [U4.51.11]):
cette norme \(\parallel .\parallel\) sur le résidu utilisant aussi la mise à l’échelle des divers termes intervenant.
Remarque:
On recommande fortement d’utiliser une valeur de tolérance inférieure ou égale à \({10}^{\text{-7}}\) (mot-clé RESI_INTE_RELA ), au moins dans les étapes les plus «délicates» du trajet de chargement.
Afin d’éviter de possibles incursions au cours des itérations de correction vers des trajets de chargements trop éloignés, pouvant conduire à diverger (en sortant du domaine de convergence locale de l’algorithme de Newton), l’algorithme d’intégration de Hujeux dit “SPECIFIQUE” impose aussi un re-découpage de l’incrément de déformation \(\Delta \boldsymbol{\varepsilon}\) si l’étape de correction aboutit à produire des incréments de facteurs de mobilisation \(\Delta {r}_{k}^{K}/{r}_{k}^{K}\ge 10\) . De plus, on impose que les valeurs prédites à ce stade des facteurs de mobilisation déviatoires vérifient: \({r}_{k}^{{m}^{+}}\le 1\) et \({r}_{k}^{{c}^{+}}\le {r}_{k}^{{m}^{+}}\) pour les mécanismes actifs. De plus, on enregistre les cas où les multiplicateurs plastiques sont négatifs. Si \(-{\eta}_{\text{tolé}}\le {\mathrm{\Delta \lambda }}_{k}^{K}/\underset{K=1,..,4}{\text{Max}}({\mathrm{\Delta \lambda }}_{k}^{K})\le 0\) , où \({\eta}_{\text{tolé}}\) désigne la valeur de RESI_INTE_RELA fourni par l’utilisateur, alors on impose: \({\mathrm{\Delta \lambda }}_{k}^{K}=0\) . Cela permet de limiter les effets «néfastes» des mécanismes en chargement «neutre», c’est-à-dire \(\Delta {\lambda}_{k}^{K}\approx 0\) .
La dernière tâche de gestion des mécanismes décrite ci-après intervient après un échec de l’algorithme de Newton local (§ [ 2.2.2]) ou bien si des critères, précisés ci-dessous, sont violés. Toute modification apportée au domaine de mécanismes potentiellement actifs \({M}^{\text{pot}}\) conduira ensuite à reprendre avec le nouvel \({M}^{\text{pot}}\) la résolution du système d’équations non linéaire. Ce processus de mise à jour est toutefois limité à 5 tentatives de reconstruction, faute de quoi on déclare une non-convergence de l’intégration locale.
Les premiers cas traités s’intéressent après chaque itération de l’algorithme de Newton local au risque de chevauchement de deux surfaces de charge déviatoires cyclique et monotone ou cyclique et cyclique dans l’espace des contraintes. Ces problèmes sont évités via des tests de proximité entre la position de l’état de contrainte dans le plan déviatoire considéré et la surface de charge considérée.
Dans le cas particulier de deux surfaces de charge déviatoires cyclique et monotone, si le critère de proximité est vérifié, le mécanisme monotone \({M}_{k}^{m}\) est ajouté directement au domaine des mécanismes \({M}^{\text{pot}}\) tout en retirant le mécanisme déviatoire cyclique \({M}_{k}^{c}\) considéré actif jusqu’à présent dans \({M}^{\text{pot}}\) .
Pour les cas de proximité de deux surfaces de charge de mécanismes déviatoires cycliques «père» et «fils», le critère porte sur la proximité de l’état de contraintes dans le plan déviatoire considéré par rapport à la valeur enregistrée \({\boldsymbol{X}}_{(k)}^{H}\) comme point de tangence entre les surfaces de charge «père» et «fils». Dès lors que le critère de proximité est atteint, le mécanisme déviatoire «fils» \({M}_{k}^{c}\) est retiré du domaine potentiel de mécanismes actifs \({M}^{\text{pot}}\) .
Après ces tests de proximité de surfaces de charge déviatoires, il est nécessaire de traiter les états de contraintes de traction causant très souvent la non-convergence du système non linéaire local à résoudre.
Si un état de traction est détecté après échec de la méthode de résolution, pour tout multiplicateur plastique associé à un mécanisme déviatoire considéré dans \({M}^{\text{pot}}\) négatif après le tir explicite d’Euler, \({(\Delta {\lambda}_{k}^{K})}_{0}\) , ce mécanisme déviatoire sera retiré de \({M}^{\text{pot}}\) . Cette solution pour relancer la résolution du système non linéaire local n’impose aucune restriction sur la valeur du prédicteur de contraintes. Ce prédicteur est établi en considérant un incrément de déformation purement élastique.
Si tous les multiplicateurs plastiques sont positifs, i-e \({(\Delta {\lambda}_{k}^{K})}_{0}\ge 0\) , alors tous les mécanismes déviatoires sont retirés de \({M}^{\text{pot}}\) : les mécanismes de consolidation intégrés seuls permettent de préparer un nouveau prédicteur pour l’étape d’intégration locale suivante après actualisation de \({M}^{\text{act}}\) .
Pour les cas différents des deux grandes classes décrites ci-dessus (proximité et traction), si la valeur maximum du vecteur résidu \(\boldsymbol{R}(\Delta {\boldsymbol{Y}}_{j})\) est porté par une loi d’évolution d’un facteur de mobilisation de mécanisme déviatoire cyclique \({M}_{k}^{c}\) , (soit l’équation \(\mathrm{LR.1}\) de [éq 6-1]), alors le mécanisme en question est retiré de \({M}^{\text{pot}}\) .
Le cas suivant s’intéresse uniquement aux mécanismes déviatoires cycliques jamais écrouis précédemment. Ces mécanismes sont alors retirés de \({M}^{\text{pot}}\) s’il y a eu échec lors de la résolution du système d’équations non linéaires.
Les derniers cas portent sur la valeur des multiplicateurs plastiques après le tir d’Euler, \({(\Delta {\lambda}_{k}^{K})}_{0}\) .
Les mécanismes associés à des valeurs négatives sont retirés de \({M}^{\text{pot}}\) .
Si tous les multiplicateurs plastiques sont positifs, on retire le mécanisme présentant la valeur la plus faible pour \({(\Delta {\lambda}_{k}^{K})}_{0}\) après le tir d’Euler.
Si aucun des traitements spécifiques décrits ci-dessus n’a été sollicité, la non-convergence du modèle à l’échelle locale est assumée et cette information est renvoyée à l’algorithme d’intégration globale de l’équilibre non linéaire.
Phase de mise à jour#
Après résolution, on met à jour le vecteur solution:
ce qui achève l’étape d’intégration locale, et on cale \({M}^{\text{act}}={M}^{\text{pot}}\) pour commencer.
Mécanismes plastiques réellement activés \({M}^{\mathrm{act}}\)#
On décrit ci-après la gestion de l’activation réelle des mécanismes plastiques. Après convergence du système d’équations locales non linéaire, on est amené à vérifier l’ensemble des mécanismes réellement activés \({M}^{\mathrm{act}}\) selon la procédure suivante, en partant de la prédiction \({M}^{\text{pot}}\) , cf. [fig. 2.2-c]. Cet arbre de décision est un élément constitutif du modèle de Hujeux, puisqu’il décrit l’enchaînement des mécanismes cycliques et les enregistrements des variables mémoratrices discrètes et les éventuelles restaurations des variables d’écrouissage, cf. aussi [fig. 1.1] et [fig. 1.1 -c].
Remarque:
Avec une intégration explicite, cf. [bib5], on procèderait plutôt avec des re-découpages (sous-incrémention) pour gérer les transitions de mécanismes.
On passe donc en revue tous les mécanismes, toute modification de \({M}^{\mathrm{act}}\) conduisant à rependre une nouvelle résolution locale non linéaire, cf. [§ 2.2.2.1]:
cas d’un mécanisme déviatorique monotone actif \({M}_{k}^{m}\) avec un multiplicateur plastique \(\Delta {\lambda}_{k}^{K}\) négatif:
si \({r}_{k}{\text{=r}}_{\text{éla}}^{d}\) , alors ce mécanisme \({M}_{k}^{m}\) est désactivé: \({M}_{\text{act}}:={M}_{\text{act}}\setminus \lbrace {M}_{k}^{m}\rbrace\) ;
sinon (remise en cause de l’estimation faite dans \({M}_{\text{pot}}\) ):
s’il existait un mécanisme cyclique \({M}_{k}^{c}\) préalablement à l’étape de \({M}_{\text{pot}}\) , on garde les variables de mémoire de ce mécanisme, on le déclare créé et on enlève le mécanisme monotone \({M}_{k}^{m}\) : \({M}_{\text{act}}:={M}_{\text{act}}\setminus \lbrace {M}_{k}^{m}\rbrace\) ,
s’il n’existait pas de mécanisme cyclique, on mémorise les variables de mémoire et on crée un mécanisme cyclique \({M}_{k}^{c}\) et on enlève \({M}_{k}^{m}\) : \({M}_{\text{act}}:{\text{=M}}_{\text{act}}\setminus \lbrace {M}_{k}^{m}\rbrace\) .
cas d’un mécanisme isotrope monotone actif \({M}_{4}^{m}\) avec un multiplicateur plastique négatif:
si \({r}_{4}{\text{=r}}_{\text{éla}}^{s}\) , alors ce mécanisme \({M}_{4}^{m}\) est désactivé: \({M}_{\text{act}}:={M}_{\text{act}}\setminus \lbrace {M}_{4}^{m}\rbrace\) ,
sinon :
s’il existait un mécanisme cyclique \({M}_{4}^{c}\) préalablement à l’étape de \({M}_{\text{pot}}\) , on garde lesvariables de mémoire de ce mécanisme, on le déclare créé et on enlève \({M}_{4}^{m}\) : \({M}_{\text{act}}:={M}_{\text{act}}\setminus \lbrace {M}_{4}^{m}\rbrace\) ,
s’il n’existe pas de mécanisme cyclique, on ne mémorise les variables de mémoire et on ne crée un mécanisme cyclique \({M}_{4}^{c}\) que si \(\dot{\overline{\left(\frac{\mid p\mid }{d\mid {p}_{c}\mid }\right)}}<0\) (on rappelle que \({p}_{c}\) dépend de \({\varepsilon}_{v}^{p}\) , qui lui-même varie en fonction de tous les mécanismes), et dans ce cas, on enlève \({M}_{4}^{m}\) .
cas d’un mécanisme déviatorique monotone \({M}_{k}^{m}\) créé mais inactif : si le seuil \({f}_{k}^{m}({t}^{+})\ge 0\) alors ce mécanisme devient actif: \({M}_{\text{act}}:{\text{=M}}_{\text{act}}\cup \left\lbrace {M}_{k}^{m}\right\rbrace\) , et les variables d’écrouissage du mécanisme cyclique \({M}_{k}^{c}\) associé prennent les valeurs initiales vierges.
cas d’un mécanisme isotrope monotone \({M}_{4}^{m}\) créé mais inactif: si le seuil \({f}_{4}^{m}({t}^{+})\ge 0\) alors ce mécanisme devient actif: \({M}_{\text{act}}:{\text{=M}}_{\text{act}}\cup \left\lbrace {M}_{4}^{m}\right\rbrace\) ; si le seuil \({f}_{4}^{m}({t}^{+})<0\) alors on ne mémorise les variables de mémoire et on ne crée un mécanisme cyclique et on n’enlève \({M}_{4}^{m}\) que si \(\dot{\overline{\left(\frac{\mid p\mid }{d\mid {p}_{c}\mid }\right)}}<0\) .
cas d’un mécanisme cyclique \({M}_{k}^{c}\) créé auparavant dans \({M}^{\text{pot}}\) à partir d’un mécanisme monotone:
si le multiplicateur plastique est négatif alors le mécanisme est désactivé: \({M}_{\text{act}}:{\text{=M}}_{\text{act}}\setminus \lbrace {M}_{k}^{c}\rbrace\) ; sinon, on passe à :
si le seuil du mécanisme monotone \({M}_{k}^{m}\) associé \({f}_{k}^{m}({t}^{+})\ge 0\) alors :
pour un mécanisme déviatorique, on passe en monotone: \({M}_{\text{act}}:={M}_{\text{act}}\setminus \lbrace {M}_{k}^{c}\rbrace \cup \left\lbrace {M}_{k}^{m}\right\rbrace\) et les variables d’écrouissage du mécanisme cyclique \({M}_{k}^{c}\) reprennent les valeurs initiales vierges,
pour un mécanisme isotrope: si la variable d’histoire de direction change de signe entre à \({t}^{-}\) et \({t}^{+}\) suite à \({M}^{\text{pot}}\) , alors on reste en cyclique \({M}_{4}^{c}\) ; sinon on passe en monotone: \({M}_{\text{act}}:={M}_{\text{act}}\setminus \lbrace {M}_{4}^{c}\rbrace \cup \left\lbrace {M}_{4}^{m}\right\rbrace\) et les variables d’écrouissage du mécanisme cyclique \({M}_{k}^{c}\) reprennent les valeurs initiales vierges.
s’il s’agit d’un mécanisme déviatorique \({M}_{k}^{c}\) déclaré actif à \({t}^{+}\) , et si \(\left(\parallel {\boldsymbol{X}}_{(k)}^{H}+\frac{{\boldsymbol{S}}_{(k)H}^{c}}{{\parallel {\boldsymbol{S}}_{(k)H}^{c}\parallel }_{\text{VM}}^{\text{2D}}}.({r}_{k}^{c}+{r}_{\text{éla}}^{\text{dc}})\parallel+ {r}_{k}^{c}-{r}_{k}^{m}\right)/{r}_{k}^{m}\ge 0\) (calculé à \({t}^{+}\) ), alors les variables mémoratrices de ce mécanisme \({M}_{k}^{c}\) sont mises à jour, et les variables d’écrouissage du mécanisme déviatorique monotone «père» sont réinitialisées à l’état vierge. Cela empêche le franchissement du seuil du mécanisme monotone par le seuil du mécanisme cyclique.
si le mécanisme cyclique \({M}_{k}^{c}\) est déclaré créé mais inactif dans \({M}^{\text{pot}}\) alors :
s’il s’agit d’un mécanisme déviatorique \({M}_{k}^{c}\) , alors :
si le seuil \({f}_{k}^{c}({t}^{+})\ge 0\) alors \({M}_{k}^{c}\) est déclaré actif et :
si \({\boldsymbol{X}}_{(k)\text{père}}^{H}\ne 0\) (s’il existe un mécanisme «père»), alors si à \({t}^{-}\) \(\frac{{\boldsymbol{S}}_{(k)\text{pèreH}}^{c}}{{\parallel {\boldsymbol{S}}_{(k)\text{pèreH}}^{c}\parallel }_{\text{VM}}^{\text{2D}}}.\left({\boldsymbol{X}}_{(k)\text{père}}^{H}-\frac{{\boldsymbol{S}}_{(k)}^{c}}{{p}_{k}.F({p}_{k}{\varepsilon}_{v}^{p})}\right)<0\) et le seuil vérifie \({\mathrm{f}}_{k}^{c}({t}^{-})\ge 0\) ce mécanisme \({M}_{k}^{c}\) reprend les variables mémoratrices du mécanisme «père» –qui redevient actif– ainsi que sa valeur de \({r}_{k}^{c}\) , puis le mécanisme «père» est réinitialisé à l’état vierge;
s’il n’existe pas de mécanisme «père» et si \(\frac{{\boldsymbol{S}}_{(k)\text{filsH}}^{c}}{{\parallel {\boldsymbol{S}}_{(k)\text{filsH}}^{c}\parallel }_{\text{VM}}^{\text{2D}}}\cdot\left({\boldsymbol{X}}_{(k)\text{fils}}^{H}-\frac{{\boldsymbol{S}}_{(k)}^{c}}{{p}_{k}.F({p}_{k}{\varepsilon}_{v}^{p})}\right)\ge 0\) à \({t}^{+}\) , alors les variables mémoratrices à \({t}^{+}\) de ce mécanisme \({M}_{k}^{c}\) sont stockées dans celles d’un mécanisme «père» dont on doit garder les valeurs, et les variables mémoratrices de ce mécanisme «fils» \({M}_{k}^{c}\) – qui devient actif – deviennent : \({\boldsymbol{X}}_{(k)\text{fils}}^{H}={\boldsymbol{X}}_{(k)\text{père}}^{H}+{\text{2r}}_{k}^{c}\frac{{\boldsymbol{S}}_{(k)\text{pèreH}}^{c}}{{\parallel {\boldsymbol{S}}_{(k)\text{pèreH}}^{c}\parallel }_{\text{VM}}^{\text{2D}}}\) et \(\frac{{{\boldsymbol{S}}_{(k)}^{c}}_{\mathrm{filsH}}}{{\parallel{{\boldsymbol{S}}_{(k)}^{c}}_{\mathrm{filsH}}\parallel}_{\mathrm{VM}}^{\mathrm{2D}}}=\frac{-{{\boldsymbol{S}}_{(k)}^{c}}_{\mathrm{pèreH}}}{{\parallel{{\boldsymbol{S}}_{(k)}^{c}}_{\mathrm{pèreH}}\parallel}_{\mathrm{VM}}^{\mathrm{2D}}}\) ;
si les variables mémoratrices estimées avant l’étape de \({M}^{\text{pot}}\) vérifient \({\boldsymbol{X}}_{(k)}^{H}\ne {\boldsymbol{X}}_{(k)}^{H}({t}^{-})\) et \({\boldsymbol{X}}_{(k)}^{H}({t}^{-})\ne 0\) , alors si à \({t}^{+}\) \(\frac{{\boldsymbol{S}}_{(k)H}^{c}}{{\parallel {\boldsymbol{S}}_{(k)H}^{c}\parallel }_{\text{VM}}^{\text{2D}}}\cdot\left({\boldsymbol{X}}_{(k)}^{H}-\frac{{\boldsymbol{S}}_{(k)}^{c}}{{p}_{k}.F({p}_{k},{\varepsilon}_{v}^{p})}\right)\ge 0\) on égalise les variables mémoratrices à \({t}^{+}\) avec ces valeurs estimées, et de plus si le seuil \({\mathrm{f}}_{k}^{c}({t}^{+})\ge 0\) , alors \({M}_{k}^{c}\) est activé et on égalise aussi les variables mémoratrices à \({t}^{-}\) et à \({t}^{+}\) , de même les valeurs de \({r}_{k}^{c}\) , et sinon, il est inutile de reprendre une résolution non linéaire locale;
si les variables mémoratrices estimées avant l’étape de \({M}^{\text{pot}}\) vérifient \({\boldsymbol{X}}_{(k)}^{H}\ne {\boldsymbol{X}}_{(k)}^{H}({t}^{-})\) et \({\boldsymbol{X}}_{(k)}^{H}({t}^{-})=0\) (le «père» est un mécanisme monotone), alors si \(\frac{{\boldsymbol{S}}_{(k)H}^{c}}{{\parallel {\boldsymbol{S}}_{(k)H}^{c}\parallel }_{\text{VM}}^{\text{2D}}}\cdot\left({\boldsymbol{X}}_{(k)}^{H}-\frac{{\boldsymbol{S}}_{(k)}^{c}}{{p}_{k}.F({p}_{k},{\varepsilon}_{v}^{p})}\right)\ge 0\) à \({t}^{+}\) les variables mémoratrices de ce mécanisme \({M}_{k}^{c}\) sont réinitialisées à l’état vierge et \({M}_{k}^{m}\) est activé;
si à \({t}^{-}\) \({r}_{k}^{c}\ne {r}_{\text{éla}}^{\text{dc}}\) , alors on enregistre les variables mémoratrices de ce mécanisme \({M}_{k}^{c}\) à l’instant \({t}^{+}\) (un nouveau «père» est créé) et on réinitialise à l’état vierge ce mécanisme \({r}_{k}^{\text{c+}}={r}_{\text{éla}}^{\text{dc}}\) , qui est créé inactif. Dans ce cas, il est inutile de reprendre une résolution non linéaire locale.
s’il s’agit d’un mécanisme isotrope \({M}_{4}^{c}\) et :
si \({\left(\frac{\mid p\mid }{d\mid {p}_{c}\mid }\right)}^{+}\le {\left(\frac{\mid p\mid }{d\mid {p}_{c}\mid }\right)}^{-}\) et \({p}_{H}^{c}/\mid {p}_{H}^{c}\mid =1\) à \({t}^{-}\) (réduction de la compression) : si le seuil \({\mathrm{f}}_{4}^{c}({t}^{+})\ge 0\) alors \({M}_{4}^{c}\) est déclaré actif sinon les variables mémoratrices de ce mécanisme isotrope à \({t}^{+}\) reprennent les valeurs avant l’étape de \({M}^{\text{pot}}\) . De plus, si \({p}_{H}^{c}=0\) à l’étape \({M}^{\text{pot}}\) , alors les mécanismes associés \({M}_{4}^{m}\) et \({M}_{4}^{c}\) sont inactivés;
si \({\left(\frac{\mid p\mid }{d\mid {p}_{c}\mid }\right)}^{+}>{\left(\frac{\mid p\mid }{d\mid {p}_{c}\mid }\right)}^{-}\) et \({p}_{H}^{c}/\mid {p}_{H}^{c}\mid =-1\) à \({t}^{-}\) (augmentation de la compression) : si le seuil \({\mathrm{f}}_{4}^{c}({t}^{+})\ge 0\) alors \({M}_{4}^{c}\) est déclaré actif sinon les variables mémoratrices de ce mécanisme isotrope à \({t}^{+}\) prennent les valeurs avant l’étape de \({M}^{\text{pot}}\) ;
si \({\left(\frac{\mid p\mid }{d\mid {p}_{c}\mid }\right)}^{+}>{\left(\frac{\mid p\mid }{d\mid {p}_{c}\mid }\right)}^{-}\) et \({p}_{H}^{c}/\mid {p}_{H}^{c}\mid =1\) à \({t}^{-}\) (réduction de la compression) :
si \({p}_{H}^{c}/\mid {p}_{H}^{c}\mid\) change de signe à cette étape, alors les variables mémoratrices de ce mécanisme isotrope cyclique à \({t}^{+}\) prennent les valeurs obtenues avant l’étape de \({M}^{\text{pot}}\) ; de plus s’il n’existe pas de \({p}_{H}^{c}\) à \({t}^{+}\) , alors les mécanismes \({M}_{4}^{m}\) et \({M}_{4}^{c}\) sont déclarés inactifs, et sinon, si le seuil \({f}_{4}^{c}({t}^{+})\ge 0\) alors \({M}_{4}^{c}\) – déclaré actif à \({t}^{-}\) – reste créé à \({t}^{+}\) et les variables mémoratrices à \({t}^{-}\) prennent les valeurs obtenues avant l’étape de \({M}^{\text{pot}}\) à \({t}^{-}\) ;
si \({p}_{H}^{c}/\mid {p}_{H}^{c}\mid\) ne change pas de signe à cette étape, alors on réactualise les variables mémoratrices de ce mécanisme isotrope (création d’un «fils») et si le seuil \({\mathrm{f}}_{4}^{c}\) calculé avec les contraintes à \({t}^{+}\) dépasse le critère initial (\({r}_{\text{éla}}^{c}\) ), alors \({M}_{4}^{c}\) est déclaré actif à \({t}^{-}\) et les variables mémoratrices à \({t}^{-}\) prennent les valeurs obtenues à \({t}^{+}\) ;
si \({\left(\frac{\mid p\mid }{d\mid {p}_{c}\mid }\right)}^{+}\le {\left(\frac{\mid p\mid }{d\mid {p}_{c}\mid }\right)}^{-}\) et \({p}_{H}^{c}/\mid {p}_{H}^{c}\mid =-1\) à \({t}^{-}\) (augmentation de la compression):
si \({p}_{H}^{c}/\mid {p}_{H}^{c}\mid\) change de signe à cette étape, alors les variables mémoratrices de ce mécanisme isotrope à \({t}^{+}\) de ce mécanisme prennent les valeurs obtenues lors de l’estimation de \({M}^{\text{pot}}\) , et si le seuil \({f}_{c}^{4}\) calculé avec les contraintes à \({t}^{+}\) , alors les variables mémoratrices à \({t}^{-}\) prennent les valeurs obtenues avant l’étape de \({M}^{\text{pot}}\) et \({M}_{4}^{c}\) est déclaré actif à \({t}^{-}\) et reste créé à \({t}^{+}\) ;
si \({p}_{H}^{c}/\mid {p}_{H}^{c}\mid\) ne change pas de signe depuis l’étape de \({M}^{\text{pot}}\) , alors on réactualise les variables mémoratrices de ce mécanisme isotrope (création d’un «fils») et si le seuil \({\mathrm{f}}_{4}^{c}\) calculé avec les contraintes à \({t}^{+}\) dépasse le critère initial (\({r}_{\text{éla}}^{c}\) ), alors \({M}_{4}^{c}\) est déclaré actif à \({t}^{-}\) et reste créé à \({t}^{+}\) , et les variables mémoratrices à \({t}^{-}\) prennent les valeurs obtenues à \({t}^{+}\) .
Remarque:
Dans le cas où \({M}^{\text{act}}\ne {M}^{\text{pot}}\) , on redémarre l’algorithme de Newton local en utilisant non pas le prédicteur élastique des contraintes \({\boldsymbol{\sigma}}^{-}+\boldsymbol{C}(\stackrel{ˆ}{{p}_{é}}).\Delta \boldsymbol{\varepsilon}\) , cf. [éq 2.2.1-1], mais le résultat obtenu à la fin de l’étape échouée précédente \({\boldsymbol{\sigma}}_{\text{échoué}}^{+}\) , cf. [éq 2.2.3-1], associé aux valeurs obtenues des variables internes \({\boldsymbol{\alpha}}^{-}\) . On espère ainsi éviter des bouclages intempestifs et en démarrant de manière plus «proche» de la solution cherchée, accélérer la convergence.
Calcul de la matrice de raideur tangente incrémentale#
On établit enfin (sur demande) la matrice tangente locale \({\boldsymbol{C}}^{T}\) du comportement incrémental, à l’itération globale courante, qui relie la variation de contrainte totale à la variation de déformation totale :
Pour cela, on exploite les conditions de «cohérence»: \({\dot{\mathrm{f}}}_{k}^{K}({\boldsymbol{\sigma}}^{+},{\varepsilon}_{v}^{\text{p+}},{r}_{k}^{\text{K+}})=0\) pour les mécanismes actifs, monotone ou cycliques \(\text{K=m,c}\) , les quantités intervenant dans ces expressions étant calculées à l’instant \({t}^{+}{\text{=t}}_{i}\) actuel. En combinant avec [éq 1.1.1-9], on tire:
cf. [éq 1.1.2-7], [éq 1.1.2-17], [éq 1.1.3-5], [éq 1.1.3-13], les diverses dérivées étant données par [éq 7-13] à [éq 7-24]. Le signe de \(\frac{\partial {\mathrm{f}}_{k}^{K}}{\partial \boldsymbol{\sigma} }{\mid}_{{\boldsymbol{Y}}^{+}}\cdot\boldsymbol{C}({\stackrel{ˆ}{p}}_{é}).\Delta \dot{\boldsymbol{\varepsilon}}\) est comparé à une tolérance de référence (R8PREM).
On en déduit l’expression de la matrice tangente incrémentale (option FULL_MECA), sur les mécanismes de \({M}^{\mathrm{act}}\) :
La matrice tangente incrémentale est ainsi établie analytiquement, cf. [§ 6], les diverses dérivées présentes dans [éq 2.2.2-2] étant calculés au [§ 6]. Le déterminant de la matrice tangente incrémentale \(\det({C}_{\text{ijrs}}^{T})\) , en option FULL_MECA, est stockée parmi les variables internes : VARI_33, cf. [§ 3.1].
Remarque:
L’évaluation de la matrice tangente incrémentale à l’aide d’une technique de perturbations (par différences finies, à titre de vérification) n’est pas possible à cause du caractère multi-mécanisme du modèle de Hujeux : l’activation effective des mécanismes n’étant pas un processus différentiable.
Fig. 218 Schéma d’intégration locale de la loi de comportement de Hujeux dans code_aster.#
Fig. 219 Algorithme d’évolution du domaine \({M}^{\text{pot}}\) des mécanismes potentiellement actifs.#
Fig. 220 Algorithme de définition du domaine \({M}^{\text{act}}\) des mécanismes activés.#
Opérateur tangent en vitesse : option RIGI_MECA_TANG#
Pour l’option RIGI_MECA_TANG, qui est utilisée lors de la prédiction globale, appelée à la première itération d’un nouvel incrément de charge \(\Delta \stackrel{ˆ}{\boldsymbol{L}}({t}_{i})\) , l’opérateur tangent global, noté \({\boldsymbol{K}}_{i-1}\) , est calculé à partir des résultats connus à l’instant \({t}^{-}{\text{=t}}_{i-1}\) [bib3].
L’opérateur tangent global est assemblé à partir des contributions de la matrice tangente en chaque point de Gauss, dite «en vitesses» :
Comme la loi de comportement est élastique non linéaire [éq 1.1.1-11], on construit un opérateur tangent élastique non linéaire.
Si à l’état précédent le tenseur des contraintes n’est à la frontière d’aucun seuil plastique, la prédiction élastique \({{\sigma}_{ij}^{+}}_{é}\) s’écrit selon [éq 2.2.1-1], à l’aide de la pression de confinement estimée \(\stackrel{ˆ}{{p}_{é}}\) [éq 2.2.1-2].
Dans ce cas, on construit la matrice \({\boldsymbol{K}}_{i-1}\) à l’aide du tenseur d’élasticité calculé pour la pression de confinement \(\stackrel{ˆ}{{p}_{é}}\) . Ce choix permet de limiter le risque d’avoir une matrice produisant des prédictions de contraintes trop éloignées des valeurs espérées, du fait de la forme non linéaire exponentielle de la loi élastique [éq 1.1.1-11].
Si à l’état précédent le tenseur des contraintes est sur la frontière d’un seuil plastique, on exploite les conditions de «cohérence» : \({\dot{\mathrm{f}}}_{k}^{K}({\boldsymbol{\sigma}}^{-},{\varepsilon}_{v}^{p-},{r}_{k}^{K-})=0\) , pour les mécanismes actifs, monotone ou cycliques \(\text{K=m,c}\) , les quantités intervenant dans ces expressions étant calculées à l’instant \({t}^{-}{\text{=t}}_{i-1}\) précédent. En combinant avec [éq 1.1.1-9], on obtient l’expression des multiplicateurs plastiques \(\dot{{\lambda}_{k}^{K}}\) :
cf. [éq 1.1.2-7], [éq 1.1.2-17], [éq 1.1.3-6], [éq 1.1.3-13], les diverses dérivées étant données en annexe par [éq 7-13] à [éq 7-24].
On en déduit l’expression de la matrice tangente «en vitesses», sur les mécanismes de \({M}^{\mathrm{act}}\) :
Remarque:
On doit noter que le caractère non associé des lois d’écoulement [éq 1.1.2-3] et [éq 1.1.2-14] des mécanismes déviatoriques fait perdre les symétries majeures dans la matrice tangente. De même, ces lois d’écoulement introduisent des couplages entre les composantes \({C}_{\mathrm{iijj}}^{\mathrm{elp}}\) et les composantes \({C}_{\mathrm{klkl}}^{\mathrm{elp}}\) pour \(k\ne l\) , ce qui est étroitement lié à la dilatance.
Remarque:
On doit noter aussi que, comme le tenseur d’élasticité \(\boldsymbol{C}(\stackrel{ˆ}{{p}_{é}})\) est isotrope, les situations de chargement où seules les directions diagonales sont sollicitées, alors les tenseurs \({({\Psi}_{\text{lt}})}_{(k)}^{K}\) ne comportent que des termes diagonaux, et donc le tenseur d’élasticité tangente \({\boldsymbol{C}}^{\mathrm{elp}}\) ne comporte pas de couplage entre termes diagonaux et termes extra-diagonaux.
Implantation dans Code_Aster#
Variables internes#
Un certain nombre de «variables internes» du modèle au sens de code_aster sont créées et stockées. Le tenseur des déformations plastiques \({\boldsymbol{\varepsilon}}^{p}\) n’est pas stocké car il est obtenu par calcul en post-traitement à partir des contraintes et des déformations totales. Les valeurs des variables internes aux points de Gauss (VARI_ELGA) sont pour toutes les modélisations:
VARI_1 |
\({r}_{1}^{m}\) (facteur de mobilisation du mécanisme déviatoire monotone, \(\text{k=}1\) ) |
VARI_2 |
\({r}_{2}^{m}\) (facteur de mobilisation du mécanisme déviatoire monotone, \(\text{k=}2\) ) |
VARI_3 |
\({r}_{3}^{m}\) (facteur de mobilisation du mécanisme déviatoire monotone, \(\text{k=}3\) ) |
VARI_4 |
\({r}_{4}^{m}\) (facteur de mobilisation du mécanisme de consolidation monotone) |
VARI_5 |
\({r}_{1}^{c}\) (facteur de mobilisation du mécanisme déviatoire cyclique, \(\text{k=}1\) ) |
VARI_6 |
\({r}_{2}^{c}\) (facteur de mobilisation du mécanisme déviatoire cyclique, \(\text{k=}2\) ) |
VARI_7 |
\({r}_{3}^{c}\) (facteur de mobilisation du mécanisme déviatoire cyclique, \(\text{k=}3\) ) |
VARI_8 |
\({r}_{4}^{c}\) (facteur de mobilisation du mécanisme de consolidation cyclique) |
VARI_9 à VARI_12 |
variables mémoratrices \({({X}_{(1)}^{H})}_{1111}\) , \(\sqrt{2}{({X}_{(1)}^{H})}_{1212}\) , \(-{({S}_{(1)H}^{c})}_{1111}/{\parallel {S}_{(1)H}^{c}\parallel }_{\text{VM}}^{\text{2D}}\) , \(-\sqrt{2}{({S}_{(1)H}^{c})}_{1212}/{\parallel {S}_{(1)H}^{c}\parallel }_{\text{VM}}^{\text{2D}}\) (normale entrante), pour le mécanisme déviatoire cyclique du plan \(\text{k=}1\) |
VARI_13 à VARI_16 |
idem pour le mécanisme déviatoire cyclique du plan \(k=2\) |
VARI_17 à VARI_20 |
idem pour le mécanisme déviatoire cyclique du plan \(k=3\) |
VARI_21 |
variable mémoratrice discontinue \({p}_{H}/(d.{P}_{\text{c0}}.{e}^{-\beta {\varepsilon}_{\text{vH}}^{p}})\) du mécanisme de consolidation |
VARI_22 |
variable mémoratrice discontinue \({p}_{H}^{c}/\mid {p}_{H}^{c}\mid\) de normale à la surface de charge du mécanisme de consolidation |
VARI_23 |
déformation plastique volumique \({\varepsilon}_{v}^{p}\) |
VARI_24 à VARI_27 |
indicateurs d’activation (1) ou non (0) des mécanismes monotones ou de passage au cyclique (-1) |
VARI_28 à VARI_31 |
indicateurs d’activation (1) ou non (0) des mécanismes cycliques |
VARI_32 |
densité normalisée pour le travail du second ordre \(\dot{\boldsymbol{\sigma}}.\dot{\boldsymbol{\varepsilon}}\) (obtenu par discrétisation: \(\Delta \boldsymbol{\sigma} .\Delta \boldsymbol{\varepsilon} /\mathrm{\parallel }\Delta \boldsymbol{\sigma} \mathrm{\parallel }.\mathrm{\parallel }\Delta \boldsymbol{\varepsilon} \mathrm{\parallel }\) ) |
VARI_33 |
Déterminant de la matrice tangente \(\frac{d\boldsymbol{\sigma} }{d\boldsymbol{\varepsilon} }\) |
VARI_34 |
Pour les schémas SPECIFIQUE et SEMI_EXPLICITE : indicateur d’état des mécanismes actifs. Pour le schéma BASCULE_EXPLICIT : indicateur d’erreur de la solution \(\sum_{i=1}^{\text{ITER_INTE_PAS}}\sum_{k\in \mathit{Actif}}{\mathrm{f}}_{k}^{m,c}\left({\boldsymbol{\sigma}}^{+},{{\varepsilon}}_{v}^{p},{r}_{k}^{+}\right)\) |
VARI_35 |
Pour le schéma BASCULE_EXPLICIT : indicateur d’erreur de la solution explicite \(\frac{{N}_{\text{échec}}}{\text{ITER_INTE_PAS}}\) (en %) |
VARI_36 à VARI_41 |
coordonnées vectorielles du point de tangence \({({X}_{(1)}^{H})}_{1111}\) , \(\sqrt{2}{({X}_{(1)}^{H})}_{1212}\) du mécanisme «père» à la surface de charge du mécanisme déviatoire du plan \(\text{k=}1\) ; coordonnées vectorielles de la normale entrante à la surface de charge du mécanisme «père» du mécanisme déviatoire; rayon du seuil déviatoire atteint par la surface de charge avant la décharge de ce mécanisme déviatoire |
VARI_42 à VARI_47 |
idem pour le mécanisme déviatoire cyclique du plan \(k=2\) |
VARI_48 à VARI_50 |
idem pour le mécanisme déviatoire cyclique du plan \(k=3\) |
Pour les schémas SPECIFIQUE et SEMI_EXPLICITE, l’indicateur d’état VARI_34 des mécanismes actifs après convergence est un octet (8 digits) obtenu par système binaire. Chaque digit est associé à un mécanismes plastique, en affectant la valeur 1 pour un mécanisme actif et la valeur 0 pour un mécanisme non sollicité.
Ceci permet de visualiser facilement avec un simple scalaire les mécanismes actifs. Par exemple, si VARI_34 = 10000010, alors les critères actifs sont: déviatoire \(\mathrm{XZ}\) monotone (\({10}^{1}\) ) et isotrope cyclique (\({10}^{7}\)).
Fonctionnalités et vérification#
La loi de HUJEUX (comportement HUJEUX pour le mot-clé COMPORTEMENT) est utilisable dans code_aster avec différentes modélisations :
version classique: 3D, D_PLAN,
éléments finis sous-intégrés: 3D_SI, D_PLAN_SI,
couplée avec les modèles de THM (cf. [R7.01.11]): 3D_HM_SI, D_PLAN_HM_SI…
Le modèle de HUJEUX est également utilisable dans Code_Aster avec différents types de schémas d’intégration (ALGO_INTE): “SPECIFIQUE”, “SEMI_EXPLICITE”, “BASCULE_EXPLICIT”.
Attention:
La validité de la formulation 3D de la loi Hujeux parait discutable (cf. *[§ 6 :ref:`6 <r7.01.23-chap-4>` ]). Il convient donc d’adopter une certaine prudence quand on l’utilise dans ce cadre.*
Voici la liste des cas de vérification disponibles :
ssnv197abcd |
essai triaxial drainé (mécanique pure) avec la loi de Hujeux; 2 pressions de pré-consolidation : 50 et 200 kPa et introduction de paramètres pour permettre une comparaison avec Gefdyn . Le cas-test SSNV197D est traité avec une rotation du repère local de \(45°\) par rapport à la verticale (direction de chargement). On compare la solution obtenue à un calcul où le maillage est tourné. |
|
ssnv204abc |
essai de consolidation drainée (mécanique pure) avec la loi de Hujeux; intérêt : tester la consolidation puis le cyclique en conditions «mécanique pure» représentatives. SSNV204b : compression isotrope d’un matériau orthotrope élastique. On dégénère la loi de Hujeux en loi élastique orthotrope, et on valide les termes sphériques du tenseur de Hooke orthotrope par rapport à un calcul élastique vrai. SSNV204c : on réalise un chargement purement isotrope conduisant à la mise en traction du sol. On valide ainsi l’activation des mécanismes plastiques de traction. |
|
ssnv205a |
essai de cisaillement cyclique drainé contrôlé en déformation. Intérêt : tester le cyclique en conditions «mécanique pure» représentatives. Comparaison avec Gefdyn. |
|
ssnv207a |
cisaillement cyclique contrôlé en contraintes avec micro-décharge. |
|
ssnv208a |
essai biaxial en conditions drainées sur sable d’Hostun dense (D_PLAN). Calcul lors du post-traitement du critère de Rice (option INDL_ELGA). |
|
ssnv210ab |
Essai de cisaillement monotone drainé en D_PLAN et 3D afin de valider les termes en cisaillement du tenseur de Hooke orthotrope. |
|
wtnv132abcde |
Test de construction par couches en conditions drainées «mécanique pure» avec la loi de Hujeux. WTNV132b : Simulation identique à la précédente surchargée d’une modélisation D_PLAN_DIL. WTNV132c : Simulation identique à la modélisation A avec une résolution globale en matrice sécante. WTNV132d : Simulation identique à la modélisation A avec modélisation D_PLAN_HM_SI. WTNV132e : Simulation identique à la modélisation C avec modélisation D_PLAN_HM_SI. |
|
wtnv133abc |
essai triaxial en conditions non drainées avec la loi de Hujeux. Intérêt : tester la consolidation puis le cyclique en conditions hydromécaniques représentatives. |
|
wtnv134ab |
essai triaxial en conditions non drainées cyclique hydromécanique. Comparaison avec Gefdyn. WTNV134b : chargement appliqué identique utilisant l’opérateur SIMU_POINT_MAT. |
|
WDNP101abcd |
[V7.34.101] |
sollicitation sismique d’une construction par couches avec la loi de Hujeux, modélisation D_PLAN_HM. WDNP101b : Simulation identique sur base physique et non modale comme le cas précédent. WDNP101c,d : Problème adimensionné et utilisant les modélisations D_PLAN_HM_SI et 3D_HM_SI. |
Bibliographie#
«Algorithme non linéaire quasi-statique». Documentation de référence de Code_Aster [R5.03.01].
LEMAITRE, J.L. CHABOCHE: Mécanique des matériaux solides, Dunod 1985.
MODARESSI: Modélisation des milieux poreux sous chargements complexes, Habilitation à Diriger des Recherches, INPG, 2003.
J.C. Hujeux: Une loi de comportement pour le chargement cyclique des sols. Génie Parasismique, Presses ENPC, Davidovici V. & al., pages 287-302, 1985.
Aubry D., J.C. Hujeux, F. Lassoudière & Y. Meimon: « A double memory model with multiple mechanisms for cyclic soil behaviours », Int. Symp. Num. Models in Geomechanics, Zürich, vol. 1, pp 3-13, (1982).
«Intégration explicite des relations de comportement non linéaire». Documentation de référence de Code_Aster [R5.03.14].
ROSCOE K.H., et al. : On the yielding of soils. Géotechnique , vol.8, pp.22-52, 1958.
«Loi de comportement CAM-CLAY». Documentation de référence de Code_Aster [R7.01.14].
COSTA D’AGUIAR. Numerical Modelling of Soil-Pile Load Transfer Mechanisms . Thèse de Doctorat, ECP/LMSSMat, 2008.
S.SICA, L.PAGANO, A. MODARESSI: Influence of past loading history on the seismic response of earth dams. Computers and Geotechnics 35, pp. 61–85, 2008.
LOPEZ-CABALLERO. Influence du comportement non linéaire du sol sur les mouvements sismiques induits dans les géo-structures . Thèse de Doctorat, ECP/LMSSMat.
PICCUEZU. Lois de comportement en géomécanique – Modélisation, mise en œuvre, identification . Thèse de Doctorat, ECP/LMSSMat, 1991.
Annexe 1 : calcul analytique de la matrice tangente d’intégration locale#
On doit résoudre par une méthode de Newton le système d’équations non linéaires locales suivantes: \(\boldsymbol{R}(\Delta \boldsymbol{Y})=0\) , écrites à la fin de l’incrément étudié, où \(\Delta \boldsymbol{Y}=(\Delta {\sigma}_{ij},\Delta {\varepsilon}_{v}^{p},\Delta {r}_{k}^{K},\Delta (\Delta {\lambda}_{k}^{K}))\) . Ceci s’écrit:
avec: \({\gamma}_{k}^{m}=1\) et \({\gamma}_{k}^{c}=\frac{2{q}_{k}^{c}.{\parallel {\boldsymbol{S}}_{(k)H}^{c}\parallel }_{\text{VM}}^{\text{2D}}}{2{q}_{k}^{c}.{\parallel {\boldsymbol{S}}_{(k)H}^{c}\parallel }_{\text{VM}}^{\text{2D}}-{\boldsymbol{S}}_{(k)H}^{c}.{\boldsymbol{S}}_{(k)}^{c}}\) , cf. [éq 1.1.2-5] et [éq 1.1.2-15].
Il faut donc établir le calcul de la matrice tangente \({\frac{d\boldsymbol{R}}{d(\Delta \boldsymbol{Y})}\mid }_{\Delta \boldsymbol{Y}}\) , ce qui nécessite le calcul analytique des différentes dérivées suivantes par rapport à \(\boldsymbol{Y}=({\sigma}_{ij},{\varepsilon}_{v}^{p},{r}_{k}^{K},\Delta {\lambda}_{k}^{K})\) :
Ces dérivées sont déterminées dans les [§ 6.1] à [§ 6.6]. Elles exploitent notamment les résultats établis par [éq 7-13] à [éq 7-24].
Remarque:
Tous les mécanismes apparaissent dans la sommation dans [éq 6-1], mais en pratique seuls figurent les mécanismes actifs.
Les directions d’écoulement plastique \({(\boldsymbol{\Psi} )}_{k}^{K}\) sont exprimées après injection dans le repère global cartésien sous la forme d’un tenseur d’ordre 2.
Dérivées de l’équation d’état#
Elles proviennent de l’équation d’état [éq 1.1.1-9] et de la déformation plastique [éq 1.1.5-1]:
avec :
Calcul des composantes \(\frac{\partial {\text{LE}}_{ij}}{\partial {\sigma}_{\text{rs}}}\) (tenseur d’ordre 4) :
La contribution des mécanismes de consolidation étant nulle dans le dernier terme. L’expression des différents termes mis en jeu dans le calcul de cette dérivée est, en exploitant les résultats [éq 7-6] et [éq 7-7] :
pour chaque mécanisme déviatoire monotone \(\text{t=}1,..,3\) , cf. [éq 1.1.2-3] et :
pour chaque mécanisme déviatoire cyclique \(\text{t=}1,...,3\) , cf. [éq 1.1.2-14].
D’où pour les mécanismes déviatoires monotones :
et pour les mécanismes déviatoires cycliques :
les divers termes étant donnés par [éq 7-7], [éq 7-8], [éq 7-9], [éq 7-10] et [éq 7-11].
Calcul des composantes \(\frac{\partial {\text{LE}}_{ij}}{\partial {\varepsilon}_{v}^{p}}\) (tenseur d’ordre 2) :
Calcul des composantes \(\frac{\partial {\text{LE}}_{ij}}{\partial {r}_{k}^{K}}\) (tenseurs d’ordre 2) :
On constate que, pour ; \(K=m,c;\text{k=}1,...,3\) , cf. [éq 6.1-1]:
D’où, la fonction \(\zeta (r)\) étant définie en [éq 1.1.2-4] et \({S}_{(k)}^{K}\) en [éq 1.1.1-5] et [éq 1.1.2-12]:
Mécanismes déviatoires monotones
Mécanismes déviatoires cycliques
la fonction \({F}_{k}\text{=F}({p}_{k},{\varepsilon}_{v}^{p})\) étant définie en [éq 1.1.2-2].
De plus, pour les mécanismes de consolidation: \(\frac{\partial {\text{LE}}_{ij}}{\partial {r}_{4}^{K}}=0\) , pour \(K=m,c\) .
Calcul des composantes \(\frac{\partial {\text{LE}}_{ij}}{\partial \Delta {\lambda}_{k}^{K}}\) (tenseurs d’ordre 2) :
Pour \(K=m,c;\text{k=}1,...,3\) :
Dérivées de l’équation d’évolution de la déformation plastique#
Elles proviennent de l’évolution de la déformation plastique volumique, cf. [éq 6.2-1] :
Calcul des composantes \(\frac{\partial \mathrm{LEVP}}{\partial {\sigma}_{ij}}\) (tenseur d’ordre 2) :
en utilisant les résultats [éq 7-7], les interventions des mécanismes de consolidation étant nulles.
Calcul de la composante \(\frac{\partial \text{LEVP}}{\partial {\varepsilon}_{v}^{p}}\) (scalaire) :
(cf. [éq 1.1.2-3] et [éq 1.1.2-14]).
Calcul des composantes \(\frac{\partial \text{LEVP}}{\partial {r}_{k}^{K}}\) (scalaires) :
Pour \(K=m,c;\text{k=}1,...,3\) , en utilisant la fonction \(\zeta (r)\) , cf. [éq 1.1.2-4] :
Mécanismes déviatoires monotones
Mécanismes déviatoires cycliques
où \({F}_{k}\text{=F}({p}_{k},{\varepsilon}_{v}^{p})\) est définie en [éq 1.1.2-2]. De plus, \(\frac{\partial \text{LEVP}}{\partial {r}_{4}^{m}}=\frac{\partial \text{LEVP}}{\partial {r}_{4}^{c}}=0\) .
Calcul des composantes \(\frac{\partial \text{LEVP}}{\partial \Delta {\lambda}_{k}^{K}}\) (scalaires) :
Pour ; \(K=m,c;\text{k=}1,..,3\) :
Mécanismes déviatoires monotones :
Mécanismes déviatoires cycliques :
Mécanismes de consolidation :
Dérivées de l’équation d’évolution de l’écrouissage déviatoire#
Elles proviennent de l’évolution de l’écrouissage des mécanismes déviatoires cf. [éq 1.1.2-5] et [éq 1.1.2-15] :
Calcul des composantes \(\frac{\partial \text{LR}.1}{\partial {\sigma}_{ij}}\) (tenseur d’ordre 2) :
Pour les mécanismes déviatoires monotones, \(\text{k=}1,...,3\) : \({\gamma}_{k}^{m}=1\) , donc :
Pour les mécanismes déviatoires cycliques, \(\text{k=}1,...,3\) :
avec l’opérateur \(\frac{\partial {({S}_{(k)}^{c})}_{\mathrm{rs}}}{\partial {\sigma}_{ij}}\) défini par [éq 6-8].
Calcul des composantes \(\frac{\partial \text{LR}.1}{\partial {\varepsilon}_{v}^{p}}\) (scalaires) :
Pour les mécanismes déviatoires monotones, \(\text{k=}1,...,3\) : \({\gamma}_{k}^{m}=1\) , donc :
Pour les mécanismes déviatoires cycliques, \(\text{k=}1,...,3\) :
avec \(\frac{\partial {({S}_{(k)}^{c})}_{ij}}{\partial {\varepsilon}_{v}^{p}}=\beta M{b}_{h}.{p}_{k}.{\left({\boldsymbol{X}}_{(k)}^{H}+\frac{{\boldsymbol{S}}_{(k)H}^{c}}{{\parallel {\boldsymbol{S}}_{(k)H}^{c}\parallel }_{\text{VM}}^{\text{2D}}}.({r}_{k}^{c}+{r}_{\text{éla}}^{\text{dc}})\right)}_{ij}\) , cf. [éq 1.1.2-12] et [éq 1.1.2-2].
Calcul des composantes \(\frac{\partial \text{LR}.1}{\partial {r}_{k}^{K}}\) (scalaires) :
Pour les mécanismes déviatoires monotones, \(\text{k=}1,...,3:{\gamma}_{k}^{m}=1\) , donc :
soit :
Pour les mécanismes déviatoires cycliques, \(\text{k=}1,...,3\) :
avec \(\frac{\partial {({S}_{(k)}^{c})}_{ij}}{\partial {r}_{k}^{c}}=-{p}_{k}(\boldsymbol{\sigma} ).F({p}_{k}(\boldsymbol{\sigma} ),{\varepsilon}_{v}^{p}).\frac{{({S}_{(k)H}^{c})}_{ij}}{{\parallel {\boldsymbol{S}}_{(k)H}^{c}\parallel }_{\text{VM}}^{\text{2D}}}\) , cf. [éq 1.1.2-12].
De plus: \(\frac{\partial \text{LR}.1}{\partial {r}_{4}^{K}}=0\) .
Calcul des composantes \(\frac{\partial \text{LR}.1}{\partial \Delta {\lambda}_{k}^{K}}\) (scalaires) :
Pour \(K=m,c;\text{ k=}1,...,3\) :
avec: \({\gamma}_{k}^{m}=1\) et \({\gamma}_{k}^{c}=\frac{2{q}_{k}^{c}.{\parallel {\boldsymbol{S}}_{(k)H}^{c}\parallel }_{\text{VM}}^{\text{2D}}}{2{q}_{k}^{c}.{\parallel {\boldsymbol{S}}_{(k)H}^{c}\parallel }_{\text{VM}}^{\text{2D}}-{\boldsymbol{S}}_{(k)H}^{c}.{\boldsymbol{S}}_{(k)}^{c}}\) .
De plus \(\frac{\partial \text{LR}.1}{\partial \Delta {\lambda}_{4}^{K}}=0\) .
Dérivées de l’équation d’évolution de l’écrouissage sphérique#
Elles proviennent de l’évolution de l’écrouissage des mécanismes de consolidation sphériques, cf. [éq 1.1.3-3] et [éq 1.1.3-11]:
Calcul des composantes \(\frac{\partial \text{LR}.2}{\partial {\sigma}_{ij}}\) (tenseur d’ordre 2) :
Calcul des composantes \(\frac{\partial \text{LR}.2}{\partial {\varepsilon}_{v}^{p}}\) (scalaires) :
Pour \(K=m,c\) :
Calcul des composantes \(\frac{\partial \text{LR}.2}{\partial {r}_{4}^{K}}\) (scalaires) :
Pour \(K=m,c\) :
Calcul des composantes \(\frac{\partial \text{LR}.2}{\partial \Delta {\lambda}_{4}^{K}}\) (scalaires) :
Pour \(K=m,c\) :
De plus: \(\frac{\partial \text{LR}.2}{\partial {r}_{k}^{K}}=0\) , pour les mécanismes déviatoires \(\text{k=}1,...,3\) .
Dérivées des critères des mécanismes déviatoriques#
Elles proviennent des critères [éq 1.1.2-1] et [éq 1.1.2-8] :
Calcul des composantes \(\frac{\partial \text{LF}.1}{\partial {\sigma}_{ij}}\) (tenseur d’ordre 2) :
Pour les mécanismes monotones, cf. [éq 7-13]:
Pour les mécanismes déviatoires, cf. [éq 7-16]:
Calcul des composantes \(\frac{\partial \text{LF}.1}{\partial {\varepsilon}_{v}^{p}}\) (scalaires) :
Pour \(\text{k=}1,...,3\) , cf. [éq 1.1.2-1] et [éq 1.1.2-10] et [éq 7-14]:
Mécanismes déviatoires monotones :
Mécanismes déviatoirescycliques :
avec \({\boldsymbol{X}}_{(k)}^{H}\) défini en [éq 1.1.2-13].
Calcul des composantes \(\frac{\partial \text{LF}.1}{\partial {r}_{k}^{K}}\) (scalaires) :
Pour \(\text{k=}1,...,3\) :
Mécanismes déviatoires monotones :
Mécanismes déviatoires cycliques :
Calcul des composantes \(\frac{\partial \text{LF}.1}{\partial \Delta {\lambda}_{k}^{K}}\) (scalaires):
Dérivées des critères des mécanismes de consolidation sphériques#
Elles proviennent des critères [éq 1.1.3-1] et [éq 1.1.3-8] :
avec dans le cas cyclique: \({p}^{c}(\boldsymbol{\sigma} ,{\varepsilon}_{v}^{p},{p}_{H},{\varepsilon}_{\text{vH}}^{p})=\mid p(\boldsymbol{\sigma} )\mid {\text{+p}}_{H}.{e}^{-\beta ({\varepsilon}_{v}^{p}-{\varepsilon}_{\text{vH}}^{p})}\) , cf. [éq 1.1.3-9].
Calcul des composantes \(\frac{\partial \text{LF}.2}{\partial {\sigma}_{ij}}\) (tenseur d’ordre 2) :
Mécanisme sphérique monotone :
Mécanisme sphérique cyclique :
Calcul des composantes \(\frac{\partial \text{LF}.2}{\partial {\varepsilon}_{v}^{p}}\) (scalaires):
Mécanisme sphérique monotone :
Mécanisme sphérique cyclique :
Calcul des composantes \(\frac{\partial \text{LF}.2}{\partial {r}_{4}^{K}}\) (scalaires):
Mécanisme sphérique monotone :
Mécanisme sphérique cyclique :
cf. [éq 7-21] et [éq 7-24], avec \({p}^{c}\) défini en [éq 1.1.3-9].
Calcul des composantes \(\frac{\partial \text{LF}.2}{\partial \Delta {\lambda}_{4}^{K}}\) (scalaires) :
Annexe 2 : notation des tenseurs, de leurs invariants et expressions de diverses dérivées#
On se place pour plus de facilité en base orthonormée \(({\mathrm{e}}_{1,}{\mathrm{e}}_{2,}{\mathrm{e}}_{3})\) de l’espace à 3 dimensions. On définit le tenseur identité d’ordre 2, ici en utilisant la notation de Walpole-Cowin :
On note, pour chaque mécanisme dans le plan \(k\) :
avec
et
\({\boldsymbol{\sigma}}_{(k)}\) étant défini par [éq 1.1.1-3], pour \({i}_{k}=1+\mathrm{mod}(k,3)\) et \({j}_{k}=1+\mathrm{mod}(k+1,3)\).
On a la relation suivante entre les termes du tenseur \({\boldsymbol{S}}_{(k)}\) :
On définit de plus les tenseurs identité d’ordre 2 (pour chaque plan) et 4 :
On donne ci-dessous les expressions de plusieurs dérivées qui interviennent plusieurs fois dans les expressions analytiques du modèle.
On remarque que: \({\boldsymbol{\mathrm{Dev}}}_{(k)}\underline{\underline{\otimes}}{\boldsymbol{{S}}}_{(k)}={\boldsymbol{S}}_{(k)}\). De plus :
On a aussi : \({\boldsymbol{\mathrm{Dev}}}_{(k)}\underline{\underline{\otimes}}{\boldsymbol{S}}_{(k)}^{c}={\boldsymbol{S}}_{(k)}^{c}\) . De plus, en \({q}_{k}=0\) (ou \({q}_{k}^{c}=0\) ), la dérivée \(\frac{\partial {q}_{k}}{\partial {\boldsymbol{S}}_{(k)}}\) (ou \(\frac{\partial {q}_{k}^{c}}{\partial {\boldsymbol{S}}_{(k)}^{c}}\) ) n’est pas définie. En pratique, on prendra la valeur \(\boldsymbol{0}\).
On utilise, de manière à alléger les écritures, la notation de Walpole-Cowin des tenseurs d’ordre 2 et 4.
car \({S}_{{j}_{k}{j}_{k}}=-{S}_{{i}_{k}{i}_{k}}\) , ce tenseur étant déviateur dans le plan \(k\) . On obtient notamment :
Par ailleurs, pour le calcul des multiplicateurs plastiques et des opérateurs tangents, on a besoin de :
En \({q}_{k}=0\) , on retire le premier terme de cette expression, cf. [éq 7 -7].
Enfin, on note que :
Annexe 3 : Validité de la formulation multi-mécanisme de la loi de Hujeux#
Représentation d’un état de contrainte par les cercles de Mohr#
L’état de contrainte quelconque d’un point M peut être représenté par un tenseur du type:
Ce tenseur symétrique est diagonalisable, avec \({\sigma}_{1}\ge {\sigma}_{2}\ge {\sigma}_{3}\) les valeurs propres (ou contraintes principales) de :math:`sigma ` telles que:
La représentation de Mohr permet de représenter l’action de cet état de contrainte sur n’importe quelle facette de normale \(n\) tournant autour du point M (Figure 1). En effet, le vecteur contrainte s’exerçant sur cette facette s’écrit :
Fig. 221 Décomposition du vecteur contrainte s’exerçant sur une facette de normale n#
La représentation de Mohr permet de représenter dans un même diagramme les composantes \(({\sigma}_{n},\tau )\) :
Ceci est illustré sur la Figure 2. Cette représentation aboutit au résultat très important qu’on met en exergue par la proposition suivante:
Proposition 1
Tout état de contrainte en un point M s’exerçant sur une facette de normale quelconque tournant autour de M coïncide exactement , dans la représentation de Mohr, à toute la zone en vert comprise entre les trois cercles de Mohr.
On tient ce résultat pour suffisamment connu pour qu’on n’ait pas à l’expliquer plus que ça. Il est très important car il définit le domaine de validité physique pour tout état de contrainte , et on va montrer que dans la loi de Hujeux, cette contrainte n’est pas respectée en 3D !
Fig. 222 Domaine de définition d’un état quelconque de contrainte dans le plan de Mohr#
Critère de plasticité de type Mohr-Coulomb#
Avant d’en arriver à Hujeux, il convient d’expliquer la relation entre la représentation de Mohr et un critère de rupture de type Mohr-Coulomb. L’expression du critère de Mohr-Coulomb est la suivante :
avec les contraintes principales \({\sigma}_{1}\ge {\sigma}_{2}\ge {\sigma}_{3}\) .
Sa représentation dans le plan de Mohr est illustrée sur la Figure 3. La signification de ce critère est de comparer l’état de contraintes au point M, représenté par le cercle de centre \(C=\frac{{\sigma}_{1}+{\sigma}_{3}}{2}\) et de rayon \(R=\frac{{\sigma}_{1}-{\sigma}_{3}}{2}\) , à la droite de rupture d’angle \(\mathrm{\phi }\) . Il convient de noter deux points importants :
L’état de contrainte au point M est représenté par le cercle ce centre \(C\) et de rayon \(R\) dans le plan de Mohr, et non par un point, car ce cercle est le lieu des couples \(({\sigma}_{n},\tau )\) pour l’ensemble des plans de normales \(n\perp {n}_{2}\) tournant autour de M ;
La tangence du cercle et de la droite de rupture est obtenue grâce au terme \(\sin\left(\mathrm{\phi }\right)\) , comme on peut le déduire à partir d’une analyse géométrique relativement simple.
Fig. 223 Représentation du critère de Mohr-Coulomb dans le plan de Mohr#
On résume les acquis de cette analyse par la proposition suivante :
Proposition 2
Dans un critère de type Mohr-Coulomb, on compare l’état de contraintes au point M représenté par un cercle de centre \(C\) et de rayon \(R\) dans le plan de Mohr à une droite de rupture d’angle \(\mathrm{\phi }\) . L’expression mathématique d’un tel critère est de la forme : \(R-\sin\left(\mathrm{\phi }\right)C\le 0\).
Expression du critère de plasticité de Hujeux#
Le critère de plasticité de Hujeux est un critère de type Mohr-Coulomb. On commence par décomposer le tenseur des contraintes par projection sur les trois plans fixes du repère cartésien, en procédant comme suit :
où \({\boldsymbol{p}}_{\vec{{e}_{i}}}\left(\boldsymbol{\sigma} \right)\) est la projection du tenseur \(\boldsymbol{\sigma}\) sur le plan de normale \({\vec{e}}_{i\in \lbrace x,y,z\rbrace }\) . Réciproquement, on aura l’inverse de l’opérateur de projection :
\({\boldsymbol{p}}_{\vec{{e}_{i}}}\left(\boldsymbol{\sigma} \right)\) est appelé tenseur réduit et peut être exprimé dans chaque plan du repère cartésien suivant la normale considérée :
\({e}_{x}\) |
\({e}_{y}\) |
\({e}_{z}\) |
|
\({\boldsymbol{p}}_{\vec{{e}_{i}}}\left(\boldsymbol{\sigma} \right)\) |
\(\left[\begin{array}{cc}{\sigma}_{yy}& {\sigma}_{yz}\\ {\sigma}_{\mathit{zy}}& {\sigma}_{zz}\end{array}\right]\) |
\(\left[\begin{array}{cc}{\sigma}_{xx}& {\sigma}_{xz}\\ {\sigma}_{\mathit{zx}}& {\sigma}_{zz}\end{array}\right]\) |
\(\left[\begin{array}{cc}{\sigma}_{xx}& {\sigma}_{xy}\\ {\sigma}_{yx}& {\sigma}_{yy}\end{array}\right]\) |
Fig. 224 Représentation des tenseurs de contraintes réduits dans chaque plan orthotrope du repère cartésien#
Fig. 225 Représentation des tenseurs de contraintes réduits dans le plan de Mohr#
Le critère de Hujeux s’écrit alors comme suit :
\(R-\sin(\mathrm{\phi })C\le 0\)
avec :
Il est clair d’après la proposition 2 que le critère de Hujeux est un critère de type Mohr-Coulomb. En conséquence, il est possible de représenter le tenseur des contraintes réduit par un cercle de Mohr dans le plan de Mohr. La signification du critère est alors de comparer le cercle de centre \(C\) et de rayon \(R\) à la droite de rupture d’angle \(\mathrm{\phi }\) .
Remarques :
Comment se justifie cette formulation? D’après notre compréhension, il semble que :
Le critère de Hujeux est dérivé par analogie du critère de Mohr-Coulomb appliqué à des chargements triaxiaux ;
La décomposition du tenseur des contraintes étant bijective et inversible est justifiée mathématiquement, et donc agréable à l’esprit.
Validité du critère de Hujeux#
Pour montrer que la décomposition des contraintes utilisée dans la loi de Hujeux n’est pas valable physiquement, il nous suffit de trouver un contre-exemple à la proposition 1. Or, en trouver un en 3D est relativement trivial : il suffit de considérer un état de contrainte 3D quelconque.
Soit par exemple l’état de contrainte suivant :
\(\sigma =\left[\begin{array}{ccc}0.6& -0.1& 0.05\\ -0.1& 0.4& 0.1\\ 0.05& 0.1& 0.3\end{array}\right]\) (en MPa)
On obtient les contraintes principales suivantes :
\(\left\lbrace \begin{array}{c}{\sigma}_{1}=0.6416\\ {\sigma}_{2}=0.4445\\ {\sigma}_{3}=0.2139\end{array}\right.\) (en MPa)
D’où on déduit les cercles de Mohr.
On calcule ensuite les cercles de Hujeux dans le plan de Mohr. Pour ce faire, on calcule pour le premier cercle :
\(\left\lbrace \begin{array}{c}R=\sqrt{{\left(\frac{{\sigma}_{xx}-{\sigma}_{yy}}{2}\right)}^{2}+{\sigma}_{xy}^{2}}\\ C=\frac{{\sigma}_{xx}+{\sigma}_{yy}}{2}\end{array}\right.\)
Et ainsi de suite pour les deux autres.
Enfin, il est intéressant de caractériser le lieu des contraintes réelles \(({\sigma}_{n},\tau )\) dans le plan de Mohr donné par une normale \({n}_{i\in \lbrace x,y,z\rbrace }\) tournant autour du point M dans le plan orthotrope de normale \({\vec{e}}_{i\in \lbrace x,y,z\rbrace }\) . L’expression des normales est la suivante:
\({n}_{x}=\left(\begin{array}{c}0\\ \cos\left(\theta \right)\\ \sin\left(\theta \right)\end{array}\right);\phantom{\rule{4em}{0ex}}{n}_{y}=\left(\begin{array}{c}\cos\left(\theta \right)\\ 0\\ \sin\left(\theta \right)\end{array}\right);\phantom{\rule{4em}{0ex}}{n}_{z}=\left(\begin{array}{c}\cos\left(\theta \right)\\ \sin\left(\theta \right)\\ 0\end{array}\right)\)
On calcule alors :
\(\left\lbrace \begin{array}{c}{\sigma}_{n,i}={\sigma}_{n}.{n}_{i}\\ {\tau}_{i}={\sigma}_{n,i}-\left({\sigma}_{n,i}.{n}_{i}\right){n}_{i}\end{array}\right.\)
La contrainte réelle s’écrit :
\(\left\lbrace \begin{array}{c}{\sigma}_{n,i}={\sigma}_{n,i}.{n}_{i}\\ {\tau}_{i}=\sqrt{{\tau}_{i}.{\tau}_{i}}\end{array}\right.\)
Les solutions sont présentées dans la Figure 5 :
Les cercles de Mohr sont représentés en traits noirs pointillés ;
Les cercles de Hujeux sont représentés en traits pleins colorés ;
Le lieu des contraintes réelles dans les plans orthotropes est représenté en traits pointillés colorés de la même couleur que les cercles de Hujeux correspondants.
D’après la proposition 1, il est clair que les cercles de Hujeux ne correspondent pas à un état de contrainte réel autour du point M. L’état de contrainte réel correspondant devrait plutôt être celui donné en pointillés, lequel vérifie bien la proposition 1. On en conclut que le critère de Hujeux, qui compare le cercle de Hujeux à une droite dans le plan de Mohr, ne représente pas non plus l’action phénoménologique d’un mécanisme réel.
Fig. 226 Représentation de l’état de contrainte au point M dans le plan de Mohr : les cercles de Mohr sont représentés en pointillés noirs ; les cercles de Hujeux en traits pleins colorés ; le lieu des contraintes réelles dans les 3 plans orthotropes cartésiens en traits pointillés colorés avec la même couleur que les cercles de Hujeux correspondants.#
Remarque
Il est intéressant de noter que l’on retrouve les cercles de Hujeux en procédant comme suit :
\({y}_{x}=\left(\begin{array}{c}0\\ -\sin\left(\theta \right)\\ \cos\left(\theta \right)\end{array}\right);\phantom{\rule{4em}{0ex}}{t}_{y}=\left(\begin{array}{c}-\sin\left(\theta \right)\\ 0\\ \cos\left(\theta \right)\end{array}\right);\phantom{\rule{4em}{0ex}}{t}_{z}=\left(\begin{array}{c}-\sin\left(\theta \right)\\ \cos\left(\theta \right)\\ 0\end{array}\right)\)
On calcule alors :
\(\left\lbrace \begin{array}{c}{\sigma}_{n,i}={\sigma}_{n,i}.{n}_{i}\\ {\tau}_{i}={\sigma}_{n,i}.{t}_{i}\end{array}\right.\)
avec \({\sigma}_{n,i}=\sigma .{n}_{i}\)
Conclusions#
L’état de contrainte utilisé pour former le critère de Hujeux est purement théorique, mais non existant . En ce sens, il ne traduit pas un mécanisme réel pour un état de contrainte tridimensionnel ;
Par contre, pour un état de contraintes « plan », c’est-à-dire celui contenu dans un plan perpendiculaire à une contrainte principale, le cercle de Hujeux coïncide parfaitement avec le cercle de Mohr. C’est le cas par exemple des problèmes en déformations ou en contraintes planes, pour lesquels le critère de Hujeux reste donc valable ;
La formulation d’un critère 3D valide passe par l’utilisation des invariants de contrainte (\({I}_{1},{I}_{2},{I}_{3}\) ou les contraintes principales).
Annexe 4 : Stratégie de redécoupage interne de la loi de Hujeux#
Le redécoupage interne programmé dans la loi de Hujeux se produit à un double niveau :
au niveau de la subroutine générique REDECE, point d’entrée de toutes les lois de comportements (LdC) dans code_aster ;
au niveau de la loi de Hujeux elle-même.
La routine REDECE gère de façon générique le redécoupage interne pour toutes les lois de comportement. Le synoptique de la routine est donnée sur la page suivante. L’idée est basée sur l’interception d’un code retour envoyé par la LdC :
IRET = 0 si l’intégration de la LdC a réussi
IRET = 1 si elle a échoué
IRET = 2 si elle a échoué
Lorsque l’intégration de la LdC a échoué, et si IRET= 2, on sort directement de la routine. Lorsque l’intégration de la LdC a échoué, et si IRET= 1, on tente un redécoupage interne en fonction de la valeur de ITER_INTE_PAS :
ITER_INTE_PAS = -1, 0, 1 : on sort de la routine
ITER_INTE_PAS > 1 : on redécoupe par ITER_INTE_PAS sur 3 niveaux
ITER_INTE_PAS < -1 : on redécoupe par –ITER_INTE_PAS sur 3 niveaux
Notons que le nombre maximal de redécoupages est ainsi de :
\(\text{NB_MAX}=\text{ITER_INTE_PAS}\times ({2}^{0}+{2}^{1}+{2}^{2})=7\times \text{ITER_INTE_PAS}\)
De plus, comme dans la routine de la loi de Hujeux, il y a aussi un redécoupage par ITER_INTE_MAXI , le nombre maximal de redécoupages pour la loi de Hujeux est ainsi de :
\(\text{NB_MAX}=7\times \text{ITER_INTE_PAS}\times \text{ITER_INTE_MAXI}\)
Enfin, la matrice tangente consistante récupérée est moyennée sur l’ensemble des pas de redécoupages. Cela n’a de sens que si celle-ci présente des propriétés de continuité suffisantes d’un sous-pas à l’autre, ce qui n’est généralement pas le cas des lois multi-mécanismes. Le synoptique de cet algorithme est donné dans la figure ci-dessous.
Fig. 227 Algorithme d’intégration du comportement au niveau de redece.F90#
La figure ci-dessous donne le synoptique du nouvel algorithme d’intégration du comportement avec les fléchages suite à un échec (IRET = 1) selon la valeur de ALGO_INTE.
Fig. 228 Algorithme d’intégration du comportement au niveau de la loi Hujeux suivant les 3 options du mot-clé ALGO_INTE (IPAL = ITER_INTE_PAS – IMAXI = ITER_INTE_MAXI)#