r7.01.38 Loi d’Iwan pour le comportement cyclique de matériaux granulaires

Résumé:

Le modèle de comportement d’Iwan [ 1 ], [ 2 ] est un modèle élastoplastique multimécanisme classique pour la description du comportement déviatorique cyclique de matériaux et particulièrement de géomatériaux. L’écrouissage est cinématique linéaire pour chaque mécanisme et la règle d’écoulement considérée est associée. Par construction, le modèle vérifie automatiquement les règles de Masing.

Le modèle permet une description aisée du comportement via des paramètres directement obtenus des courbes de variation du module de cisaillement sécant. L’utilisateur doit simplement fournir les paramètres d’une fonction hyperbolique décrivant la variation du module sécant en fonction de la déformation de cisaillement.

Il est implémenté dans l’environnement MFront, ce qui facilite l’implémentation et la maintenabilité du modèle.

Intégration numérique de la loi de comportement

Le choix dans MFront est porté sur l’intégration des systèmes d’équations différentiels écrits de la manière suivante:

(3884)\[\dot{Y}=G\left(Y,t\right)\]

Où \(G\) est une fonction a priori non linéaire et supposée a minima continûment dérivable. Une intégration implicite de type thêta-méthode est choisie pour l’intégration numérique du système non-linéaire. Dans ce cas le système s’écrit de la manière suivante:

(3885)\[Y{|}_{t+\Delta t}=Y{|}_{t}+\Delta tG\left(\left\lbrace 1-\theta \right\rbrace Y{|}_{t}+\theta Y{|}_{t+\Delta t},t+\theta \Delta t\right)\]

Ou de manière équivalente:

(3886)\[F\left(\Delta Y\right)=\Delta Y-\Delta tG\left(\left\lbrace 1-\theta \right\rbrace Y{|}_{t}+\theta Y{|}_{t+\Delta t},t+\theta \Delta t\right)=0\]

La méthode de Newton-Raphson est choisie pour le calcul du zéro de la fonction \(F\left(\Delta Y\right)\) . Soit:

(3887)\[F\left(\Delta Y\right)=F\left(\Delta {Y}^{-}\right)+{\left(\frac{dF}{d\Delta Y}\right)}_{\Delta Y=\Delta {Y}^{-}}\left(\Delta Y+\Delta {Y}^{-}\right)\]

Où \(J=\frac{dF}{d\Delta Y}\) est la matrice jacobienne de la fonction \(F\) . L’approximation \(\Delta {Y}^{+}\) est obtenue en considérant \(F\left(\Delta {Y}^{+}\right)=0\) , soit:

(3888)\[\Delta {Y}^{+}=\Delta {Y}^{-}-\left({J}^{-1}\left(\Delta {Y}^{-}\right)\right)F\left(\Delta {Y}^{-}\right)\]

La difficulté de la méthode de Newton-Raphson repose sur le calcul de la matrice jacobienne, dont le calcul analytique est fourni dans le section 2.4 . Un schéma purement implicite est choisi pour l’intégration des équations (\(\theta =1\) ), afin de donner des résultats physiquement satisfaisants.

Suivant la démarche présentée dans le document [ 4 ], pour un écrouissage cinématique linéaire il est possible de réduire le système tensoriel à un système scalaire, grâce à la colinéarité entre l’incrément de déformation plastique \({\dot{\epsilon}}_{n}^{p}\) et la différence \(S-{X}_{n}\) en absence de déformation plastique, \({S}_{n}^{e}\) . Cela peut s’écrire de la manière suivante pour un écrouissage cinématique multiple:

(3889)\[{S}_{n}^{e}={S}^{-}+\dot{S}-{X}_{n}^{-}={S}^{-}+2\mu {\dot{\epsilon}}_{n}^{e}-{X}_{n}^{-}\]

Donc:

(3890)\[S-{X}_{n}={S}_{n}^{e}-{\dot{X}}_{n}={S}_{n}^{e}-{C}_{n}{\dot{\epsilon}}_{n}^{p}\]

