r7.01.25 Lois de comportement des joints des barrages : JOINT_MECA_RUPT, JOINT_MECA_FROT etJOINT_MECA_ENDO

Résumé:

Ce document décrit des lois surfaciques permettant de modéliser la rupture et le frottement entre les lèvres d’une fissure ou d’un joint. La loi JOINT_MECA_RUPT ne représente que le phénomène d’endommagement, elle est basée sur une formulation cohésive de la rupture. La loi JOINT_MECA_FROT est une version élastoplastique d’une loi de frottement de Mohr-Coulomb. La loi JOINT_MECA_ENDO est une loi plus complète, elle combine des deux mécanismes: l’endommagement et la plasticité sont couplés à travers de l’énergie libre. Cette dernière loi est écrite dans le formalisme standard généralisé [bib18].

En mécanique ces loiss’appuient sur les modélisations des joints standards XXX_JOINT. Les lois sont dédiées à la modélisation des barrages, plus précisément des joints béton/rocher ou des joints entre les plots d’un barrage. En fonction du type de chargement l’utilisation de l’une ou de l’autre loi peut être choisie pour différentes parties de l’ouvrage.

Pour pouvoir simuler le comportement des barrages réels certaines spécificités de la construction ont été introduites, notamment le clavage et le sciage. Le clavage est une procédure industrielled’injection du béton sous pression entre les plots d’un barrage. C’est une étape intermédiaire de la construction d’un barrage-voûte, il sert à renforcer son étanchéité après la phase de construction de plots verticaux. Le sciage est une procédure similaire pendant laquelle le barrage est scié afin de relâcher les contraintes de compression. Chaque procédure est définie via les mots-clef PRES_CLAVAGE et SCIAGE dans DEFI_MATERIAU. Le clavage et le sciage ne sont pas activés pour l’instant pour la loi JOINT_MECA_ENDO.

L’effet de pression hydrostatique dû à la présence de fluide dans le joint est pris en compte via la contrainte de Terzaghi. Toutes les lois admettent soit la modélisation sans couplage (pression imposée dans l’option PRES_FLUIDE), soit une modélisation couplée hydromécanique (XXX_JOINT_HYME) ce qui active la propagation des sous-pressions à l’interface barrage-rocher.

Formulation théorique de JOINT_MECA_RUPT

La loi JOINT_MECA_RUPT accepte une modélisation couplée hydromécanique, mais ces deux phénomènes peuvent être traités séparément. Dans un premier temps on décrira la partie mécanique de la loi, qui englobe rupture, contact, procédure de clavage et pression imposée. Elle s’appuie sur les modélisations XXX_JOINT (R3.06.09).

Remarque

Les deux prochaines sections présentent le concept théorique général. En première lecture on peut les omettre et passer directement à la r7.01.25-vecteur-contrainte , qui fournit suffisamment d’éléments pour comprendre l’implémentation numérique de la loi.

Loi cohésive en mécanique

Le cadre théorique choisi pour modéliser la rupture des joints de barrages en béton est basé sur les modèles cohésifs (r7.02.11, [bib1], [bib6]), car en rupture fragile, pour éviter le problème de contraintes infinies en fond de fissure, on peut introduire des forces de cohésion qui imposent un critère d’amorçage en contrainte. Les forces à caractère évanescent s’exercent alors entre les particules de part et d’autre du plan de séparation de la fissure pendant son ouverture (voir Fig. 229).

Fig. 229 Schéma d’une fissure cohésive

Du point de vue physique on considère que l’ouverture de la fissure coûte une énergie proportionnelle à sa longueur en 2D et à sa surface en 3D. On l’appelle énergie de surface que l’on exprime à l’aide de la densité d’énergie \(\Psi ={\int}_{\Gamma}\psi (\vec{\delta})d\Gamma\) [1]_. Le champ de déplacement à l’équilibre \(\mathrm{u}\) est obtenu en minimisant la somme de l’énergie élastique \(\Phi\) , de l’énergie de surface, et du travail des efforts extérieurs \({W}^{\mathrm{ext.}}\) . La solution est obtenus en utilisant une approche variationnelle de la rupture. L’état qui réalise le minimum de l’énergie totale correspond à un état d’équilibre mécanique:

(3382)\[\underset{\mathrm{u}}{\min}\left(\Phi (u)+\Psi (\vec{\delta}(u))+{W}^{\mathit{ext}}\right)\]

La surface de discontinuité est discrétisée en 2D ou 3D par des éléments finis de joint (voir documentation r3.06.09). Le saut de déplacement dans l’élément \(\vec{\delta}=({\delta}_{n},{\delta}_{\mathit{t1}},{\delta}_{\mathit{t2}})\) est une fonction linéaire des déplacements nodaux. La force [2]_ de cohésion qui s’exerce sur les lèvres de la fissure est notée \(\vec{\sigma}\) , elle est définie par la dérivée de la densité d’énergie de surface par rapport au saut de déplacement. On appelle loi cohésive une relation entre \(\vec{\sigma}\) et le saut de déplacement \(\vec{\delta}\) . Les paramètres matériaux les plus pertinents, qui décrivent le joint d’un barrage sont:

• les deux rigidités en sollicitations normale \({K}_{n}\) et tangentielle \({K}_{t}\) , qui caractérisent la surface et les matériaux de remplissage de la fissure;

• et la contrainte critique à la rupture \({\sigma}_{max}\)

