r7.01.43 Loi élastoplastique avec critère de Mohr-Coulomb lissé#
Résumé:
Ce document présente la méthode de résolution de la loi élastoplastique non associée avec surface de charge lissée de Mohr-Coulomb dédiée particulièrement aux géo-matériaux granulaires, accessible via l’environnement MFront . Le critère de Mohr-Coulomb d’origine est une pyramide: les arêtes et le sommet constituent des difficultés pour l’algorithme de résolution, cf. [R7.01.28]. On adopte ici la formulation en invariants et la solution de lissage (régularisation) proposée par Abbo et Sloan dans [ 2 ], avec une méthode d’intégration implicite directe en temps. Elle est implémentée dans l’environnement MFront , cf. [ 4 ], en utilisant l’algorithme développé par [ 3 ] au sein de la communauté http://www.opengeosys.org ` <http://www.opengeosys.org/>`_.
Formulation théorique#
Variables d’état, loi d’état#
Les variables d’état sont: la déformation totale \(\epsilon \left(\mathrm{u}\right)\) , la déformation plastique \({\epsilon}^{p}\) .
La densité d’énergie libre s’écrit:
où \(\mathrm{A}\) désigne le tenseur d’élasticité linéaire; on le suppose isotrope (2 coefficients élastiques indépendants en 3D). On ne considère pas d’écrouissage.
La loi de comportement tridimensionnelle élastoplastique s’écrit alors (loi d’état):
Critère (surface de charge)lissé#
Le critère délimitant le domaine élastique (surface de charge) est formulé avec les 3 invariants tensoriels:
avec \({I}_{1}=\text{tr}\left(\sigma \right)=3p\) le triple de la pression moyenne, \({J}_{2}=\frac{1}{2}\left({\sigma}^{D}:{\sigma}^{D}\right)=\frac{1}{3}{\left({\Vert \sigma \Vert }_{\mathit{VM}}\right)}^{2}\) ; \({J}_{3}=\det\left({\sigma}^{D}\right)\) , ayant noté \({\sigma}^{D}=\sigma -p\mathrm{Id}\) le déviateur du tenseur des contraintes, et l’angle de Lode modifié défini par: \({\theta}_{L}=\frac{1}{3}\arcsin\left(-\frac{3\sqrt{3}{J}_{3}}{2\sqrt{{J}_{2}^{3}}}\right)\phantom{\rule{2em}{0ex}}\in \phantom{\rule{2em}{0ex}}\left[-\frac{\pi}{6},\frac{\pi}{6}\right]\) , et les trois paramètres: la cohésion \(c\) , l’angle de frottement \({\varphi}_{\mathit{pp}}\) et la «troncature de traction» \({a}_{\mathit{tt}}\) , qui sert à l’approximation hyperbolique utilisée pour le lissage au sommet de la pyramide de Mohr-Coulomb. cf. [ 2 ].
Remarque 1 :
Selon [ 2 ], si \({a}_{\mathit{tt}}=0\) , le critère (3) retrouve le sommet de la pyramide de Mohr-Coulomb, donc il n’est plus lissé. Si \({a}_{\mathit{tt}}=\frac{1}{2}c.\mathit{cotan}{\varphi}_{\mathit{pp}}\) , alors au sommet du critère lissé, on vérifie: \({I}_{1}=0\) .
On sait que les contraintes principales peuvent être écrites selon:
La fonction \(K\left({\theta}_{L}\right)\) dans (3)est définie par:
avec un paramètre supplémentaire \({\theta}_{T}\in \phantom{\rule{2em}{0ex}}\left[0,\frac{\pi}{6}\right]\) appelé angle de transition et lesfonctions \(A\left({\theta}_{L}\right),B\left({\theta}_{L}\right),C\left({\theta}_{L}\right),\) définiespar:
Remarque 2 :
Selon [ 2 ], on ne doit pas prendre une valeur de \({\theta}_{T}\) trop proche de \(\frac{\pi}{6}\) , au risque d’un mauvais conditionnement; et une valeur typique de \({\theta}_{T}\) peut être prise égale à \(\frac{5\pi }{36}\) *.* Avec cette valeur de \({\theta}_{T}\) et avec \({a}_{\mathit{tt}}=\frac{1}{20}c.\mathit{cotan}{\varphi}_{\mathit{pp}}\) , le lissage hyperbolique proposé conduit à une erreur au maximum égale à 0,13% dans le plan déviatorique par rapport à la pyramidede Mohr-Coulomb d’origine.
Les et sont une représentation graphique de la surface de charge de Mohr-Coulomb originale dans l’espace des contraintes principales. La montre l’effet du lissage [ 2 ]sur les arêtes de la surface de charge de Mohr-Coulomb.
Figure 2.2-1: Représentation de la surface de charge de Mohr-Coulomb dans l’espace tridimensionnel des contraintes principales.
Figure 2.2-2: Représentation de la surface de charge de Mohr-Coulomb dans le plan :math:`pi ` des déviateurs des contraintes (tout vecteur représenté dans ce plan correspond à une contrainte déviatorique).
Figure 2.2-3: Représentation de la surface de charge de Mohr-Coulomb lissée dans le plan \(\pi\) des déviateurs des contraintes, pour l’ angle de frottement \({\varphi}_{\mathit{pp}}=25°\) et deux valeurs de \({\theta}_{T}\) (tirée de [ 3 ]).
Un enrichissement du critère () est apporté en ajoutant un terme d’écrouissage à la cohésion, sous la forme suivante:
où le paramètre \({r}_{\mathit{eg}}\) (exprimé sans unités) permet d’introduire un écrouissage linéaire de cohésion.
Potentiel d’écoulement plastique#
L’écoulement plastique de cette loi n’est passtandard (loi non associée), donc non orienté selon la normale à la surface de charge. Les auteurs [ 2 ]proposent une expression de l’écoulement plastique à l’aide d’un autrepotentiel de même forme analytique que celle du critère (), en faisant intervenir l’angle de dilatance \(\psi\) :
Lorsque \(\psi ={\varphi}_{\mathit{pp}}\) , la loi d’écoulement plastique devient associée.
La fonction \({K}_{g}\left({\theta}_{L}\right)\) dans (7)est définie par:
lesfonctions \({A}_{\text{g}}\left({\theta}_{L}\right),{B}_{\text{g}}\left({\theta}_{L}\right),{C}_{\text{g}}\left({\theta}_{L}\right),\) étantdéfiniespar:
Lorsque l’état de contrainte se trouve sur la surface de charge \(\text{f}\left(\sigma \right)=0\) , le taux de déformation plastique est normal au potentiel \(\text{g}\left(\sigma \right)\) et s’écrit avec lemultiplicateur plastique \(d\lambda\) :
où le tenseur d’ordre deux normal extérieur à la surface est:
La seule vraie variable interne thermodynamique est la déformation plastique équivalente \({\Vert {\epsilon}^{p}\Vert }_{\mathit{VM}}\) .
On a naturellement \(\frac{\partial \text{g}\left(\sigma \right)}{\partial {I}_{1}\left(\sigma \right)}=\frac{1}{3}\cdot \sin\psi\) et on rappelle que: \(\frac{\partial {J}_{2}\left(\sigma \right)}{\partial \sigma }={\sigma}^{D}\) .
Remarque 3 :
L e potentiel d’écoulement choisi, de même forme que le critère, donc le tenseur d’écoulement \({\mathrm{n}}_{\text{g}}\) , n’est pas affecté par la possible introduction d’un écrouissage linéaire de cohésion, cf. () *.* En effet, l’évolution de la cohésion ne modifie pas la direction normale aux faces de la pyramide \(\text{g}\left(\sigma \right)=0\) .
Dissipationd’énergie par écoulement plastique#
La densité de puissance mécanique dissipée par le mécanisme de plasticité est:
La dissipation d’énergie mécanique s’obtient par cumul au cours du trajet de chargement.
Paramètres de la relation de comportement#
Les paramètres de la relation de comportement sont:
YoungModulus |
Module de Young (positif) |
\(E\) |
PoissonRatio |
Coefficient de Poisson |
\(\nu\) |
Cohesion |
Cohésion (positif ou nul) |
\(c\) |
FrictionAngle |
Angle de frottement (fourni en °) |
\({\varphi}_{\mathit{pp}}\) |
DilatancyAngle |
Angle de dilatance (fourni en °) |
\(\psi\) |
TransitionAngle |
Angle de Lode de transition (fourni en °, inférieur à 30 °) |
\({\theta}_{T}\) |
TensionCutOff |
Troncature de traction (positive ou nulle) |
\({a}_{\mathit{tt}}\) |
HardeningCoef |
Écrouissage de cohésion (positif ou nul) |
\({r}_{\mathit{eg}}\) |
Intégration numérique de la loi de comportement#
Le choix dans MFront est porté sur l’intégration locale des systèmes d’équations différentiels non linéaires du premier ordre écrits de la manière suivante:
Remarque 4 :
On doit noter que l’interface nécessite la fourniture des composantes des déformations élastiques à l’état initial. Par défaut elles sont supposées nulles.
Intégration implicite directe#
Une intégration implicite de type \(\theta\) -méthode est choisie dans [ 4 ] pour l’intégration numérique du système non linéaire; de fait on prend \(\theta =1\) (schéma totalement implicite de manière à respecter exactement le critère à la fin du pas de temps).
La prédiction est élastique linéaire:
puis on calcule les valeurs de prédiction des troisinvariants tensoriels \({I}_{1}\) , \({J}_{2}\) , \({J}_{3}\) , et onévalue le critère (), cf. § 2.2 .
Si l’on n’est pas en charge plastique, alors \(d\lambda =0\) et :math:`{sigma}^{+}={sigma}^{-}+mathrm{A}.depsilon ` . Sinon, on entame la résolution du zéro des équations non linéaires.
La méthode itérative de Newton-Raphson est choisie pour le calcul du zéro des équations non linéaires exprimant le respect de l’équation de cohérence \(\dot{\text{f}}=0\) , pour le critère () cf. § 2.2 et l’expression de la relation de comportement () avec la direction d’écoulement plastique, cf. § 2.3 , fournissant le taux de multiplicateur plastique \(d\lambda\) et le taux de contraintes en utilisant la prédiction élastique \({\sigma}^{e}\) ; la matrice jacobienne réactualisée est utilisée à chaque itération locale.
La contrainte finale \({\sigma}^{+}\) est calculée par:
La tolérance choisie de l’intégration locale non linéaire est fournie par le mot-clé RESI_INTE_RELA dans le mot-clé facteur RELATION des opérateurs de calcul non linéaire. Cette valeur est transmise dans la variable @Epsilonde MFront , dont la valeur par défaut est: \({10}^{-14}\) . Il est conseillé de ne pas adopter de valeur supérieure à \({10}^{-10}\) . Le mot-clé ITER_INTE_MAXI dans le mot-clé facteur RELATION est également affecté à la variable @iterMax.
Leparamètre numérique spécifique positifnoté \({r}_{\mathit{eg}}\) (exprimé sansunité) est fourni avec les paramètres de lissage: l’angle de Lode de transition \({\theta}_{T}\) et la troncature de traction \({a}_{\mathit{tt}}\) . Il ajoute une contribution de type écrouissage au critère, qui sera choisie faible devant \(1\) . Le terme \(-c.\cos{\varphi}_{\mathit{pp}}\) dans le critère ( ) devient, cf.( ) :
Après convergence de l’intégration locale la mise à jour des variables est réalisée.
Mise à jour à la fin du pas d’intégration#
Une fois les équations non linéaires résolues on a directement la mise à jour du tenseur contrainte finale \({\sigma}^{+}\) .
On calcule aussi la déformation plastique équivalente EquPlasticStrain, égale à \({\lambda}^{+}\) , la déformation plastique volumique VolPlasticStrain, la densité volumique de dissipation mécanique cumulée Dissipation.
Expression de la matrice tangente cohérente#
La matrice tangente cohérente analytique du système non linéaire à résoudre par la méthode de Newton utilise les dérivées des diverses grandeurs définies pour exprimer l’équation de cohérence \(\dot{\text{f}}=0\) pour le critère (), cf.§ 2.2 et pour exprimer la relation de comportement () avec la direction d’écoulement plastique, cf. § 2.3 , par rapport à \(\lambda\) (en tenant compte aussi de ()), et par rapport à \(\sigma\) . La matrice tangente cohérente du système non linéaire est ainsi définie par:
La dimensionnalisation de la ligne associée à l’équation de cohérence utilise le module de Young afin d’avoir une dimension de l’ordre des contraintes et ainsi mieux conditionner cette matrice.
On voit aussi par là l’avantage de choisir le paramètre d’écrouissage \({r}_{\mathit{eg}}\) non nul qui permet d’ajouter une contribution non nulle sur la diagonale de cette matrice.
Étapes MFront#
Le comportement est défini dans le fichier MohrCoulombAS.mfront .
Nommage |
nom = “MohrCoulombAS”, symbol_mfront = “MohrCoulombAS”, |
Parser/DSL |
Implicit |
Algorithme |
NewtonRaphson algo_inte = (“NEWTON”,) |
Paramètres schéma |
@Theta 1. @IterMax 30 @Epsilon 1.E-14 |
Modélisations |
“3D”, “AXIS”, “D_PLAN” |
Déformations |
“PETIT” “PETIT_REAC” “GDEF_LOG” |
Matrice tangente |
syme_matr_tang = (“No”,) |
Variables localesMFront#
eel |
Tenseur de déformations élastiques |
eto |
Tenseur de déformations totales |
sig |
Tenseur des contraintes |
sig_el s_el |
Tenseur des contraintes élastiques et son déviateur |
dt |
incrément de pas de temps |
@StateVariable:lam |
Equivalent Plastic Strain |
@MaterialProperty |
Young: YoungModulus nu: PoissonRatio C:Cohesion phi: FrictionAngle psi: DilatancyAngle lodeT: TransitionAngle a: TensionCutOff reg: HardeningCoef |
F, Fy |
critère |
n |
Stensor: direction d’écoulement plastique |
nf |
Dérivée du critère par rapport à sig |
dfeel_ddeel dfeel_ddlam dflam_ddeel dflam_ddlam |
Termes de la matrice Jacobienne |
dlam |
Incrément de multiplicateur plastique |
epv du ip |
VolPlasticStrain Dissipation PlasticIndex |
Variables internes de la loi de comportement#
Voici la liste des variables internes stockées. On indique également leur nommage.
\(V1\) à \(V6\) en 3D \(V1\) à \(V4\) en 2D |
Composantes du tenseur des déformations élastiques (demandées par MFront ) |
\(V7\) en 3D (\(V5\) en 2D) |
Déformation plastique équivalente \({\Vert {\epsilon}^{p}\Vert }_{\mathit{VM}}\) : vraie variable interne thermodynamique: EquPlasticStrain |
\(V8\) en 3D (\(V6\) en 2D) |
Déformation plastique volumique \(\text{tr}{\epsilon}^{p}\) : VolPlasticStrain |
\(V9\) en 3D ( \(V7\) en 2D) |
Densité d’énergie mécanique dissipée \(D(t)={\int}_{0}^{t}\sigma :{d\epsilon }^{p}\) : Dissipation |
\(V10\) en 3D (\(V8\) en 2D) |
Indicateur d’activation de la plasticité (1) ou non (0): PlasticIndex |
Vérification et validationde la loi de comportement#
Voici la liste des cas de vérification et de validationdisponibles:
intitulé |
documentation |
Titre et objet |
SSNV232 |
essai triaxial drainé monotone (point matériel). La solution est comparée à une solution analytique. |
|
SSNV233 |
essai de torsiondrainé monotone (point matériel). |
|
WTNV142 |
essai triaxial non drainé monotone (point matériel). |
|
COMP012 |
essai triaxial drainé monotone à déplacement imposé (point matériel). La solution est comparée à une solution analytique. Avec la macro-commande CALC_ESSAI_GEOMECA. |
|
COMP001 |
[V6.07.1 01 ` <https://www.code-aster.org/V2/doc/default/fr/man_v/v6/v6.07.112.pdf>`_ ] |
Vérification de comportements élastoplastiques sur point matériepour unchangement d’unités et après rotation du référentiel. |
SSNP104 |
en chargement monotone: semelle rigide posée sur un sol homogène. Comparaison avec les résultats obtenus avec la loi élastoplastique de Mohr-Coulomb d’origine et des résu.ltats analytiques. Versions vérification et validation |
Bibliographie#
[R7.01.28] Loi de Mohr-Coulomb.
Abbo, A. J. and Sloan, S. W. A smooth hyperbolic approximation to the Mohr-Coulomb yield criterion . Computers & Structures. February 1995. Vol. 54, no. 3, p. 427–441. DOI 10.1016/0045-7949(94)00339-5 .
Nagel, T., Minkley, W., Böttcher, N., Naumov, D., Görke, U.-J. and Kolditz, O. Implicit numerical integration and consistent linearization of inelastic constitutive models of rock salt . Computers & Structures. April 2017. Vol. 182, p. 87–103. DOI 10.1016/j.compstruc.2016.11.010 . Voir aussi: http://tfel.sourceforge.net/MohrCoulomb.html.
Invariant-based implementation of the Mohr-Coulomb elasto-plastic model in OpenGeoSys using MFront : http://tfel.sourceforge.net/MohrCoulomb.html
NGUYEN Q.S. – Stabilité et mécanique non linéaire – Hermès, Paris, 2000.