r5.04.03 Modélisations second gradient

Résumé

On présente ici les modélisations second gradient issue des travaux de Chambon et al ([1], [2], [3]) et second gradient de dilatation issue des travaux de Fernandes et al ([4], [5]), et de Gantier ([14]).

Ces modélisations s’inscrivent dans le cadre des milieux à microstructure et ont comme objectif numérique de donner des résultats convergents par rapport à la finesse de la discrétisation spatiale pour éviter d’obtenir des solutions localisées. Elles doivent être utilisées dès lors que les lois rhéologiques mécaniques considérées présentent un adoucissement du comportement traduisant l’endommagement ou la dégradation d’un matériau avant fissuration.

La modélisation second gradient de dilatation est, en fait, une approche simplifiée du second gradient restreinte aux matériaux assujettis au phénomène de dilatance. L’objectif numérique est de réduire de façon significative les temps de calcul. Elle se prête particulièrement bien aux géomatériaux.

Remarque: La modélisation second gradient a été retirée du code. Seule la modélisation gradient de dilatation est disponible.

Table des matières

Le modèle second gradient de dilatation

Le principe des milieux à microstructure est basé sur la prise en compte des déformations microscopiques pour introduire dans l’expression du modèle une longueur interne. Si l’efficacité de la régularisation apportée par le modèle second gradient ne fait aucun doute, les temps de calcul des simulations peuvent devenir prohibitifs. Pour les diminuer dans le cadre particulier des matériaux dilatants (pour lesquels la déformation volumique évolue en fonction du chargement) on restreint la généralité de l’effet régularisant aux milieux à microstructure dilatants.

Les milieux à microstructure dilatants

La cinématique de ces milieux est définie par le champ de déplacement habituel \({u}_{i}\) , la variation volumique microscopique notée \(\chi\) et ses gradients. Par dualité à cette cinématique enrichie, sont introduites les variables statiques des contraintes macroscopiques classiques \({\sigma}_{ij}\) , la contrainte microscopique de dilatation \(\kappa\) et les doubles contraintes vectorielles de dilatation \({S}_{j}\) . La variable \(\kappa\) est la composante conjuguée de la déformation volumique relative (du champ macroscopique par rapport au microscopique) \({\epsilon}_{V}-\chi\) , tandis que les composantes \({S}_{j}\) définissent un vecteur qui est le conjugué du gradient de la dilatation microscopique \(\frac{\partial \chi }{\partial {x}_{j}}\) .

De façon analogue à l’expression de la puissance virtuelle des milieux à microstructure du chapitre 1.1 , et en négligeant de nouveau l’expression des efforts extérieurs de volume pour simplification analytique, on trouve pour les milieux à microstructure dilatants, pour tout champ cinématiquement admissible \(({u}_{i}^{\text{*}},{\chi }^{\text{*}})\)

(2888)\[\underset{\Omega}{\int}({\sigma}_{ij}\frac{\partial {u}_{i}^{\text{*}}}{\partial {x}_{j}}+\kappa ({\epsilon}_{V}^{\text{*}}-{\chi }^{\text{*}})+{S}_{j}\frac{\partial {\chi }^{\text{*}}}{\partial {x}_{j}})\text{dv}=\underset{\partial \Omega }{\int}({t}_{i}{u}_{i}^{\text{*}}+{\mathrm{m\chi }}^{\text{*}})\text{ds}\]

Pour laquelle

(2889)\[{t}_{i}=\left({\sigma}_{ij}+{\text{κδ}}_{ij}\right){n}_{j}\]

(2890)\[m={S}_{j}{n}_{j}\]

sont les conditions aux limites, exprimées sur la frontière \(\partial \Omega ` , conjuguées par dualité aux variables cinématiques :math:`{u}_{i}\) et \(\chi\) respectivement.

Les relations d’équilibre de ce problème s’écrivent

(2891)\[\frac{\partial {\sigma}_{ij}}{\partial {x}_{j}}+\frac{\partial \kappa }{\partial {x}_{i}}=0\]

(2892)\[\kappa +\frac{\partial {S}_{j}}{\partial {x}_{j}}=0\]

Le système d’équations composé de [(), (), () et ()] est obtenu classiquement par application du théorème de la divergence et par une intégration par partie de ().

Le modèle second gradient de dilatation

Suivant un principe analogue à celui mis en œuvre pour le modèle second gradient du chapitre 1.2 , on introduit une contrainte mathématique pour forcer l’égalité entre les déformations volumiques macroscopique \({\epsilon}_{V}\) et microscopique \(\chi\)

(2893)\[\chi ={\epsilon}_{V}\]

L’expression des puissances virtuelles est alors la suivante, pour tout champ cinématiquement admissible \({u}_{i}^{}\)