L’incrément de déformation plastique du mécanisme \(>\) s’écrit donc:

(3891)\[{\dot{\epsilon}}_{n}^{p}=\frac{3}{2}\frac{{\dot{\lambda}}_{n}}{{q}_{n}}\left({S}_{n}^{e}-{C}_{n}{\dot{\epsilon}}_{n}^{p}\right)=\frac{3}{2}\frac{{\dot{\lambda}}_{n}}{{q}_{n}^{e}}{S}_{n}^{e}\]

Avec \({q}_{n}^{e}\) la norme de \({S}_{n}^{e}\) . En remplaçant \({\dot{\epsilon}}_{n}^{p}\) dans l’expression (), on obtient:

(3892)\[S-{X}_{n}={S}_{n}^{e}\left(1-\frac{\frac{3}{2}{C}_{n}\dot{{\lambda}_{n}}}{{q}_{n}^{e}}\right)\]

En injectant l’expression de \(S-{X}_{n}\) dans la condition \(f=0\) , on obtient:

(3893)\[\left({q}_{n}^{e}-\frac{3}{2}{C}_{n}\dot{{\lambda}_{n}}\right)-{Y}_{n}=0\]

Dans la mesure où \({Y}_{n}>0\) et \({q}_{n}^{e}>0\) , cette équation admet la solution suivante pour le multiplicateur plastique:

(3894)\[\dot{{\lambda}_{n}}=\frac{{q}_{n}^{e}-{Y}_{n}}{\frac{3}{2}{C}_{n}}\]

Donc, forcément:

(3895)\[0<\dot{{\lambda}_{n}}<\frac{{q}_{n}^{e}}{\frac{3}{2}{C}_{n}}\]

Le système à résoudre est donc directement basé sur le multiplicateur plastique et non sur le tenseur d’écrouissage, qui est une variable auxiliaire (@AuxiliaryStateVariable) dans la routine MFront.

Dans la mesure où la mise à jour du tenseur de contraintes dans le bloc @ComputeStresss’effectue avant l’itération de Newton de l’algorithme de résolution, dans la routine MFront \({q}_{n}^{e}\) est obtenu directement à partir de la valeur mise à jour du tenseur déviatorique \(\mathrm{S}\) , \({q}_{n}\) .

Étape de prédiction

L’étape de prédiction est gérée à partir du mot-clé facteur PREDICTION du mot-clé NEWTON dans l’opérateur de résolution non linéaire, le bloc @Predictor sur MFront n’ayant pas besoin d’être renseigné pour une résolution avec code_aster .

Étape de correction par algorithme de Newton

Le vecteur d’inconnues étant donné par \(Y=\left({\epsilon}^{e},{\dot{\lambda}}_{n}\right)\) , les itérations de l’algorithme de Newton consistent à résoudre le système suivant:

(3896)\[\begin{split}\lbrace \begin{array}{c}{f}_{\epsilon}^{e}={\dot{\epsilon}}^{e}-\dot{\epsilon}+\sum_{n=1}^{N}{\dot{\lambda}}_{n}{\Psi}_{n}=0\\ {f}_{\dot{{\lambda}_{n}}}=\dot{{\lambda}_{n}}-\frac{{q}_{n}-{Y}_{n}}{\frac{3}{2}{C}_{n}}\end{array}\end{split}\]

L’écriture des équations en forme scalaire permet de réduire la taille du système d’équations à résoudre et de cette manière d’améliorer la performance du modèle.

Étape de mise à jour

Une fois la convergence de la résolution locale achevée, le vecteur des inconnues est mis à jour classiquement de la manière suivante:

(3897)\[{Y}^{+}={Y}_{\text{j+}1}={Y}_{j}+\delta {Y}_{\text{j+}1}\]

On profite également pour mettre à jour les variables auxiliaires \({X}_{n}\) et \({\epsilon}^{p}\) dans le bloc @UpdateAuxiliaryStateVars.

Calcul des termes de la matrice jacobienne

