r7.02.15 Modélisation des fissures avec couplage hydro-mécanique en milieu poreux saturé

Résumé:

On propose un nouvel élément de joint permettant de modéliser en 2D les discontinuités avec couplage-hydromécanique. Cet élément est utilisé pour mailler une fissure ou un trajet de fissure dans un maillage d’éléments THM volumiques classiques dans le cas d’un matériau poreux saturé. Il permet de modéliser le comportement d’un joint hydraulique ou d’une fissure sous pression de fluide.

On présente ici les équations constitutives de la mécanique et de l’hydraulique dans la fissure ainsi que les lois de comportement correspondantes. On détaille également la discrétisation adoptée et le vecteur force interne et l’opérateur tangent introduits dans l’opérateur STAT_NON_LINE.

Table des matières

1 Description_des_versions_du_document

Présentation du problème

Hypothèses

On considère un milieu poreux, noté \(\Omega\) à l’instant actuel et dont la frontière est notée \(\partial \Omega\) . Ce volume est séparé en deux parties \({\Omega}^{+}\) et \({\Omega}^{-}\) par une interface \(\Gamma\) . Toutes les grandeurs dans \({\Omega}^{+}\) (resp. \({\Omega}^{-}\) ) sont indicées d’un + (reps. d’un -).

On se limite à l’étude d’un milieu poreux saturé en deux dimensions.

Notations

Grandeurs mécaniques

Le champ de déplacement dans le massif et la discontinuité de déplacement à travers la fissure sont notés respectivement \(u=(\begin{array}{c}{u}_{x}\\ {u}_{y}\end{array})\) et \(〚u〛=(\begin{array}{c}〚{u}_{x}〛\\ 〚{u}_{y}〛\end{array})\) .

La contrainte totale et la contrainte effective sont notées respectivement \(\sigma\) et \(\sigma '\) dans le massif et \(\tau\) et \(\tau '\) sur l’interface.

Grandeurs hydrauliques

La pression de pore du fluide intersticiel est noté \({p}^{+}\) et \({p}^{-}\) dans les massifs environnant la fissure. La pression de fluide dans la fissure est notée \(p\) . Le gradient de pression dans les massifs est noté \(\nabla {p}^{+}\) ou \(\nabla {p}^{-}\) . Le gradient longitudinal de la pression de la fissure est noté \({\nabla}_{l}p\) et est défini par

(4495)\[{\nabla}_{l}p=\nabla p-\frac{\partial p}{\partial n}n\]

où \(n\) est la normale à l’interface.

Les apports massiques de fluide sont notés \(m\) dans le massif (unité : \({\mathit{M.L}}^{-3}\) ) . Dans la fissure, ils sont notés \(w\) (unité : \({\mathit{M.L}}^{-2}\) ) et sont intégrés sur l’épaisseur \(\varepsilon\) de la fissure :

Le flux hydraulique massique dans le massif est noté \(M\) (unité : \({\mathit{M.T}}^{-1}.{L}^{-2}\) ) . Dans la fissure, le débit hydraulique massique est noté \(W\) (unité : \({\mathit{M.T}}^{-1}.{L}^{-1}\) ).

\(\varepsilon\) est l’ouverture de la fissure et est donc relié au saut de déplacement normal par

(4495)\[\varepsilon ={\varepsilon}_{0}+〚u〛.n\]

Où \({\varepsilon}_{0}\) est l’épaisseur initiale.

Équations continues

Mécanique

Équations constitutives

La formulation des modèles à zone cohésive et leur intégration numérique sont présentées dans la documentation de référence [R7.02.11]. On en fait ici un très bref rappel.

On cherche le champ de déplacement \(u\) à l’équilibre par minimisation de l’énergie

(4513)\[E(u)={\int}_{\Omega /\Gamma }\phi (\varepsilon (u))d\Omega -{W}^{\text{ext}}(u)+{\int}_{\Gamma}\Psi (〚u〛)d\Gamma\]

où \(\phi\) est la densité d’énergie mécanique volumique, \({W}^{\text{ext}}\) est le travail des forces extérieurs et \(\Psi\) est la densité d’énergie de surface.