(2894)\[\underset{\Omega}{\int}({\sigma}_{ij}\frac{\partial {u}_{i}^{\text{*}}}{\partial {x}_{j}}+{S}_{j}\frac{{\partial}^{2}{u}_{i}^{\text{*}}}{\partial {x}_{i}\partial {x}_{j}})\text{dv}=\underset{\partial \Omega }{\int}({p}_{i}{u}_{i}^{\text{*}}+{\text{Pn}}_{i}{\text{Du}}_{i}^{\text{*}})\text{ds}\]

où \({p}_{i}\) et \(P\) sont les conditions aux limites définies par

(2895)\[{p}_{i}={\sigma}_{ij}{n}_{j}-{n}_{i}{n}_{j}{\text{DS}}_{j}-\frac{{\text{DS}}_{j}}{{\text{Dx}}_{j}}{n}_{i}-\frac{{\text{DS}}_{j}{n}_{j}}{{\text{Dx}}_{i}}+\frac{{\text{Dn}}_{p}}{{\text{Dx}}_{p}}{S}_{j}{n}_{j}{n}_{i}\]

(2896)\[P={S}_{j}{n}_{j}\]

Comme pour le modèle second gradient, l’hypothèse sur l’égalité entre champ de déformations volumiques microscopique et macroscopique () a un impact direct sur l’expression des conditions aux limites car les variables \({u}_{i}^{\text{*}}\) et \({\chi }^{\text{*}}\) ne sont plus indépendantes.

Pour des raisons de simplicité, on suppose que \(P=0\) . La conséquence de cette hypothèse est que \({S}_{j}{n}_{j}=0\) sur la frontière, ce qui réduit l’expression () à

(2897)\[{p}_{i}={\sigma}_{ij}{n}_{j}-\frac{\partial {S}_{j}}{\partial {x}_{j}}{n}_{i}\]

Une propriété remarquable de cette simplification vient du fait que l’expression des conditions aux limites se décompose alors en la partie classique \({\sigma}_{ij}{n}_{j}\) et un second terme qui n’induit pas de composantes de cisaillement.

L’équation d’équilibre s’écrit:

(2898)\[\frac{\partial {\sigma}_{ij}}{\partial {x}_{j}}-\frac{{\partial}^{2}{S}_{j}}{\partial {x}_{i}\partial {x}_{j}}=0\]

Discrétisation spatiale par éléments finis

Discrétisation historique du modèle de second gradient de dilatation

Le but de la démarche du second gradient de dilatation est de définir un modèle régularisant, permettant d’assurer l’indépendance des résultats par rapport aux discrétisations spatiales, en introduisant un minimum d’inconnues nodales dans son approche par éléments finis. Pour cela on se limite aux applications considérant des matériaux dilatants. Or, discrétiser par la méthode des éléments finis l’expression () a comme conséquence d’imposer que le champ des inconnues de déplacement ainsi que sa divergence soient continus et dérivables. Cela revient a prendre en compte des éléments finis C1-continus.

On propose alors d’introduire la contrainte mathématique () dans l’expression des milieux à microstructure dilatants au moyen d’un couplage de multiplicateurs de Lagrange et de pénalisation. On obtient alors pour tout champ cinématiquement admissible \(\left({u}_{i}^{\text{*}},{\chi }^{\text{*}},{\lambda}^{\text{*}}\right)\)

(2899)\[\begin{split}\begin{array}{c}\underset{\Omega}{\int}\left({\sigma}_{ij}\frac{\partial {u}_{i}^{\text{*}}}{\partial {x}_{j}}+{S}_{j}\frac{\partial {\chi }^{\text{*}}}{\partial {x}_{j}}-\lambda \left({\epsilon}_{V}^{\text{*}}-{\chi }^{\text{*}}\right)+{\lambda}^{\text{*}}\left({\epsilon}_{V}-\chi \right)+r\left({\epsilon}_{V}-\chi \right)\left({\epsilon}_{V}^{\text{*}}-{\chi }^{\text{*}}\right)\right)\text{dv}\\ \phantom{\rule{2em}{0ex}}=\underset{\partial \Omega }{\int}\left({p}_{i}{u}_{i}^{\text{*}}+{\text{Pn}}_{i}{\text{Du}}_{i}^{\text{*}}\right)\text{ds}\end{array}\end{split}\]

cette formulation historique est disponible sous le nom de DIL. La mise en œuvre numérique proposée est détaillée au chapitre 3 . On précise ici uniquement les espaces d’approximation des variables définissant le champ cinématique. Les interpolations polynomiales sont P2-P1-P1(en déformations planes et en 3D):

1. des fonctions de formes du second ordre pour les variables de déplacements \(u\)

2. des fonctions de formes du premier ordre pour le tenseur des déformations volumiques microscopique \(\chi\) et les multiplicateurs de Lagrange \(\lambda\)

Figure 2.3-a: Discrétisation élément fini du modèle second gradient de dilatation pour le 3D.