On introduit, par ailleurs, trois paramètres numériques adimensionnels: \({P}_{\mathit{rupt}}\) , \({P}_{\mathit{cont}}\) et \(\alpha ` . Le premier pilote la régularisation de la pente d'adoucissement en rupture, le second la pénalisation du contact et le troisième assure une reprise progressive d’efforts tangentiels en fonction de l'ouverture normale. Ce dernier peut être associé à la taille relative des aspérités des surfaces en contact: :math:\)alpha in [0,2]` .

Dans les lois cohésives standards (R7.02.11) l’énergie de surface dépend du vecteur de déplacement [3] et les contraintes sont définies en tant que les dérivées premières de l’énergie:

(3383)\[ \begin{align}\begin{aligned}{\sigma}_{n}=\frac{\partial {\psi}_{n}({\delta}_{n})}{\partial {\delta}_{n}}\\\textrm{et}\\\vec{{\sigma}_{t}}=\frac{\partial {\psi}_{t}(\vec{{\delta}_{t}})}{\partial \vec{{\delta}_{t}}}\end{aligned}\end{align} \]

Avec l’énergie:

(3384)\[\Psi (\vec{\delta})\equiv {\int}_{\Gamma}\psi ({\delta}_{n},\vec{{\delta}_{t}})d\Gamma\]

À la différence de ces lois cohésives standards, dans ce modèle, uniquement la partie normale de la loi est dérivée à partir de l’énergie de surface, alors que la composante tangentielle de la loi est donnée d’une manière explicite [4].

Dans les deux cas, l’irréversibilité de la fissuration est prise en compte via une condition d’accroissement d’ouverture normale maximale de joint.

Énergie de surface pour le comportement normal

La densité d’énergie de surface \(\psi ` , en un point donné de la fissure, dépend explicitement du saut de déplacement normal entre les lèvres de la fissure :math:`{\delta}_{n}\) . Elle varie aussi en fonction de l’état du joint, ce qui est pris en compte via une variable interne seuil \(\kappa ⩾0\) , qui gère l’irréversibilité de la fissuration. Cette dernière mémorise la plus grande norme du saut atteinte au cours de l’ouverture. Sa loi d’évolution entre deux incréments de chargement successifs \(-\) et \(+\) s’écrit:

(3385)\[{\kappa}^{+}=\max({\kappa}^{-},{\delta}_{n}^{+})\]

Suivant la valeur d’ouverture du joint, on pourra se retrouver dans une de ces trois situations:

1. Le joint comprimé se trouve dans le régime de contact;

2. Pour une ouverture normale positive, si cette dernière dépasse le seuil on parle de régime dissipatif(dissipation d’énergie au cours de la fissuration);

3. Dans le cas intermédiaire le joint est dans un régime linéaire (décharge ou recharge linéaire sans dissipation d’énergie).

La partie normale de l’énergie s’écrit:

(3386)\[\begin{split}{\psi}_{n}({\delta}_{n},\kappa )=\left\{ \begin{array}{ccc}{\psi}_{n}^{\mathit{con}}({\delta}_{n},\kappa )\hfill & \text{si}& {\delta}_{n}<0\\ {\psi}_{n}^{\mathit{lin}}({\delta}_{n},\kappa )\hfill & \text{si}& 0⩽{\delta}_{n}<\kappa \\ {\psi}_{n}^{\mathit{dis}}({\delta}_{n})\hfill & \text{si}& {\delta}_{n}⩾\kappa \end{array}\right.\end{split}\]

Et la partie tangentielle:

(3387)\[\begin{split}{\psi}_{t}(\vec{{\delta}_{t}})=\left\{ \begin{array}{ccc}{\psi}_{t}^{\mathit{fer}}(\vec{{\delta}_{t}})& \text{si}& {\delta}_{n}<0\\ {\psi}_{t}^{\mathit{ouv}}(\vec{{\delta}_{t}})& \text{si}& {\delta}_{n}⩾0\end{array}\right.\end{split}\]

D’une manière synthétique l’énergie de surface s’écrit de la manière suivante:

(3388)\[\psi ({\delta}_{n},\kappa )=H({\delta}_{n}-\kappa ){\psi}_{n}^{\mathit{dis}}({\delta}_{n})+H({\delta}_{n})\cdot H(\kappa -{\delta}_{n}){\psi}_{n}^{\mathit{lin}}({\delta}_{n},\kappa )+H(-{\delta}_{n})\cdot {\psi}_{n}^{\mathit{con}}({\delta}_{n})\]

Où \(H(x)\) est la fonction de Heaviside (\(H(x)=0\text{si}x<0,H(x)=1\text{si}x⩾0\) ).

Une des particularités de la loi JOINT_MECA_RUPT est que, quelque soit le régime, le comportement est toujours linéaire au niveau des contraintes. Au niveau des énergies, nous obtenons des fonctions quadratiques, qui sont données dans la suite à une constante additive près, qui elle dépend du seuil.

Le comportement normal de joint est séparé en trois régimes: contact, linéaire et dissipatif. La fonction \({\psi}_{n}^{\mathit{con}}({\delta}_{n})=\frac{{P}_{\mathit{con}}{K}_{n}{\delta}_{n}^{2}}{2}\) assure la condition de contact (non interpénétration des lèvres de la fissure) elle est régularisée afin d’obtenir une meilleure convergence numérique. En faisant varier le paramètre \({P}_{\mathrm{con}}\) , on peut changer la rigidité normale pour le joint fermé.

En régime linéaire, dans le cas où une fissure existante évolue sans dissiper d’énergie, la densité d’énergie correspondante ne dépend que du paramètre de rigidité normale: \({\psi}_{n}^{\mathit{lin}}({\delta}_{n},\kappa )=\frac{{K}_{a}(\kappa ){\delta}_{n}^{2}}{2}\) , où \({K}_{a}(\kappa )=({P}_{\mathit{rupt}}^{-1}+1){\sigma}_{max}/\kappa -{K}_{n}{P}_{\mathit{rupt}}^{-1}\) est une fonction décroissante du seuil de rupture \(\kappa\) .

Dans le régime dissipatif, afin d’obtenir un adoucissement linéaire, la densité de l’énergie correspondante a une forme quadratique en fonction de l’ouverture:

(3389)\[\begin{split}{\psi}_{n}^{\mathit{dis}}({\delta}_{n})=\left\{ \begin{array}{ccc}{\sigma}_{max}(1+{P}_{\mathit{rupt}}^{-1}){\delta}_{n}-{K}_{n}{P}_{\mathit{rupt}}^{-1}{\delta}_{n}^{2}/2& \text{si}& {\delta}_{n}<{\sigma}_{max}(1+{P}_{\mathit{rupt}})/{K}_{n}\\ {\sigma}_{max}^{2}{(1+{P}_{\mathit{rupt}})}^{2}/({\mathrm{2P}}_{\mathit{rupt}}{K}_{n})& \text{si}& {\delta}_{n}⩾{\sigma}_{max}(1+{P}_{\mathit{rupt}})/{K}_{n}\end{array}\right.\end{split}\]

Fig. 230 Densité d’énergie de surface en fonction du saut de déplacement pour différentes valeurs du seuil d’endommagement \(\kappa\)

La constante additive fait translater les énergies élastiques \({\psi}_{n}^{\mathit{lin}}\) et \({\psi}_{n}^{\mathit{con}}\) définies précédemment de sorte que peu importe l’état d’endommagement du joint on obtient toujours le même taux de restitution de l’énergie (constante de Griffith \({G}_{f}\) , voir figure Fig. 230 et (61) . Elle ne dépend que du seuil et n’affecte pas l’expression des contraintes.

\({P}_{\mathrm{rupt}}\) est introduit de telle sorte que plus il augmente plus l’ouverture critique à la rupture est importante et par conséquence plus cela fait accroître l’énergie dissipative de Griffith \({G}_{f}\) .

En résumé : le comportement normal de la loi JOINT_MECA_RUPT est piloté par l’évolution de la densité d’énergie surfacique, celle ci se présente sous la forme d’un puits potentiel. Elle comporte trois régimes principaux: contact, traction linéaire élastique, endommagement/adoucissement, dont les profils correspondant sont approximés par des fonctions quadratiques en ouverture. Le joint commence à s’endommager quand la contrainte normale atteint la valeur critique \({\sigma}_{n}={\sigma}_{max}\) *. Plus le joint est endommagé plus le puits énergétique s’aplatit voir (Figure* Fig. 232 ). Le paramètre d’endommagement (le seuil) lui évolue seulement dans ce dernier régime en partant de sa valeur initiale pour le joint sain \({\kappa}_{0}={\sigma}_{max}/{K}_{n}\) jusqu’à la valeur ultime pour le joint complètement endommagé \({\kappa}_{\mathit{rupt}}={\sigma}_{max}(1+{P}_{\mathit{rupt}})/{K}_{n}\) *.* \({P}_{\mathit{con}}\) est une constante définie par l’utilisateur qui change le niveau de pénalisation en contact (voir Figure Fig. 232 ). \({P}_{\mathrm{rupt}}\) change l’énergie dissipée par unité de surface unitaire \({G}_{f}={\sigma}_{\max}^{2}(1+{P}_{\mathrm{rupt}})/({\mathrm{2K}}_{n})\) . Au final, on obtient l’expression suivante pour l’énergie normale:

(3390)\[\begin{split}{\psi}_{n}(\textcolor[rgb]{0,0,1}{{\delta}_{n}},\kappa )=\left\{ \begin{array}{ccc}{\sigma}_{max}(1+{P}_{\mathit{rupt}}^{-1})(\kappa -{\kappa}_{0})/2+{P}_{\mathit{con}}{K}_{n}\textcolor[rgb]{0,0,1}{{\delta}_{n}^{2}}/2& \text{si}& {\delta}_{n}<0\\ {\sigma}_{max}(1+{P}_{\mathit{rupt}}^{-1})(\kappa -{\kappa}_{0})/2+\left[{\sigma}_{max}(1+{P}_{\mathit{rupt}}^{-1})/\kappa -{K}_{n}{P}_{\mathit{rupt}}^{-1}\right]\textcolor[rgb]{0,0,1}{{\delta}_{n}^{2}}/2& \text{si}& 0⩽{\delta}_{n}<\kappa \\ {\sigma}_{max}(1+{P}_{\mathit{rupt}}^{-1})(\textcolor[rgb]{0,0,1}{{\delta}_{n}}-{\kappa}_{0}/2)-{K}_{n}{P}_{\mathit{rupt}}^{-1}\textcolor[rgb]{0,0,1}{{\delta}_{n}^{2}}/2& \text{si}& \kappa ⩽{\delta}_{n}<{\kappa}_{\mathit{rupt}}\\ {G}_{f}& \text{si}& {\delta}_{n}⩾{\kappa}_{\mathit{rupt}}\end{array}\right.\end{split}\]

Cette fonction est continue et dérivable, ce qui assure la continuité des contraintes.

Vecteur contrainte

Le vecteur contrainte dans l’élément est noté \(\vec{\sigma}=({\sigma}_{n},\vec{{\sigma}_{t}})\) [5], il peut être séparé en plusieurs régimes. Pour la partie normale le vecteur-contrainte est égal à la somme des dérivées de la densité d’énergie de surface et de la densité d’énergie de pénalisation en contact par rapport au saut.

(3391)\[{\sigma}_{n}=H({\delta}_{n}-\kappa ){\sigma}_{n}^{\mathit{dis}}+H(\kappa -{\delta}_{n})H({\delta}_{n}){\sigma}_{n}^{\mathit{lin}}+H(-{\delta}_{n}){\sigma}_{n}^{\mathit{con}}\]

Il suffit doncde dériver les expressions données dans la section précédente r7.01.25-energie-surface pour obtenir la composante normale des contraintes [6]. La partie tangentielle \({\overrightarrow{\sigma}}_{t}^{\mathrm{fer}}=f({K}_{t},\overrightarrow{{\delta}_{t}})\) s’écrit:

(3392)\[{\vec{\sigma}}_{t}=H({\delta}_{n}){\vec{\sigma}}_{t}^{\mathit{ouv}}+H(-{\delta}_{n}){\vec{\sigma}}_{t}^{\mathit{fer}}\]

C’est une fonction de l’ouverture tangentielle, elle met en jeu la rigidité tangentielle pour le joint fermé. En fonction du profil de l’interface le comportement tangentiel pour le joint ouvert peut être varié \({\vec{{\sigma}_{t}}}^{\mathit{ouv}}={\vec{{\sigma}_{t}}}^{\mathit{fer}}\) pour les surfaces en créneau (comme par exemple sur la r7.01.25-fig-joint-lisse), ou bien \({\overrightarrow{{\sigma}_{t}}}^{\mathrm{ouv}}\equiv 0\) pour les surfaces très lisses (r7.01.25-fig-joint-rougeux).

Contraintes normales

Considérons l’expression la contrainte normale:

(3393)\[\begin{split}{\sigma}_{n}({\delta}_{n},\kappa )=\left\{ \begin{array}{ccc}{P}_{\mathit{con}}{K}_{n}{\delta}_{n}& \text{si}& {\delta}_{n}<0\\ \left[{\sigma}_{max}(1+{P}_{\mathit{rupt}}^{-1})/\kappa -{K}_{n}{P}_{\mathit{rupt}}^{-1}\right]{\delta}_{n}& \text{si}& 0⩽{\delta}_{n}<\kappa \\ {\sigma}_{max}(1+{P}_{\mathit{rupt}}^{-1})-{K}_{n}{P}_{\mathit{rupt}}^{-1}{\delta}_{n}& \text{si}& \kappa ⩽{\delta}_{n}<{\kappa}_{\mathit{rupt}}\\ 0& \text{si}& {\delta}_{n}⩾{\kappa}_{\mathit{rupt}}\end{array}\right.\end{split}\]

L’évolution de cette contrainte dans la zone de traction en fonction du saut est représenté sur la Fig. 231. Les flèches représentent le sens d’évolution possible de la contrainte suivant que le processus d’ouverture est réversible (régime linéaire) ou non (régime dissipatif). A l’amorçage, le joint se comporte d’abord de façon élastique linéaire, puis dès que la contrainte normale atteint la valeur critique \({\sigma}_{n}={\sigma}_{\max}\) , il a un comportement adoucissant: il perd progressivement sa rigidité, ce qui donne le régime linéaire mais pas élastique, il est caractérisé par une pente d’adoucissement de la rupture \(-{K}_{n}/{P}_{\mathrm{rupt}}\) . La rigidité élastique pour le joint sain définit la valeur initiale du seuil d’endommagement \({\kappa}_{0}={\sigma}_{\max}/{K}_{n}\) . Le seuil de la rupture est donné par \({\kappa}_{\mathrm{rupt}}={\sigma}_{\max}(1+{P}_{\mathrm{rupt}})/{K}_{n}\) . Plus \({P}_{\mathit{rupt}}\) est important plus l’énergie de dissipation accroît \({G}_{f}={\sigma}_{\max}^{2}(1+{P}_{\mathrm{rupt}})/({\mathrm{2K}}_{n})\) .

Fig. 231 Dépendance des contraintes en fonction de l’ouverture

Contrainte de pénalisation du contact

La valeur de la pente de pénalisation en contact est donnée par la relation suivante:

(3394)\[ \begin{align}\begin{aligned}{\sigma}_{n}({\delta}_{n})={P}_{\mathit{con}}{K}_{n}{\delta}_{n}\\\textrm{si}\\{\delta}_{n}<0\end{aligned}\end{align} \]

Fig. 232 Contrainte cohésive normale en fonction du saut normal pour le joint partiellement endommagé

Le paramètre numérique PENA_CONTACT, donné par l’utilisateur, permet de jouer sur la pente de la pénalisation du contact (voir figure Fig. 232). Il vaut par défaut \(1\) , ce qui correspond au cas où la pente du contact est identique à celle de la rigidité en ouverture. Si on choisit une valeur supérieure à \(1\) , on augmente la pénalisation. Ceci permet de modéliser, par exemple, la reprise des efforts par le béton partiellement endommagé en traction. Pour une valeur inférieure à \(1\) , on diminue la pénalisation, ce qui permet de simuler le béton partiellement endommagé en compression.

Contrainte tangentielle

Pour les joints de barrage à faible ouverture, on observe que indépendamment du régime de chargement normal, la contrainte tangentielle varie toujours linéairement, la rigidité tangentielle est fonction de l’ouverture normale. Dans le cas extrême d’une surface de contact parfaitement lisse, la rigidité tangentielle chute brutalement à zéro à l’ouverture normale positive. En conséquence, l’énergie surfacique de la loi n’est plus continue, ce qui génère en principe un pic dans la contrainte normale (fonction delta \(\textcolor[rgb]{1,0,0}{\delta}(x)\) ) à l’ouverture. C’est pour cette raison que nous renonçons dans la version actuelle de la loi de garder le formalisme énergétique complet. La loi tangentielle est alors postulée d’une manière empirique sous forme linéaire avec la rigidité tangentielle dépendante de l’ouverture normale:

(3395)\[\begin{split}\vec{{\sigma}_{t}}({\delta}_{t},{\delta}_{n})=\left\{ \begin{array}{ccc}{K}_{t}\cdot (\vec{{\delta}_{t}}-{\vec{\delta}}_{\mathit{shift}})& \text{si}& {\delta}_{n}<0\\ (1-{\delta}_{n}/{\kappa}_{\mathit{rupt}}^{\tan}){K}_{t}\cdot (\vec{{\delta}_{t}}-{\vec{\delta}}_{\mathit{shift}})& \text{si}& 0⩽{\delta}_{n}<{\kappa}_{\mathit{rupt}}^{\tan}\\ 0& \text{si}& {\delta}_{n}⩾{\kappa}_{\mathit{rupt}}^{\tan}\end{array}\right.\end{split}\]

La variable d’histoire \(\vec{{\delta}_{\mathit{shift}}}\) représente le décalage du point d’équilibre en glissement tangentiel à joint ouvert: \(\vec{{\delta}_{\mathit{shift}}}\) est égale à la dernière valeur de \(\vec{{\delta}_{t}}\) pour laquelle le joint a été ouvert \({\delta}_{n}\ge {\kappa}_{\mathit{rupt}}^{\tan}\) . Nous introduisons le seuil de rupture tangentielle \({\kappa}_{\mathit{rupt}}^{\tan}={\kappa}_{\mathit{rupt}}\tan(\alpha \pi /4)\) , dont la valeur peut être modifiée par l’utilisateur avec \(\alpha \in [0,2]\) et le mot-clé ALPHA. Pour une valeur zéro la pente tangentielle change brutalement à l’ouverture, pour la valeur \(\alpha =2\) la rigidité tangentielle n’évolue pas. Par souci de compatibilité avec les lois de comportement développées dans le code Gefdyn, nous ne faisons pas de correction de la composante normale des contraintesdans la phase de transition, celle ci est toujours donnée par (4846), ce qui donne une matrice tangente non symétrique et par conséquence le régime à dissipation non contrôlée. L’évolution de la contrainte tangentielle se sépare en trois régimes: joint élastique en compression; joint élastique partiellement ouvert avec une rigidité diminuée; joint rompu complètement (voir équation (4848) et figure Fig. 233.)

Fig. 233 Illustration du couplage entre le cisaillement et l’ouverture normale de joint : I dans le régime en compression; II dans le régime d’ouverture partiel; III dans le régime d’ouverture complète. On note \({K}_{t}^{a}\equiv (1-{\delta}_{n}/{\kappa}_{\mathrm{rupt}}^{\tan}){K}_{t}\)

Opérateur tangent

Comme le comportement dans chacun des régimes est linéaire, le calcul de la matrice tangente est aisé:

(3396)\[\begin{split}\frac{\partial {\sigma}_{n}({\delta}_{n})}{\partial {\delta}_{n}}=\left\{ \begin{array}{ccc}{P}_{\mathit{con}}{K}_{n}& \text{si}& {\delta}_{n}<0\\ {\sigma}_{max}(1+{P}_{\mathit{rupt}}^{-1})/\kappa -{K}_{n}{P}_{\mathit{rupt}}^{-1}& \text{si}& 0⩽{\delta}_{n}<\kappa \\ -{K}_{n}{P}_{\mathit{rupt}}^{-1}& \text{si}& \kappa ⩽{\delta}_{n}<{\kappa}_{\mathit{rupt}}\\ 0& \text{si}& {\delta}_{n}⩾{\kappa}_{\mathit{rupt}}\end{array}\right.\end{split}\]

(3397)\[ \begin{align}\begin{aligned}\frac{\partial {\sigma}_{n}({\delta}_{n})}{\partial {\delta}_{t}}=0\\\begin{split}\frac{\partial \vec{{\sigma}_{t}}(\vec{{\delta}_{t}},{\delta}_{n})}{\partial \vec{{\delta}_{t}}}=\left\{ \begin{array}{ccc}{K}_{t}\mathit{Id}& \text{si}& {\delta}_{n}<0\\ (1-{\delta}_{n}/{k}_{\mathit{rupt}}^{\tan}){K}_{t}\mathit{Id}& \text{si}& 0⩽{\delta}_{n}<{\kappa}_{\mathit{rupt}}^{\tan}\\ 0& \text{si}& {\delta}_{n}⩾{\kappa}_{\mathit{rupt}}^{\tan}\end{array}\right.\end{split}\end{aligned}\end{align} \]

(3398)\[\begin{split}\frac{\partial \vec{{\sigma}_{t}}(\vec{{\delta}_{t}},{\delta}_{n})}{\partial \vec{{\delta}_{t}}}=\left\{ \begin{array}{ccc}{K}_{t}\mathit{Id}& \text{si}& {\delta}_{n}<0\\ (1-{\delta}_{n}/{k}_{\mathit{rupt}}^{\tan}){K}_{t}\mathit{Id}& \text{si}& 0⩽{\delta}_{n}<{\kappa}_{\mathit{rupt}}^{\tan}\\ 0& \text{si}& {\delta}_{n}⩾{\kappa}_{\mathit{rupt}}^{\tan}\end{array}\right.\end{split}\]

(3399)\[\begin{split}\frac{\partial \vec{{\sigma}_{t}}(\vec{{\delta}_{t}},{\delta}_{n})}{\partial {\delta}_{n}}=\left\{ \begin{array}{ccc}0& \text{si}& {\delta}_{n}<0\\ -{K}_{t}(\vec{{\delta}_{t}}-{\vec{\delta}}_{\mathit{shift}})/{k}_{\mathit{rupt}}^{\tan}& \text{si}& 0⩽{\delta}_{n}<{\kappa}_{\mathit{rupt}}^{\tan}\\ 0& \text{si}& {\delta}_{n}⩾{\kappa}_{\mathit{rupt}}^{\tan}\end{array}\right.\end{split}\]

Notons que la matrice tangente n’est pas symétrique. Ceci résulte de la non répercussion sur les contraintes normales de la régularisation de l’évolution de la contrainte tangentielle à l’ouverture du joint (les termes singuliers ne sont pas pris en compte).

Réalisation numérique du clavage

Le clavage correspond physiquement à une procédure d’injection de coulis de béton sous pression entre les plots du barrage, il est caractérisé par un simple paramètre la pression locale de coulis injecté \({\sigma}_{\mathit{nc}}⩾0\) . Pour mettre en place le clavage l’utilisateur doit définir une fonction de pression de clavage, mot clef PRES_CLAVAGE, qui dépend à la foi de l’espace (clavage aux différents endroits à différente pression) et du temps (plusieurs clavages successifs). Les endroits où la pression de clavage est négative ne sont pas clavés.

La procédure est modélisée par la modification de l’épaisseur des joints concernés. Si le joint est en forte compression initialement, le clavage ne l’influence pas. Si par contre le joint est ouvert ou pas suffisamment comprimé (c’est-à-dire si \({\sigma}_{n}>-{\sigma}_{\mathrm{nc}}\) ), le clavage se traduira par le changement du paramètre d’épaisseur totale du joint noté \({\delta}_{\mathit{nf}}^{+}={\delta}_{\mathit{nf}}^{-}+{\delta}_{n}^{+}+{\sigma}_{\mathit{nc}}/({P}_{\mathit{con}}{K}_{n})\) .

Le clavage préserve alors l’ouverture normale de chaque joint tout en translatant la loi de comportement selon l’axe de l’abscisse. La nouvelle position d’équilibre local possède une contrainte normale égale à la pression du béton injecté \({\sigma}_{n}=-{\sigma}_{\mathrm{nc}}\) . Cette procédure est appliquée localement. Bien que chaque joint soumis au clavage se retrouve dans le régime linéaire en compression, la sélection elle même des joints à claver est «non-linéaire». La procédure globale devient, de ce fait, aussi non-linéaire et il ne suffit pas de modifier juste les contraintes normales des joints clavés pour obtenir à nouveau l’équilibre mécanique après le clavage. En conséquence un calcul d’équilibre mécanique est réalisé après ce clavage «local» pour remettre le système dans son état d’équilibre. Si le jeu du joint à claver est modifié, la mise à jour de l’épaisseur des joints est effectuée et ainsi de suite tant qu’il existe les points où \({\sigma}_{n}>-{\sigma}_{\mathrm{nc}}\) . Les joints concernés se ferment et se mettent progressivement en compression tout en suivant la courbe de comportement normal. Cette procédure n’est faite qu’une seule fois pendant la phase de «construction» numérique de barrage. Une illustration est donnée sur la Fig. 236. Les points «in» et «out» correspondent respectivement aux valeurs de la contrainte avant et après le clavage.

Le comportement normal du joint après le clavage peut être modifié. Ainsi, la résistance à la traction peut être restaurée soit partiellement, soit complètement (c’est le cas sur la Fig. 236) selon l’endommagement des joints avant la procédure de clavage. Dans la procédure telle qu’elle a été développée nous avons fait le choix de ne pas restaurer la résistance à la traction, celle-ci garde sa valeur courante.

Fig. 234 L’évolution de la contrainte normale pendant le clavage : joint partiellement endommagé

Fig. 235 L’évolution de la contrainte normale pendant le clavage : joint complètement endommagé

Afin d’identifier la matrice tangente pour le clavage il est utile de représenter cette procédure comme une modification temporaire de la loi de comportement initiale par une loi incrémentale. Pour le régime tangentiel la procédure de clavage n’importe pas de modification: la contrainte associée dépend toujours d’une manière linéaire du déplacement tangentiel et la pente de chargement varie en fonction de l’ouverture de joint. Le couplage normal mérite d’être analysé plus en détail. Supposons qu’on se trouve à un état de chargement donné \({\delta}_{n+}\) et on essaie de retrouver l’état de contrainte infinitésimalement proche. Regardons d’abord le cas où \({\delta}_{n}\) diminue. Vu que l’épaisseur de joint était déjà mise à jour à l’itération de Newton précédente le joint rentre dans le domaine de compression et l’activation de clavage ne joue aucun rôle à cette étape. On retrouve donc la pente de compression de la loi initiale. Par contre si l’on charge le joint en traction la procédure de clavage s’active et on procède à l’adaptation de l’épaisseur de joint en sorte de faire plafonner la contrainte normale par la pression de clavage. De facto le comportement devient bilinéaire avec une pente de compression de la loi initiale en compression et la pente dans la zone de traction, qui s’annule complètement (voir Fig. 236).

Fig. 236 Modification de la loi de comportement pendant la phase de clavage

Comme le comportement dans chacun des régimes reste linéaire, le calcul de la matrice tangente est toujours aisé:

(3400)\[\begin{split}\frac{\partial {\sigma}_{n}({\delta}_{n})}{\partial {\delta}_{n}}=\left\{ \begin{array}{ccc}{P}_{\mathit{con}}{K}_{n}& \text{si}& {\delta}_{n}<{\delta}_{\mathit{nf}}\\ 0& \text{si}& {\delta}_{n}⩾{\delta}_{\mathit{nf}}\end{array}\right.\end{split}\]

(3401)\[\frac{\partial {\sigma}_{n}({\delta}_{n})}{\partial {\delta}_{t}}=0\]

(3402)\[\begin{split}\frac{\partial \vec{{\sigma}_{t}}(\vec{{\delta}_{t}},{\delta}_{n})}{\partial \vec{{\delta}_{t}}}=\left\{ \begin{array}{ccc}{K}_{t}{I}_{d}& \text{si}& {\delta}_{n}<{\delta}_{\mathit{nf}}\\ (1-{\delta}_{n}/{k}_{\mathit{rupt}}^{\tan}){K}_{t}{I}_{d}& \text{si}& {\delta}_{\mathit{nf}}⩽{\delta}_{n}<{\kappa}_{\mathit{rupt}}^{\tan}\\ 0& \text{si}& {\delta}_{n}⩾{\kappa}_{\mathit{rupt}}^{\tan}\end{array}\right.\end{split}\]

(3403)\[\begin{split}\frac{\partial \vec{{\sigma}_{t}}(\vec{{\delta}_{t}},{\delta}_{n})}{\partial {\delta}_{n}}=\left\{ \begin{array}{ccc}0& \text{si}& {\delta}_{n}<{\delta}_{\mathit{nf}}\\ -{K}_{t}(\vec{{\delta}_{t}}-{\vec{\delta}}_{\mathit{shift}})/{k}_{\mathit{rupt}}^{\tan}& \text{si}& {\delta}_{\mathit{nf}}⩽{\delta}_{n}<{\kappa}_{\mathit{rupt}}^{\tan}\\ 0& \text{si}& {\delta}_{n}⩾{\kappa}_{\mathit{rupt}}^{\tan}\end{array}\right.\end{split}\]

Variables internes

La loi JOINT_MECA_RUPT possède dix-huit variables internes. Du point de vue de la loi de comportement, seules la première et la dixième sont stricto sensu des variables internes. Les autres fournissent des indications sur l’état hydromécanique du joint à un instant donné.

\(V1\)	Seuil \(\kappa\) en saut (plus grande norme atteinte).
\(V2\)	Indicateur de dissipation \(=0\) si régime linéaire, \(=1\) si régime dissipatif.
\(V3\)	Indicateur d’endommagement normal \(=0\) sain, \(=1\) endommagé, \(=2\) cassé
\(V4\)	Pourcentage d’endommagement normal (dans la zone adoucissante) \(V4\in [0,1]\)
\(V5\)	Indicateur d’endommagement tangentiel \(=0\) sain, \(=1\) endommagé, \(=2\) cassé
\(V6\)	Pourcentage d’endommagement tangentiel \(V6\in [0,1]\)
\(V7\)	Saut normal \({\delta}_{n}\) .
\(V8\) \(V9\)	Sauts tangentiels, respectivement \({\delta}_{t1}\) et \({\delta}_{t2}\) (\(V9\) nul en 2D).
\(V10\)	Épaisseur \({\delta}_{\mathit{nf}}\) du joint clavé.
\(V11\)	Contrainte mécanique normale \({\sigma}_{n}\) (sans pression de fluide)
\(V12\) \(V13\) \(V14\)	Composantes du gradient de pression dans le repère global (uniquement pour les modélisationsxxx_JOINT_HYME). Respectivement \({\partial}_{x}p\) , \({\partial}_{y}p\) et \({\partial}_{z}p\) .
\(V15\) \(V16\) \(V17\)	Composantes du flux hydraulique dans le repère global (uniquement pour les modélisations xxx_JOINT_HYME). Respectivement \({w}_{x}\) , \({w}_{y}\) et \({w}_{z}\) .
\(V18\)	Pression de fluide \(p\) imposée par l’utilisateur (PRES_FLUIDE) dans le cas des modélisations xxx_JOINTou pression de fluide interpolée à partir de celle calculée aux nœuds milieux des éléments de joint (degré de liberté du problème) pour lesmodélisations : xxx_JOINT_HYME.
\(V19\) \(V20\)	Deux composantesdu vecteur de glissement tangentiel \(\vec{{\delta}_{\mathit{shift}}}\)

[1]

\(\Gamma\) représente le contour de la fissure

[2]

force par unité de surface, homogène à une contrainte.

[3]

Le vecteur à deux composantes \(\vec{{\delta}_{t}}=({\delta}_{t1},{\delta}_{t2})\) désigne le saut tangentiel.

[4]

Notre formulation sort donc du formalisme énergétique de la rupture. Nous pouvons imaginer une amélioration future de la loi en y introduisant une fonction de régularisation dérivable.

[5]

Dans notre cas \({\sigma}_{n}=\partial {\psi}_{n}({\delta}_{n})/\partial {\delta}_{n};\)

[6]

Les expressions des contraintes peuvent être postulées directement sans passer par la formulation énergétique

Formulation théorique de JOINT_MECA_FROT

On considère la loi de frottement de Coulomb, qui ne dépend que d’un seul paramètre \(\mu \in \left(0,\infty \right]\) . Elle réalise la condition de non-interpénétration des lèvres en contact (condition de Signorini) en établissant un lien local entre la contrainte tangentielle et normale dans la phase de glissement: \(∣∣{\vec{\sigma}}_{t}∣∣=\mu {\sigma}_{n}\) . Plusieurs régularisations de cette loi sont faites afin de faciliter son implémentation numérique. Premièrement la condition de Signorini doit être rendue dérivable, ce qui est aisé si l’on suppose que le comportement des surfaces en contact suit une loi élastique. De même pour la pente de changement de direction de glissement dans la comportement tangentiel. De plus pour la modélisation des barrages en béton on observe expérimentalement une résistance à la traction non négligeable entre les joints. Toutes ces considérations nous ramènent vers une loi de Mohr-Coulomb dont la représentation dans le plan de Mohr est donnée sur la figure Fig. 237, elle décrit la phase de glissement de joints entre la fondation et le barrage ou entre les plots d’un barrage (de façon simplifiée) tout en prenant en compte les effets les plus pertinents.

Fig. 237 Loi de frottement de Coulomb en 2D

La loi JOINT_MECA_FROT est une variante élastoplastique de la loi de Mohr-Coulomb, elle dépend de quatre paramètres : la rigidité normale \({K}_{n}\) , la rigidité tangentielle \({K}_{t}\) , l’adhésion \(c\) (qui est liée à la résistance à la traction maximale \({R}_{t}=c/\mu\) ) et le coefficient de frottement du joint \(\mu\) . Par ailleurs nous introduisons un paramètre d’écrouissage isotrope, qui permet de régulariser la pente tangentielle dans la phase de glissement, on le note \({\rm K}\) . Le modèle élastoplastique introduit ne porte que sur la partie tangente de la loi de comportement. Il n’y a pas de partie plastique du déplacement pour la partie normale: celle-ci est toujours élastique. Le saut de déplacement tangentiel est décomposé en une partie élastique \({\vec{\delta}}_{t}^{\mathit{el}}\) et une partie plastique \({\vec{\delta}}_{t}^{\mathit{pl}}\) , on désigne par \(\lambda\) le saut de déplacement tangentiel cumulé. La loi d’écoulement est orthogonale au plan de coupure du cône de glissement \({\sigma}_{n}=\text{const}\) (cercle 2D pour un cône 3D). Ce qui donne, strictement parlant, une loi d’écoulement globale non-associée. La formulation mécanique en vitesse d’une telle loi donne les équations suivantes:

(3404)\[\begin{split}\left\{ \begin{array}{c}{\vec{\delta}}_{t}={\vec{\delta}}_{t}^{\mathit{el}}+{\vec{\delta}}_{t}^{\mathit{pl}}\\ \vec{{\sigma}_{t}}={K}_{t}{\vec{\delta}}_{t}^{\mathit{el}}\equiv {K}_{t}({\vec{\delta}}_{t}-{\vec{\delta}}_{t}^{\mathit{pl}})\\ {\sigma}_{n}=min({K}_{n}{\delta}_{n},{R}_{t})\end{array}\right.\end{split}\]

Avec:

(3405)\[\begin{split}\left\{ \begin{array}{c}f(\vec{\sigma},\lambda )=\parallel {\vec{\sigma}}_{t}\parallel +\mu {\sigma}_{n}-c-{\rm K}\lambda ⩽0\\ f\cdot \dot{\lambda}=0;\dot{\lambda}⩾0\\ {\dot{\vec{\delta}}}_{t}^{\mathit{pl}}=\dot{\lambda}\frac{{\vec{\sigma}}_{t}}{\parallel {\vec{\sigma}}_{t}\parallel }\end{array}\right.\end{split}\]

Tant qu’on est dans la zone élastique \(f(\vec{\sigma},\lambda )<0\) les relations entre les sauts d’ouverture de joint et les contraintes sont linéaires et le paramètre du saut tangentiel plastique n’évolue pas \({\vec{\delta}}_{t}^{\mathit{pl}}=\mathit{const}\) . Dès qu’on touche les bords du cône de glissement définis par \(f(\vec{\sigma},\lambda )=0\) , l’évolution de saut tangentiel plastique est régi par la loi d’écoulement non-associée (4796). La régularisation de la fonction seuil d’écoulement avec le terme d’écrouissage [1]_ \({\rm K}\lambda >0\) est nécessaire afin de rendre inversible la matrice tangente de la loi et d’éviter donc le problème de solutions multiples dans le cas de chargement en forces imposées. La résistance à la traction du joint varie dans l’intervalle \((0,{R}_{t})\) , c’est une fonction de la contrainte tangentielle, elle est nulle pour une contrainte tangentielle supérieure au paramètre d’adhésion \(c\) et vaut le maximum si la contrainte tangentielle est nulle (voir Fig. ). Dans la version actuelle de la loi, la résistance à la traction maximale n’est pas affectée par le phénomène de glissement, elle n’évolue pas à cause du terme d’écrouissage (Figure ). On suppose aussi qu’une fois atteinte la valeur de la résistance à la traction maximale, la contrainte normale n’évolue plus. Cette dernière hypothèse reste valable tant que le nombre de joints rompus par le cisaillement plastique et par la suite sollicités en traction est négligeable. Afin de pouvoir prendre en compte ce phénomène plus rigoureusement, il est nécessaire d’introduire une régularisation de rupture-frottement dans la phase de glissement en traction, ce qui représente une difficulté théorique assez importante [bib3].

Fig. 238 Evolution du cône de glissement due à l’écrouissage

Discrétisation implicite de la loi de frottement

La version élastoplastique de la loi de frottement est formulée en vitesse, ce qui facilite sa discrétisation numérique. La version incrémentale de la loi reste strictement équivalente à la version continue à condition d’avoir des pas de changement infinitésimaux. Afin de limiter le nombre de pas de chargement, nous adaptons la version continue de la loi à des incréments finis, de façon implicite, c.-a-d. que les conditions de glissement sont écrites dans l’état d’équilibre final. L’algorithme utilisé est celui du retour radial avec prédiction élastique. Par convention nous notons par un signe «-» les variables à l’état d’équilibre précédent, l’état courant est noté par des variables habituelles sans signe supplémentaire (voir [bib8]).

Fig. 239 Loi de frottement de Coulomb

Les équations continues de la loi s’écrivent d’une façon discrétisée:

(3406)\[\begin{split}\left\{ \begin{array}{c}\lambda ={\lambda}^{-}+\Delta \lambda \\ {\vec{\delta}}_{t}^{\mathit{pl}}={\vec{\delta}}_{t}^{\mathit{pl}-}+\Delta {\vec{\delta}}_{t}^{\mathit{pl}}\\ \vec{{\sigma}_{t}}={K}_{t}({\vec{\delta}}_{t}-{\vec{\delta}}_{t}^{\mathit{pl}})\\ {\sigma}_{n}=min({K}_{n}{\delta}_{n},{R}_{t})\end{array}\right.\end{split}\]

Avec la loi suivante:

(3407)\[\begin{split}\left\{ \begin{array}{c}f(\vec{\sigma},\lambda )=\parallel {\vec{\sigma}}_{t}\parallel +\mu {\sigma}_{n}-c-{\rm K}\lambda ⩽0\\ f(\vec{\sigma},\lambda )\cdot \Delta \lambda =0;\Delta \lambda ⩾0\\ {\Delta \vec{\delta}}_{t}^{\mathit{pl}}=\Delta \lambda \frac{{\vec{\sigma}}_{t}}{\parallel {\vec{\sigma}}_{t}\parallel }\end{array}\right.\end{split}\]

Dans un algorithme de Newton les sauts de déplacement \(\vec{\delta}=({\delta}_{n},{\vec{\delta}}_{t})\) ainsi que les variables internes à l’instant précédent \({\vec{\delta}}_{t}^{\mathit{pl}-}\) et \({\lambda}^{-}\) étant connues, pour résoudre la loi il suffit d’obtenir les valeurs de contraintes \(\vec{\sigma}=({\sigma}_{n},{\vec{\sigma}}_{t})\) et toutes les variables internes à l’instant courant (dans notre cas \({\vec{\delta}}_{t}^{\mathit{pl}}\) et \(\lambda\) ). La loi l’équation d’évolution pour la composante normale est complètement découplée du mouvement tangentiel, on peut donc la résoudre immédiatement:

(3408)\[{\sigma}_{n}=min({K}_{n}{\delta}_{n},{R}_{t})\]

Nous obtenons donc un jeu de cinq d’équations et une inégalité, pour obtenir cinq inconnues scalaires:

(3409)\[\vec{{\sigma}_{t}}={K}_{t}({\vec{\delta}}_{t}-{\vec{\delta}}_{t}^{\mathit{pl}-}-\Delta {\vec{\delta}}_{t}^{\mathit{pl}})\]

Et:

(3410)\[\begin{split}\left\{ \begin{array}{c}f(\vec{\sigma},\lambda )\equiv \parallel {\vec{\sigma}}_{t}\parallel +\mu {\sigma}_{n}-c-{\rm K}{\lambda}^{-}-{\rm K}\Delta \lambda ⩽0\\ f(\vec{\sigma},\lambda )\cdot \Delta \lambda =0;\Delta \lambda ⩾0\\ {\Delta \vec{\delta}}_{t}^{\mathit{pl}}=\Delta \lambda \frac{{\vec{\sigma}}_{t}}{\parallel {\vec{\sigma}}_{t}\parallel }\end{array}\right.\end{split}\]

Pour simplifier ce problème mathématique regardons plus en détail la condition de Karush-Kuhn-Tucker \(f(\vec{\sigma},\lambda )\cdot \Delta \lambda =0\) . On a deux possibilités: soit on glisse \(\Delta \lambda >0\) , soit on est dans le domaine élastique \(\Delta \lambda =0\) . Si l’on est dans le domaine élastique, alors \(\Delta {\vec{\delta}}_{t}^{\mathit{pl}}=0\) et on obtient la solution élastique si, et seulement si:

(3411)\[{f}_{\mathit{el}}({\sigma}_{n},{\vec{\sigma}}_{\tau}^{-},{\lambda}^{-})\equiv {K}_{t}\parallel {\vec{\delta}}_{\tau}-{\vec{\delta}}_{\tau}^{\mathit{pl}-}\parallel +\mu {\sigma}_{n}-c-{\rm K}{\lambda}^{-}⩽0\]

Dans la pratique si (49) est satisfaite alors la prédiction élastique est la solution du problème:

(3412)\[\begin{split}\left\{ \begin{array}{c}\lambda ={\lambda}^{-}\\ {\vec{\delta}}_{t}^{\mathit{pl}}={\vec{\delta}}_{t}^{\mathit{pl}-}\\ \vec{{\sigma}_{t}}={K}_{t}({\vec{\delta}}_{t}-{\vec{\delta}}_{t}^{\mathit{pl}-})\\ {\sigma}_{n}=min({K}_{n}{\delta}_{n},{R}_{t})\end{array}\right.\end{split}\]

Si la condition (49) n’est pas satisfaite, alors \(\Delta \lambda >0\) et on est dans la phase de glissement. On obtient alors un système de trois équations non-linéaires avec trois inconnues \(\Delta {\vec{\delta}}_{t}^{\mathit{pl}}\) et \(\Delta \lambda\) :

(3413)\[\begin{split}\left\{ \begin{array}{c}f(\vec{\sigma},\lambda )\equiv \parallel {K}_{t}({\vec{\delta}}_{t}-{\vec{\delta}}_{t}^{\mathit{pl}-}-\Delta {\vec{\delta}}_{t}^{\mathit{pl}})\parallel +\mu {\sigma}_{n}-c-{\rm K}{\lambda}^{-}-{\rm K}\Delta \lambda =0\\ {\Delta \vec{\delta}}_{t}^{\mathit{pl}}=\Delta \lambda \frac{({\vec{\delta}}_{t}-{\vec{\delta}}_{t}^{\mathit{pl}-}-\Delta {\vec{\delta}}_{t}^{\mathit{pl}})}{\parallel {\vec{\delta}}_{t}-{\vec{\delta}}_{t}^{\mathit{pl}-}-\Delta {\vec{\delta}}_{t}^{\mathit{pl}}\parallel }\end{array}\right.\end{split}\]

Cette équation peut-être résolue en éliminant \(\Delta {\overrightarrow{\delta}}_{t}^{\mathrm{pl}}\) , suivant la procédure couramment utilisée pour les calculs en plasticité:

(3414)\[{\Delta \vec{\delta}}_{t}^{\mathit{pl}}\lbrace \parallel {\vec{\delta}}_{t}-{\vec{\delta}}_{t}^{\mathit{pl}-}-\Delta {\vec{\delta}}_{t}^{\mathit{pl}}\parallel +\Delta \lambda \rbrace =\Delta \lambda ({\vec{\delta}}_{t}-{\vec{\delta}}_{t}^{\mathit{pl}-})\]

En prenant la norme de cette dernière équation et en notant que \(\parallel \Delta {\overrightarrow{\delta}}_{t}^{\mathrm{pl}}\parallel =\Delta \lambda\) , on obtient la norme du vecteur de la contrainte tangentielle actualisée, qu’on insère dans l’équation (51) afin d’obtenir une équation scalaire pour \(\Delta \lambda\) :

(3415)\[\parallel {\vec{\delta}}_{t}-{\vec{\delta}}_{t}^{\mathit{pl}-}-\Delta {\vec{\delta}}_{t}^{\mathit{pl}}\parallel =\parallel {\vec{\delta}}_{t}-{\vec{\delta}}_{t}^{\mathit{pl}-}\parallel -\Delta \lambda\]

Et:

(3416)\[\parallel {K}_{t}({\vec{\delta}}_{t}-{\vec{\delta}}_{t}^{\mathit{pl}-})\parallel -{K}_{t}\Delta \lambda +\mu {\sigma}_{n}-c-{\rm K}{\lambda}^{-}-{\rm K}\Delta \lambda =0\]

Une fois \(\Delta \lambda\) connu, il suffit de remarquer la colinéarité des vecteurs suivant \(\Delta {\overrightarrow{\delta}}_{t}^{\mathrm{pl}}\uparrow \uparrow {\overrightarrow{\delta}}_{t}-{\overrightarrow{\delta}}_{t}^{\mathrm{pl}-}-\Delta {\overrightarrow{\delta}}_{t}^{\mathrm{pl}}\) , d’où \(\Delta {\overrightarrow{\delta}}_{t}^{\mathrm{pl}}\uparrow \uparrow {\overrightarrow{\delta}}_{t}-{\overrightarrow{\delta}}_{t}^{\mathrm{pl}-}\) . Ce qui permet de réécrire l’équation (54) sous une forme simplifiée, qui donne la valeur de la deuxième inconnue \(\Delta {\overrightarrow{\delta}}_{t}^{\mathrm{pl}}\) :

(3417)\[{\Delta \vec{\delta}}_{t}^{\mathit{pl}}=\Delta \lambda \frac{{\vec{\delta}}_{t}-{\vec{\delta}}_{t}^{\mathit{pl}-}-\Delta {\vec{\delta}}_{t}^{\mathit{pl}}}{\parallel {\vec{\delta}}_{t}-{\vec{\delta}}_{t}^{\mathit{pl}-}-\Delta {\vec{\delta}}_{t}^{\mathit{pl}}\parallel }=\Delta \lambda \frac{{\vec{\delta}}_{t}-{\vec{\delta}}_{t}^{\mathit{pl}-}}{\parallel {\vec{\delta}}_{t}-{\vec{\delta}}_{t}^{\mathit{pl}-}\parallel }\]

Cette solution correspond en fait aux glissements dans la direction de contrainte tangentielle en prédiction élastique. Ceci implique que le changement de la direction de glissement se fera essentiellement dans la zone élastique à condition que les pas de chargement soient petits. La solution finale en cas de glissement, obtenue à partir de l’équation (55), s’écrit comme:

(3418)\[\begin{split}\left\{ \begin{array}{ccc}{\sigma}_{n}& =& min({K}_{n}{\delta}_{n},{R}_{t})\\ {f}_{\mathit{el}}({\sigma}_{n},{\vec{\sigma}}_{t}^{-},{\lambda}^{-})& =& {K}_{t}\parallel {\vec{\delta}}_{t}-{\vec{\delta}}_{t}^{\mathit{pl}-}\parallel +\mu {\sigma}_{n}-c-{\rm K}{\lambda}^{-}\\ \lambda & =& {\lambda}^{-}+\frac{{f}_{\mathit{el}}({\sigma}_{n},{\vec{\sigma}}_{t}^{-},{\lambda}^{-})}{{K}_{t}+{\rm K}}\\ {\vec{\delta}}_{t}^{\mathit{pl}}& =& {\vec{\delta}}_{t}^{\mathit{pl}-}+\frac{{f}_{\mathit{el}}({\sigma}_{n},{\vec{\sigma}}_{t}^{-},{\lambda}^{-})}{{K}_{t}+{\rm K}}\cdot \frac{{\vec{\delta}}_{t}-{\vec{\delta}}_{t}^{\mathit{pl}-}}{\parallel {\vec{\delta}}_{t}-{\vec{\delta}}_{t}^{\mathit{pl}-}\parallel }\\ \vec{{\sigma}_{t}}& =& {K}_{t}({\vec{\delta}}_{t}-{\vec{\delta}}_{t}^{\mathit{pl}})\end{array}\right.\end{split}\]

En résumé, à partir de la prédiction élastique on vérifie d’abord l’inégalité (49), si elle est satisfaite, alors la solution est donnée par (50), sinon la solution est donnée par l’équation (56).

Matrice tangente

Pour la loi JOINT_MECA_FROT la matrice tangente est calculée en implicite, ce qui renforce la robustesse des calculs [2]_. Comme il est démontré dans la référence [bib7], un tel schéma numérique est inconditionnellement stable pour les lois à écrouissage positif \(K\ge 0\) . Dans le cas du régime élastique (l’inégalité (58) est satisfaite), la matrice tangente prend une forme simple, elle est diagonale:

(3419)\[\begin{split}\left(\begin{array}{ccc}{K}_{n}& 0& 0\\ 0& {K}_{t}& 0\\ 0& 0& {K}_{t}\end{array}\right)\end{split}\]

Dans le cas de glissement (l’inégalité (58) n’est pas satisfaite) la matrice tangente est obtenue par dérivation des équations (56). Les dérivées par rapport à l’ouverture normale dépendent de l’état du joint. Pour le joint fermé cela donne:

(3420)\[\begin{split}\text{si}{\delta}_{n}<\frac{c}{\mu {K}_{n}}\text{}\Rightarrow \text{}\left\{ \begin{array}{ccc}\frac{\partial {\sigma}_{n}}{\partial {\delta}_{n}}& =& {K}_{n}\\ \frac{\partial \vec{{\sigma}_{t}}}{\partial {\delta}_{n}}& =& -\mu \frac{{\vec{\delta}}_{t}-{\vec{\delta}}_{t}^{\mathit{pl}-}}{\parallel {\vec{\delta}}_{t}-{\vec{\delta}}_{t}^{\mathit{pl}-}\parallel }\cdot \frac{{K}_{n}{K}_{t}}{{K}_{t}+{\rm K}}\end{array}\right.\end{split}\]

Pour le joint ouvert toutes les dérivées correspondantes sont nulles:

(3421)\[\begin{split}\text{si}{\delta}_{n}⩾\frac{c}{\mu {K}_{n}}\text{}\Rightarrow \text{}\left\{ \begin{array}{ccc}\frac{\partial {\sigma}_{n}}{\partial {\delta}_{n}}& =& 0\\ \frac{\partial \vec{{\sigma}_{t}}}{\partial {\delta}_{n}}& =& 0\end{array}\right.\end{split}\]

Les dérivées par rapport à l’ouverture tangentielle ne dépendent pas de l’ouverture de joint:

(3422)\[\begin{split}\left\{ \begin{array}{ccc}\frac{\partial {\sigma}_{n}}{\partial {\vec{\delta}}_{t}}& =& 0\\ \frac{\partial \vec{{\sigma}_{t}}}{\partial {\vec{\delta}}_{t}}& =& \frac{{\rm K}{K}_{t}}{{K}_{t}+{\rm K}}\mathrm{Id}+\frac{-\mu {\sigma}_{n}+c+{\rm K}{\lambda}^{-}}{{K}_{t}+{\rm K}}\cdot \frac{{K}_{t}}{\parallel {\vec{\delta}}_{t}-{\vec{\delta}}_{t}^{\mathit{pl}-}\parallel }\cdot \left(\mathrm{Id}-\frac{({\vec{\delta}}_{t}-{\vec{\delta}}_{t}^{\mathit{pl}-})\otimes ({\vec{\delta}}_{t}-{\vec{\delta}}_{t}^{\mathit{pl}-})}{{\parallel {\vec{\delta}}_{t}-{\vec{\delta}}_{t}^{\mathit{pl}-}\parallel }^{2}}\right)\end{array}\right.\end{split}\]

(3423)\[\begin{split}\left\{ \begin{array}{ccc}\frac{\partial {\sigma}_{n}}{\partial {\vec{\delta}}_{t}}& =& 0\\ \frac{\partial \vec{{\sigma}_{t}}}{\partial {\vec{\delta}}_{t}}& =& \frac{{\rm K}{K}_{t}}{{K}_{t}+{\rm K}}\mathrm{Id}+\frac{-\mu {\sigma}_{n}+c+{\rm K}{\lambda}^{-}}{{K}_{t}+{\rm K}}\cdot \frac{{K}_{t}}{\parallel {\vec{\delta}}_{t}-{\vec{\delta}}_{t}^{\mathit{pl}-}\parallel }\cdot \left(\mathrm{Id}-\frac{({\vec{\delta}}_{t}-{\vec{\delta}}_{t}^{\mathit{pl}-})\otimes ({\vec{\delta}}_{t}-{\vec{\delta}}_{t}^{\mathit{pl}-})}{{\parallel {\vec{\delta}}_{t}-{\vec{\delta}}_{t}^{\mathit{pl}-}\parallel }^{2}}\right)\end{array}\right.\end{split}\]

Remarque

La matrice tangente dans la phase plastique (en glissement) est non-symétrique, elle est dégénérée si l’écrouissage est nul (\({\rm K}=0\) ). On peut exhiber un vecteur propre associé à la valeur propre nulle:

(3424)\[\begin{split}\left( \begin{array}{ccc} K_n & & 0 \\ \mu K_n \dfrac{\overrightarrow{\delta}_t - \overrightarrow{\delta}_t^{\mathrm{pl}-}}{\parallel \overrightarrow{\delta}_t - \overrightarrow{\delta}_t^{\mathrm{pl}-} \parallel} & & \dfrac{-\mu \sigma_n + c}{\parallel \overrightarrow{\delta}_t - \overrightarrow{\delta}_t^{\mathrm{pl}-} \parallel} \left( \mathrm{Id} - \dfrac{ (\overrightarrow{\delta}_t - \overrightarrow{\delta}_t^{\mathrm{pl}-}) \otimes (\overrightarrow{\delta}_t - \overrightarrow{\delta}_t^{\mathrm{pl}-}) }{ \parallel \overrightarrow{\delta}_t - \overrightarrow{\delta}_t^{\mathrm{pl}-} \parallel^2 } \right) \end{array} \right) \cdot \left( \begin{array}{c} 0 \\ \\ \overrightarrow{\delta}_t - \overrightarrow{\delta}_t^{\mathrm{pl}-} \end{array} \right) = \left( \begin{array}{c} 0 \\ 0 \\ 0 \end{array} \right)\end{split}\]

C’est pour cette raison que le paramètre d’écrouissage isotrope est introduit.

Remarque

Dans les études on constate une difficulté de convergence liée au fait qu’on passe de l’état de plastification (par exemple: écrouissage dû au chargement en force imposés) vers la solution élastique. La matrice tangente dans la phase de prédiction n’est pas «bonne» pour récupérer la phase élastique et le calcul ne converge que lorsque le pas est très découpé. Pour avoirla convergence il faut juste faire un pas avec la matrice tangente élastique. Il suffit d’utiliser le mot clef PAS_MINI_ELAS dans le mot-clef facteur NEWTON de STAT_NON_LINE ou DYNA_NON_LINE.

Variables internes

La loi JOINT_MECA_FROT possède dix-huit variables internes. Du point de vue de la loi de comportement, seules la première, la troisième et la quatrième sont stricto sensu des variables internes. Les autres fournissent des indications sur l’état de hydromécanique du joint à un instant donné.

\(V1\)	Paramètre croissant indiquant le déplacement tangentiel plastiquecumulée \(\lambda\) .
\(V2\)	Indicateur de glissement \(=0\) si régime linéaire élastique, \(=1\) si régime plastique.
\(V3\) \(V4\)	Vecteur de plasticité tangentielle \({\vec{\delta}}^{\mathit{pl}}\) (indique la position d’équilibre ). \(\mathrm{V4}\) est mis à zéro en 2D.
\(V5\)	Indicateur d’ouverture complète \(=0\) fermé \(({\sigma}_{n}<c/\mu )\) , \(=1\) ouvert \(({\sigma}_{n}=c/\mu )\) .
\(V6\)	Norme de la contrainte tangente \(\Vert {\vec{\sigma}}_{\tau}\Vert\) .
\(V7\) \(V8\) \(V9\)	Valeur du saut dans le repère local : saut normal \({\delta}_{n}\) , saut tangentiel 1 \({\delta}_{t2}\) et saut tangentiel 2 \({\delta}_{t2}\) (\(V9\) nul en 2D).
\(V10\)	Variable non utilisée.
\(V11\)	Contrainte mécanique normale \({\sigma}_{n}\) (sans pression de fluide).
\(V12\) \(V13\) \(V14\)	Composantes du gradient de pression dans le repère global (uniquement pour les modélisations xxx_JOINT_HYME). Respectivement \({\partial}_{x}p\) , \({\partial}_{y}p\) et \({\partial}_{z}p\) .
\(V15\) \(V16\) \(V17\)	Composantes du flux hydraulique dans le repère global (uniquement pour les modélisations xxx_JOINT_HYME). Respectivement \({w}_{x}\) , \({w}_{y}\) et \({w}_{z}\) .
\(V18\)	Pression de fluide \(p\) imposée par l’utilisateur (PRES_FLUIDE) dans le cas des modélisations xxx_JOINT ou pression de fluide interpolée à partir de celle calculée aux nœuds milieux des éléments de joint (degré de liberté du problème) pour les modélisations : xxx_JOINT_HYME.

Prise en compte de l’amortissement en dynamique

L’amortissement dynamique dans chaque élément de joint 3D_JOINT est représenté par une matrice d’amortissement \(\mathrm{C}\) , calculée à chaque pas de calcul par l’option AMOR_MECA pour ce type de modélisation. L’opérateur assemblé \(\mathrm{C}\) intervient alors dans le premier membre de l’équation d’équilibre du problème dynamique et de ce fait elle est invariante des schémas d’intégration.

Fig. 240 Prise en compte de l’amortissement

Les termes de cette matrice sont obtenus par l’intégration de densités surfaciques d’amortissement normal \({A}_{n}\) (renseigné par le mot clé AMOR_NOR du comportement JOINT_MECA_FROT) ou tangentiel \({A}_{t}\) (renseigné par le mot clé AMOR_TAN du comportement JOINT_MECA_FROT) intégrées sur la surface \(\mathit{Surf}\) d’une face d’élément 3D_JOINT.

Ces termes sont équivalents à la répartition en parallèle de caractéristiques de discret de type DIS_T sur chaque segment joignant chaque couple de nœuds sommets en vis-à-vis d’une face à l’autre de l’élément (voir la figure Fig. 240). Ces caractéristiques sont affectées avec leur pleine valeur seulement si l’élément de joint est en compression: soit si la septième composante de variable interne du comportement JOINT_MECA_FROT, valeur du saut dans le repère local, \(\mathit{V7}={\delta}_{n}\) est négative.

Les éléments de la matrice \(\mathrm{C}\) sont directement calculées pour les degrés de liberté des nœuds sommets. En conséquence, si le nombre de nœuds sommets d’une face d’élément de joint est \({n}_{s}\) , on va affecter pour la matrice d’amortissement mécanique les caractéristiques locales \({A}_{n}.\mathit{Surf}/{n}_{s}\) ou \({A}_{t}.\mathit{Surf}/{n}_{s}\) réparties dans le repère local de l’élément sur chacun des termes diagonaux. On affecte leurs valeurs opposées aux termes de couplage. C’est donc comme si on avait \({n}_{s}\) discrets d’amortissement de type DIS_T affectés sur des mailles SEG2.

Dans le cas où l’élément de joint n’est pas en compression: soit si la septième composante de variable interne du comportement JOINT_MECA_FROT, valeur du saut dans le repère local, \(\mathit{V7}={\delta}_{n}\) n’est pas négative, alors on affecte les valeurs calculées précédentes d’un coefficient renseigné par le mot clé COEF_AMOR (nul par défaut) dans la loi de comportement JOINT_MECA_FROT.

[1]

Écrouissage isotrope sur le plan tangent

[2]

Dans les tentatives précédentes [bib5], [bib2], d’introduction en version explicite [bib4], il est apparu que ces modélisations s’avèrent généralement peu performantes numériquement, leur extension en 3D est aussi assez délicate [bib4]

Formulation théorique de JOINT_MECA_ENDO

L’un des points clés dans la description précise des interfaces entre les roches est l’identification de leur domaine élastique (nommé également surface de charge ou le critère). En raison leur nature adoucissante, non seulement les contacts rocheux mais aussi toute sorte d’interface analogue faisant intervenir les matériaux similaires comme le béton, l’argile ou le sol peuvent généralement être classés en fonction de leur résistance à divers chargements en mode mixte. Cette approche est étroitement liée aux règles de sécurité de base dans les applications industrielles, où il est admis en première approximation que les géomatériaux pourraient présenter une défaillance instable une fois la charge critique atteinte.

Tout au long du siècle dernier, des machines d’essais spécifiques et des protocoles de mesures correspondants ont été mis en place afin d’harmoniser et de standardiser la caractérisation expérimentale des interfaces fragiles cités ci-dessus. Pour n’en citer que quelques-unes, le test de traction simple ou encore le cisaillement à compression fixée sont toutes les deux des expériences bien documentées qui sont exécutées en pratique pour cataloguer la résistance des matériaux en repérant seulement quelques points de leur surface d’élastique multidimensionnelle. Cette vision réduite « en un seul point » de la résistance des matériaux est de plus en plus scrutée ces derniers temps, et de multiples évolutions ont été adoptées. Par exemple, les améliorations constantes des logiciels d’éléments finis permettent de nouveaux types de modélisation, avec des chargements allant au-delà du domaine élastique pour explorer un comportement post-pic plus subtil. Alors que la complexité de la description post-pic pourrait être atteinte à travers de divers formalismes théoriques, la plupart d’entre eux reposent encore sur la définition initiale de la surface de charge, et la question de cette forme de domaine élastique reste la pierre angulaire de toute identification de modèle non linéaire. Même si des progrès considérables ont été réalisés au cours des dernières décennies, cette identification précise de la surface de charge reste aujourd’hui une tâche assez difficile.

Une caractéristique commune des géomatériaux est leur forte résistance aux charges de compression. Dans des applications industrielle courante c’est précisément cette propriété de résistance à la compression plus élevée qui est exploitée dans le but d’augmenter la robustesse globale de la structure. Comme les interfaces considérés sont composées des géomatériaux il est naturel de supposer une similitude entre les surfaces de charges correspondantes. Pour les géomatériaux, selon leur origine physique, une grande variété de critères définissant le domaine élastique sont employés afin de modéliser leur comportement mécanique. La surface la plus simple, admettant une résistance infinie en compression hydrostatique, est celle en forme de cône linéaire. Il a été introduit pour la première fois il y a plus d’un siècle [bib9] comme la combinaison de l’hypothèse de frottement de Coulomb [bib10] avec le principe de coupure de tension de Galileo-Rankine. La description initiale de la contrainte principale, communément appelée la forme de Mohr-Coulomb, a ensuite été généralisée dans les travaux de Drucker et Prager à une relation déviateur-trace plus lisse [bib11]. La relation directe entre les chargements de cisaillement et de compression, qui est la signature principale de ces critères linéaires, a permis le développement de diverses relations constitutives basées sur la même dépendance simplifiée [bib12]. Nous retenons à ce stade que la surface conique linéaire est couramment utilisé à ce jour pour la modélisation de géomatériaux.

Avec de nombreux grands projets d’infrastructure en cours, des critères plus complexes ont émergé à mesure que de nouvelles données expérimentales devenaient disponibles. Pour des gammes plus larges de chargements, la dépendance de cisaillement en compression semblait s’écarter de la courbe linéaire. En 1980 Evert Hoek et Edwin T. Brown [bib13] proposent une nouvelle forme particulière de critère d’une forme quadratique. Il reproduit la singularité de type Mohr-Coulomb pour un faible chargement en traction en tenant compte simultanément de la réduction de la résistance au cisaillement pour une forte compression. Ce critère purement empirique a été largement utilisé dans la conception des excavations souterraines [bib14]. Les travaux plus récents sur les joints rocheux semblent montrer la même modification de la surface de charge.

La loi JOINT_MECA_ENDO tente de capter partiellement ce phénomène, en partant de la surface de charge élastique linéaire cette dernière évolue grâce à l’écrouissage cinématique, la résistance à la traction est captée par le couplage endommageant de l’amplitude d’écrouissage. La forme apparente de la surface de charge n’est plus linéaire et de ce fait elle s’approche des surfaces de charge des joints réels. Nous avons également l’avantage de pouvoir représenter le modèle des zones cohésives et le modèle élastoplastique au sein de même et unique formalisme, que nous décrivons rapidement ci-dessous.

Comme les lois précédente la loi JOINT_MECA_ENDO accepte également une modélisation couplée hydromécanique, le couplage HM y est identique et ces deux phénomènes peuvent être traités séparément. Nous décrivons ici la partie mécanique de la loi, qui unifie les deux phénomènes mécanique: la rupture et le frottement. Elle s’appuie sur les modélisations XXX_JOINT (R3.06.09).

Remarque :Les deux prochaines sections présentent le concept théorique général. En première lecture on peut les omettre et passer directement à la section suivante, qui fournit suffisamment d’éléments pour comprendre l’implémentation numérique de la loi.

Modèle standard généralisé pour les joints

Le cadre théorique choisi pour modéliser le couplage entre la rupture et le frottement des joints de barrages en béton est basé sur le modèle standard généralisé [bib18], [bib12] , [bib20]. Ce modèle permet de garantir des bonnes propriétés thermodynamique en opérant exclusivement avec les notions énergétique. Initialement développé pour représenter les phénomènes tridimensionnels nous l’adaptons pour le modèle des joints.

En mécanique des interfaces, deux grandeurs physiques sont typiquement introduites. Le premier est le vecteur de la densité des contraintes de Cauchy, qui caractérise l’état des efforts locaux, tandis que le second est le vecteur de saut de déplacement, qui permet de mesurer la position relatives des lèvres de part et d’autre de l’interface à une certaine configuration de référence près.

Si la configuration de référence est connue, les deux grandeurs sont objectives et mesurables expérimentalement. Par conséquent un des buts de la mécanique des interfaces est d’établir un lien entre ces deux grandeurs, que l’on appelle relation constitutive matérielle ou la loi de comportement.

En suivant les idées originale de Cauchy et de la manière la plus simple nous pouvons postuler que la contrainte est une fonction bijective du saut de déplacement. Pas toutes les relations contrainte-saut bijective ne sont toutefois admissibles d’un point de vue thermodynamique. Par exemple, le travail des efforts internes pour les relations constitutives avec un opérateur tangent non symétrique dépend du chemin de chargement, ce qui est physiquement inadmissible pour un processus élastique entièrement réversible. L’hyperélasticité restreint cette relation contrainte-saut en admettant l’existence d’une fonction scalaire de densité d’énergie de déformation à partir de laquelle la contrainte est dérivée. Le chargement mécanique d’un matériau hyperélastique est considéré comme un processus réversible à dissipation nulle, la seule variable d’état étant le saut de déplacement. En ce sens, le formalisme hyper-élastique est conservatif et compatible avec une description thermodynamique générale [bib15]. Les idées sont similaires pour éteindre la notion des lois admissibles thermodynamiquement pour les comportements en présence de variables enregistrant l’historique de chargement.

Afin de bien prendre en compte le processus de dissipation d’énergie, bien présent en cas de frottement de lèvres, il est nécessaire d’introduire au moins une variable d’état supplémentaire. Cela se fait généralement via le vecteur de glissement, appelé le saut de déplacement plastique. Ce dernier est censé décomposer le saut de déplacement total en contributions élastiques et plastiques non réversible. Bien que cette variable plastique ait une signification physique claire de déformation résiduelle pour le système totalement déchargé, elle n’est généralement pas directement observable, c’est pourquoi elle est également appelée variable d’état interne (cachée). Afin de capter non seulement les glissements irréversibles, mais également prendre en compte l’atténuation de la rugosité de l’interface en phase de glissement une autre variable interne doit être introduit. Nous faisons ici le choix d’introduire une variable scalaire d’endommagement, dont la particularité est d’être contenue dans le borne [0,1]. Zéro correspond à l’état sain et un - à l’état complètement endommagé.

La description thermodynamique de la mécanique des interface repose soit sur la définition de la densité d’énergie libre de Helmholtz [bib15], soit sur l’énergie libre de Gibbs [bib16], [bib17].

Comme il a été expliqué ci-dessus, pour l’évolution isotherme en présence de plasticité et d’endommagement, cette fonction dépend de trois variables d’état, dont deux sont cachées. Nous faisons un choix de description de Helmholtz avec le saut de déplacement total. La projection sur l’interface du tenseur des contraintes de Cauchy est alors le conjugué dual du saut de déplacement total \(\vec{\sigma}=\partial \psi /\partial \vec{\delta}\) , tandis que la réponse énergétique du système sur l’évolution de la plasticité et d’endommagement est capturée par les efforts généralisées correspondantes \(\vec{X}=-\frac{\partial \psi }{\partial \vec{p}};Y=-\frac{\partial \psi }{\partial \alpha }\) .

La plupart des modèles élasto-plastiques admettent une séparation additive des déformations élastiques et plastiques. Cela conduit à la forme quadratique la plus simple de l’énergie libre décrivant l’état plastique résiduel. Nous adaptons la même hypothèse pour les joints :

(3425)\[\psi (\vec{\delta},\vec{p})=K\frac{{(\vec{\delta}-\vec{p})}^{2}}{2}\]

où \(K\) est le module d’élasticité tensoriel supposé constant. Notez que nous adoptons des notations abrégées pour simplifier la description. En particulier, dans ce cas élastique linéaire, l’effort plastique généralisé et la contrainte de Cauchy sont égales : \(\vec{\sigma}=\vec{X}\) .

En ajoutant le terme quadratique plastique dans l’énergie libre on introduit de facto un écrouissage cinématique, car

(3426)\[ \begin{align}\begin{aligned}\psi (\vec{\delta},\vec{p},\alpha )=K\frac{{(\vec{\delta}-\vec{p})}^{2}}{2}+H(\alpha )\frac{{\vec{p}}^{2}}{2}\\\textrm{et}\\\vec{X}=\vec{\sigma}-H(\alpha )\vec{p}\end{aligned}\end{align} \]

\(H\) joue le rôle d’amplitude d’écrouissage cinématique. Afin de reproduire le phénomène d’adoucissement nous supposons que l’amplitude d’écrouissage diminue avec le développement de plasticité. Le domaine élastique réversible, caractérisé par l’absence d’évolution de la plasticité et d’endommagement, est défini dans l’espace des contraintes observables :

\(\vec{p}=\text{const et}\alpha =\text{const}\to \vec{\sigma}\in \text{domaine élastique}\) .

Thermodynamiquement, il faut s’assurer que le taux de dissipation d’énergie soit positif :

(3427)\[\text{D}=\vec{\sigma}\cdot \dot{\vec{\delta}}-\dot{\psi}\ge 0\]

où \(\dot{}\) correspond à la dérivée temporelle: \(\dot{\psi}=\partial \psi /\partial t\) . En considérant d’abord la décharge élastique, nous obtenons la relation \(\vec{\sigma}=\frac{\partial \psi }{\partial \vec{\delta}}\) comme conséquence du processus à dissipation nulle \(\text{D}\left(\vec{p}=\text{const};\alpha =\text{const}\right)=0\) , réduisant alors le second principe thermodynamique eq. à l’inégalité directe suivante : \(\text{D}=\vec{X}\cdot \dot{\vec{p}}+Y\cdot \dot{\alpha}\ge 0\) .

Pour obtenir une description complète du matériau, il faut introduire les lois d’évolutions des variables cachées: \(\dot{\vec{p}}\text{et}\dot{\alpha}\) . Cela se fait soit via un potentiel de dissipation [bib18], soit via un potentiel plastique [bib19] ou directement via une règle d’écoulement explicite. Pour un écoulement arbitraire, la preuve de la positivité du taux de dissipation eq. est une tâche significativement difficile. C’est pourquoi des hypothèses restrictives supplémentaires sont couramment introduites. Par exemple, si l’évolutionde variables cachées dépendent exclusivement des efforts généralisés \((\dot{\vec{p}};\dot{\alpha})\iff (\vec{X};Y)\) , le taux de dissipation du matériau devient une fonction exclusive de ces efforts \(\text{D}=\text{D}(\vec{X};Y)\) . Le principal avantage d’une telle description est que la positivité de la dissipation est indépendante du chemin de chargement : celle-ci dépend exclusivement du choix de fonctions qui pilotent l’écoulement. Ainsi l’ensemble du modèle peut être analysé dans l’espace des forces généralisées [bib18].

Nous nous inspirons des travaux récents de [bib12] en introduisant l’énergie libre de Helmholtz pour les interfaces de la forme suivante :

(3428)\[\psi (\vec{\delta},\vec{p},\alpha )={K}_{n}\frac{{({\delta}_{n}-{p}_{n})}^{2}}{2}+{K}_{t}\frac{{\Vert \vec{{\delta}_{t}}-\vec{{p}_{t}}\Vert }^{2}}{2}+{A}_{n}(\alpha )\frac{{{p}_{n}}^{2}}{2}+{A}_{t}(\alpha )\frac{{\vec{{p}_{t}}}^{2}}{2}\]

avec \({A}_{n}(\alpha )={B}_{n}\cdot {(1-\alpha )}^{3}/\sqrt{\alpha}\) et \({A}_{t}(\alpha )={B}_{t}\cdot {(1-\alpha )}^{3}/\sqrt{\alpha}\) . Nous utilisons la décomposition du saut de déplacement en partie normale et tangentielle \(\vec{\delta}=({\delta}_{n},{\delta}_{t1},{\delta}_{t2})=({\delta}_{n},{\vec{\delta}}_{t})\) . Les deux rigidités en sollicitations normale \({K}_{n}\) et tangentielle \({K}_{t}\) , qui caractérisent la surface et les matériaux de remplissage de la fissure (la même signification que pour les lois JOINT_MECA_RUPT et JOINT_MECA_FROT). Les deux amplitudes de couplage cinématique (normale \({B}_{n}\) et tangentielle \({B}_{t}\) ) sont également les paramètres matériaux.

Les expressions des contraintes et d’efforts généralisés donnent un jeux d’équations suivant:

(3429)\[ \begin{align}\begin{aligned}\begin{split}\begin{array}{c}{\sigma}_{n}=\frac{\partial \psi (\vec{\delta},\vec{p},\alpha )}{\partial {\delta}_{n}}={K}_{n}({\delta}_{n}-{p}_{n})\\ {X}_{n}=\frac{\partial \psi (\vec{\delta},\vec{p},\alpha )}{\partial {p}_{n}}={K}_{n}({\delta}_{n}-{p}_{n})-{A}_{n}\cdot {p}_{n}\end{array}\end{split}\\\textrm{et}\\\begin{split}\begin{array}{c}\vec{{\sigma}_{t}}=\frac{\partial \psi (\vec{\delta},\vec{p},\alpha )}{\partial {\delta}_{t}}={K}_{t}(\vec{{\delta}_{t}}-\vec{{p}_{t}})\\ \vec{{X}_{t}}=\frac{\partial \psi (\vec{\delta},\vec{p},\alpha )}{\partial \vec{{p}_{t}}}={K}_{t}(\vec{{\delta}_{t}}-\vec{{p}_{t}})-{A}_{t}\cdot \vec{{p}_{t}}\end{array}\end{split}\\Y=-{A}_{n}'{p}_{n}^{2}/2-{A}_{t}'{\Vert \vec{{p}_{t}}\Vert }^{2}/2\end{aligned}\end{align} \]

Où \({A}_{n}'=\frac{d{A}_{n}(\alpha )}{d\alpha }\) . En choisissant le domaine élastique comme l’intersection de deux domaines définit dans les espaces de contrainte plastique généralisé (cône de révolution autours de \({X}_{n}\) ) et d’effort généralisé dual à la variable d’endommagement nous obtenons la loi d’écoulement suivante:

(3430)\[ \begin{align}\begin{aligned}\begin{split}\begin{array}{c}{f}_{X}(\vec{X})=\Vert \vec{{X}_{t}}\Vert +\mu {X}_{n}-\overline{c}\le 0\\ \dot{\vec{p}}=\dot{\lambda}\frac{\partial {f}_{X}}{\partial \vec{X}}\\ {f}_{X}(\vec{X})\dot{\lambda}=0;\dot{\lambda}\ge 0\end{array}\end{split}\\\textrm{et}\\\begin{split}\begin{array}{c}Y=-{A}_{n}^{I}{p}_{n}^{2}/2-{A}_{t}'{\Vert \vec{{p}_{t}}\Vert }^{2}/2\le {D}_{1}\\ \dot{\alpha}\ge 0\\ (Y-{D}_{1})\cdot \dot{\alpha}=0\end{array}\end{split}\end{aligned}\end{align} \]

\({D}_{1}\) est un paramètre matériau qui jour le rôle de pénalisation en rupture, comme \({P}_{\mathit{rupt}}\) pour la loi JOINT_MECA_RUPT. Il peut être démontré, qu’avec notre choix de la fonction d’adoucissement l’évolution de variable d’endommagement et de plasticité est simultanée (voir thèse Fontana, [bib20]). Ceci simplifie la description pour la version discrétisée de la loi :

(3431)\[\begin{split}\begin{array}{c}{\sigma}_{n}={K}_{n}({\delta}_{n}-{p}_{n}^{-}-\Delta {p}_{n})\\ \vec{{\sigma}_{t}}={K}_{t}(\vec{{\delta}_{t}}-{\vec{p}}_{t}^{-}-\Delta \vec{{p}_{t}})\\ \vec{{X}_{t}}={K}_{t}(\vec{{\delta}_{t}}-{\vec{p}}_{t}^{-}-\Delta \vec{{p}_{t}})-{A}_{t}({\alpha}^{-}+\Delta \alpha )\cdot ({\vec{p}}_{t}^{-}+\Delta \vec{{p}_{t}})\\ {X}_{n}={K}_{n}({\delta}_{n}-{p}_{n}^{-}-\Delta {p}_{n})-{A}_{n}({\alpha}^{-}+\Delta \alpha )\cdot ({p}_{n}^{-}+\Delta {p}_{n})\\ Y=-{B}_{n}\cdot S'({\alpha}^{-}+\Delta \alpha )\cdot {({p}_{n}^{-}+\Delta {p}_{n})}^{2}/2-{B}_{t}\cdot S'({\alpha}^{-}+\Delta \alpha )\cdot {\Vert {\vec{p}}_{t}^{-}+\Delta \vec{{p}_{t}}\Vert }^{2}/2\end{array}\end{split}\]

où nous avons introduit la fonction scalaire \(S(\alpha )={(1-\alpha )}^{3}/\sqrt{\alpha}\) . Les seuils gardent s’expriment facilement en version discrétisée implicite:

(3432)\[\begin{split}\begin{array}{c}{f}_{X}(\vec{X})=\Vert \vec{{X}_{t}}\Vert +\mu {X}_{n}-\overline{c}\le 0\\ \Delta {p}_{n}=\mu \Delta \lambda \\ \Delta \vec{{p}_{t}}=\Delta \lambda \frac{\vec{{X}_{t}}}{\Vert \vec{X}\Vert }\\ {f}_{X}(\vec{X})\Delta \lambda =0;\Delta \lambda \ge 0\\ Y=-{B}_{n}\cdot S'({\alpha}^{-}+\Delta \alpha )\cdot {({p}_{n}^{-}+\Delta {p}_{n})}^{2}/2-{B}_{t}\cdot S'({\alpha}^{-}+\Delta \alpha )\cdot {\Vert {\vec{p}}_{t}^{-}+\Delta \vec{{p}_{t}}\Vert }^{2}/2\le {D}_{1}\\ (Y-{D}_{1})\cdot \Delta \alpha =0\end{array}\end{split}\]

La résolution de ce système d’équation se fait par la procédure standard de la prédiction élastique. En phase plastique endommageante, nous obtenons un système d’équation implicite avec quatre inconnus suivants \(\Delta {p}_{n};\Delta \vec{{p}_{t}};\Delta \lambda ;\Delta \alpha ` . En utilisant le fait que :math:\)Delta vec{{p}_{t}}uparrow uparrow {K}_{t}(vec{{delta}_{t}}-{vec{p}}_{t}^{-})-{A}_{t}({alpha}^{+})cdot {vec{p}}_{t}^{-}` le système peut être simplifié est ramené à une équation pour la seule inconnue scalaire \({\alpha}^{+}={\alpha}^{-}+\Delta \alpha ` . Pour cela l’inconnu :math:\)Delta lambda ` est exprimé à partir des deux fonctions seuil. A partir du seuil plastique nous obtenons:

(3433)\[\Delta \lambda =\frac{\Vert {K}_{t}(\vec{{\delta}_{t}}-{\vec{p}}_{t}^{-})-{A}_{t}({\alpha}^{+})\cdot {\vec{p}}_{t}^{-}\Vert +\mu {K}_{n}({\delta}_{n}-{p}_{n}^{-})-\mu {A}_{n}({\alpha}^{+})\cdot {p}_{n}^{-}-\overline{c}}{{\mu}^{2}{K}_{n}+{K}_{t}+{\mu}^{2}{A}_{n}({\alpha}^{+})+{A}_{t}({\alpha}^{+})}\]

En parallèle à partir du seuil d’endommagement nous avons:

(3434)\[\begin{split}\begin{array}{c}F({\alpha}^{+})=-({\mu}^{2}{B}_{n}+{B}_{t})\Delta {\lambda}^{2}-2\left(\mu {B}_{n}+{B}_{t}{\vec{p}}_{t}^{-}\frac{{K}_{t}(\vec{{\delta}_{t}}-{\vec{p}}_{t}^{-})-{A}_{t}({\alpha}^{+})\cdot {\vec{p}}_{t}^{-}}{\Vert {K}_{t}(\vec{{\delta}_{t}}-{\vec{p}}_{t}^{-})-{A}_{t}({\alpha}^{+})\cdot {\vec{p}}_{t}^{-}\Vert }\right)\Delta \lambda -\\ -{B}_{n}{p}_{n}^{-2}-{B}_{t}{\vec{p}}_{t}^{-2}-\frac{2{D}_{1}}{S'({\alpha}^{+})}=0\end{array}\end{split}\]

En remplaçant \(\Delta \lambda ` de l’éq. dans l’éq. , nous obtenons une équation scalaire non-linéaire pour :math:`{\alpha}^{+}\) . Cette équation non-linéaire est résolue par la méthode de dessication alternée avec la méthode de corde.

La même jeux d’équation est résolution dans la pointe du cône en adoptant la règle d’écoulement selon les normales sortante:

(3435)\[G({\alpha}^{+})=-{B}_{n}{\left(\frac{{K}_{n}{\delta}_{n}-\overline{c}/\mu }{{K}_{n}+{A}_{n}({\alpha}^{+})}\right)}^{2}-{B}_{t}{\left(\frac{{K}_{t}\vec{{\delta}_{t}}}{{K}_{t}+{A}_{t}({\alpha}^{+})}\right)}^{2}-\frac{2{D}_{1}}{S'({\alpha}^{+})}=0\]

L’algorithme de résolution se résume schématiquement ainsi:

Variables internes

La loi JOINT_MECA_ENDO possède dix-huit variables internes. Du point de vue de la loi de comportement, seules la première, la troisième et la quatrième sont stricto sensu des variables internes. Les autres fournissent des indications sur l’état de hydromécanique du joint à un instant donné.

\(V1\)	Paramètre croissant indiquant la plasticitécumulée au cours de glissement \(\lambda\) .
\(V2\)	Indicateur de glissement \(=0\) si le régime estlinéaire élastique, \(=1\) si le régime est plastique en cisaillement sous compression, \(=2\) si le régime est plastique en traction.
\(V3\)	Plastification normale à l’orientation du joint \({p}_{n}\) [1]_
\(V4\) \(V5\)	Vecteur de la plastification tangentielle \(\vec{{p}_{t}}\) . \(V5\) est mis à zéro en 2D.
\(V6\)	Variable d’endommagement \(\alpha :0\le \alpha \le 1\) .
\(V7\)	Saut normal \(\mathit{V7}={\delta}_{n}\) dans le repère local.
\(V8\) \(V9\)	Saut tangentiel \({\delta}_{t1}\) et \({\delta}_{t2}\) dans le repère local (\(V9\) nul en 2D).
\(V10\)	Variable non utilisée (réservé pour les modifications futures de la loi).
\(V11\)	Contrainte mécanique normale \({\sigma}_{n}\) (sans pression de fluide).
\(V12\) \(V13\) \(V14\)	Composantes du gradient de pression dans le repère global (uniquement pour les modélisationsxxx_JOINT_HYME). Respectivement \({\partial}_{x}p\) , \({\partial}_{y}p\) et \({\partial}_{z}p\) .
\(V15\) \(V16\) \(V17\)	Composantes du flux hydraulique dans le repère global (uniquement pour les modélisations xxx_JOINT_HYME). Respectivement \({w}_{x}\) , \({w}_{y}\) et \({w}_{z}\) .
\(V18\)	Pression de fluide \(p\) imposée par l’utilisateur (PRES_FLUIDE) dans le cas des modélisations xxx_JOINTou pression de fluide interpolée à partir de celle calculée aux nœuds milieux des éléments de joint (degré de liberté du problème) pour lesmodélisations : xxx_JOINT_HYME.
\(V19\) \(V20\)	Contraintes mécanique tangentielles \({\sigma}_{t1}\) et \({\sigma}_{t2}\) (\(V20\) nul en 2D).

[1]

Attention: nous utilisons \(({p}_{n},\vec{{p}_{t}})`avec indices pour la plasticité et sans indice pour la pression :math:`p\)

Prise en compte de la pression hydrostatique sans couplage

Fig. 241 Illustration d’un calcul possible de stabilité d’un barrage avec le profil de pression imposé

Bien que la modélisation XXX_JOINT ne couple pas la mécanique et l’hydraulique, on peut toutefois introduire explicitement l’influence d’un fluide sur la mécanique via une pression imposée. La présence du fluide dans le joint modifie la contrainte mécanique normale \({\sigma}_{n}\to {\sigma}_{n}-p\) . En mettant une pression importante on est capable de faire rompre le joint par un simple effet hydraulique. Pour prendre en compte les effets hydrostatiques la loi mécanique est décalée vers le bas (Figure ) en fonction de la valeur de pression \(p\) en chaque point d’intégration.

L’implémentation numérique est aisée en cas d’écriture complète des lois mécaniques sous forme explicite non incrémentale en fonction des déplacements et des variables internes (il faut exclure la dépendance des contraintes à l’instant plus en fonction de contraintes à l’instant précédent). Dans ce cas la seule modification de courbe normale est suffisante pour introduire le couplage:

(3436)\[{\sigma}_{n}={\sigma}_{n}^{\mathit{meca}}({\delta}_{n},{\delta}_{t})-p\]

En se limitant à ce type de phénomène physique, il est possible de faire des études où le profil de pression est imposé par utilisateur, par exemple une étude de stabilité de barrage sous hypothèse conservative (Figure ), c’est-à-dire en présence de sous-pression, dont la forme est très pénalisante. Afin de faire un calcul avec une pression imposée l’utilisateur doit définir une fonction, par le mot-clef PRES_FLUIDE, qui dépend à la fois de l’espace (profil de pression non-homogène) et du temps (évolution du profil de pression).

En se limitant à ce type de phénomène physique, il est possible de faire des études où le profil de pression est imposé par utilisateur, par exemple une étude de stabilité de barrage sous hypothèse conservative (Figure ), c’est-à-dire en présence de sous-pression, dont la forme est très pénalisante. Afin de faire un calcul avec une pression imposée l’utilisateur doit définir une fonction, par le mot-clef PRES_FLUIDE, qui dépend à la fois de l’espace (profil de pression non-homogène) et du temps (évolution du profil de pression).

Formulation théorique du couplage hydromécanique

Les lois introduites peuvent s’appuyer sur une modélisation hydromécanique couplée, par les éléments XXX_JOINT_HYME. Dans cette partie on parlera de la partie hydraulique de la loi, ainsi que du couplage lui même; tous les détails sur la partie mécanique de la loi ont été décrits précédemment.

Modélisation hydraulique

Le fluide s’écoule des zones de haute pression vers celles de basse pression. Une manière théorique de prendre en compte l’écoulement stationnaire est d’associer à l’état hydraulique donné une énergie [1]_ \(H(p(x))\) dépendant de la distribution de pression. La première hypothèse consiste à supposer que l’énergie dépend explicitement de la variation de pression et pas de la pression elle-même \(H=H(\nabla p(x))\) . En prenant la forme convexe la plus simple possible de cette dépendance en gradient, on obtient ainsi l’énergie \(H=C{(\overrightarrow{\nabla}p)}^{2}/2\) où \(C\) est un paramètre de la loi, qui ne dépend pas de la pression.

En calculant les efforts généralisés correspondant au champ de gradient de pression on obtient la première loi de Fick. Le flux hydraulique est proportionnel au gradient de pression:

(3437)\[\overrightarrow{w}=\frac{\partial H}{\partial \overrightarrow{\nabla}p}=C\overrightarrow{\nabla}p\]

Dans ce formalisme énergétique on cherche le champ de pression à l’équilibre par minimisation de l’énergie hydraulique \(\underset{p(\overrightarrow{x})}{\min}\underset{\Omega}{\int}H(\overrightarrow{\nabla}p(\overrightarrow{x}))d\Omega\) . Ce qui donne une équation d’équilibre ressemblant à celle de la mécanique \(div\overrightarrow{w}=0\) . Dans le cadre de ce modèle la résolution d’équation d’équilibre hydraulique est équivalente à une résolution de problème mécanique en quasi-statique, où le flux hydraulique est équivalente aux contrainte \(\overrightarrow{w}\iff \sigma\) , le champ de pression correspond au champ de déplacement \(p(\vec{x})\iff u(\vec{x})\) et enfin le gradient de pression s’apparente au champ de déformation \(\overrightarrow{\nabla}p\mathrm{\iff }\varepsilon\) .

Influence de l’hydraulique sur la mécanique

La présence du fluide dans le joint rajoute une contrainte hydrostatique et ce fait modifie la contrainte mécanique normale \({\sigma}_{n}\to {\sigma}_{n}-p\) . En mettant une pression importante on est capable de faire rompre le joint par un simple effet hydraulique. On peut décaler la loi mécanique vers le bas en fonction de la valeur de pression \(p\) en chaque point pour prendre en compte les effets de la pression, voir les figures r7.01.25-fig-joint-rup et r7.01.25-fig-joint-frot.

Influence de la mécanique sur l’hydraulique

Dans le cas d’écoulement de fluide à travers une fissure le flux hydraulique doit augmenter avec l’ouverture (\({\delta}_{n}\) ) de cette dernière (\(\overrightarrow{w}~O({\delta}_{n})\overrightarrow{\nabla}p\) ). Dans la loi de Poiseuille, qui a été trouvée empiriquement pour l’écoulement laminaire d’un fluide visqueux et incompressible, la dépendance de flux en ouverture est cubique (la loi est souvent appelée la loi cubique). La partie hydraulique de la loi utilise ce type de couplage. Les équations à résoudre s’écrivent de la manière suivante:

(3438)\[div\overrightarrow{w}=0;\overrightarrow{w}=\frac{\rho}{12\stackrel{ˉ}{\mu}}{\delta}_{n}^{3}\overrightarrow{\nabla}p\]

Dans le cas d’un écoulement de fluide à travers des jonctions d’un barrage, on constate des flux non-négligeables même pour les joints fermés. Il est nécessaire alors de définir une épaisseur minimale \({\epsilon}_{\min}\) , mot clef OUV_MIN, en dessous de laquelle le flux atteint sa valeur minimale. Nous régularisons les équations d’écoulement de la manière suivante:

(3439)\[\vec{w}=\frac{\rho}{12\stackrel{̄}{\mu}}\max{({ϵ}_{min},{ϵ}_{min}+{\delta}_{n})}^{3}\vec{\nabla}p\]

Pour un gradient de pression non-nul le flux n’atteint jamais la valeur nulle, \(\min\overrightarrow{w}~{\epsilon}_{\min}^{3}\overrightarrow{\nabla}p\) , ce qui correspond à l’écoulement à travers des parois perméables du joint fermé.

Couplage hydromécanique

Le couplage hydromécanique fait intervenir les deux mécanismes décrits précédemment: d’un coté le fluide agit par pression sur les lèvres de joint, de l’autre coté plus la fissure est ouverte plus l’écoulement de fluide est aisé. En absence des forces extérieures le calcul hydromécanique se présente schématiquement sous cette forme:

(3440)\[\begin{split}\left\{ \begin{array}{cc}\vec{w}=\vec{w}(\vec{\delta}(u),\vec{\nabla}p);& div\vec{w}=0\\ \vec{\sigma}=\vec{\sigma}(\vec{\delta}(u),p);& div\vec{\sigma}=0\end{array}\text{}\Rightarrow \text{}\vec{Y}=\vec{Y}(\vec{X});div\vec{Y}=0\right.\end{split}\]

La résolution des équations d’équilibre hydromécanique est équivalente à la résolution du problème mécanique en quasi-statique, où on introduit les contraintes généralisées \(\overrightarrow{Y}=(\overrightarrow{w},\overrightarrow{\sigma})\) , et le champ vectoriel des inconnues \(\overrightarrow{X}=(p,u)\) .

Matrice tangente

Vu que les efforts généralisés ne dépendent de \(u\) qu’à travers de \(\overrightarrow{\delta}(u)\) , pour calculer la matrice tangente de couplage hydromécanique, il est nécessaire de ne connaître que les quatre termes suivants:

(3441)\[ \begin{align}\begin{aligned}\frac{\partial \overrightarrow{\sigma}}{\partial \overrightarrow{\delta}}\\\textrm{,}\\\frac{\partial \overrightarrow{\sigma}}{\partial p}` , :math:`\frac{\partial \overrightarrow{w}}{\partial \overrightarrow{\nabla}p}` et :math:`\frac{\partial \overrightarrow{w}}{\partial \overrightarrow{\delta}}\end{aligned}\end{align} \]

Le premier terme est le même qu’en mécanique pure, il est donné dans l’équation (4849). Le deuxième terme est trivial, car la seule composante non-nulle est égale à \(\partial {\sigma}_{n}/\partial p=-1\) . Le terme hydraulique diagonal prend une forme simple car le flux hydraulique ne dépend que du gradient de pression:

(3442)\[\frac{\partial \vec{w}}{\partial \vec{\nabla}p}=\frac{\rho}{12\stackrel{̄}{\mu}}\max{({ϵ}_{min},{ϵ}_{min}+{\delta}_{n})}^{3}\]

Dans le dernier terme seule la dérivée par rapport à l’ouverture normale n’est pas nulle:

(3443)\[\frac{\partial \overrightarrow{w}}{\partial {\delta}_{n}}=\frac{\rho}{4\stackrel{ˉ}{\mu}}{({\epsilon}_{\min}+{\delta}_{n})}^{2}\overrightarrow{\nabla}p\]

Cette quantité est égale à zéro pour une fissure fermée \({\delta}_{n}<0\) . La matrice tangente ainsi formulée n’est pas symétrique.

[1]

Nous utilisons une notation simplifiée, le terme exact serait : taux de densité de l’énergie.

Fonctionnalités et validation

Deux lois de comportement JOINT_MECA_RUPT et JOINT_MECA_FROT sont introduites. Elles sont validées sur les cas tests élémentaires ssnp162 et le pseudo barrage-poids ssnp142 . La procédure du clavage est validée sur la simulation d’injection du coulis entre deux blocs rectangulaires encastrés au sol ssnp143 . La procédures du sciage est validées par le sciage de deux blocs rectangulaires avec différents types d’encastrement ssnp143 .

Validation en mécanique pure, modélisations type XXX_JOINT
Loi: JOINT_MECA_RUPT	Loi: JOINT_MECA_FROT
Tests: ssnp162a/b/c; ssnp142a/b; Clavage:ssnp143a/b Sciage: ssnp143c/d/e/f	Tests: ssnp162d/e/f; ssnp142c/d Sciage: ssnp143c/d/g/h

Validation hydromécanique couplée, modélisations type XXX_JOINT_HYME
Loi: JOINT_MECA_RUPT	Loi: JOINT_MECA_FROT
Tests: ssnp162g/h/i; ssnp142e/f	Tests: ssnp162j/k/l; ssnp142g/h

Bibliographie

[bib1]

BARENBLATT G. I.; The mathematical theory of equilibrium cracks in brittle fracture,Adv. Appl. Mech ., 55-129, 1962

[bib2]

K.Kazymyrenko, J. Laverne; Avancement de la modélisation mécanique des joints-plots, CR-AMA.09.039, 2009

[bib3]

C.Kazymyrenko, J.Laverne; État des lieux et perspectives pour la modélisation des joints dans les barrages, CR-AMA.10.357, 2010

[bib4] (1,2)

16. Divoux, Modélisation du comportement hydromécanique des discontinuités dans les structures et les fondations rocheuses : application aux barrages béton , 1997

[bib5]

16. Kolmayer; Modélisation du comportement hydromécanique des joints dans les barrages en béton, IH/CODHY/GCED/00001/AO,

[bib6]

LAVERNE J., Formulation énergétique de la rupture par des modèles de forces cohésives: considérations théoriques et implantations numériques, 2004

[bib7]

Nguen Quoc Sun; ,International Journal for numerical methods in engineering, 817-832, 1977

[bib8]

Intégration des relations de comportement élasto-plastique de Von Mises, Manuel de Référence du Code_Aster,

[bib9]

15. Mohr. Welche umst ̈ande bedingen die elastizit ̈atsgrenze und den bruch eines materials? Zeitschrift des Vereines deutscher Ingenieure, 44:1524–1530, 1900.

[bib10]

3. 1. Coulomb. Essai sur une application des règles des maximis et minimis à quelquels problèmes de statique, relatifs à l’architecture. Académie Royale Des Sciences Paris Mem. Math. Phys., 7:343–382, 1776.

[bib11]

4. 3. Drucker and W. Prager. Soil mechanics and plastic analysis or limit design. Quarterly of Applied Mathematics, 10:157–165, 1952.

[bib12] (1,2,3)

Jean-Jacques Marigo and Kyrylo Kazymyrenko. A micromechanical inspired model for the coupled to damage elastoplastic behavior of geomaterials under compression. Mechanics & Industry, 20(105), 2019

[bib13]

Evert Hoek and Edwin T. Brown. Empirical strength criterion for rock masses. Journal of the Geotechnical Engineering Division, 106:1013–1035, 1980.

[bib14]

5. Hoek and E. T. Brown. Underground Excavation in Rock. E & FN Spon, 1982.

[bib15] (1,2)

16. Germain, Q. S. Nguyen, and P. Suquet. Continuum thermodynamics. Journal of Applied Mechanics, 50(4b):1010–1020, 1983. doi: 10.1115/1.3167184.

[bib16]

1. 6. Collins and G. T. Houlsby. Application of thermomechanical principles to the modelling of geotechnical materials. Proceedings of the Royal Society A: Mathematical, Physical, and Engineering Sciences, 453(1964):1975–2001, 1997. doi: 10.1098/rspa.1997.0107.

[bib17]

7. 20. Houlsby and A. M. Puzrin. Principles of Hyperplasticity: An Approach to Plasticity Theory Based on Thermodynamic Principles. Springer Verlag London Limited, 2007. doi: 10.1007/978-1-84628-240-9.

[bib18] (1,2,3,4)

Bernard Halphen and Quoc Son Nguyen. Sur les matériaux standard généralisés. Journal de Mécanique, 14(1):39–63, January 1975.

[bib19]

16. 1. Veermer and R. de Borst. Non-associated plasticity for soils, concrete and rock. HERON, 29(3):1–64, 1984.

[bib20] (1,2)

«Problèmes d’interfaces pour les ouvrages hydrauliques» thèse de doctorat Ilaria Fontana, Université de Montpellier et EDF R&D, 2022