La donnée de \(\Psi\) caractérise le comportement mécanique de l’interface. Elle est exprimée par une loi cohésive qui relie la force de cohésion \(\tau\) qui s’exerce sur les lèvres de la fissure au saut de déplacement par :

(4513)\[\tau =\frac{\partial \Psi }{\partial 〚u〛}\]

Par ailleurs, on fait l’hypothèse des contraintes effectives. Dans le massif, le tenseur contrainte totale est décomposé en :

(4513)\[\sigma =\sigma '+{\sigma}_{p}I\]

où \(\sigma '\) désigne la contrainte effective et \({\sigma}_{p}\) la contrainte hydraulique.

Sur la discontinuité, le vecteur contrainte totale est décomposé en :

(4513)\[\tau =\tau '+{\tau}_{p}n\]

où \(\tau '\) désigne la contrainte effective.

Lois de comportement

Contraintes effectives

Dans le massif, on fait l’hypothèse des contraintes effectives de Biot pour les milieux saturés, en introduisant la contrainte hydraulique \({\sigma}_{p}\) fonction de la pression \(p\) (en l’occurence \({p}^{+}\) ou \({p}^{-}\) dans le massif) telle que :

(4513)\[{\sigma}_{p}=-bp\]

où \(b\) est le coefficient de Biot du matériau.

Sur la discontinuité, dans le cadre d’un modèle de zone cohésive, on distingue trois zones :

• une zone ouverte où la contrainte totale sur les lèvres est égale à \(pn\) ;

• une zone de transition entre le milieu ouvert et le milieu sain. Dans cette zone, on observe l’apparition de micro-fissuration et de plasticité. On fait ainsi l’hypothèse d’incompressibilité plastique de la matrice et c’est la contrainte effective de Terzaghi qui contrôle l’ouverture de la zone. La contrainte totale s’écrit alors

(4513)\[\tau =\tau '-pn\]

• une zone saine où il y a adhérence et non interpénétration des lèvres. Tant que la valeur de la contrainte reste inférieure à la contrainte critique, l’ouverture reste nulle.

On fait donc l’hypothèse, en cohérence avec le comportement de ces trois zones, que

(4513)\[{\tau}_{p}=-p\]

Lois cohésives

Le modèle présenté ici est compatible avec les lois cohésives dont l’énergie est régularisée en zéro :

• CZM_LIN_REG,

• CZM_EXP_REG.

Ces lois permettent de prendre en compte :

• La non interpénétration des lèvres de fissure par pénalisation ;

• La propagation de fissure par une loi adoucissante ;

• L’irréversibilité de la fissuration.

Les paramètres matériau de l’interface sont la contrainte critique à la rupture \({\sigma}_{c}\) et l’énergie de rupture \({G}_{c}\) . Les paramètres numériques sont le paramètre de régularisation de l’énergie PENA_ADHERENCE et le paramètre de pénalisation de l’interpénétration PENA_CONTACT.

Loi de Bandis

Dans ce cas, le modèle de joint ne prend pas en compte la propagation de fissure. Le joint sépare deux massifs en étant ouvert tout du long et est simplement doté d’une loi qui relie l’ouverture à la contrainte effective.

Dans le cas de la loi empirique de Bandis, la relation entre la contrainte effective normale \(\tau {'}_{n}\) et la fermeture de fissure normale \(U={U}_{\max}-\varepsilon\) (où \(U\) désigne \(〚u〛.n\) ) est donnée par

(4513)\[d\tau {'}_{n}=-{K}_{ni}\frac{\mathrm{dU}}{{(1-\frac{U}{{U}_{\max}})}^{\gamma}}\]

où \({U}_{\max}\) est la fermeture asymptotique de la fissure, \({K}_{ni}\) la rigidité initiale normale de la discontinuité et \(\gamma\) un coefficient empirique qui varie entre 2 et 6 et qui dépend de la rugosité des parois du joint. La loi de comportement () est représentée en figure .

Illustration 3.1.1: Loi de Goodman : courbe contrainte-fermeture

Dans la direction tangentielle, le comportement est supposé élastique.

(4513)\[\tau {'}_{t}={\tau}_{t}={K}_{t}〚{u}_{t}〛\]

Le mot clé de Code_Aster correspondant à cette loi est JOINT_BANDIS. Ce mot clé est à renseigner dans DEFI_MATERIAU ainsi que dans le comportement (dans RELATION_KIT).

Hydraulique

Équations constitutives

Illustration 3.2.1: conservation de la masse dans la fissure (à gauche) et profil typique de pression à travers l’interface (à droite)

L’équation de conservation de la masse dans la fissure s’écrit

(4513)\[\frac{\partial w}{\partial t}+\frac{\partial W}{\partial s}=0\]

L’épaisseur de fissure étant faible, on considère que le champ de pression est continu à travers l’interface. En tout point de l’interface, on a donc

(4513)\[{p}_{\Gamma}^{+}=p={p}_{\Gamma}^{-}\]

En revanche, les discontinuités de flux sont autorisées à travers l’interface (voir figure droite)

Lois de comportement

Les apports massiques de fluide dans la fissure ouverte s’écrivent

(4513)\[w=\rho \varepsilon\]

Où \(\rho\) est la masse volumique du fluide.

L’évolution des apports massiques est donnée par :

(4513)\[\frac{\mathrm{dw}}{\rho}=(\frac{1}{N}+\frac{\varepsilon}{{K}_{f}})\mathrm{dp}+d\varepsilon\]

où \({K}_{f}\) est le module de compressibilité du fluide et \(N\) le module de Biot de la zone cohésive.

On considère que l’écoulement dans la fissure est darcéen. Le flux massique peut donc s’écrire

(4513)\[W=-\varepsilon \rho {\lambda}_{H}(\frac{\partial p}{\partial s}-\rho g\frac{\partial z}{\partial s})\]

où \({\lambda}_{H}\) est la conductivité hydraulique de la fissure et z désigne la coordonnée suivant l’axe vertical. Elle est donnée par

(4513)\[{\lambda}_{H}=\frac{K(\varepsilon )}{\mu (T)}\]

La relation entre le flux et le gradient de pression est donné par la loi cubique. Le flux volumique à travers la fissure s’écrit alors

(4513)\[\frac{W}{\rho}=-\frac{{\varepsilon}^{3}}{12\mu }(\frac{\partial p}{\partial s}-\rho g\frac{\partial z}{\partial s})\]

et la perméabilité intrinsèque \(K\) en fonction de l’ouverture

(4513)\[K(\varepsilon )=\frac{{\varepsilon}^{2}}{12}\]

On note \(T\) la transmissivité de la fissure, définie par

(4513)\[T=\varepsilon {\lambda}_{H}\rho =\frac{\rho {\varepsilon}^{3}}{12\mu }\]

Cette loi est retrouvée analytiquement par la loi de Poiseuille lorsqu’on étudie un écoulement laminaire entre deux plaques lisses séparées par une distance \(\varepsilon\) petite devant les autres dimensions. Sa validité a également été mise en évidence sur des fissures dans les roches pour une large gamme de paramètres.

Lorsque les parois de fissure sont imperméable, on peut modéliser le massif avec des éléments finis mécaniques classiques. Dans ce cas, l’équation d’écoulement devient singulière dans la partie de fissure fermée et non endommagée. On introduit donc un paramètre de régularisation OUV_FICT qui permet d’avoir un écoulement fictif dans la partie de fissure fermée. Sur les éléments fermés, on prend \(\varepsilon\) égal à OUV_FICT.

Formulation variationnelle

Mécanique

On note \({U}_{\mathrm{ad}}\) l’ensemble des champs de déplacements admissibles, c’est-à-dire les éléments de \({({H}^{1}(\Omega ))}^{d}\) vérifiant les conditions aux limites en déplacement sur la partie de \(\partial \Omega\) supportant de telles conditions.

Les conditions d’optimalité pour l’énergie () donnent la formulation variationnelle suivante :

(4519)\[\begin{split}\lbrace \begin{array}{c}\sigma =\sigma '+{\sigma}_{p}I\\ \tau =\tau '-pn\\ {\int}_{\Omega \setminus \Gamma }\sigma :\epsilon (\stackrel{ˆ}{u})d\Omega +{\int}_{\Gamma}\tau .〚\stackrel{ˆ}{u}〛d\Gamma ={W}_{\mathrm{ext}}(\stackrel{ˆ}{u})\forall \stackrel{ˆ}{u}\in {U}_{\mathrm{ad}}\end{array}\end{split}\]

Hydraulique

On note \({P}_{\mathrm{ad}}^{+}\) (resp. \({P}_{\mathrm{ad}}^{-}\) ) l’ensemble des champs de pression admissibles sur \({\Omega}^{+}\) (resp. \({\Omega}^{-}\) ) , c’est-à-dire les éléments de \({H}^{1}({\Omega}^{+})\) (resp. \({H}^{1}({\Omega}^{-})\) ) vérifiant les conditions aux limites en pression sur \(\partial {\Omega}_{P}^{+}\) , la partie de \(\partial {\Omega}^{+}\) supportant des conditions aux limites en pression, (resp. \(\partial {\Omega}_{P}^{-}\) ). Et on note \({P}_{\mathrm{ad}}\) l’ensemble des champs de pression admissibles sur \(\Gamma\) , c’est-à-dire les éléments de \({H}^{1}(\Gamma )\) vérifiant les conditions aux limites en pression sur \(\partial {\Gamma}_{P}\) .

(4519)\[-{\int}_{{\Omega}^{+}}\frac{\partial m}{\partial t}{\stackrel{ˆ}{p}}^{+}d\Omega +{\int}_{{\Omega}^{+}}M.{\nabla \stackrel{ˆ}{p}}^{+}d\Omega ={\int}_{\partial {\Omega}_{F}^{+}}F.n{\stackrel{ˆ}{p}}^{+}d\Gamma +{\int}_{\Gamma}{q}^{+}{\stackrel{ˆ}{p}}^{+}d\Gamma \forall {\stackrel{ˆ}{p}}^{+}\in {P}_{\mathrm{ad}}^{+}\]

(4519)\[-{\int}_{{\Omega}^{-}}\frac{\partial m}{\partial t}{\stackrel{ˆ}{p}}^{-}d\Omega +{\int}_{{\Omega}^{-}}M.{\nabla \stackrel{ˆ}{p}}^{-}d\Omega ={\int}_{\partial {\Omega}_{F}^{-}}F.n{\stackrel{ˆ}{p}}^{-}d\Gamma +{\int}_{\Gamma}{q}^{-}{\stackrel{ˆ}{p}}^{-}d\Gamma \forall {\stackrel{ˆ}{p}}^{-}\in {P}_{\mathrm{ad}}^{-}\]

(4519)\[-{\int}_{\Gamma}\frac{\partial w}{\partial t}\stackrel{ˆ}{p}d\Gamma +{\int}_{\Gamma}W{\nabla}_{l}\stackrel{ˆ}{p}d\Gamma +{\int}_{\Gamma}({q}^{+}+{q}^{-})\stackrel{ˆ}{p}d\Gamma ={\int}_{\partial {\Gamma}_{F}}F.n\stackrel{ˆ}{p}\mathrm{ds}\forall \stackrel{ˆ}{p}\in {P}_{\mathrm{ad}}\]

(4519)\[{\int}_{\Gamma}({p}^{+}-p){\stackrel{ˆ}{q}}^{+}d\Gamma =0\forall {\stackrel{ˆ}{q}}^{+}\in {H}^{-1}(\Gamma )\]

(4519)\[{\int}_{\Gamma}({p}^{-}-p){\stackrel{ˆ}{q}}^{-}d\Gamma =0\forall {\stackrel{ˆ}{q}}^{-}\in {H}^{-1}(\Gamma )\]

Discrétisation temporelle

On adopte une discrétisation en temps implicite. Les notations indicées par \(n\) sont les quantités en début de pas de temps et celles indicées par \(n+1\) sont les quantités en fin de pas de temps.

Le pas de temps est noté \(\Delta t={t}_{n+1}-{t}_{n}\) .

Par la suite, en l’absence de précision, les notations non indicées désigneront les grandeurs en fin de pas de temps.

Construction d’un élément de joint

La trajectoire de fissure \(\Gamma\) est définie a priori, on peut donc la discrétiser par des éléments de joint. Les sous-domaines \({\Omega}^{+}\) et \({\Omega}^{-}\) , situés de part et d’autre de \(\Gamma\) , sont discrétisés par des éléments THM volumiques classiques de telle sorte que leurs nœuds à l’interface coïncident.

Degrés d’interpolation des degrés de liberté

L’élément de joint avec couplage hydro-mécanique reprend la formulation mécanique des éléments de joint classiques (mais en passant d’un élément linéaire à un élément quadratique) et interagit avec les éléments THM voisins. Les degrés d’interpolation de des degrés de liberté est donc choisit en cohérence avec ces éléments voisins.

Afin d’être compatibles avec les éléments THM du massif :

• les déplacements sont interpolés quadratiquement (\(\mathit{P2}\) -continus)

• les pressions sont interpolées linéairement (\(\mathit{P1}\) -continues).

Un degré de liberté qui ne préexiste pas dans aucun des deux précédents éléments apparaît dans la formulation. Il s’agit du multiplicateur de Lagrange hydraulique (variables \({q}^{+}\) et \({q}^{-}\) dans les équations variationnelles à ). Le degré d’interpolation de ces multiplicateurs de Lagrange hydrauliques doit être choisit de manière à respecter une condition LBB discrète. Ils sont donc pris constants par élément.

Description de l’élément

L’élément construit est d’épaisseur nulle et on maille l’ensemble du trajet de fissure avec cet élément. Les bords inférieurs et supérieurs sont reliés au reste de la structure. L’écoulement du fluide le long de la fissure est écrit le long du plan milieu de l’élément (en gris sur la figure ).

Illustration 5.2.1: Vue « éclatée » de l’élément de joint avec couplage hydro-mécanique et sans propagation

Algorithme de résolution

Le problème variationnel (-) peut s’écrire sous la forme :

(4535)\[F(U)={L}^{i}(U)-{L}^{\mathrm{meca}}=0\]

où \(U\) désigne les déplacements généralisés.

L’opérateur tangent associé est noté \(DF=\frac{\partial F}{\partial U}\) .

Contraintes et déformations généralisées

On note \({U}^{\mathrm{el}}\) le vecteur des inconnues nodales sur l’élément \(\mathrm{el}\) , \({E}_{pi}^{\mathrm{el}}\) le vecteur des déformations généralisées au point d’intégration \(pi\) de l’élément \(\mathrm{el}\) , \({\Sigma}_{pi}^{\mathrm{el}}\) le vecteur des contraintes généralisées

(4535)\[\begin{split}{\Sigma}_{pi}^{\mathrm{el}}={\left\lbrace \begin{array}{c}\tau '\\ -p\\ w\\ W\\ {q}^{+}\\ {q}^{-}\\ {p}^{+}-p\\ {p}^{-}-p\end{array}\right\rbrace }^{\mathrm{el}}\end{split}\]

(4535)\[\begin{split}{E}_{pi}^{\mathrm{el}}={\left\lbrace \begin{array}{c}〚u〛\\ p\\ \nabla p\\ {p}^{+}\\ {p}^{-}\\ {q}^{+}\\ {q}^{-}\end{array}\right\rbrace }^{\mathrm{el}}\end{split}\]

Les déformations généralisées sont obtenus par la relation

(4535)\[{E}_{pi}^{\mathrm{el}}={Q}_{pi}^{\mathrm{el}}.{U}^{\mathrm{el}}\]

où \({Q}_{pi}^{\mathrm{el}}\) est la matrice de passage des degrés de libertés nodaux aux déformations généralisées aux point d’intégration \(pi\) . On la calcule à partir des fonctions de forme et de l’orientation de l’élément. En effet, les déplacements mécaniques de chaque nœud sont définis dans le repère global. Afin d’en extraire les composantes normale et tangentielles à l’élément de joint, on applique une matrice de rotation \({\Theta}_{pi}^{\mathit{el}}\) aux vecteurs déplacements nodaux.

Intégration

Pour intégrer sur l’élément les termes du vecteur, on choisit une méthode d’intégration sélective. Celle-ci permet d’éviter les oscillations numériques pour les problèmes où la structure est soumise à un choc et où les phénomènes mécaniques prédominent . Cette méthode consiste à intégrer aux sommets de l’élément les termes faisant intervenir une dérivée par rapport au temps et à intégrer aux points de Gauss les termes permanents.

Le vecteur force interne s’écrit

(4535)\[{L}^{i}(U)={\sum}_{\mathrm{el}}({\sum}_{g}{\omega}_{g}{R}_{g}^{\mathrm{el}}({U}^{\mathrm{el}})+{\sum}_{s}{\omega}_{s}{R}_{s}^{\mathrm{el}}({U}^{\mathrm{el}}))\]

en notant \({R}_{g}^{\mathrm{el}}\) et \({R}_{s}^{\mathrm{el}}\) les valeurs respectivement au point de Gauss \(g\) et au sommet \(s\) des forces nodales et \({\omega}_{g}\) et \({\omega}_{s}\) les poids d’intégration respectivement de \(g\) et \(s\) .

L’opérateur tangent s’écrit

(4535)\[DF(U)={\sum}_{\mathrm{el}}({\sum}_{g}{\omega}_{g}D{F}_{g}^{\mathrm{el}}(U)+{\sum}_{s}{\omega}_{s}D{F}_{s}^{\mathrm{el}}(U))\]

en notant \(D{F}_{g}^{\mathrm{el}}\) et \(D{F}_{s}^{\mathrm{el}}\) les valeurs respectivement au point de Gauss \(g\) et au sommet \(s\) de l’opérateur tangent.

L’opérateur tangent et les forces nodales sont donc calculées différemment aux points de Gauss et aux sommets. En revanche, toutes les composantes des contraintes généralisées et toutes les variables internes sont calculées à la fois aux points de Gauss et aux points d’intégration.

Vecteur forces internes : options RAPH_MECA et FULL_MECA

On adopte la décomposition suivante

(4535)\[{{\stackrel{ˆ}{E}}_{pi}^{\mathrm{el}}}^{T}.{\stackrel{ˉ}{\Sigma}}_{pi}^{\mathrm{el}}={\stackrel{ˉ}{\Sigma}}_{u}〚{\stackrel{ˆ}{u}}_{g}〛+{\stackrel{ˉ}{\Sigma}}_{p}\stackrel{ˆ}{p}+{\stackrel{ˉ}{\Sigma}}_{\nabla p}\nabla \stackrel{ˆ}{p}+{\stackrel{ˉ}{\Sigma}}_{{p}^{+}}{\stackrel{ˆ}{p}}^{+}+{\stackrel{ˉ}{\Sigma}}_{{p}^{-}}{\stackrel{ˆ}{p}}^{-}+{\stackrel{ˉ}{\Sigma}}_{{q}^{+}}{\stackrel{ˆ}{q}}^{+}+{\stackrel{ˉ}{\Sigma}}_{{q}^{-}}{\stackrel{ˆ}{q}}^{-}\]

où \({\stackrel{ˆ}{E}}_{pi}^{\mathrm{el}}=(〚{\stackrel{ˆ}{u}}_{g}〛,\stackrel{ˆ}{p},\nabla \stackrel{ˆ}{p},{\stackrel{ˆ}{p}}^{+},{\stackrel{ˆ}{p}}^{-},{\stackrel{ˆ}{q}}^{+},{\stackrel{ˆ}{q}}^{-})\) est une déformation généralisée virtuelle calculée à partir d’un vecteur déplacement généralisé virtuel.

À partir des formulations variationnelles discrètes et en répartissant les termes stationnaires aux points de Gauss et les termes instationnaires aux sommets de l’élément, on obtient :

• aux points de Gauss

(4535)\[\begin{split}\begin{array}{ccc}{\stackrel{ˉ}{\Sigma}}_{u}^{\mathrm{el}}& =& \tau '-{p}_{g}n\\ {\stackrel{ˉ}{\Sigma}}_{p}^{\mathrm{el}}& =& \Delta t({q}_{g}^{+}+{q}_{g}^{-})\\ {\stackrel{ˉ}{\Sigma}}_{\nabla p}^{\mathrm{el}}& =& \Delta tW\\ {\stackrel{ˉ}{\Sigma}}_{{p}^{+}}^{\mathrm{el}}& =& -\Delta t{q}^{+}\\ {\stackrel{ˉ}{\Sigma}}_{{p}^{-}}^{\mathrm{el}}& =& -\Delta t{q}^{-}\\ {\stackrel{ˉ}{\Sigma}}_{{q}^{+}}^{\mathrm{el}}& =& {p}^{+}-p\\ {\stackrel{ˉ}{\Sigma}}_{{q}^{-}}^{\mathrm{el}}& =& {p}^{-}-p\end{array}\end{split}\]

• aux sommets

(4535)\[\begin{split}\begin{array}{ccc}{\stackrel{ˉ}{\Sigma}}_{u}^{\mathrm{el}}& =& 0\\ {\stackrel{ˉ}{\Sigma}}_{p}^{\mathrm{el}}& =& {w}_{n}-{w}_{n+1}\\ {\stackrel{ˉ}{\Sigma}}_{\nabla p}^{\mathrm{el}}& =& 0\\ {\stackrel{ˉ}{\Sigma}}_{{p}^{+}}^{\mathrm{el}}& =& 0\\ {\stackrel{ˉ}{\Sigma}}_{{p}^{-}}^{\mathrm{el}}& =& 0\\ {\stackrel{ˉ}{\Sigma}}_{{q}^{+}}^{\mathrm{el}}& =& 0\\ {\stackrel{ˉ}{\Sigma}}_{{q}^{-}}^{\mathrm{el}}& =& 0\end{array}\end{split}\]

On a par ailleurs :

(4535)\[{{\stackrel{ˆ}{U}}^{\mathrm{el}}}^{T}.{R}_{pi}^{\mathrm{el}}={{\stackrel{ˆ}{E}}_{pi}^{\mathrm{el}}}^{T}.{\stackrel{ˉ}{\Sigma}}_{pi}^{\mathrm{el}}\]

ce qui donne

(4535)\[{R}_{pi}^{\mathrm{el}}={{Q}_{pi}^{\mathrm{el}}}^{T}.{\stackrel{ˉ}{\Sigma}}_{pi}^{\mathrm{el}}\]

Opérateur tangent : options FULL_MECA et RIGI_MECA_TANG

L’opérateur tangent au point de Gauss \(g\) est donné par :

(4535)\[\begin{split}\frac{\partial {\Sigma}_{g}^{\mathrm{el}}}{\partial {E}_{g}^{\mathrm{el}}}=\left[\begin{array}{cccccccc}\frac{\partial {\tau}_{n}'}{\partial 〚{u}_{n}〛}& \frac{\partial {\tau}_{n}'}{\partial 〚{u}_{t}〛}& -1& 0& 0& 0& 0& 0\\ \frac{\partial {\tau}_{t}}{\partial 〚{u}_{n}〛}& \frac{\partial {\tau}_{t}}{\partial 〚{u}_{t}〛}& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& \Delta t& \Delta t\\ \Delta t\frac{\partial W}{\partial {u}_{n}}& 0& \Delta t\frac{\partial W}{\partial p}& -\Delta tT& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& -\Delta t& 0\\ 0& 0& 0& 0& 0& 0& 0& -\Delta t\\ 0& 0& -1& 0& 1& 0& 0& 0\\ 0& 0& -1& 0& 0& 1& 0& 0\end{array}\right]\end{split}\]

L’opérateur tangent au sommet \(s\) est donné par :

(4535)\[\begin{split}\frac{\partial {\Sigma}_{s}^{\mathrm{el}}}{\partial {E}_{s}^{\mathrm{el}}}=\left[\begin{array}{cccccccc}0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ -\frac{\partial w}{\partial 〚{u}_{n}〛}& 0& -\frac{\partial w}{\partial p}& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0\end{array}\right]\end{split}\]

L’opérateur tangent au point d’intégration \(pi\) est finalement donné par :

(4535)\[D{F}_{pi}^{\mathrm{el}}={{Q}_{pi}^{\mathrm{el}}}^{T}.\frac{\partial {\Sigma}_{pi}^{\mathrm{el}}}{\partial {E}_{pi}^{\mathrm{el}}}.{Q}_{pi}^{\mathrm{el}}\]

Vecteur forces nodales : option FORC_NODA

L’option FORC_NODA est utilisée par STAT_NON_LINE lors de la phase de prédiction.

Au point d’intégration \(pi\) , le vecteur forces nodales au temps \(n+1\) est définie par :

(4535)\[{F}_{pi,n+1}^{\mathrm{el}}={{Q}_{pi}^{\mathrm{el}}}^{T}.{\stackrel{ˉ}{S}}_{pi}^{\mathrm{el}}\]

où \({\stackrel{ˉ}{S}}_{pi}^{\mathrm{el}}\) est nul aux sommets et est donné aux points de Gauss par :

(4535)\[\begin{split}\begin{array}{ccc}{\stackrel{ˉ}{S}}_{u,n+1}^{\mathrm{el}}& =& {\tau}_{g,n}-{p}_{g,n}n\\ {\stackrel{ˉ}{S}}_{p,n+1}^{\mathrm{el}}& =& \Delta t({q}_{g,n}^{+}+{q}_{g,n}^{-})\\ {\stackrel{ˉ}{S}}_{\nabla p,n+1}^{\mathrm{el}}& =& \Delta t{W}_{g,n}\\ {\stackrel{ˉ}{S}}_{{p}^{+},n+1}^{\mathrm{el}}& =& -\Delta t{q}_{g,n}^{+}\\ {\stackrel{ˉ}{S}}_{{p}^{-},n+1}^{\mathrm{el}}& =& -\Delta t{q}_{g,n}^{-}\\ {\stackrel{ˉ}{S}}_{{q}^{+},n+1}^{\mathrm{el}}& =& {p}_{g,n}^{+}-{p}_{g,n}\\ {\stackrel{ˉ}{S}}_{{q}^{-},n+1}^{\mathrm{el}}& =& {p}_{g,n}^{-}-{p}_{g,n}\end{array}\end{split}\]

Variables internes

Comme pour les modèles de THM classiques, les variables de \(1\) à \(\mathrm{NVIM}\) ( \(\mathrm{NVIM}\) dépendant de la loi de comportement mécanique utilisée) concernent la loi de mécanique utilisée. Les variables internes suivantes sont :

• \({V}_{\mathrm{NVIM}+1}\) : \(\rho -{\rho}_{0}\) , variation de la masse volumique,

• \({V}_{\mathrm{NVIM}+2}\) : \(\varepsilon\) , ouverture de la fissure.

Validation

Le tableau suivant récapitule les cas tests de validation des modèles présentés dans cette note.

Titre	Nom du test	Documentation	Modèle	Loi
Modélisation axisymétrique d’un joint avec couplage hydromécanique.	wtna111	V7.33.11	AXIS_JHMS	JOINT_BANDIS
Déplétion d’un réservoir	wtnp125	V7.32.125	PLAN_JHMS	JOINT_BANDIS
Injection de gaz dans un milieu poreux fracturé	wtnp126	V7.32.126	PLAN_JHMS	JOINT_BANDIS
Essais de fendage par coin du béton sous pression de fluide	wtnp128	V7.32.128	PLAN_JHMS	JOINT_BANDIS CZM_LIN_REG CZM_EXP_REG

Bibliographie