Modélisation incompressible du modèle de second gradient de dilatation

Cette modélisation proposée par Gantier [14] est une variation de la modélisation présentée au paragraphe adaptée aux comportements où la dilatance tend vers zéro avec l’écrouissage. Cette modélisation s’inspire des éléments incompressibles à trois champs (voir R3.06.08). Elle n’est pas disponible pour les modélisations de THM avec second gradient.

(2900)\[\begin{split}\begin{array}{c}\underset{\Omega}{\int}\left({\sigma}_{ij}\left({{\varepsilon}}_{ij}^{D}(u\text{*})+\frac{1}{3}\chi \text{*}{\delta}_{ij}\right)+{S}_{j}\frac{\partial {\chi }^{\text{*}}}{\partial {x}_{j}}+\lambda \left({\epsilon}_{V}^{\text{*}}-{\chi }^{\text{*}}\right)+{\lambda}^{\text{*}}\left({\epsilon}_{V}-\chi \right)+r\left({\epsilon}_{V}-\chi \right)\left({\epsilon}_{V}^{\text{*}}-{\chi }^{\text{*}}\right)\right)\text{dv}\\ \phantom{\rule{2em}{0ex}}=\underset{\partial \Omega }{\int}\left({p}_{i}{u}_{i}^{\text{*}}+{\text{Pn}}_{i}{\text{Du}}_{i}^{\text{*}}\right)\text{ds}\end{array}\end{split}\]

Dans cette expression, \({{\varepsilon}}^{D}\) représente la déformation déviatorique et \({\delta}_{ij}\) le symbole de Kronecker.

Cette formulation est donc inspirée des éléments incompressibles à trois champs: c’est la variable \(\chi\) qui intervient dans la loi de comportement de premier gradient.

Par ailleurs, il a été démontré (voir Gantier [14]) que l’utilisation du terme de pénalisation est optionnelle: l’utilisation d’un multiplicateur de Lagrange est suffisante pour imposer l’égalité entre déformations volumiques microscopique et macroscopique.

Les fonctions de forme utilisées sont identiques à celles de la formulation historique du second gradient. Elles sont de type P2-P1-P1:

1. des fonctions de formes du second ordre pour les variables de déplacements \(u\)

2. des fonctions de formes du premier ordre pour le tenseur des déformations volumiques microscopique \(\chi\) et les multiplicateurs de Lagrange \(\lambda\)

Intégration numérique des modèles second gradient

On détaille ci-dessous les champs de déformations, de contraintes associés et la matrice tangente en mécanique avec la formulation historique du modèle de second gradient de dilatation. Il n’y a aucune variable interne spécifique au second gradient. Pour rappel, en mécanique, les variables nodales définissant les degrés de liberté sont les suivantes: \(\left({u}_{i},\chi ,\lambda \right)\) .

Pour les modélisations de THM avec second gradient, il n’existe pas de couplages entre les déformations et contraintes généralisées de second gradient et celles de thermo-hydrique. Ainsi, les champs de déformations et de contraintes généralisées, ainsi que les matrices de rigidités élémentaires associés au second gradient sont identiques à ceux des modélisations demécanique.

\(E=\left(\begin{array}{c}\epsilon \\ {\epsilon}_{V}\\ \chi \\ \frac{\partial \chi }{\partial {x}_{j}}\\ \lambda \end{array}\right)\)

\(\Sigma =\left(\begin{array}{c}\sigma \\ r\left({\epsilon}_{V}-\chi \right)+\lambda \\ -r\left({\epsilon}_{V}-\chi \right)-\lambda \\ {S}_{j}\\ \left({\epsilon}_{V}-\chi \right)\end{array}\right)\)

La matrice tangente élémentaire de la modélisation second gradient est composée, entre autres:

1. de la matrice tangente de rigidité élémentaire \({D}^{1d}\) associée à la loi de comportement de premier gradient

2. de la matrice tangente de rigidité élémentaire \({D}^{2d}\) associée à la loi de comportement second gradient de dilatation qui lie les doubles contraintes \({S}_{j}\) aux gradients des déformations volumiques microscopiques \(\frac{\partial \chi }{\partial {x}_{j}}\) .

\({K}^{\text{el}}=\left(\begin{array}{ccccc}{\left({D}^{1d}\right)}_{\text{dim}\left(j\right)}& 0& 0& 0& 0\\ 0& r& -r& 0& 1\\ 0& -r& r& 0& -1\\ 0& 0& 0& {\left({D}^{2d}\right)}_{\text{dim}\left(j\right)}& 0\\ 0& 1& -1& 0& 0\end{array}\right)\)

Conseils/Procédure pour l’utilisation des modèles second gradient