Afin de calculer l’opérateur tangent cohérent (directive @TangentOperatorde MFront), il faut disposer des termes de la matrice jacobienne \(J\) du système d’équations à résoudre. Dans le cas présent, la matrice jacobienne s’écrit de la manière suivante:

(3898)\[\begin{split}J=\left[\begin{array}{cc}\frac{\partial {f}_{\dot{{\lambda}_{n}}}}{\partial \dot{{\lambda}_{n}}}{I}_{4}& \frac{\partial {f}_{\epsilon}^{e}}{\partial \dot{{\lambda}_{n}}}\otimes I\\ \frac{\partial {f}_{\dot{{\lambda}_{n}}}}{\partial {\dot{\epsilon}}^{e}}\otimes I& \frac{\partial {f}_{\epsilon}^{e}}{\partial {\dot{\epsilon}}^{e}}\end{array}\right]\end{split}\]

MFront permet d’obtenir ces termes numériquement par l’utilisation de l’algorithme NewtonRaphson_NumericalJacobian. Mais la version développée ici utilise une matrice analytique, plus performante.Dans ce qui suit, on propose une description analytique des dérivées partielles en considérant \(\theta =1\) .

Les termes à calculer sont les suivants: \(\frac{\partial {f}_{{\dot{\lambda}}_{n}}}{\partial {\dot{\lambda}}_{n}}\) (scalaire), \(\frac{\partial {f}_{{\dot{\lambda}}_{n}}}{\partial {\dot{\epsilon}}^{e}}\) (tenseur d’ordre deux), \(\frac{\partial {f}_{\epsilon}^{e}}{\partial \dot{{\lambda}_{n}}}\) (tenseur d’ordre deux) et \(\frac{\partial {f}_{\epsilon}^{e}}{\partial {\dot{\epsilon}}^{e}}\) (tenseur d’ordre quatre). Les dérivées partielles par rapport au multiplicateur plastique \({\dot{\lambda}}_{n}\) pour chaque mécanisme \(n\) sont immédiates:

(3899)\[\frac{\partial {f}_{\dot{{\lambda}_{n}}}}{\partial \dot{{\lambda}_{n}}}=1\]

Et:

(3900)\[\frac{\partial {f}_{\epsilon}^{e}}{\partial \dot{{\lambda}_{n}}}={\Psi}_{n}\]

Il reste à calculer les dérivées par rapport à la déformation élastique \({\dot{\epsilon}}^{e}\) .On obtient l’expression suivante pour \(\frac{\partial {f}_{\epsilon}^{e}}{\partial {\dot{\epsilon}}^{e}}\) :

(3901)\[\frac{\partial {f}_{\epsilon}^{e}}{\partial {\dot{\epsilon}}^{e}}=\frac{2\mu \dot{{\lambda}_{n}}}{{q}_{n}}\left(M-\Psi \otimes \Psi \right)\]

Avec l’expression suivantepour le tenseur d’ordre 4 \(M\) :

(3902)\[M=\frac{3}{2}\left({I}_{4}-\frac{1}{3}I\otimes I\right)\]

\({I}_{4}\) étant le tenseur identité d’ordre quatreet \(\otimes\) le produit tensoriel.Le calcul de \(\frac{\partial {f}_{\dot{{\lambda}_{n}}}}{\partial \dot{{\epsilon}^{e}}}\) pour chaque mécanisme \(n\) est plus délicat, car il dépend si le mécanisme est actif ou pas pendant les itérations de Newton.

Si le mécanisme n’est pas actif (\({f}_{\dot{{\lambda}_{n}}}<0\) ):

(3903)\[\frac{\partial {f}_{\dot{{\lambda}_{n}}}}{\partial {\dot{\epsilon}}^{e}}=0\]

Si le mécanisme estactif (\({f}_{\dot{{\lambda}_{n}}}>0\) ):

(3904)\[\frac{\partial {f}_{\dot{{\lambda}_{n}}}}{\partial {\dot{\epsilon}}^{e}}=-\frac{4}{3}\frac{\mu}{{C}_{n}}{\Psi}_{n}:{I}_{4}\]