La démarche ci-dessous est la même quelle que soit la modélisation locale (3D, D_PLAN, en mécanique ou en THM). Des exemples sont disponibles dans la base des cas-tests. On peut citer, entre autres, les essais biaxiaux en élasticité linéaire en compression référencés ssll117, la construction d’une colonne de sol en non-linéaire (loi de comportement Hujeux en milieu poreux couplé hydromécanique) référencée wtnv132 ou l’essai triaxial en élasticité linéaire référencé sslv117.

L’élasticité linéaire second gradient

Pour utiliser les modélisations second gradient il faut définir deux lois de comportements pour décrire respectivement les relations macroscopiques entre contraintes \({\sigma}_{ij}\) et déformations \({\epsilon}_{ij}\) et les relations microscopiques entre doubles contraintes (\({\Sigma}_{\text{ijk}}\) pour le second gradient ou \({S}_{j}\) pour le second gradient de dilatation) et gradient de déformations \(\frac{\partial {f}_{ij}}{\partial {x}_{k}}\) pour le second gradient ou gradient de déformation volumique \(\frac{\partial \chi }{\partial {x}_{j}}\) pour le second gradient de dilatation.

Pour le choix du comportement de type premier gradient – celui qui intervient dans la modélisation locale décrite en section 3 – il n’y a aucune restriction: toutes les lois de comportements sont théoriquementpossibles. En pratique, il est cependant nécessaire de se restreindre à des lois de comportements pour lesquelles la variation de volume est non nulle.

En revanche, on dispose actuellement d’une seule relation de comportement de type second gradient. Il s’agit de l’élasticité linéaire proposée par Mindlin ([6], [7]) dont voici l’écriture en 2D

(2901)\[\begin{split}\left[\begin{array}{c}{\Sigma}_{111}\\ {\Sigma}_{112}\\ {\Sigma}_{121}\\ {\Sigma}_{122}\\ {\Sigma}_{211}\\ {\Sigma}_{212}\\ {\Sigma}_{221}\\ {\Sigma}_{222}\end{array}\right]=\left[\begin{array}{cccccccc}{a}^{12345}& 0& 0& {a}^{23}& 0& {a}^{12}& {a}^{12}& 0\\ 0& {a}^{145}& {a}^{145}& 0& {a}^{25}& 0& 0& {a}^{12}\\ 0& {a}^{145}& {a}^{145}& 0& {a}^{25}& 0& 0& {a}^{12}\\ {a}^{23}& 0& 0& {a}^{34}& 0& {a}^{25}& {a}^{25}& 0\\ 0& {a}^{25}& {a}^{24}& 0& {a}^{34}& 0& 0& {a}^{23}\\ {a}^{12}& 0& 0& {a}^{25}& 0& {a}^{145}& {a}^{145}& 0\\ {a}^{12}& 0& 0& {a}^{25}& 0& {a}^{145}& {a}^{145}& 0\\ 0& {a}^{12}& {a}^{12}& 0& {a}^{23}& 0& 0& {a}^{12345}\end{array}\right].\left[\begin{array}{c}{\chi }_{111}\\ {\chi }_{112}\\ {\chi }_{121}\\ {\chi }_{122}\\ {\chi }_{211}\\ {\chi }_{212}\\ {\chi }_{221}\\ {\chi }_{222}\end{array}\right]\end{split}\]

où \({\chi }_{\text{ijk}}=\frac{\partial {f}_{ij}}{\partial {x}_{k}}\) et tous les termes de la matrice dépendent de cinq constantes selon la relation suivante:

\(\lbrace \begin{array}{c}{a}^{12345}=2({a}^{1}+{a}^{2}+{a}^{3}+{a}^{4}+{a}^{5})\\ {a}^{23}={a}^{2}+{\mathrm{2a}}^{3}\\ {a}^{12}={a}^{1}+\frac{{a}^{2}}{2}\\ {a}^{145}=\frac{{a}^{1}}{2}+{a}^{4}+\frac{{a}^{5}}{2}\\ {a}^{25}=\frac{{a}^{2}}{2}+{a}^{5}\\ {a}^{34}=2({a}^{3}+{a}^{4})\end{array}\)

Fernandes et al [5] ont appliqué ce modèle pour le second gradient de dilatation en simplifiant l’expression ()) en 2D par

(2902)\[ \begin{align}\begin{aligned}\begin{split}\left[\begin{array}{c}{S}_{1}\\ {S}_{2}\end{array}\right]=\left[\begin{array}{cc}3{a}^{1}& 0\\ 0& 3{a}^{1}\end{array}\right].\left[\begin{array}{c}\frac{\partial \chi }{\partial {x}_{1}}\\ \frac{\partial \chi }{\partial {x}_{2}}\end{array}\right]\end{split}\\ Attention, en 3D, cette relation devient:\end{aligned}\end{align} \]

(2903)\[\begin{split}\left[\begin{array}{c}{S}_{1}\\ {S}_{2}\\ {S}_{3}\end{array}\right]=\left[\begin{array}{ccc}4{a}^{1}& 0& 0\\ 0& 4{a}^{1}& 0\\ 0& 0& 4{a}^{1}\end{array}\right].\left[\begin{array}{c}\frac{\partial \chi }{\partial {x}_{1}}\\ \frac{\partial \chi }{\partial {x}_{2}}\\ \frac{\partial \chi }{\partial {x}_{3}}\end{array}\right]\end{split}\]

Modélisations disponibles dans Code_Aster

En mécanique, les modélisations disponibles sontD_PLAN_INCO_DIL (déformations planes) et 3D_DIL. Pour ces modélisations, deux formulations sont disponibles: DIL pour la formulation historique, DIL_INCO pour la formulation incompressible.

En THM, les modélisations de second gradient disponibles sont détaillées dans les documentations U3.13.08 et U3.14.07.

Estimation du paramètre de régularisation A1 (2nd gradient de dilatation)

Dans le but de caractériser au mieux la valeur à affecter au paramètre de régularisation, A1, du 2nd gradient de dilatation, il est possible de déterminer une borne supérieure à ce paramètre en fonction de la matrice tangente en vitesse de tout modèle de comportement adoucissant.

De cette façon, l’utilisateur pourra déterminer en fonction de la taille des mailles du problème qu’il traite la valeur du paramètre A1 la mieux adaptée a son problème.

Ce calcul est basé sur un problème analytique \(\mathrm{2D}\) d’une bande de cisaillement [bib2 et bib11].

La rupture dans les géomatériaux se caractérise souvent par la formation de zones à déformations localisées, relatant un passage de zones de déformations homogènes vers des modes de déformations non homogènes. L’apparition de ce phénomène peut être vu sur le plan théorique comme le changement spontané du mode de déformation, assimilé à une bifurcation d’une branche de l’équilibre, c’est à dire l’intersection de deux branches de solutions fonctions des paramètres de contrôle. Dans le cadre des milieux continus, il est possible de dégager sous des conditions restreintes des critères de bifurcation permettant d’identifier le paramètre de contrôle provoquant l’apparition de la bande de cisaillement, ainsi que les orientations potentielles de celle-ci. En revanche, pour les milieux continus classiques, le régime de «post-localisation», notamment la largeur de la bande de cisaillement, ne peut pas être caractérisée. Les résultats expérimentaux montrent néanmoins que cette largeur est un élément intrinsèque aux propriétés du matériau, lié à sa microstructure (forme et taille des grains, Desrues et Viggiani [bib12]) et son état initial (indice des vides et état de contraintes). L’utilisation de la théorie du second gradient local de dilatation permet d’enrichir la cinématique du milieu traduisant ainsi les effets de la microstructure à l’échelle globale. Dans cette section, nous allons extraire d’un problème analytique bidimensionnel les éléments permettant de caractériser la largeur de la bande de cisaillement durant le régime de «post-localisation».

La formulation du problème en vitesse s’écrit alors sous la forme suivante, en considérant les forces de volume constantes:

(2904)\[\underset{\Omega}{\int}{\dot{\sigma}}_{ij}{\dot{\epsilon}}_{ij}({u}^{\text{*}})+{\dot{S}}_{j}{\eta}_{j}({u}^{\text{*}})\mathrm{dv}=\underset{\partial \Omega }{\int}{\dot{F}}_{i}{u}_{i}^{\text{*}}\mathrm{ds}\forall {u}^{\text{*}}\in {\text{V}}_{0}\]

\({\dot{\sigma}}_{ij}\) , \({\dot{S}}_{j}\) et \({\dot{F}}_{i}\) sont les dérivéestemporelles des termes introduits dans l’équation .

Résolution analytique

Le problème analytique bidimensionnel abordé ci-dessous a été résolu par [bib2] pour un modèle de second gradient complet avec le modèle local élastoplastique de Mohr-Coulomb décrit par [bib13]. L’exemple consiste à appliquer le modèle de second gradient de dilatation dans une couche cisaillée et à déterminer les solutions du problème en vitesse.