En vérifiant les valeurs obtenues pour la jacobienne numérique (option @CompareToNumericalJacobian), il a été observée que pour des valeurs de \({f}_{\dot{{\lambda}_{n}}}\) proches de zéro (positives ou négatives), \(\frac{\partial {f}_{\dot{{\lambda}_{n}}}}{\partial {\dot{\epsilon}}^{e}}\) vaut la moitié de l’expression obtenue analytiquement. Le traitement dans ce cas est fait en considérant une valeur critique de la valeur de la surface de charge de \({10}^{-6}\) .

Implantation dans code_aster

Le nombre de surfaces d’écrouissage cinématique du modèle est fixé à douze. Les premières onze surfaces permettent de couvrir la gamme de déformations de cisaillement entre \({10}^{-5}\) et \(2,0\times {10}^{-2}\) . On considère un comportement linéaire élastique pour des déformations plus petites que \({10}^{-5}\) et une déformation de cisaillement maximale de \(2,0\times {10}^{-2}\) . Les valeurs d’interpolation utilisées sont les suivantes:

γ1=1.00000000e-05

γ2=2.15443469e-05

γ3=4.64158883e-05

γ4=1.00000000e-04

γ5=2.15443469e-04

γ6=4.64158883e-04

γ7=1.00000000e-03

γ8=2.15443469e-03

γ9=4.64158883e-03

γ10=1.00000000e-02

γ11=2.00000000e-02

Au-delà d’une déformation de cisaillement de \(2,0\times {10}^{-2}\) , les résultats du modèle seront approximatifs car la dernière surface, positionnée à une déformation de cisaillement de \(1,0\times {10}^{-1}\) , permet simplement de faire en sorte que le modèle fournisse une réponse entre \(2,0\times {10}^{-2}\) et \(1,0\times {10}^{-1}\) .

Variables internes

Les variables internes (StateVariableet AuxiliaryStateVariabledans le langage MFront) suivantes sont disponibles lors d’un calcul avec la loi de comportement de Iwan:

Eel(1-6)	Tenseur de déformations élastiques
pp (7-18)	Vecteur de multiplicateurs plastiques scalaires \(>\)
X (19-91)	Vecteur de tenseurs d’écrouissage cinématique \({\text{X}}_{\text{n}}\)
fn(92-103)	Vecteur de valeurs de la surface de charge

Descriptiondu modèle sous MFront

Le comportement est défini dans le fichier Iwan.mfront.

Parser/DSL	Implicit
Algorithm	NewtonRaphson
	@Theta 1. @IterMax 50 @Epsilon 1.E-12
Variables internes (@StateVariable)	real pp[1:12]
Variables internes auxilliaires (@AuxiliaryStateVariable)	Stensor X[1:12] real fn[1:12]
Variables de commandes (@ExternalStateVariable)	aucune
Modélisations	“3D” “AXIS” “D_PLAN”
Déformations	“PETIT” “PETIT_REAC” “GDEF_LOG”

Fonctionnalités et vérification

Voici la liste des cas de vérification disponibles:

ssnv205b	[V6.04.205]	essai de cisaillement cyclique contrôlé en déformation.
ssnv207b	[V6.04.207]	cisaillement cyclique contrôlé en déformationavec micro-décharge.
COMP012 d	[V6.07.112]	Essai de cisaillement cyclique avec la commande CALC_ESSAI_GEOMECA

Bibliographie

[bib1]

IWAN W. D., A distributed element model for hysteresis and its steady state dynamic response, Journal of Applied Mechanics , 88 (14), pp. 893-900, 1966.

[bib2]

IWAN W. D., On a class of models for the yielding behaviour of continuous and composite systems, Journal of Applied Mechanics , 89 (13), pp. 612-617, 1967.

[bib3]

GONDOMZADEH A.: Interaction dynamique sol-structure: influence de non-linéarités de comportement du sol, 2011.

[bib4]

«Relation de comportement élastoplastique à écrouissage cinématique linéaire et isotrope non linéaire. Modélisations 3D et contraintes planes». Documentation de référence de Code_Aster [R5.03.16].