Figure 4.3-a: Domaine étudié On considère une couche définie dans le plan \((x,y)\) avec \(z\) l’axe perpendiculaire au plan. Chambon et al. [bib2] considèrent une couche infinie comprise entre deux plans définis par \(x=0\) et \(x=l\) (Figure ). Les champs de déplacement sont notés \(u\) dans la direction \(x\) et \(v\) dans la direction \(y\) . \(\dot{u}\) et \(\dot{v}\) sont supposés ne dépendre que de \(x\) et leurs dérivées par rapport à \(x\) seront notées \('\) .

L’état initial du matériau est défini homogène. Les évolutions des conditions aux limites appliquées au domaine sont les suivantes:

• pour \(x=0\) , \(\dot{u}=0\) , \(\dot{v}=0\)

• pour \(x=l\) , \({\dot{F}}_{i}\) sont connues.

Les forces de volume sont supposées constantes.

Les conditions de symétrie appliquées à la frontière \(x=0\) indiquent que le domaine étudié à une largeur totale de \(\mathrm{2l}\) .

Ce problème peut être vu comme l’analyse du comportement d’une bande de localisation, où l’orientation de la bande est supposée et l’état de contraintes est laissé libre de toute restriction. La frontière du domaine en \(x=l\) définit la zone de transition entre la zone de localisation (domaine étudié) et le domaine d’un solide quelconque où les solutions sont supposées régulières et stables.

Le gradient du champ de déplacement dans le domaine prend la forme suivante:

(2905)\[\begin{split}{\dot{u}}_{i,j}=\left[\begin{array}{cc}\dot{u}'& 0\\ \dot{v}'& 0\end{array}\right]\end{split}\]

Le champ de déformation macroscopique s’écrit explicitement:

(2906)\[\begin{split}{\dot{\epsilon}}_{ij}=\frac{1}{2}({\dot{u}}_{i,j}+{\dot{u}}_{j,i})=\left[\begin{array}{cc}\dot{u}'& \frac{1}{2}\cdot \dot{v}'\\ \frac{1}{2}\cdot \dot{v}'& 0\end{array}\right]\end{split}\]

Le seul terme non nul du gradient des déformations volumiques est:

(2907)\[\frac{\partial {\dot{\epsilon}}_{v}}{\partial {x}_{1}}=\dot{u}\text{''}\]

L’application du principe des travaux virtuels au champ de déplacement imposé après 2 intégrations par parties donne les équations d’équilibre suivantes pour le problème en vitesse:

(2908)\[\begin{split}\frac{\partial {\dot{\sigma}}_{ij}}{\partial {x}_{j}}-\frac{{\partial}^{2}{\dot{S}}_{j}}{\partial {x}_{i}\partial {x}_{j}}=0\phantom{a}\Rightarrow \phantom{a}\lbrace \begin{array}{ccc}{\dot{\sigma}}_{11}^{'}-{\dot{S}}_{1}^{\text{''}}& =& 0\\ {\dot{\sigma}}_{12}^{'}& =& 0\end{array}\end{split}\]

Les conditions aux limites fournissent en \(x=l\) :

(2909)\[\begin{split}\lbrace \begin{array}{ccc}{\dot{\sigma}}_{11}-{\dot{S}}_{1}^{'}& =& {\dot{F}}_{1}\\ {\dot{\sigma}}_{12}& =& {\dot{F}}_{2}\\ 3{a}^{1}{\dot{u}}^{'''}& =& 0\end{array}\end{split}\]

Après intégration spatiale des équations d’équilibre en tenant compte des conditions aux limites pour un problème homogène, les solutions du problème respectent les équations suivantes:

(2910)\[\begin{split}\lbrace \begin{array}{ccc}{\dot{\sigma}}_{11}-{\dot{S}}_{1}^{'}& =& {\dot{F}}_{1}\\ {\dot{\sigma}}_{12}& =& {\dot{F}}_{2}\end{array}\end{split}\]

En utilisant les lois constitutives des premier et second gradients, les équations deviennent:

(2911)\[\begin{split}\lbrace \begin{array}{ccc}{H}_{1111}{\dot{u}}^{'}+{H}_{1112}{\dot{v}}^{'}-{\mathrm{3a}}^{1}{\dot{u}}^{\text{'''}}& =& {\dot{F}}_{1}\\ {H}_{1211}{\dot{u}}^{'}+{H}_{1212}{\dot{v}}^{'}& =& {\dot{F}}_{2}\end{array}\end{split}\]

Le système d’équations différentielles couplées peut se réduire à une équation, la seconde équation exprimant une relation linéaire entre les premiers gradients des déplacements verticaux et horizontaux dans la bande.

(2912)\[\begin{split}\lbrace \begin{array}{ccc}{H}_{1111}{\dot{u}}^{'}+{H}_{1112}\frac{{\dot{F}}_{2}-{H}_{1211}{\dot{u}}^{'}}{{H}_{1212}}-{\mathrm{3a}}^{1}{\dot{u}}^{\text{'''}}& =& {\dot{F}}_{1}\\ \frac{{\dot{F}}_{2}-{H}_{1211}{\dot{u}}^{'}}{{H}_{1212}}& =& {\dot{v}}^{'}\end{array}\end{split}\]

Les solutions de ce système d’équations différentielles linéaires d’ordre 2 s’expriment comme la somme d’une solution particulière répondant aux conditions aux limites imposées et d’une solution partielle établie à partir des racines du polynôme caractéristique de ce système: \({\dot{u}}^{'}={\dot{u}}_{0}^{'}+{\dot{u}}_{p}^{'}\)

On suppose pour \({\dot{u}}_{p}^{'}\) une solution du type \({\dot{u}}_{p}^{'}={e}^{\eta x}\) qui permet d’établir l’équation caractéristique suivante:

(2913)\[{H}_{1111}-{H}_{1112}\frac{{H}_{1211}}{{H}_{1212}}-{\mathrm{3a}}^{1}{\eta}^{2}=0\]

On peut alors donner l’expression de la solution partielle en fonction des paramètres matériau du premier et second gradient:

(2914)\[ \begin{align}\begin{aligned}\textrm{soit}\\{\eta}^{2}=\frac{{H}_{1111}{H}_{1212}-{H}_{1112}{H}_{1211}}{{\mathrm{3a}}^{1}{H}_{1212}}\end{aligned}\end{align} \]

La loi constitutive de premier gradient étant élasto-plastique, les solutions seront différentes selon l’état de la région considérée. Dans une région à comportement élastique, la loi constitutive s’écrit \(\dot{\sigma}=A\dot{\epsilon}\) et dans une région à comportement plastique \(\dot{\sigma}={H}^{\mathrm{ep}}\dot{\epsilon}\) .

Finalement les solutions du problème prennent les formes suivantes en fonction du signe de

(2915)\[{\eta}^{2}=\Delta =\frac{{H}_{1111}{H}_{1212}-{H}_{1112}{H}_{1211}}{{\mathrm{3a}}^{1}{H}_{1212}}\]

• Si \(\Delta >0\) , alors

(2916)\[{\dot{u}}^{'}={\dot{u}}_{0}^{'}+({C}_{11}\exp({\eta}_{1}x)+{C}_{12}\exp({\eta}_{2}x))\]

avec \({\eta}_{i}=\pm \sqrt{\Delta}\) et \({\dot{u}}_{0}^{'}\) obtenue à partir de l’équation de la solution particulière:

(2917)\[ \begin{align}\begin{aligned}{\dot{u}}_{0}^{'}=\frac{{H}_{1212}{\dot{F}}_{1}-{H}_{1112}{\dot{F}}_{2}}{{H}_{1111}{H}_{1212}-{H}_{1112}{H}_{1211}}\\\textrm{et}\\{\dot{v}}_{0}^{'}=\frac{{H}_{1111}{\dot{F}}_{2}-{H}_{1211}{\dot{F}}_{1}}{{H}_{1111}{H}_{1212}-{H}_{1112}{H}_{1211}}\end{aligned}\end{align} \]

De plus, les solutions partielles \({\dot{u}}_{p}^{'}\) vérifient les équations suivantes:

(2918)\[\frac{{H}_{1111}{H}_{1212}-{H}_{1112}{H}_{1211}}{{\mathrm{3a}}^{1}{H}_{1212}}{\dot{u}}_{p}^{'}-{\dot{u}}_{p}^{\text{'''}}=0\]

A partir des conditions de symétrie du domaine étudié et \({\dot{u}}_{p}^{'}(l)=0\) , on en déduit que les coefficients \({C}_{11}\) et \({C}_{12}\) sont nuls.

• Si \(\Delta <0\) , alors

(2919)\[{\dot{u}}^{'}={\dot{u}}_{0}^{'}+({C}_{21}\cos(\omega x)+{C}_{22}\sin(\omega x))\]

avec \(\omega =\sqrt{\frac{{H}_{1211}{H}_{1112}-{H}_{1212}{H}_{1111}}{{\mathrm{3a}}^{1}{H}_{1212}}}\) . De même que pour les solutions précédemment établies pour les zones où \(\Delta >0\) , la solution partielle \({\dot{u}}_{p}^{'}\) vérifie l’équation .

L’expression de la solution dans la zone où \(\Delta <0\) fait apparaître des fonctions trigonométriques. Il est donc clair qu’une zone de localisation peut apparaître dans la structure. On peut alors admettreque la structure privilégie une longueur interne de \({l}_{c}=2\pi /\omega\) , s’exprimant selon le premier mode de plus basse énergie.

On constate aussi que la dépendance de la longueur interne \({l}_{c}\) est en \(\sqrt{{a}^{1}}\) . Il est intéressant de remarquer également que la condition d’apparition des solutions de bifurcations est contrôlée uniquement par les termes du modèle de premier gradient. Autrement dit, la prise en compte d’une cinématique enrichie restreinte aux milieux de second gradient ne modifie pas la valeur du paramètre de contrôle provoquant l’apparition d’une solution bifurquée en bandes de cisaillement.

Résolution numérique

L’option de calcul élémentaire PDIL_ELGA, détermine, pour un état initial donné des contraintes et variables internes, la valeur de A1_LC2 établie selon la formule suivante:

(2920)\[\text{A1\_LC2}=\frac{{a}^{1}}{{{l}_{c}}^{2}}=\frac{\left[{H}_{1211}{H}_{1112}-{H}_{1111}{H}_{1212}\right]}{3{H}_{1212}{(2\pi )}^{2}}\]

pour différentes orientations de la bande de cisaillement. L’écriture explicite des composantes du tenseur de rigidité du modèle local classique est obtenue via les routines de calcul des matrices tangentes en vitesse.

Une rotation de la bande de cisaillement d’un angle \(\theta\) implique une rotation locale appliquée aux composantes du tenseur des contraintes. Elle permet d’estimer les nouvelles composantes du tenseur constitutif local du modèle de comportement classique, portant sur le premier gradient des déplacements.

La discrétisation angulaire utilisée est tout d’abord de \(5°\) . Autour du premier maximum relevé, on effectue une discrétisation à \(1°\) puis \(0.2°\) . De cette façon, on s’assure d’obtenir une valeur précise du paramètre A1_LC2.

L’utilisateur peut ensuite définir la valeur de A1 adaptée à la discrétisation spatiale de la structure étudiée, sachant qu’un minimum de 6 éléments sur LC apparaît nécessaire pour garantir l’indépendance aux maillages des résultats.

La validation numérique de l’option de calcul PDIL_ELGA est portée dans les cas-tests SSNV208A, SSNP125A et WTNV132C.

Vérification

La liste des tests de validation pour le second gradient :

Cas-test	description
ssll117 (modélisations a et b)	Test de D_PLAN_DIL (voir V3.01.117)
sslv117(modélisations a à d)	Test de 3D_DIL (voir V3.01.117)
wtvn132 (modélisations b et f)	Test de D_PLAN_HM_SI_DIL etD_PLAN_HMS_DIL (voir V7.31.32)
wtnl101c	Test de D_PLAN_THMS_DIL (voir V7.30.101)
wtnp119f	Test de D_PLAN_HH2MS_DIL(voir V7.32.119)
wtnl100i	Test de 3D_HMS_DIL(voir V7.30.100)
wtnv142b	Test de 3D_HM_SI_DIL(voir V7.31.142)
wtnv113k	Test de 3D_THMS_DIL(voir V7.31.113)
wtnv136e	Test de 3D_HH2MS_DIL(voir V7.31.136)

Bibliographie

1. Chambon R, Caillerie D, El Hassan N : “One-dimensional localisation studied with a second grade model”, Eur. J. Mech. A/Solids (1998), vol 17, pp 637-656.

2. Chambon R, Caillerie D, Matsushima T: “Plastic continuum with microstructure, local second gradient theories for geomaterials: localization studies”, Int. J. Solids and Struct. (2001), vol 38, pp 8503-8527.

3. Matsushima T, Chambon R, Caillerie D: “Large strain finite element analysis of local second gradient models, application to localization”, Int. J. Num. Meth. Eng. (2002), vol 54, pp 499-521.

4. Fernandes R: “Modélisation numérique objective des problèmes hydromécaniques couplés dans le cas des géomatériaux”, Thèse de doctorat de l’Université Joseph Fourier – Grenoble, 23 janvier 2009.

5. Fernandes R, Chavant C, Chambon R: “A simplified second gradient model for dilatant materials : theory and numerical implementation”, Int. J. Solids and Struct. (2008), vol 45, pp 5289-5307.

6. Mindlin RD: “Microstructure in linear elasticity”, Arch. Ration. Mech. Anal. (1964), vol 16, pp 51-78.

7. Mindlin RD: “Second gradient of strain and surface-tension in linear elasticity”, Int. J. Solids and Struct. (1965), vol 1, pp 417-738.

8. Germain P: “La méthode des puissances virtuelles en mécanique des milieux continues. Première partie: théorie du second gradient”, Journal de Mécanique (1973), vol 12, pp 235-274.

9. Germain P: “The method of virtual power in continuum mechanics. Part 2 : Microstructure”, SIAM J. Appl. Math. (1973), vol 25, pp 556-575.

10. Zervos A, Papanastasiou P, Vardoulakis I: “A finite element displacement formulation for gradient elastoplasticity”, Int. J. Num. Meth. Eng. (2001), vol 50, pp 1369-1388.

11. Foucault A.: “Modélisation du comportement cyclique des ouvrages en terre intégrant des techniques de régularisation ”, Thèse de doctorat de l’Ecole Centrale Paris – Laboratoire MSS-Mat, 21 juin 2010.

12. Desrues J., Viggiani G. “Strain localization in sand: an overview of the experimental results obtained in Grenoble using stereophotogrammetry”. Int. J. Num. Anal. Meth. Geom. (2004), vol. 28: pp. 279-321.

13. Vardoulakis I., Sulem J. “ Bifurcation Analysis in Geomechanics” . London: Blackie, 1995.

14. Gantier M.: «Modélisation numérique robuste et fiable de la fissuration des roches et des interfaces», Thèse de doctorat de l’Université Grenoble Alpes – Grenoble, 17 novembre 2021.