r7.01.10 Modélisations THHM. Généralités et algorithmes

Résumé

Les modules THMsont ceux qui traitent les équations de la mécanique des milieux continus en utilisant la théorie des milieux poreux éventuellement non saturés et en considérant que les phénomènes mécaniques, thermiques et hydrauliques sont complètement couplés. Nous présentons ici les équations d’équilibre, ou équations de conservation résolues par ces modules. Nous donnons une définition des contraintes généralisées et déformations généralisées, permettant de définir de façon assez générale ce qu’est une loi de comportement THM - du moins ce que les modules considérés considèrent ainsi - et permettant de traiter les équations non linéaires exhibées dans le cadre des algorithmes de l’opérateur STAT_NON_LINE.

Présentation du problème : hypothèses, notations

Dans ce chapitre, on s’attache principalement à présenter le milieu poreux et ses caractéristiques.

Description du milieu poreux

Le milieu poreux considéré est un volume constitué d’une matrice solide plus ou moins homogène, plus ou moins cohérente (très cohérente dans le cas du béton, peu dans le cas du sable). Entre les éléments solides, on trouve des pores. On distingue les pores fermés qui n’échangent rien avec leurs voisins et les pores connectés dans lesquels les échanges sont nombreux. Lorsque l’on parle de porosité, c’est bien des pores connectés dont on parle. A l’intérieur de ces pores se trouvent au plus deux constituants présents au plus sous deux phases. Le système est considéré comme fermé.

Notations

Les grandeurs associées à un constituant \(c\) présent sous une phase \(p\) sont notées \({X}_{c}^{p}\) . L’indice du composant \(c\) peut varier de \(1\) à \(2\) et celui de la phase également. Ces constituants peuvent être (et seront indicés comme tels en cas de besoin):

\(w\) pour l’eau liquide,

\(\text{ad}\) pour l’air dissous,

\(\text{as}\) pour l’air sec,

\(\text{vp}\) pour la vapeur d’eau.

Le milieu poreux à l’instant actuel est noté \(\Omega ` , sa frontière :math:\)partial Omega ` , et il est noté \({\Omega}_{0},\partial {\Omega}_{0}\) à l’instant initial.

\(n\) désigne la normale en un point de \(\partial \Omega\) , image de la normale \({n}_{0}\) à \(\partial {\Omega}_{0}\) . Nous noterons \(d(\partial \Omega )\) (respectivement \(d(\partial {\Omega}_{0})\) ) l’élément de surface de \(\partial \Omega\) (respectivement \(\partial {\Omega}_{0}\) ).

Le milieu est défini par:

• des paramètres (vecteur position \(x\) , temps \(t\) ),

• des variables (déplacements, pressions, température),

• des grandeurs intrinsèques (contraintes et déformations, apports massiques, chaleur, enthalpies, flux hydrauliques, thermique…).

Pour la phase solide, on fait l’hypothèse des petits déplacements.

Les différentes notations sont explicitées ci-après.

Variables descriptives du milieu

Ce sont les variables dont la connaissance en fonction du temps et du lieu permettent de connaître complètement l’état du milieu. Ces variables se décomposent en deux catégories:

• variables géométriques,

• variables d’état thermodynamique.

Variables géométriques

Dans tout ce qui suit, on adopte une représentation lagrangienne par rapport au squelette (au sens de [bib1]) et les coordonnées \(\mathrm{x}={\mathrm{x}}_{s}(t)\) sont celles d’un point matériel attaché au squelette. Tous les opérateurs de dérivation spatiaux sont définis par rapport à ces coordonnées.

Les déplacements du squelette sont notés \(u(x,t)=(\begin{array}{c}{u}_{x}\\ {u}_{y}\\ {u}_{z}\end{array})\) .

Variables d’état thermodynamiques

Les variables thermodynamiques sont:

• les pressions des constituants: étant donné que nous considérons qu’il y a au plus deux constituants, il y aura au plus deux équations de conservation de la masse, et donc par dualité au plus deux variables de pression,

• la température du milieu \(T(x,t)\) .

Champs descriptifs du milieu

Les inconnues principales, qui sont aussi les inconnues nodales (notées \(U(x,t)\) dans ce document) sont:

• 2 ou 3 (selon la dimension d’espace) déplacements \({u}_{x}(x,t),{u}_{y}(x,t),{u}_{z}(x,t)\) pour les modélisations KIT_HM, KIT_HHM, KIT_THM, KIT_THHM,

• la température \(T(x,t)\) pour les modélisations KIT_THH, KIT_THM, KIT_THHM,

• deux pressions \({p}_{1}(x,t),{p}_{2}(x,t)\) pour les modélisations KIT_HHM, KIT_THH, KIT_THHM,

• une pression \({p}_{1}(x,t)\) pour les modélisations KIT_HM, KIT_THM.

Grandeurs

Les équations d’équilibre sont:

• la conservation de la quantité de mouvement pour la mécanique,

• la conservation des masses de fluide pour l’hydraulique,

• la conservation de l’énergie pour la thermique.

Les équations d’équilibre font intervenir directement les contraintes généralisées. Les contraintes généralisées sont reliées aux déformations généralisées par les lois de comportement. Les déformations généralisées se calculent directement à partir des variables d’état et de leurs gradients spatiaux temporels.

Les lois de comportement peuvent utiliser des quantités annexes, souvent rangées dans les variables internes. Ces quantités ne sont pas décrites dans ce document qui ne traite pas des lois de comportement à proprement parler.

Grandeurs caractéristiques du milieu hétérogène

• La porosité eulérienne: \(\phi\) .

Si on note \({\Omega}_{\varphi}\) la partie du volume \(\Omega\) occupée par les vides dans la configuration courante, on a:

\(\varphi =\frac{{\Omega}_{\varphi}}{\Omega}\) éq 2.2.2.1-1

La définition de la porosité \(\varphi\) est donc celle de la porosité eulérienne.

• La saturation de la phase \(p\) : \({S}^{p}\) .

Si on note \({\Omega}^{p}\) le volume total occupé par la phase \(p\) , dans la configuration courante, on a par définition:

\({S}^{p}=\frac{{\Omega}^{p}}{{\Omega}_{j}}\) éq 2.2.2.1-2

Cette saturation est donc finalement une proportion variant entre \(0\) et \(1\) . Dans les équations de bilan, il est clair que c’est le produit de la porosité par la saturation \(\phi {S}^{p}\) qui va intervenir. On peut donc légitimement se demander pourquoi ce n’est pas cette quantité qui est prise comme inconnue. La réponse vient de ce que c’est la saturation \({S}^{p}\) qui intervient plus simplement dans les lois de comportement.

• La masse volumique eulérienne du constituant \(c\) dans la phase \(p\) : \({\rho}_{c}^{p}\) .

Si on note \({M}_{c}^{p}\) la masse de la phase \(p\) du constituant \(c\) , dans de volume \(\Omega\) du squelette dans la configuration courante, on a par définition:

\({M}_{c}^{p}={\int}_{{\Omega}^{p}}{{\rho}_{c}^{p}d\Omega }^{p}={\int}_{{\Omega}_{\phi }}{\rho}_{c}^{p}{S}^{p}d{\Omega}_{\phi }={\int}_{\Omega}{\rho}_{c}^{p}{S}^{p}\phid \Omega\) éq 2.2.2.1-3

La masse volumique de la phase \(p\) est simplement la somme des masses volumiques de ses constituants:

\({\rho}^{p}=\sum_{c}{\rho}_{c}^{p}\)

• La masse volumique homogénéisée lagrangienne: \(r\) .

A l’instant courant, la masse du volume \(\Omega\) , \({M}_{\Omega}\) est donnée par:

\({M}_{\Omega}={\int}_{{\Omega}_{0}}rd{\Omega}_{0}\) éq 2.2.2.1-4

Grandeurs mécaniques

• Le tenseur des déformations \(\varepsilon (u)(x,t)=\frac{1}{2}(\nabla u{}^{T}+\nabla u)\) ,

• Le tenseur des contraintes qui s’exercent sur le milieu poreux: \(\sigma\) .

Ce tenseur se décompose en un tenseur des contraintes effectives plus un tenseur de contraintes de pression \(\sigma =\sigma '+{\sigma}_{p}I\) . \(\sigma '\text{et}{\sigma}_{p}\) sont des composantes des contraintes généralisées. Ce découpage est finalement assez arbitraire, mais correspond tout de même à une hypothèse assez communément admise, au moins pour les milieux saturés en liquide.

Grandeurs hydrauliques

• Les apports massiques en constituants \({m}_{c}^{p}\) (unité: kilogramme par mètre cube). Ils représentent la masse de fluide apportée entre les instants initiaux et actuels. Ils font partie des contraintes généralisées.

\({m}_{c}^{p}=J{\rho}_{c}^{p}\phi {S}^{p}-{\rho}_{c}^{{p}_{0}}{\phi }_{0}{S}_{0}^{p}\) éq 2.2.2.3-1

Les apports massiques permettent de définir la masse volumique globale vue par rapport à la configuration de référence: \(r(x,t)={r}_{0}+{m}_{\mathrm{lq}}(x,t)+{m}_{\mathrm{vp}}(x,t)+{m}_{\mathrm{as}}(x,t)\) , où \({r}_{0}\) désigne la masse volumique homogénéisée à l’état initial.

• Les flux hydrauliques:

\({w}_{c}^{p}\) (unité: kilogramme/seconde/mètre carré) en représentation eulérienne

\({M}_{c}^{p}\) (unité: kilogramme/seconde/mètre carré) en représentation lagrangienne

On note \({v}_{c}^{p}\) la vitesse du constituant \(c\) dans la phase \(p\) , \(J\) le Jacobien de la transformation matérielle et \({v}_{S}=\frac{du}{\mathrm{dt}}\) la vitesse du squelette. \({\rho}_{c0}^{p},{\phi }_{0},{S}_{0}^{p}\) désignent les masses volumiques, la porosité et les saturations à l’instant initial. Par définition:

\({w}_{c}^{p}={\rho}_{c}^{p}\phi {S}^{p}({v}_{c}^{p}-{v}_{s})\) éq 2.2.2.3-2

La forme lagrangienne de

notée

est obtenue en écrivant:

\({M}_{c}^{p}.{n}_{0}d(\partial {\Omega}_{0})={w}_{c}^{p}.nd(\partial \Omega )\) éq 2.2.2.3-3

Les variables \({m}_{1},{M}_{1}\) et \({m}_{2},{M}_{2}\) se rapportent chacune à un constituant de masse conservative .

On pose par principe:

\(\begin{array}{}{m}_{1}={m}_{1}^{1}+{m}_{1}^{2};{M}_{1}={M}_{1}^{1}+{M}_{1}^{2}\\ {m}_{2}={m}_{2}^{1}+{m}_{2}^{2};{M}_{2}={M}_{2}^{1}+{M}_{2}^{2}\end{array}\)

Ce que nous écrirons:

\(\begin{array}{}{m}_{\mathrm{constituant}}=\sum_{\begin{array}{}\mathrm{nb}\mathrm{phase}\mathrm{du}\\ \mathrm{constituant}\end{array}}{m}_{\mathrm{constituant}}^{\mathrm{phase}}\\ {M}_{\mathrm{constituant}}=\sum_{\begin{array}{}\mathrm{nb}\mathrm{phase}\mathrm{du}\\ \mathrm{constituant}\end{array}}{M}_{\mathrm{constituant}}^{\mathrm{phase}}\end{array}\)

Dans les applications, on pourrait par exemple avoir:

• 2 constituants: air et eau,

• 2 phases pour l’eau,

• 1 phase pour l’air.

On aurait alors: \({m}_{1}^{1}\text{et}{M}_{1}^{1}\) : apport de masse et flux d’eau liquide

\({m}_{1}^{2}\text{et}{M}_{1}^{2}\) : apport de masse et flux de vapeur

\({m}_{2}^{1}\text{et}{M}_{2}^{1}\) : apport de masse et flux d’air sec

\({m}_{2}^{2}\text{et}{M}_{2}^{2}\) : inexistants

• Les pressions:

Puisque nous considérons qu’il peut y avoir deux constituants autres que le solide, il y a deux équations de conservation de la masse, et donc deux multiplicateurs associés, c’est-à-dire deux pressions \({p}_{1}\text{et}{p}_{2}\) . Aucune hypothèse n’est faite sur ce que signifient ces deux pressions \({p}_{1}\text{et}{p}_{2}\) . Cela dépendra des lois de comportement et de la façon de les écrire. Par exemple on peut choisir:

\(\begin{array}{}{p}_{1}=\text{pression capillaire}(\text{p}(\text{gaz})-\text{p}(\text{liquide}))\\ {p}_{2}=\text{pression de gaz}(\text{vapeur + air})\end{array}\)

Grandeurs thermiques

• La chaleur non convectée \(Q'\) (voir plus loin) (unité: Joule),

• Les enthalpies massiques des constituants \({{h}^{m}}_{c}^{p}\) (unité: Joule/Kelvin/kilogramme),

• Le flux de chaleur : \(q\) (unité: J/s/mètre carré).

Données extérieures

• La force massique \({\mathrm{F}}^{m}\) (en pratique la gravité),

• Les sources de chaleur :math:`Theta ` ,

• Les conditions aux limites portant soit sur des variables imposées, soit sur des flux imposés.

Dérivées particulaires, densités volumiques et massiques

La description que nous faisons du milieu est lagrangienne par rapport au squelette. On trouvera dans [bib1] une définition de la notion de squelette: «la matrice (partie solide+porosité occluse) constituant la partie matérielle du squelette et l’espace poreux connecté du volume élémentaire en question constituent le point matériel du squelette ou particule du squelette».

Soit \(a\) un champ quelconque sur \(\Omega\) , soit \({x}_{s}(t)\) la coordonnée d’un point attaché au squelette que nous suivons dans son mouvement et soit \({x}_{\mathrm{fl}}(t)\) la coordonnée d’un point attaché au fluide. On note \(\dot{a}=\frac{{d}^{s}a}{\mathrm{dt}}\) la dérivée temporelle dans le mouvement du squelette:

\(\dot{a}=\frac{{d}^{s}a}{\mathrm{dt}}=\underset{\deltat \to 0}{\lim}\frac{a({x}_{s}(t+\deltat ),t+\deltat )-a({x}_{s}(t),t)}{\deltat }\)

\(\dot{a}\) est appelée dérivée particulaire et souvent notée \(\frac{\mathrm{da}}{\mathrm{dt}}\) (par exemple dans [bib1]). Nous préférons utiliser une notation qui rappelle que la configuration utilisée pour repérer une particule est celle du squelette par rapport auquel une particule de fluide a une vitesse relative. Pour une particule de fluide le repérage \({x}_{s}(t)\) est quelconque, c’est à dire que la particule de fluide qui occupe la position \({x}_{s}(t)\) à l’instant \(t\) n’est pas la même que celle qui occupe la position \({x}_{S}(t')\) à un autre instant \(t'\) .

Soit alors \(A=\underset{\Omega}{\int}ad\Omega\) une quantité liée à une densité volumique \(a\) , laquelle densité est elle même portée en partie par les grains solides et par les fluides . Soit \({{a}^{m}}_{c}^{p}\) la densité massique de \(a\) portée par la phase fluide \(p\) du constituant \(c\) et soit \({a}_{s}\) la densité volumique de \(a\) liée aux grains solides. Toutes ces définitions reviennent finalement à écrire:

\(A={\int}_{\Omega}\mathrm{ad}\Omega ={A}_{s}+{A}_{\mathrm{fl}}={\int}_{\Omega}{a}_{s}d\Omega +{\int}_{\Omega}{a}_{\mathrm{fl}}d\Omega ={\int}_{\Omega}({a}_{s}+\sum_{p,c}{\rho}_{c}^{p}{\mathrm{jS}}^{p}{{a}^{m}}_{c}^{p})d\Omega\) éq 2.3-1

En suivant [bib1], nous notons \(\frac{{d}^{\mathrm{fl}}{A}_{\mathrm{fl}}}{\mathrm{dt}}\) la dérivée de \({A}_{\mathrm{fl}}\) si nous suivons \(\Omega\) dans le mouvement du fluide et \(\frac{{d}^{s}{A}_{s}}{\mathrm{dt}}\) la dérivée de \({A}_{s}\) si nous suivons \(\Omega\) dans le mouvement du squelette.

Nous définissons alors:

\(\frac{\mathrm{DA}}{\mathrm{Dt}}=\frac{{d}^{s}{A}_{s}}{\mathrm{dt}}+\frac{{d}^{\mathrm{fl}}{A}_{\mathrm{fl}}}{\mathrm{dt}}=\frac{{d}^{s}}{\mathrm{dt}}\int{}_{\Omega}\text{}{a}_{s}d\Omega +\frac{{d}^{\mathrm{fl}}}{\mathrm{dt}}\int{}_{\Omega}\text{}\sum_{p,c}{\rho}_{c}^{p}\phi {S}^{p}{a}^{{m}_{c}^{p}}d\Omega\) éq 2.3-2

La densité \({{a}^{m}}_{c}^{p}\) est transportée avec une vitesse relative de \(({v}_{c}^{p}-{v}_{s})\) par rapport au squelette. Compte tenu de la définition de \(\dot{a}=\frac{{d}^{s}a}{\mathrm{dt}}\) , et de la définition \({w}_{c}^{p}={\rho}_{c}^{p}\phi {S}^{p}({v}_{c}^{p}-{v}_{s})\) , on voit facilement que la dérivée totale de \(A\) par rapport au temps s’écritfinalement:

\(\frac{\mathrm{DA}}{\mathrm{Dt}}=\int{}_{\Omega}\text{}(\dot{a}+\sum_{p,c}\text{Div}({{a}^{m}}_{c}^{p}{w}_{c}^{p}))d\Omega\) éq 2.3-3

Remarque:

Dans la mesure où nous avons fait l’hypothèse des petits déplacements du squelette, \(\dot{a}=\frac{{d}^{s}a}{\mathrm{dt}}\) peut se confondre avec la dérivée partielle par rapport au temps \(\frac{\partial a}{\partial t}\) et \({v}_{s}\) peut être considérée comme nulle. De même, dans la suite de la note nous confondrons les représentations lagrangiennes et eulériennes des flux, \({M}_{c}^{p}\) et \({w}_{c}^{p}\) .

Équations continues

Mécanique : conservation de la quantité de mouvement

Nous notons \(\sigma\) le tenseur des contraintes de Cauchy et \(s\) le second tenseur (symétrique) de Piola-Kirchhoff.

Nous notons \(P\) le gradient de la transformation \({x}_{0}={x}_{S}(0)\to {x}_{S}({x}_{0},t)\) .

\(P=\frac{\partial {x}_{S}({x}_{0},t)}{\partial {x}_{0}}\)

On a: \(s=\detP.{P}^{-1}.\sigma .{P}^{-T}\) .

Les équations d’équilibre mécanique s’écrivent dans la configuration \({\Omega}_{0}\) :

\({\text{Div}}_{0}(P.s)+r{F}^{m}=0\)

Nous avons noté \({\text{Div}}_{0}\) l’opérateur de divergence par rapport aux variables d’espace \({x}_{0}\) de la configuration \({\Omega}_{0}\) .

Dans la mesure où nous faisons l’hypothèse des petits déplacements et des petites déformations du squelette, cette équation peut être approchée par:

\(\text{Div}\sigma +r{F}^{m}=0\) éq 3.1-1

Nous verrons plus loin que nous adoptons toujours la décomposition \(\sigma =\sigma '+{\sigma}_{p}I\) , où \(\sigma '\) désigne la contrainte effective. C’est donc à la charge du module d’intégration des équations d’équilibre de faire la somme: \(\sigma =\sigma '+{\sigma}_{p}I\) .

Hydraulique : conservation de la masse

L’écriture eulérienne de la conservation de la masse fluide pour le constituant \(c\) s’écrit:

\(\frac{{d}^{\mathrm{fl}}}{\mathrm{dt}}\int{}_{\Omega}\text{}\sum_{p}{\rho}_{c}^{p}\phi {S}^{p}d\Omega =0\)

On peut alors appliquer [éq 2.3-1] en prenant: \({a}_{s}=0\) et \({{a}^{m}}_{c}^{p}=1\) et [éq 2.3-3] donnera:

\(\sum_{p}\frac{{d}^{s}{\rho}_{c}^{p}\phi {S}^{p}}{\mathrm{dt}}+\sum_{p}\text{Div}({w}_{c}^{p})=0\)

En utilisant la définition des apports massiques [éq 2.2.2.3-3] , la définition des flux lagrangiens [éq2.2.2.3-2] on trouve la forme lagrangienne de la conservation de la masse fluide:

\(\lbrace \begin{array}{}\dot{{m}_{1}}+{\text{Div}}_{0}({M}_{1})=0\\ \dot{{m}_{2}}+{\text{Div}}_{0}({M}_{2})=0\end{array}\) éq 3.2-1

Équation de l’énergie

Pour les fonctions thermodynamiques, nous adoptons systématiquement une décomposition du type [éq2.3-1]. Cela correspond au fait que les différentes énergies ont toutes une partie portée par le solide et une partie portée par les fluides. La partie portée par le solide est caractérisée par une densité volumique alors que les parties portées par le fluide sont caractérisées par des densités massiques, comme nous l’avons montré au paragraphe [§2.3].

Énergie internetotale : \(E=\int{}_{\Omega}\text{}({e}_{s}+\sum_{p,c}{\rho}_{c}^{p}\phi {S}^{p}{{e}^{m}}_{c}^{p})d\Omega\) éq 3.3.1

Entropie totale: \(S=\int{}_{\Omega}\text{}({s}_{s}+\sum_{p,c}{\rho}_{c}^{p}\phi {S}^{p}{{s}^{m}}_{c}^{p})d\Omega\) éq 3.3.2

Enthalpie totale: \(H=\int{}_{\Omega}\text{}({h}_{s}+\sum_{p,c}{\rho}_{c}^{p}\phi {S}^{p}{{h}^{m}}_{c}^{p})d\Omega\) éq 3.3.3

Énergie libre: \(\lbrace \begin{array}{}\Psi =E-TS\\ {\Psi}_{s}={e}_{s}-{\text{Ts}}_{s}\\ {{\Psi}^{m}}_{c}^{p}={{e}^{m}}_{c}^{p}-T{{s}^{m}}_{c}^{p}\end{array}\) éq 3.3.4

Enthalpie libre: \(\lbrace \begin{array}{}G=H-TS\\ {g}_{s}={h}_{s}-T{s}_{s}\\ {{g}^{m}}_{c}^{p}={{h}^{m}}_{c}^{p}-T{{s}^{m}}_{c}^{p}\end{array}\) éq 3.3.5

Enfin, en notant \(\dot{Q}(\Omega )\) le taux de chaleur reçu par un volume \(\Omega\) , on a par définition:

\(\dot{Q}(\Omega )=\underset{\partial \Omega }{\int}q.nd\Gamma +\underset{\Omega}{\int}\Thetad \Omega\) éq 3.3.6

On rappelle enfin que l’enthalpie des fluides se calcule par la formule:

\(h=e+\frac{p}{\rho}\) éq 3.3.7

Le premier principe

Avec les définitions données plus haut, il s’écrit:

\(-\sum_{p,c}\text{Div}\left({{h}^{m}}_{c}^{p}{\mathrm{M}}_{c}^{p}\right)+\sigma :\dot{\epsilon}+\sum_{p,c}{\mathrm{M}}_{c}^{p}.{\mathrm{F}}^{m}+\Theta -\text{Div}q=0\) éq 3.3.1-1

Cette écriture correspond à l’équation (22) du chapitre III-2-3 de [bib1], dans laquelle nous avons négligé les termes d’inertie. Pour les milieux homogènes, elle correspond à l’équation (31) du paragraphe IV-3-2 de [bib3].

Le deuxième principe

Sa forme assez bien connue est:

\(\dot{s}+\sum_{p,c}\text{Div}({{s}^{m}}_{c}^{p}{M}_{c}^{p})+\text{Div}(\frac{q}{T})-\frac{\Theta}{T}\ge 0\) éq 3.3.2-1

En utilisant les considérations thermodynamiques classiques [bib1] liées à l’introduction de l’enthalpie libre [éq 3.3.5], on montre que l’on doit nécessairement avoir:

\(\sigma -\frac{\partial \Psi }{\partial \varepsilon }=0\) éq 3.3.2-2

\({{g}^{m}}_{c}^{p}-\frac{\partial \Psi }{\partial {m}_{c}^{p}}=0\) éq 3.3.2-3

\(s+\frac{\partial \Psi }{\partial T}=0\) éq 3.3.2-4

Équation de l’énergie

Assez souvent, on considère que, les transformations étant réversibles, le second principe fournit finalement une égalité. De plus, on remplace dans [éq 3.3.2-1] la température inconnue \(T\) par une valeur constante dite température de référence. Il s’agit finalement d’une linéarisation de [éq 3.3.2-1] justifiée si les variations de températures sont «petites». Notons que le terme de transport \(\sum_{p,c}\text{Div}({{s}^{m}}_{c}^{p}{M}_{c}^{p})\) complique le traitement de la non linéarité due à la présence de la température en dénominateur des autres termes de [éq 3.3.2-1].

Nous travaillons en enthalpie de façon à lever cette difficulté. On part de l’équation du premier principe [éq 3.3.1-1] dans laquelle on injecte les équations [éq 3.3.2-2], [éq 3.3.2-3], [éq 3.3.2-4], et la définition de l’enthalpie libre [éq 3.3-5] et l’on obtient:

\(T\dot{s}+\sum_{p,c}({{h}^{m}}_{c}^{p}\dot{{m}_{c}^{p}}-T{{s}^{m}}_{c}^{p}\dot{{m}_{c}^{p}})=-\sum_{p,c}\text{Div}({{h}^{m}}_{c}^{p}{M}_{c}^{p})+\sum_{p,c}{M}_{c}^{p}.{F}^{m}+\Theta -\text{Div}q\) éq 3.3.3-1

On pose alors:

\(\deltaq '=T\deltas -T\sum_{p,c}{s}^{{m}_{c}^{p}}\delta {m}_{c}^{p}\) éq 3.3.3-2

La quantité \(Q'\) a la dimension d’une énergie par unité de volume. Elle représente la chaleur reçue par le système dans une transformation pour laquelle il n’y a pas d’apports de chaleur par entrée de fluide ayant une enthalpie. Bien que \(\deltaq '\) ne soit pas une différentielle exacte, nous prenons cette quantité comme variable d’état.

Finalement, l’équation d’énergie retenue a la forme suivante:

\(\sum_{p,c}{{h}^{m}}_{c}^{p}\dot{{m}_{c}^{p}}+\dot{Q'}+\sum_{p,c}\text{Div}({{h}^{m}}_{c}^{p}{M}_{c}^{p})+\text{Div}q-\sum_{p,c}{M}_{c}^{p}.{F}^{m}=\Theta\) éq 3.3.3-3

Écriture variationnelle des équations d’équilibre

Mécanique

Nous notons \({U}_{ad}\) l’ensemble des champs de déplacement cinématiquement admissibles, c’est-à-dire les éléments de \({({H}^{1}(\Omega ))}^{3}\) vérifiant les conditions aux limites en déplacement sur la partie de \(\partial \Omega\) supportant de telles conditions [bib3].

La forme variationnelle de [éq 3.1-1] est:

\(\lbrace \begin{array}{c}\sigma =\sigma '+{\sigma}_{p}\mathrm{I}\\ \int{}_{\Omega}\text{}\sigma .\epsilon \left(\mathrm{v}\right)=\int{}_{\Omega}\text{}r{\mathrm{F}}^{m}.\mathrm{v}d\Omega +\int{}_{\partial \Omega }\text{}{\mathrm{f}}^{\mathit{ext}}.\mathrm{v}d\Gamma \text{}\forall \mathrm{v}\in {U}_{\mathit{ad}}\end{array}\) éq 4.1-1

Hydraulique

Nous notons \({P}_{1\mathrm{ad}}\) (resp. \({P}_{2\mathrm{ad}}\) ) l’ensemble des champs de pression admissibles \({\pi}_{1}\) (resp. \({\pi}_{2}\) ), c’est-à-dire les éléments de \({H}^{1}(\Omega )\) vérifiant les conditions aux limites en pression \({P}_{1}\) (resp. \({P}_{2}\) ) sur la partie de \(\partial \Omega\) supportant de telles conditions [bib3]. La forme variationnelle de [éq 3.2-1] est:

\(\begin{array}{}-{\int}_{\Omega}(\dot{{m}_{1}^{1}}+\dot{{m}_{1}^{2}}).{\pi}_{1}+{\int}_{\Omega}({M}_{1}^{1}+{M}_{1}^{2}).\nabla {\pi}_{1}d\Omega =\text{}\\ {\int}_{\partial \Omega }({M}_{1\mathrm{ext}}^{1}+{M}_{1\mathrm{ext}}^{2}).{\pi}_{1}d\Gamma \forall {\pi}_{1}\in {P}_{1\mathrm{ad}}\\ -{\int}_{\Omega}(\dot{{m}_{2}^{1}}+\dot{{m}_{2}^{2}}).{\pi}_{2}+{\int}_{\Omega}({M}_{2}^{1}+{M}_{2}^{2}).\nabla {\pi}_{2}d\Omega =\text{}\\ {\int}_{\partial \Omega }({M}_{2\mathrm{ext}}^{1}+{M}_{2\mathrm{ext}}^{2}).{\pi}_{2}d\Gamma \forall {\pi}_{2}\in {P}_{2\mathrm{ad}}\end{array}\rbrace\) éq 4.2-1

qui fait apparaître les flux hydrauliques scalaires \({M}_{i\mathrm{ext}}^{j}\) sur les bords.

Thermique

Nous notons \({T}_{\mathrm{ad}}\) l’ensemble des champs de température admissibles \(\tau\) , c’est-à-dire les éléments de \({H}^{1}(\Omega )\) vérifiant les conditions aux limites en température sur la partie de \(\partial \Omega\) supportant de telles conditions. [bib3]. La forme variationnelle de [éq 3.3.3-3] est:

\(\begin{array}{}{\int}_{\Omega}\dot{Q}'.\tau d\Omega +\sum_{p,c}{\int}_{\Omega}{{h}^{m}}_{c}^{p}\dot{{m}_{c}^{p}}.\tau d\Omega -{\int}_{\Omega}(\sum_{p,c}{{h}^{m}}_{c}^{p}{M}_{c}^{p}+q).\nabla \tau d\Omega =\text{}\\ {\int}_{\Omega}(\Theta +\sum_{p,c}{M}_{c}^{p}.{F}^{m}).\tau d\Omega -{\int}_{\partial \Omega }(\sum_{p,c}{{h}^{m}}_{c}^{p}{M}_{c\mathrm{ext}}^{p}+{q}_{\mathrm{ext}}).\tau d\Gamma \\ \forall \tau \in {T}_{\mathrm{ad}}\end{array}\) éq 4.3-1

Notons que, contrairement à d’autres présentations, et notamment [bib8] nous n’avons pas injecté les équations de conservation de la masse, et nous avons intégré par partie le terme de transport \(\sum_{p,c}\text{Div}({{h}^{m}}_{c}^{p}{M}_{c}^{p})\) .

Ce dernier point a pour avantage de ne pas faire apparaître des dérivées d’ordre supérieur, et, au contraire de faire apparaître naturellement des conditions aux limites relatives à l’entrée de chaleur liée aux flux hydrauliques: \(\sum_{p,c}{\int}_{\partial \Omega }{{h}^{m}}_{c}^{p}{M}_{c\mathrm{ext}}^{p}.\tau d\Gamma\) .

On pourra en fait considérer que les conditions de flux thermiques définissent directement:

\(\stackrel{~}{{q}_{\mathit{ext}}}={{h}^{m}}_{c}^{p}{M}_{c\mathit{ext}}^{p}+{q}_{\mathit{ext}}\)

Discrétisation en temps

Dans ce chapitre, nous nous contentons de reprendre les formulations variationnelles en leur appliquant une discrétisation par rapport au temps de type théta-schéma pour l’hydraulique et la thermique. Il s’agit d’une méthode générale d’intégration des équations différentielles [bib12] et [bib13].

\(\theta\) est un paramètre numérique compris entre \(0\) et \(1\) . Pour les équations différentielles linéaires (ce qui n’est pas notre cas …) ce schéma est inconditionnellement stable pour \(\theta \ge 1/2\) , il est d’ordre \(1\) pour \(\theta \ne 1/2\) et d’ordre 2 pour \(\theta =1/2\) . Néanmoins, il peut être préférable d’utiliser une valeur différente de \(1/2\) , et ceci pour des raisons d’oscillations parasites [bib12].

Les quantités indicées par \(+\) sont les quantités en fin de pas de temps, et celles indicées par \(–\) sont celles du début du pas de temps. On note:

\(\Delta t={t}^{+}-{t}^{-}\)

\({a}^{\theta}=\theta {a}^{+}+\left(1-\theta \right){a}^{-}\text{}\forall a\)

Mécanique

\(\lbrace \begin{array}{c}{\sigma}^{+}=\sigma {'}^{+}+{\sigma}_{p}^{+}\mathrm{I}\\ {\int}_{\Omega}{\sigma}^{+}.\varepsilon (\mathrm{v})d\Omega ={\int}_{\Omega}{r}^{+}{\mathrm{F}}^{\mathrm{m}+}.\mathrm{v}d\Omega +{\int}_{\partial \Omega }{\mathrm{f}}^{\mathrm{ext}+}.vd\Gamma \text{}\forall \mathrm{v}\in {U}_{\mathit{ad}}\end{array}\) éq 5.1-1

Hydraulique

\(\begin{array}{}-{\int}_{\Omega}({m}_{1}^{{1}^{+}}+{m}_{1}^{{2}^{+}}).{\pi}_{1}d\Omega +\theta \Delta t{\int}_{\Omega}({M}_{1}^{{1}^{+}}+{M}_{1}^{{2}^{+}}).\nabla {\pi}_{1}d\Omega =\\ -{\int}_{\Omega}({m}_{1}^{{1}^{-}}+{m}_{1}^{{2}^{-}}).{\pi}_{1}d\Omega -(1-\theta )\Delta t{\int}_{\Omega}({M}_{1}^{{1}^{-}}+{M}_{1}^{{2}^{-}}).\nabla {\pi}_{1}d\Omega \\ +\Delta t{\int}_{\partial \Omega }({M}_{1\mathrm{ext}}^{{1}^{\theta}}+{M}_{1\mathrm{ext}}^{{2}^{\theta}}).{\pi}_{1}d\Gamma \forall {\pi}_{1}\in {P}_{1\mathrm{ad}}\\ -{\int}_{\Omega}({m}_{2}^{{1}^{+}}+{m}_{2}^{{2}^{+}}).{\pi}_{2}d\Omega +\theta \Delta t{\int}_{\Omega}({M}_{2}^{{1}^{+}}+{M}_{2}^{{2}^{+}}).\nabla {\pi}_{2}d\Omega =\\ -{\int}_{\Omega}({m}_{2}^{{1}^{-}}+{m}_{2}^{{2}^{-}}).{\pi}_{2}d\Omega -(1-\theta )\Delta t{\int}_{\Omega}({M}_{2}^{{1}^{-}}+{M}_{2}^{{2}^{-}}).\nabla {\pi}_{2}d\Omega \\ +\Delta t{\int}_{\partial \Omega }({M}_{2\mathrm{ext}}^{{1}^{\theta}}+{M}_{2\mathrm{ext}}^{{2}^{\theta}}).{\pi}_{2}d\Gamma \forall {\pi}_{2}\in {P}_{2\mathrm{ad}}\end{array}\rbrace\) éq 5.2-1

Thermique

\(\begin{array}{}{\int}_{\Omega}(Q{'}^{+}-Q{'}^{-})\tau d\Omega -\theta \Delta t{\int}_{\Omega}(\sum_{p,c}{h}_{c}^{mp+}{M}_{c}^{p+}+{q}^{+})\nabla \tau d\Omega \\ -(1-\theta )\Delta t{\int}_{\Omega}(\sum_{p,c}{h}_{c}^{mp-}{M}_{c}^{p-}+{q}^{-})\nabla \tau d\Omega +\theta {\int}_{\Omega}(\sum_{p,c}{h}_{c}^{mp+}({m}_{c}^{p+}-{m}_{c}^{p-}))\tau d\Omega \\ +(1-\theta ){\int}_{\Omega}{h}_{c}^{mp-}({m}_{c}^{p+}-{m}_{c}^{p-})\tau d\Omega =\\ +\theta \Delta t{\int}_{\Omega}\sum_{p,c}{M}_{c}^{p+}.{F}^{m}\tau d\Omega +(1-\theta )\Delta t{\int}_{\Omega}\sum_{p,c}{M}_{c}^{p-}.{F}^{m}\tau d\Omega \\ +\Delta t{\int}_{\Omega}{\Theta}^{\theta}\tau d\Omega -\Delta t{\int}_{\partial \Omega }(\sum_{p,c}{h}_{c}^{m{p}^{\theta}}{M}_{c\mathrm{ext}}^{{p}^{\theta}}+{q}_{\mathrm{ext}}^{\theta}).\tau d\Gamma \forall \tau \in {T}_{\mathrm{ad}}\end{array}\rbrace\) éq 5.3-1

On peut de nouveau considérer que les conditions de flux thermiques définissent directement:

\({\tilde{q}}_{\mathrm{ext}}^{\theta}=\sum_{p,c}{h}_{c}^{m{p}^{\theta}}{M}_{c\mathrm{ext}}^{{p}^{\theta}}+{q}_{\mathrm{ext}}^{\theta}\)

Principe des travaux virtuels, déformations et contraintes généralisés, lois de comportement

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

En se reportant aux formulations variationnelles [éq 4.1-1], [éq 4.2-1] et [éq 4.3-1], il apparaît que l’on peut choisir:

Pour les contraintes généralisées:

\(\Sigma =\left\lbrace \begin{array}{}\sigma ',{\sigma}_{p};\\ {m}_{1}^{1},{M}_{1}^{1},{h}_{1}^{m1};{m}_{1}^{2},{M}_{1}^{2},{h}_{1}^{m2};\\ {m}_{2}^{1},{M}_{2}^{1},{h}_{2}^{m1};{m}_{2}^{2},{M}_{2}^{2},{h}_{2}^{m2};\\ Q',q\end{array}\right\rbrace\) éq 6.1-1

Pour les déformations généralisées:

\(\mathrm{\Epsilon }=\left\lbrace u,\varepsilon (u);{p}_{1},\nabla {p}_{1};{p}_{2},\nabla {p}_{2};T,\nabla T\right\rbrace\) éq 6.1-2

On remarque le fait que les déformations généralisées contiennent les déplacements. Cela est dû au terme \({\int}_{\Omega}r{F}^{m}.v\) de la formulation variationnelle de l’équation de conservation de la quantité de mouvement [éq 4.1-1], lequel terme couple finalement les contraintes généralisées et les déplacements du fait de [éq 6.3.4-1]. Les déformations généralisées contiennent la pression et la température parce que les équations associées sont paraboliques.

Principe des travaux virtuels

L’ensemble des équations non linéaires à résoudre peut se mettre sous la forme:

\(\mathrm{R}\left(\mathrm{U}\right)={\mathrm{L}}^{\mathit{meca}}\) éq 6.2-1

où \(U\) désigne les déplacements généralisés, c’est-à-dire: \(U=\left\lbrace {u}_{x},{u}_{y},{u}_{z},{p}_{1},{p}_{2},T\right\rbrace\) dans le cas le plus général. Les forces internes \(R\) s’expriment à partir d’un principe de travaux virtuels généralisé. Dans le cas de la mécanique des milieux continus «classique», c’est à dire quand il n’y a pas d’autre constituant que le solide, on a l’habitude de définir les forces internes par:

\({w}^{T}.R={\int}_{\Omega}\varepsilon (w).\sigma d\Omega\) \(\forall w\) , champ de déplacement cinématiquement admissible.

Dans cette formulation, le champ de déformation \(\varepsilon\) ne dépend que du champ de déplacement et de ses dérivées spatiales, éventuellement de façon non linéaire si l’on prend en compte des déformations finies. On écrit symboliquement:

\(R={Q}^{T}\sigma\)

La loi de comportement relie les contraintes \(\sigma\) aux déformations \(\varepsilon\) .

Dans le cadre de la théorie des milieux poreux développée ici, nous essayons de nous rapprocher le plus possible de cette formulation en introduisant des contraintes généralisées \(\Sigma\) et des déformations généralisées \(E\) ; Les déformations généralisées ne dépendent que du champ de déplacements généralisés \(U\) et de ses dérivées spatiales. L’opérateur \(U\to E(U)\) est un opérateur de dérivation par rapport au champ de coordonnées.

La loi de comportement permet de calculer \(\Sigma\) en fonction de \(E\) .

Par contre, nous ne pouvons pas écrire directement \({W}^{T}.R={\int}_{\Omega}\mathrm{\Epsilon }(W).\Sigmad \Omega\) , pour les raisons suivantes:

• les équations que nous traitons sont des équations évolutives en temps et les dérivées par rapport au temps des quantités interviennent,

• les équations sont non linéaires à cause des termes de transport liés à la représentation eulérienne des fluides: ces termes non linéaires ne figurent que dans l’équation de thermique,

• le choix des inconnues fait que les termes non linéaires de transport interviennent dans les contraintes généralisées. Soit un terme de transport dans l’équation [éq 4.3-1], \({\int}_{\Omega}(\sum_{p,c}{h}_{c}^{mp}{M}_{c}^{p}+q).\nabla \taud \Omega\) , étant donné que l’on a pris comme inconnues principales pour l’hydraulique les pressions, les quantités \({M}_{c}^{p}\) liées aux vitesses des fluides appartiennent aux contraintes généralisées, de même que les enthalpies \({h}_{c}^{mp}\) , et le terme de transport donné en exemple est linéaire en déformation généralisée et quadratique en contrainte généralisée. Ceci fait une différence avec la formulation de la théorie des milieux continus classiques où les termes de grande déformation sont des quantités quadratiques des déformations.

Pour toutes ces raisons, nous introduisons un champ noté \(\overline{\Sigma}\) tel que:

\(R={Q}^{T}\overline{\Sigma}\) éq 6.2-2

\({Q}^{T}\) est défini par:

\(\begin{array}{}{W}^{T}.R={\int}_{\Omega}E(W).\overline{\Sigma}d\Omega \\ \forall W\text{champ de déplacement généralisé cinématiquement admissible}\end{array}\) éq 6.2-3

On voit facilement et très classiquement que \({Q}^{T}\) est le transposé de l’opérateur \(Q\) tel que:

\(E=\mathrm{QU}\) éq 6.1-4

Le champ \(\overline{\Sigma}\) est une fonction linéaire de \(\dot{\Sigma}\) et non linéaire de \(\Sigma\) :

\(\overline{\Sigma}=\overline{\Sigma}(\dot{\Sigma},\Sigma )\) éq 6.1-5

Après discrétisation en temps, \(\overline{{\Sigma}^{+}}\) deviendra une fonction non linéaire de \({\Sigma}^{+}\) et \({\Sigma}^{-}\) :

\(\overline{{\Sigma}^{+}}=\overline{{\Sigma}^{+}}({\Sigma}^{+},{\Sigma}^{-})\) éq 6.1-6

Notons enfin que pour des raisons algorithmiques (entre autres), on a besoin de connaître la dérivée des forces internes par rapport aux déplacements généralisés:

\(\frac{\partial R}{\partial U}=\frac{\partial R}{\partial \overline{\Sigma}}\frac{\partial \overline{\Sigma}}{\partial \Sigma }\frac{\partial \Sigma }{\partial E}\frac{\partial E}{\partial U}={Q}^{T}\frac{\partial \overline{\Sigma}}{\partial \Sigma }\frac{\partial \Sigma }{\partial \Sigma }Q\)

Il est clair que \(\frac{\partial \overline{\Sigma}}{\partial \Sigma }\) ne dépend que de la forme des équations d’équilibre.

Lois de comportement

Une loi de comportement sera définie simplement comme une relation quelconque entre contraintes généralisées et déformations généralisées. Les variables internes sont définies comme des champs nécessaires au calcul des contraintes, dont l’évolution est donnée par les lois de comportement, mais qui n’interviennent pas directement dans les équations d’équilibre.

De plus, nous considérons que les lois de comportement sont écrites sous forme incrémentale et qu’elles sont locales. En notant :math:`` les variables internes, une loi de comportement est donc une relation:

\(\Sigma ,\khi ,\dot{E}\to \dot{\Sigma},\dot{\khi }\)

Après discrétisation en temps, la loi de comportement devient une relation:

\({\Sigma}^{-},{\khi }^{-},{E}^{-},{E}^{+}\to {\Sigma}^{+},{\khi }^{+}\)

La loi de comportement devra également fournir le seul terme qui dans l’expression \(\frac{\partial R}{\partial U}={Q}^{T}\frac{\partial \overline{\Sigma}}{\partial \Sigma }\frac{\partial \Sigma }{\partial E}Q\) dépend d’elle, à savoir \(\frac{\partial \Sigma }{\partial E}\) . Finalement une relation de comportement est une relation:

\({\Sigma}^{-},{\khi }^{-},{E}^{-},{E}^{+}\to {\Sigma}^{+},{\khi }^{+},\frac{\partial {\Sigma}^{+}}{\partial {E}^{+}}\) éq 6.3-1

Dans les paragraphes suivants, nous précisons certains aspects des lois de comportement en distingunt les contributions mécaniques, hydrauliques et thermiques.

Loi de comportement mécanique

Écriture générale

Les variables internes sont notées \(\khi\) . Une loi de comportement de mécanique, dans le cadre THM s’écrit:

\(\lbrace \begin{array}{}{\sigma}^{+}={\sigma}^{+}({\varepsilon}^{+},{p}_{1}^{+},{p}_{2}^{+},{T}^{+};{\varepsilon}^{-},{p}_{1}^{-},{p}_{2}^{-},{T}^{-},{\sigma}^{-},{\khi }^{-})\\ {\khi }^{+}={\khi }^{+}({\varepsilon}^{+},{p}_{1}^{+},{p}_{2}^{+},{T}^{+};{\varepsilon}^{-},{p}_{1}^{-},{p}_{2}^{-},{T}^{-},{\sigma}^{-},{\khi }^{-})\end{array}\) éq 6.3.1.1-1

Cas des contraintes effectives

Dans le cas de l’hypothèse des contraintes effectives, on a la décomposition: \(\sigma =\sigma '+{\sigma}_{p}I\) où \(\sigma '\) est le tenseur des contraintes effectives et \({\sigma}_{p}\) est un scalaire.

Les variables internes sont séparées en deux parties: les variables internes mécaniques \({\khi }_{\sigma}\) et les variables internes hydrauliques \({\khi }_{H}\) . La loi de comportement mécanique se scinde alors en deux lois, dont la première peut être une loi déjà existante dans le cadre habituel de la thermomécanique.

\(\lbrace \begin{array}{}\sigma {'}^{+}=\sigma {'}^{+}({\varepsilon}^{+},{T}^{+};{\varepsilon}^{-},{T}^{-},\sigma {'}^{-},\khi {'}_{\sigma}^{-})\\ \khi {'}_{\sigma}^{+}=\sigma {'}_{\sigma}^{+}({\varepsilon}^{+},{T}^{+};{\varepsilon}^{-},{T}^{-},\sigma {'}^{-},\khi {'}_{\sigma}^{-})\end{array}\) éq 6.3.1.2-1

\(\lbrace \begin{array}{}{\sigma}_{p}^{+}={\sigma}_{p}^{+}({p}_{1}^{+},{p}_{2}^{+};{p}_{1}^{-},{p}_{2}^{-},{\khi }_{H}^{-})\\ {\khi }_{H}^{+}={\khi }_{H}^{+}({p}_{1}^{+},{p}_{2}^{+};{p}_{1}^{-},{p}_{2}^{-},{\khi }_{H}^{-})\end{array}\) éq 6.3.1.2-2

Les dépendances indiquées par les équations [éq 6.3.1.2-1] et [éq 6.3.1.2-2] n’ont pas de justification théorique a priori. Il s’agit simplement de montrer les dépendances les plus générales possibles du point de vue de la programmation informatique. On remarque dans cette décomposition que la dépendance par rapport à la thermique a été laissée dans les contraintes effectives; typiquement, on pense que les lois sur les contraintes effectives s’écrivent comme en thermomécanique classique:

\(\sigma {'}^{+}=\sigma {'}^{+}({\varepsilon}^{+}-{\alpha}^{+}{T}^{+}I,{\varepsilon}^{-}-{\alpha}^{-}{T}^{-}I,\sigma {'}^{-},{\khi }_{\sigma}^{-})\)

Choix des contraintes

Du fait de l’utilisation assez fréquente de l’hypothèse des contraintes effectives, on décide que le vecteur des contraintes pour la partie mécanique contient dans tous les cas le tenseur des contraintes effectives \(\sigma '\) et le scalaire \({\sigma}_{p}\) . Dans le cas général où l’hypothèse des contraintes effectives n’est pas retenue, on aura simplement: \({\sigma}_{p}=0\) . C’est donc à la charge du module d’intégration des équations d’équilibre (et non pas des lois de comportement) de faire la somme: \({\sigma}^{+}=\sigma {'}^{+}+{\sigma}_{p}^{+}I\) .

Hydraulique

La loi de comportement hydraulique fournit les relations suivantes:

\(\lbrace \begin{array}{}\begin{array}{}{m}_{c}^{p+}={m}_{c}^{p+}({\varepsilon}^{+},{p}_{1}^{+},{p}_{2}^{+},{T}^{+};{\varepsilon}^{-},{p}_{1}^{-},{p}_{2}^{-},{T}^{-},{m}_{c}^{p-},{M}_{c}^{p-},{\chi }_{H}^{-})\\ {M}_{c}^{p+}={M}_{c}^{p+}(\begin{array}{}{\varepsilon}^{+},{p}_{1}^{+},\nabla {p}_{1}^{+},{p}_{2}^{+},\nabla {p}_{2}^{+},{T}^{+},\nabla {T}^{+};\\ {\varepsilon}^{-},{p}_{1}^{-},\nabla {p}_{1}^{-},{p}_{2}^{-},\nabla {p}_{2}^{-};\\ {T}^{-},\nabla {T}^{-},{M}_{c}^{p-},{\chi }_{H}^{-};{F}^{m+}\end{array})\end{array}\rbrace \begin{array}{}\forall c\\ \forall p\end{array}\\ {\chi }_{H}^{+}={\chi }_{H}^{+}({\varepsilon}^{+},{p}_{1}^{+},{p}_{2}^{+},{T}^{+};{\varepsilon}^{-},{p}_{1}^{-},{p}_{2}^{-},{T}^{-},{m}_{c}^{p-},{\chi }_{H}^{-})\end{array}\) éq 6.3.2-1

On remarque que le champ de gravité est une donnée de la loi de comportement hydraulique parce que l’évolution du vecteur de flux suit des relations du type:

\(M={\lambda}_{H}{\rho}^{\mathrm{fl}}\left[-\nabla P+{\rho}^{\mathrm{fl}}{F}^{m}\right]\)

Loi de comportement thermique

Les lois de comportement donnent:

\(\lbrace \begin{array}{}Q{'}^{+}=Q{'}^{+}({\varepsilon}^{+},{p}_{1}^{+},{p}_{2}^{+},{T}^{+};{\varepsilon}^{-},{p}_{1}^{-},{p}_{2}^{-},{T}^{-},Q{'}^{-})\\ {h}_{c}^{mp+}={h}_{c}^{mp+}({\varepsilon}^{+},{p}_{1}^{+},{p}_{2}^{+},{T}^{+};{\varepsilon}^{-},{p}_{1}^{-},{p}_{2}^{-},{T}^{-},{h}_{c}^{mp-})\forall c\text{et}\forall p\\ {q}^{+}={q}^{+}({\varepsilon}^{+},{p}_{1}^{+},{p}_{2}^{+},{T}^{+},\nabla {T}^{+};{\varepsilon}^{-},{p}_{1}^{-},{p}_{2}^{-},{T}^{-},\nabla {T}^{-},{q}^{-})\\ {\khi }_{T}^{+}={\khi }_{T}^{+}({\varepsilon}^{+},{p}_{1}^{+},{p}_{2}^{+},{T}^{+},\nabla {T}^{+};{\varepsilon}^{-},{p}_{1}^{-},{p}_{2}^{-},{T}^{-},\nabla {T}^{-},{\khi }_{T}^{-})\end{array}\) éq 6.3.3-1

Notons que nous avons introduit d’éventuelles variables internes liées à la thermique.

Masse volumique homogénéisée

Par définition, la masse volumique homogénéisée, qui intervient dans l’équation d’équilibre de la mécanique [éq 3.1-1] est donnée par:

\({r}^{+}={r}_{0}+{m}_{1}^{1+}+{m}_{1}^{2+}+{m}_{2}^{1+}+{m}_{2}^{2+}\) éq 6.3.4-1

Cette équation n’est pas une loi de comportement, mais elle fait partie des équations de conservation. Elle est intégrée dans le module de calcul des équations d’équilibre, les modules de calcul des lois de comportement n’ayant pas à la traiter.

Opérateur de masse

Cet opérateur est appelable par l’option MASS_MECApour tous les éléments concernés par les modélisations contenant de l’hydromécanique en 3D et en D_PLAN, soit: “THM”, “THHM”, “HHM”, “HH2M”, “THH2M”. Il permet ainsi de traiter une analyse dynamique. Il comprend uniquement des termes associés aux degrés de liberté en translation, soit respectivement DX, DY en 2D, et DX , DY , DZ en 3D.

Ces termes sont évalués dans chaque élément avant assemblage à partir des fonctions de forme des éléments isoparamétriques \(N\) décrites dans le document [R3.01.00] selon l’expression:

\({M}_{ij}^{e}=\underset{{\Omega}_{e}}{\int}\rho {N}_{i}.{N}_{j}d\Omega\)

L’expression précédente contient la masse volumique homogénéisée \(\rho\) initiale (correspondant à \({r}_{0}\) au paragraphe [§6.3.4]) entrée derrière l’opérande RHOdu mot-clé THM_DIFFU de la commande DEFI_MATERIAU. En effet, on juge pour les applications dynamiques courantes de courte durée (du type chargement sismique) que la modification de la masse volumique est suffisamment faible, cf. [bib15]. Ainsi on évite une mise à jour de la matrice de masse.

Les termes élémentaires \({M}_{ij}^{e}\) sont complétés lors de leur assemblage par des termes nuls associés aux degrés de libertés supplémentaires spécifiques à la modélisation couplée. Il en résulte lors de leur expansion dans la matrice de masse assemblée un décalage de degrés de liberté liés à la présence éventuelle des composantes “PRESS1”, “PRESS2”, et “TEMPE’dans la modélisation.

Algorithme de résolution

Algorithme non linéaire de résolution des équations d’équilibre

Dans le cas général de la modélisation (coefficients variables, désaturation, convection) le problème variationnel présenté ci-dessus [éq 4.1-1] à [éq 4.3-1] est non linéaire par rapport aux champs de déplacement, pression et température. Après discrétisation par éléments finis, on obtient un système matriciel non linéaire. La matrice de l’opérateur tangent contient de plus un terme non symétrique traité comme tel. On utilise dans tous les cas de modélisation le solveur non linéaire STAT_NON_LINE reposant sur une méthode de Newton-Raphson, décrite en [bib5]. On introduit la fonctionnelle vectorielle:

\(F(U)=R(U)-{L}^{\mathrm{meca}}\) éq 7.1-1

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

Pour les modules THM, objets de la présente note, l’opérateur \({L}^{\text{meca}}\) ne dépend pas des déplacements généralisés. Tous les termes dépendant des déplacements généralisés ont été introduits dans \(R\) , et c’est justement pour cette raison que les déplacements se retrouvent dans les déformations généralisées. Notons à ce sujet le traitement très particulier du terme \({\int}_{\Omega}r{F}^{m}.v\) de l’équation [éq 4.1-1].

D’après [éq 6.3.4-1], \({\int}_{\Omega}r{F}^{m}.vd\Omega ={\int}_{\Omega}({r}_{0}+{m}_{1}^{1+}+{m}_{1}^{2+}+{m}_{2}^{1+}+{m}_{2}^{2+}){F}^{m}.vd\Omega\) .

Nous avons choisi de scinder ce terme en deux:

Le terme \({\int}_{\Omega}{r}_{0}{F}^{m}.vd\Omega\) est une contribution à \({L}^{\mathrm{meca}}\) si l’utilisateur a renseigné l’opérande PESANTEUR du chargement utilisé (défini par la commande AFFE_CHAR_MECA), alors que le terme \({\int}_{\Omega}({m}_{1}^{1+}+{m}_{1}^{2+}+{m}_{2}^{1+}+{m}_{2}^{2+}){F}^{m}.vd\Omega\) , qui dépend des contraintes généralisées est une contribution à \(\mathrm{R}\) .

Passage des valeurs nodales aux valeurs aux points de Gauss

Comme dans tous les codes d’éléments finis, les termes sont calculés par boucle sur les éléments et boucle sur les points de gauss. En notant \({R}_{g}^{\mathrm{el}}\) et \(D{F}_{g}^{\mathrm{el}}\) les valeurs au point de Gauss \(g\) de l’élément \(\mathrm{el}\) des forces nodales et de l’opérateur tangent, et \({w}_{g}^{\mathrm{el}}\) le poids d’intégration lié à ce point de Gauss, on a:

\(R(U)=\sum_{\mathrm{el}}\sum_{g}{w}_{g}^{\mathrm{el}}{R}_{g}^{\mathrm{el}}(U)\) éq 7.2-1

\(DF(U)=\sum_{\mathrm{el}}\sum_{g}{w}_{g}^{\mathrm{el}}D{F}_{g}^{\mathrm{el}}(U)\) éq 7.2-2

Notons alors \({U}^{\mathrm{el}}\) le vecteur des inconnues nodales sur un élément fini \(\mathrm{el}\) . On peut ainsi avoir:

\(\text{par exemple}{U}^{\mathrm{el}}=\begin{array}{c}\begin{array}{c}u\\ v\\ w\\ {p}_{1}\\ {p}_{2}\\ T\end{array}\rbrace \text{noeud 1}\\ \begin{array}{c}u\\ v\\ w\\ {p}_{1}\\ {p}_{2}\\ T\end{array}\rbrace \text{noeud 2}\\ \begin{array}{c}u\\ v\\ w\\ {p}_{1}\\ {p}_{2}\\ T\end{array}\rbrace \text{noeud 3}\end{array}\) .

Notons aussi \({E}_{g}^{\mathrm{el}}\) le vecteur des déformations généralisées au point de Gauss \(g\) de l’élément \(\mathrm{el}\) et \({\Sigma}_{g}^{\mathrm{el}}\) le vecteur de contraintes généralisées pour le point de Gauss \(g\) de l’élément \(\mathrm{el}\) . Dans le cas le plus complet on a ainsi:

\({E}_{g}^{\mathrm{el}}={\left\lbrace \begin{array}{c}u\\ \varepsilon (u)\\ {p}_{1}\\ \nabla {p}_{1}\\ {p}_{2}\\ \nabla {p}_{2}\\ T\\ \nabla T\end{array}\right\rbrace }_{g}^{\mathrm{el}};{\Sigma}_{g}^{\mathrm{el}}={\left\lbrace \begin{array}{c}\sigma '\\ {\sigma}_{p}\\ {m}_{1}^{1}\\ {M}_{1}^{1}\\ {h}^{{m}_{1}^{1}}\\ {m}_{1}^{2}\\ {M}_{1}^{2}\\ {h}^{{m}_{1}^{2}}\\ {m}_{2}^{1}\\ {M}_{2}^{1}\\ {h}^{{m}_{2}^{1}}\\ {m}_{2}^{2}\\ {M}_{2}^{2}\\ {h}^{{m}_{2}^{2}}\\ Q'\\ q\end{array}\right\rbrace }_{g}^{\mathrm{el}}\)

Les fonctions de forme des éléments finis permettent de calculer alors la matrice \({Q}_{g}^{\mathrm{el}}\) de passage des inconnues nodales aux déformations généralisées aux points de Gauss définie par:

\({E}_{g}^{\mathrm{el}}={Q}_{g}^{\mathrm{el}}.{U}_{g}^{\mathrm{el}}\) éq 7.2-3

Vecteurs et matrices selon les options

Les présentations des deux paragraphes suivants sont faites dans le cas le plus général où on a une équation de mécanique, deux équations d’hydraulique et une équation de thermique. Les indices \(g\) et \(\mathrm{el}\) sont désormais omis, mais il est clair que ce qui est décrit s’applique à chaque point de Gauss de chaque élément.

Résidu ou force nodale : options RAPH_MECA et FULL_MECA

On répartit les termes de la formulation variationnelle selon le principe suivant:

Si \({E}_{g}^{\text{*}\mathrm{el}}\) désigne un champ de déformation virtuel, \({E}_{g}^{\text{*}\mathrm{el}}=(v,\varepsilon (v),{\pi}_{1},\nabla {\pi}_{1},{\pi}_{2},\nabla {\pi}_{2},\tau ,\nabla \tau )\) calculé à partir d’un vecteur de déplacements nodaux virtuels \({U}^{\text{*}\mathrm{el}}\) , on peut définir: \({E}_{g}^{\text{*}\mathrm{el}T}.{\overline{\Sigma}}_{g}^{\mathrm{el}}(U)={\overline{\Sigma}}_{1}v+{\overline{\Sigma}}_{2}(v)+{\overline{\Sigma}}_{3}{\pi}_{1}+{\overline{\Sigma}}_{4}\nabla {\pi}_{1}+{\overline{\Sigma}}_{5}{\pi}_{2}+{\overline{\Sigma}}_{6}\nabla {\pi}_{2}+{\overline{\Sigma}}_{7}\tau +{\overline{\Sigma}}_{8}\nabla \tau\) . On reprend alors les formulations variationnelles discrètes [éq 5.1-1], [éq 5.2-1], [éq 5.3-1], et on y remplace les intégrales \({\int}_{\Omega}fd\Omega\) par \(\sum_{\mathrm{el}}\sum_{g}{w}_{g}^{\mathrm{el}}{f}_{g}^{\mathrm{el}}\) pour toutes les intégrantes \(f\) . On distingue les termes multipliant respectivement \(v\) , \(\varepsilon (v)\) , \({\pi}_{1}\) , \(\nabla {\pi}_{1}\) , \({\pi}_{2}\) , \(\nabla {\pi}_{2}\) , \(\tau\) et \(\nabla \tau\) , et on trouve:

Indice	\(\overline{\Sigma}\)	associé à
1	\(-({m}_{1}^{1+}+{m}_{1}^{2+}+{m}_{2}^{1+}+{m}_{2}^{2+}){F}^{m+}\)	\(v\)
2	\(\sigma {'}^{+}+{\sigma}_{p}^{+}I\)	\(\varepsilon (v)\)
3	\(-{m}_{1}^{1+}-{m}_{1}^{2+}+{m}_{1}^{1-}+{m}_{1}^{2-}\)	\({\pi}_{1}\)
4	\(\theta \Delta t ({M}_{1}^{1+}+{M}_{1}^{2+})+(1-\theta )\Delta t ({M}_{1}^{1-}+{M}_{1}^{2-})\)	\(\nabla {\pi}_{1}\)
5	\(-{m}_{2}^{1+}-{m}_{2}^{2+}+{m}_{2}^{1-}+{m}_{2}^{2-}\)	\({\pi}_{2}\)
6	\(\theta \Delta t ({M}_{2}^{1+}+{M}_{2}^{2+})+(1-\theta )\Delta t ({M}_{2}^{1-}+{M}_{2}^{2-})\)	\(\nabla {\pi}_{2}\)
7	\(\begin{array}{}-Q{'}^{+}+Q{'}^{-}\\ -(\theta {h}_{1}^{m1+}+(1-\theta ){h}_{1}^{m1-})({m}_{1}^{1+}-{m}_{1}^{1-})-(\theta {h}_{1}^{m2+}+(1-\theta ){h}_{1}^{m2-})({m}_{1}^{2+}-{m}_{1}^{2-})\\ -(\theta {h}_{2}^{m1+}+(1-\theta ){h}_{2}^{m1-})({m}_{2}^{1+}-{m}_{2}^{1-})-(\theta {h}_{2}^{m2+}+(1-\theta ){h}_{2}^{m2-})({m}_{2}^{2+}-{m}_{2}^{2-})\\ +\Delta t\theta ({M}_{1}^{1+}+{M}_{1}^{2+}+{M}_{2}^{1+}+{M}_{2}^{2+}).{F}^{m}+\Delta t(1-\theta )({M}_{1}^{1-}+{M}_{1}^{2-}+{M}_{2}^{1-}+{M}_{2}^{2-}).{F}^{m}\end{array}\)	\(\tau\)
8	\(\begin{array}{}+\theta \Delta t ({h}_{1}^{m1+}{M}_{1}^{1+}+{h}_{1}^{m2+}{M}_{1}^{2+}+{h}_{2}^{m1+}{M}_{2}^{1+}+{h}_{2}^{m2+}{M}_{2}^{2+}+{q}^{+})+\\ +(1-\theta )\Delta t ({h}_{1}^{m1-}{M}_{1}^{1-}+{h}_{1}^{m2-}{M}_{1}^{2-}+{h}_{2}^{m1-}{M}_{2}^{1-}+{h}_{2}^{m2-}{M}_{2}^{2-}+{q}^{-})\end{array}\)	\(\nabla \tau\)

Remarque:

Dans le premier terme \({\stackrel{ˉ}{\Sigma}}_{1}\) ne figure pas le terme \(-{r}_{0}{F}^{m}\) car il est mis dans le chargement extérieur \({L}^{\mathrm{meca}}\) et calculé par l’option de calcul du chargement extérieur de pesanteur.

En utilisant la définition [éq 7.2-1] de \({R}_{g}^{\mathrm{el}}\) , on a:

\({U}^{\text{*}{\mathrm{el}}^{T}}.{R}_{g}^{\mathrm{el}}={E}_{g}^{\text{*}{\mathrm{el}}^{T}}.{\overline{\Sigma}}_{g}^{\mathrm{el}}\) , ce qui nous donne encore:

\({R}_{g}^{\mathrm{el}}={Q}_{g}^{{\mathrm{el}}^{T}}.{\overline{\Sigma}}_{g}^{\mathrm{el}}\)

Cette dernière égalité n’est que la forme locale au niveau d’un point de Gauss de [éq 6.2-2].

Opérateur tangent : options FULL_MECA, RIGI_MECA_TANG

Dans ce qui suit, si \(X\) désigne un vecteur de composantes \({X}^{i}\) et \(Y\) un vecteur de composantes \({Y}^{j}\) , \(\left[\frac{\partial X}{\partial Y}\right]\) désignera une matrice dont l’élément occupant la ligne \(i\) et la colonne \(j\) est \(\frac{\partial {X}^{i}}{\partial {Y}^{j}}\) .

Pour calculer l’opérateur tangent, on calculera les quantités suivantes:

\(\left[\text{DRDE}\right]=\)

DR1U	DR1E	DR1P1	DR1GP1	DR1P2	DR1GP2	DR1T	DR1GT
DR2U	DR2E	DR2P1	DR2GP1	DR2P2	DR2GP2	DR2T	DR2GT
DR3U	DR3E	DR3P1	DR3GP1	DR3P2	DR3GP2	DR3T	DR3GT
DR4U	DR4E	DR4P1	DR4GP1	DR4P2	DR4GP2	DR4T	DR4GT
DR5U	DR5E	DR5P1	DR5GP1	DR5P2	DR5GP2	DR5T	DR5GT
DR6U	DR6E	DR6P1	DR6GP1	DR6P2	DR6GP2	DR6T	DR6GT
DR7U	DR7E	DR7P1	DR7GP1	DR7P2	DR7GP2	DR7T	DR7GT
DR8U	DR8E	DR8P1	DR8GP1	DR8P2	DR8GP2	DR8T	DR8GT

Où l’on a noté:

\(\begin{array}{cc}\text{DRiU}=\frac{\partial {F}_{i}}{\partial u}& \text{DRiGP}1=\frac{\partial {F}_{i}}{\partial \nabla {p}_{1}}\\ \text{DRiE}=\frac{\partial {F}_{i}}{\partial \varepsilon }& \text{DRiGP}2=\frac{\partial {F}_{i}}{\partial \nabla {p}_{2}}\\ \text{DRiP}1=\frac{\partial {F}_{i}}{\partial {p}_{1}}& \text{DRiT}=\frac{\partial {F}_{i}}{\partial T}\\ \text{DRiP}2=\frac{\partial {F}_{i}}{\partial {p}_{2}}& \text{DRiDT}=\frac{\partial {F}_{i}}{\partial \nabla T}\end{array}\)

Pour faire ces calculs on considère que les lois de comportement fournissent, pour les options correspondantes, toutes les dérivées suivantes:

\(D\Sigma \mathrm{DE}=\left[\begin{array}{cccccccc}\frac{\partial \sigma '}{u}& \frac{\partial \sigma '}{\partial \varepsilon }& \frac{\partial \sigma '}{\partial {p}_{1}}& \frac{\partial \sigma '}{\partial \nabla {p}_{1}}& \frac{\partial \sigma '}{\partial {p}_{2}}& \frac{\partial \sigma '}{\partial \nabla {p}_{2}}& \frac{\partial \sigma '}{\partial T}& \frac{\partial \sigma '}{\partial \nabla T}\\ \frac{\partial {\sigma}_{p}}{\partial u}& \frac{\partial {\sigma}_{p}}{\partial \varepsilon }& \frac{\partial {\sigma}_{p}}{\partial {p}_{1}}& \frac{\partial {\sigma}_{p}}{\partial \nabla {p}_{1}}& \frac{\partial {\sigma}_{p}}{\partial {p}_{2}}& \frac{\partial {\sigma}_{p}}{\partial \nabla {p}_{2}}& \frac{\partial {\sigma}_{p}}{\partial T}& \frac{\partial {\sigma}_{p}}{\partial \nabla T}\\ \frac{\partial {m}_{1}^{1}}{\partial u}& \frac{\partial {m}_{1}^{1}}{\partial \varepsilon }& \frac{\partial {m}_{1}^{1}}{\partial {p}_{1}}& \frac{\partial {m}_{1}^{1}}{\partial \nabla {p}_{1}}& \frac{\partial {m}_{1}^{1}}{\partial {p}_{2}}& \frac{\partial {m}_{1}^{1}}{\partial \nabla {p}_{2}}& \frac{\partial {m}_{1}^{1}}{\partial T}& \frac{\partial {m}_{1}^{1}}{\partial \nabla T}\\ \frac{\partial {M}_{1}^{1}}{\partial u}& \frac{\partial {M}_{1}^{1}}{\partial \varepsilon }& \frac{\partial {M}_{1}^{1}}{\partial {p}_{1}}& \frac{\partial {M}_{1}^{1}}{\partial \nabla {p}_{1}}& \frac{\partial {M}_{1}^{1}}{\partial {p}_{2}}& \frac{\partial {M}_{1}^{1}}{\partial \nabla {p}_{2}}& \frac{\partial {M}_{1}^{1}}{\partial T}& \frac{\partial {M}_{1}^{1}}{\partial \nabla T}\\ \frac{\partial {h}^{{m}_{1}^{1}}}{\partial u}& \frac{\partial {h}^{{m}_{1}^{1}}}{\partial \varepsilon }& \frac{\partial {h}^{{m}_{1}^{1}}}{\partial {p}_{1}}& \frac{\partial {h}^{{m}_{1}^{1}}}{\partial \nabla {p}_{1}}& \frac{\partial {h}^{{m}_{1}^{1}}}{\partial {p}_{2}}& \frac{\partial {h}^{{m}_{1}^{1}}}{\partial \nabla {p}_{2}}& \frac{\partial {h}^{{m}_{1}^{1}}}{\partial T}& \frac{\partial {h}^{{m}_{1}^{1}}}{\partial \nabla T}\\ \frac{\partial {m}_{1}^{2}}{\partial u}& \frac{\partial {m}_{1}^{2}}{\partial \varepsilon }& \frac{\partial {m}_{1}^{2}}{\partial {p}_{1}}& \frac{\partial {m}_{1}^{2}}{\partial \nabla {p}_{1}}& \frac{\partial {m}_{1}^{2}}{\partial {p}_{2}}& \frac{\partial {m}_{1}^{2}}{\partial \nabla {p}_{2}}& \frac{\partial {m}_{1}^{2}}{\partial T}& \frac{\partial {m}_{1}^{2}}{\partial \nabla T}\\ \frac{\partial {M}_{1}^{2}}{\partial u}& \frac{\partial {M}_{1}^{2}}{\partial \varepsilon }& \frac{\partial {M}_{1}^{2}}{\partial {p}_{1}}& \frac{\partial {M}_{1}^{2}}{\partial \nabla {p}_{1}}& \frac{\partial {M}_{1}^{2}}{\partial {p}_{2}}& \frac{\partial {M}_{1}^{2}}{\partial \nabla {p}_{2}}& \frac{\partial {M}_{1}^{2}}{\partial T}& \frac{\partial {M}_{1}^{2}}{\partial \nabla T}\\ \frac{\partial {h}^{{m}_{1}^{2}}}{\partial u}& \frac{\partial {h}^{{m}_{1}^{2}}}{\partial \varepsilon }& \frac{\partial {h}^{{m}_{1}^{2}}}{\partial {p}_{1}}& \frac{\partial {h}^{{m}_{1}^{2}}}{\partial \nabla {p}_{1}}& \frac{\partial {h}^{{m}_{1}^{2}}}{\partial {p}_{2}}& \frac{\partial {h}^{{m}_{1}^{2}}}{\partial \nabla {p}_{2}}& \frac{\partial {h}^{{m}_{1}^{2}}}{\partial T}& \frac{\partial {h}^{{m}_{1}^{2}}}{\partial \nabla T}\\ \frac{\partial {m}_{2}^{1}}{\partial u}& \frac{\partial {m}_{2}^{1}}{\partial \varepsilon }& \frac{\partial {m}_{2}^{1}}{\partial {p}_{1}}& \frac{\partial {m}_{2}^{1}}{\partial \nabla {p}_{1}}& \frac{\partial {m}_{2}^{1}}{\partial {p}_{2}}& \frac{\partial {m}_{2}^{1}}{\partial \nabla {p}_{2}}& \frac{\partial {m}_{2}^{1}}{\partial T}& \frac{\partial {m}_{2}^{1}}{\partial \nabla T}\\ \frac{\partial {M}_{2}^{1}}{\partial u}& \frac{\partial {M}_{2}^{1}}{\partial \varepsilon }& \frac{\partial {M}_{2}^{1}}{\partial {p}_{1}}& \frac{\partial {M}_{2}^{1}}{\partial \nabla {p}_{1}}& \frac{\partial {M}_{2}^{1}}{\partial {p}_{2}}& \frac{\partial {M}_{2}^{1}}{\partial \nabla {p}_{2}}& \frac{\partial {M}_{2}^{1}}{\partial T}& \frac{\partial {M}_{2}^{1}}{\partial \nabla T}\\ \frac{\partial {h}^{{m}_{2}^{1}}}{\partial u}& \frac{\partial {h}^{{m}_{2}^{1}}}{\partial \varepsilon }& \frac{\partial {h}^{{m}_{2}^{1}}}{\partial {p}_{1}}& \frac{\partial {h}^{{m}_{2}^{1}}}{\partial \nabla {p}_{1}}& \frac{\partial {h}^{{m}_{2}^{1}}}{\partial {p}_{2}}& \frac{\partial {h}^{{m}_{2}^{1}}}{\partial \nabla {p}_{2}}& \frac{\partial {h}^{{m}_{2}^{1}}}{\partial T}& \frac{\partial {h}^{{m}_{2}^{1}}}{\partial \nabla T}\\ \frac{\partial {m}_{2}^{2}}{\partial u}& \frac{\partial {m}_{2}^{2}}{\partial \varepsilon }& \frac{\partial {m}_{2}^{2}}{\partial {p}_{1}}& \frac{\partial {m}_{2}^{2}}{\partial \nabla {p}_{1}}& \frac{\partial {m}_{2}^{2}}{\partial {p}_{2}}& \frac{\partial {m}_{2}^{2}}{\partial \nabla {p}_{2}}& \frac{\partial {m}_{2}^{2}}{\partial T}& \frac{\partial {m}_{2}^{2}}{\partial \nabla T}\\ \frac{\partial {M}_{2}^{2}}{\partial u}& \frac{\partial {M}_{2}^{2}}{\partial \varepsilon }& \frac{\partial {M}_{2}^{2}}{\partial {p}_{1}}& \frac{\partial {M}_{2}^{2}}{\partial \nabla {p}_{1}}& \frac{\partial {M}_{2}^{2}}{\partial {p}_{2}}& \frac{\partial {M}_{2}^{2}}{\partial \nabla {p}_{2}}& \frac{\partial {M}_{2}^{2}}{\partial T}& \frac{\partial {M}_{2}^{2}}{\partial \nabla T}\\ \frac{\partial {h}^{{m}_{2}^{2}}}{\partial u}& \frac{\partial {h}^{{m}_{2}^{2}}}{\partial \varepsilon }& \frac{\partial {h}^{{m}_{2}^{2}}}{\partial {p}_{1}}& \frac{\partial {h}^{{m}_{2}^{2}}}{\partial \nabla {p}_{1}}& \frac{\partial {h}^{{m}_{2}^{2}}}{\partial {p}_{2}}& \frac{\partial {h}^{{m}_{2}^{2}}}{\partial \nabla {p}_{2}}& \frac{\partial {h}^{{m}_{2}^{2}}}{\partial T}& \frac{\partial {h}^{{m}_{2}^{2}}}{\partial \nabla T}\\ \frac{\partial Q'}{\partial u}& \frac{\partial Q'}{\partial \varepsilon }& \frac{\partial Q'}{\partial {p}_{1}}& \frac{\partial Q'}{\partial \nabla {p}_{1}}& \frac{\partial Q'}{\partial {p}_{2}}& \frac{\partial Q'}{\partial \nabla {p}_{2}}& \frac{\partial Q'}{\partial T}& \frac{\partial Q'}{\partial \nabla T}\\ \frac{\partial q}{\partial u}& \frac{\partial q}{\partial \varepsilon }& \frac{\partial q}{\partial {p}_{1}}& \frac{\partial q}{\partial \nabla {p}_{1}}& \frac{\partial q}{\partial {p}_{2}}& \frac{\partial q}{\partial \nabla {p}_{2}}& \frac{\partial q}{\partial T}& \frac{\partial q}{\partial \nabla T}\end{array}\right]\)

Remarque:

Dans ces expressions, les dérivées par rapport à \(u\) sont toutes nulles, mais nous gardons l’écriture compte tenu de la définition des matrices \({Q}_{g}^{\text{el}}\) que nous avons adoptée.

L’appel aux lois de comportement fournira les morceaux de la matrice \(D\Sigma \text{DE}\) selon les équations présentes:

\(\left[\text{DMECDE}\right]=\left[\begin{array}{c}\frac{\partial \sigma '}{\partial \varepsilon }\\ \frac{\partial {\sigma}_{p}}{\partial \varepsilon }\end{array}\right];\left[\text{DMECP1}\right]=\left[\begin{array}{cc}\frac{\partial \sigma '}{\partial {p}_{1}}& \frac{\partial \sigma '}{\partial \nabla {p}_{1}}\\ \frac{\partial {\sigma}_{p}}{\partial {p}_{1}}& \frac{\partial {\sigma}_{p}}{\partial \nabla {p}_{1}}\end{array}\right];\left[\text{DMECP}2\right]=\left[\begin{array}{cc}\frac{\partial \sigma '}{\partial {p}_{2}}& \frac{\partial \sigma '}{\partial \nabla {p}_{2}}\\ \frac{\partial {\sigma}_{p}}{\partial {p}_{2}}& \frac{\partial {\sigma}_{p}}{\partial \nabla {p}_{2}}\end{array}\right]\left[\text{DMECDT}\right]=\left[\begin{array}{cc}\frac{\partial \sigma }{\partial T}& \frac{\partial \sigma }{\partial \nabla T}\\ \frac{\partial {\sigma}_{p}}{\partial T}& \frac{\partial {\sigma}_{p}}{\partial \nabla T}\end{array}\right]\)

\(\left[\text{DP11DE}\right]=\left[\begin{array}{c}\frac{\partial {m}_{1}^{1}}{\partial \varepsilon }\\ \frac{\partial {M}_{1}^{1}}{\partial \varepsilon }\\ \frac{\partial {h}^{{m}_{1}^{1}}}{\partial \varepsilon }\end{array}\right];\left[\text{DP11P1}\right]=\left[\begin{array}{cc}\frac{\partial {m}_{1}^{1}}{\partial {p}_{1}}& \frac{\partial {m}_{1}^{1}}{\partial \nabla {p}_{1}}\\ \frac{\partial {M}_{1}^{1}}{\partial {p}_{1}}& \frac{\partial {M}_{1}^{1}}{\partial \nabla {p}_{1}}\\ \frac{\partial {h}^{{m}_{1}^{1}}}{\partial {p}_{1}}& \frac{\partial {h}^{{m}_{1}^{1}}}{\partial \nabla {p}_{1}}\end{array}\right];\left[\text{DP11P2}\right]=\left[\begin{array}{cc}\frac{\partial {m}_{1}^{1}}{\partial {p}_{2}}& \frac{\partial {m}_{1}^{1}}{\partial \nabla {p}_{2}}\\ \frac{\partial {M}_{1}^{1}}{\partial {p}_{2}}& \frac{\partial {M}_{1}^{1}}{\partial \nabla {p}_{2}}\\ \frac{\partial {h}^{{m}_{1}^{1}}}{\partial {p}_{2}}& \frac{\partial {h}^{{m}_{1}^{1}}}{\partial \nabla {p}_{2}}\end{array}\right]\left[\text{DP11DT}\right]=\left[\begin{array}{cc}\frac{\partial {m}_{1}^{1}}{\partial T}& \frac{\partial {m}_{1}^{1}}{\partial \nabla T}\\ \frac{\partial {M}_{1}^{1}}{\partial T}& \frac{\partial {M}_{1}^{1}}{\partial \nabla T}\\ \frac{\partial {h}^{{m}_{1}^{1}}}{\partial T}& \frac{\partial {h}^{{m}_{1}^{1}}}{\partial \nabla T}\end{array}\right]\)

\(\left[\text{DP12DE}\right]=\left[\begin{array}{c}\frac{\partial {m}_{1}^{2}}{\partial \varepsilon }\\ \frac{\partial {M}_{1}^{2}}{\partial \varepsilon }\\ \frac{\partial {h}^{{m}_{1}^{2}}}{\partial \varepsilon }\end{array}\right];\left[\text{DP12P1}\right]=\left[\begin{array}{cc}\frac{\partial {m}_{1}^{2}}{\partial {p}_{1}}& \frac{\partial {m}_{1}^{2}}{\partial \nabla {p}_{1}}\\ \frac{\partial {M}_{1}^{2}}{\partial {p}_{1}}& \frac{\partial {M}_{1}^{2}}{\partial \nabla {p}_{1}}\\ \frac{\partial {h}^{{m}_{1}^{2}}}{\partial {p}_{1}}& \frac{\partial {h}^{{m}_{1}^{2}}}{\partial \nabla {p}_{1}}\end{array}\right];\left[\text{DP12P2}\right]=\left[\begin{array}{cc}\frac{\partial {m}_{1}^{2}}{\partial {p}_{2}}& \frac{\partial {m}_{1}^{2}}{\partial \nabla {p}_{2}}\\ \frac{\partial {M}_{1}^{2}}{\partial {p}_{2}}& \frac{\partial {M}_{1}^{2}}{\partial \nabla {p}_{2}}\\ \frac{\partial {h}^{{m}_{1}^{2}}}{\partial {p}_{2}}& \frac{\partial {h}^{{m}_{1}^{2}}}{\partial \nabla {p}_{2}}\end{array}\right]\left[\text{DP12DT}\right]=\left[\begin{array}{cc}\frac{\partial {m}_{1}^{2}}{\partial T}& \frac{\partial {m}_{1}^{2}}{\partial \nabla T}\\ \frac{\partial {M}_{1}^{2}}{\partial T}& \frac{\partial {M}_{1}^{2}}{\partial \nabla T}\\ \frac{\partial {h}^{{m}_{1}^{2}}}{\partial T}& \frac{\partial {h}^{{m}_{1}^{2}}}{\partial \nabla T}\end{array}\right]\)

\(\left[\text{DP21DE}\right]=\left[\begin{array}{c}\frac{\partial {m}_{2}^{1}}{\partial \varepsilon }\\ \frac{\partial {M}_{2}^{1}}{\partial \varepsilon }\\ \frac{\partial {h}^{{m}_{2}^{1}}}{\partial \varepsilon }\end{array}\right];\left[\text{DP21P1}\right]=\left[\begin{array}{cc}\frac{\partial {m}_{2}^{1}}{\partial {p}_{1}}& \frac{\partial {m}_{2}^{1}}{\partial \nabla {p}_{1}}\\ \frac{\partial {M}_{2}^{1}}{\partial {p}_{1}}& \frac{\partial {M}_{2}^{1}}{\partial \nabla {p}_{1}}\\ \frac{\partial {h}^{{m}_{2}^{1}}}{\partial {p}_{1}}& \frac{\partial {h}^{{m}_{2}^{1}}}{\partial \nabla {p}_{1}}\end{array}\right];\left[\text{DP21P2}\right]=\left[\begin{array}{cc}\frac{\partial {m}_{2}^{1}}{\partial {p}_{2}}& \frac{\partial {m}_{2}^{1}}{\partial \nabla {p}_{2}}\\ \frac{\partial {M}_{2}^{1}}{\partial {p}_{2}}& \frac{\partial {M}_{2}^{1}}{\partial \nabla {p}_{2}}\\ \frac{\partial {h}^{{m}_{2}^{1}}}{\partial {p}_{2}}& \frac{\partial {h}^{{m}_{2}^{1}}}{\partial \nabla {p}_{2}}\end{array}\right]\left[\text{DP21DT}\right]=\left[\begin{array}{cc}\frac{\partial {m}_{2}^{1}}{\partial T}& \frac{\partial {m}_{2}^{1}}{\partial \nabla T}\\ \frac{\partial {M}_{2}^{1}}{\partial T}& \frac{\partial {M}_{2}^{1}}{\partial \nabla T}\\ \frac{\partial {h}^{{m}_{2}^{1}}}{\partial T}& \frac{\partial {h}^{{m}_{2}^{1}}}{\partial \nabla T}\end{array}\right]\)

\(\left[\text{DP22DE}\right]=\left[\begin{array}{c}\frac{\partial {m}_{2}^{2}}{\partial \varepsilon }\\ \frac{\partial {M}_{2}^{2}}{\partial \varepsilon }\\ \frac{\partial {h}^{{m}_{2}^{2}}}{\partial \varepsilon }\end{array}\right];\left[\text{DP22P1}\right]=\left[\begin{array}{cc}\frac{\partial {m}_{2}^{2}}{\partial {p}_{1}}& \frac{\partial {m}_{2}^{2}}{\partial \nabla {p}_{1}}\\ \frac{\partial {M}_{2}^{2}}{\partial {p}_{1}}& \frac{\partial {M}_{2}^{2}}{\partial \nabla {p}_{1}}\\ \frac{\partial {h}^{{m}_{2}^{2}}}{\partial {p}_{1}}& \frac{\partial {h}^{{m}_{2}^{2}}}{\partial \nabla {p}_{1}}\end{array}\right];\left[\text{DP22P2}\right]=\left[\begin{array}{cc}\frac{\partial {m}_{2}^{2}}{\partial {p}_{2}}& \frac{\partial {m}_{2}^{2}}{\partial \nabla {p}_{2}}\\ \frac{\partial {M}_{2}^{2}}{\partial {p}_{2}}& \frac{\partial {M}_{2}^{2}}{\partial \nabla {p}_{2}}\\ \frac{\partial {h}^{{m}_{2}^{2}}}{\partial {p}_{2}}& \frac{\partial {h}^{{m}_{2}^{2}}}{\partial \nabla {p}_{2}}\end{array}\right]\left[\text{DP22DT}\right]=\left[\begin{array}{cc}\frac{\partial {m}_{2}^{2}}{\partial T}& \frac{\partial {m}_{2}^{2}}{\partial \nabla T}\\ \frac{\partial {M}_{2}^{2}}{\partial T}& \frac{\partial {M}_{2}^{2}}{\partial \nabla T}\\ \frac{\partial {h}^{{m}_{2}^{2}}}{\partial T}& \frac{\partial {h}^{{m}_{2}^{2}}}{\partial \nabla T}\end{array}\right]\)

\(\left[\text{DTDE}\right]=\left[\begin{array}{c}\frac{\partial Q'}{\partial \varepsilon }\\ \frac{\partial q}{\partial \varepsilon }\end{array}\right];\left[\text{DTDP1}\right]=\left[\begin{array}{cc}\frac{\partial Q'}{\partial {p}_{1}}& \frac{\partial Q'}{\partial \nabla {p}_{1}}\\ \frac{\partial q}{\partial {p}_{1}}& \frac{\partial q}{\partial \nabla {p}_{1}}\end{array}\right];\left[\text{DTDP2}\right]=\left[\begin{array}{cc}\frac{\partial Q'}{\partial {p}_{2}}& \frac{\partial Q'}{\partial \nabla {p}_{2}}\\ \frac{\partial q}{\partial {p}_{2}}& \frac{\partial q}{\partial \nabla {p}_{2}}\end{array}\right]\left[\text{DTDT}\right]=\left[\begin{array}{cc}\frac{\partial Q'}{\partial T}& \frac{\partial Q'}{\partial T}\\ \frac{\partial q}{\partial T}& \frac{\partial q}{\partial T}\end{array}\right]\)

Par ailleurs, en dérivant l’expression du résidu par rapport aux contraintes, on définit:

\(D\overline{\Sigma}D\Sigma =\left[\begin{array}{cccccccccccccccc}\frac{\partial \overline{{\Sigma}_{1}}}{\partial \sigma '}& \frac{\partial \overline{{\Sigma}_{1}}}{\partial {\sigma}_{p}}& \frac{\partial \overline{{\Sigma}_{1}}}{\partial {m}_{1}^{1}}& \frac{\partial \overline{{\Sigma}_{1}}}{\partial {M}_{1}^{1}}& \frac{\partial \overline{{\Sigma}_{1}}}{\partial {h}^{{m}_{1}^{1}}}& \frac{\partial \overline{{\Sigma}_{1}}}{\partial {m}_{1}^{2}}& \frac{\partial \overline{{\Sigma}_{1}}}{\partial {M}_{1}^{2}}& \frac{\partial \overline{{\Sigma}_{1}}}{\partial {h}^{{m}_{1}^{2}}}& \frac{\partial \overline{{\Sigma}_{1}}}{\partial {m}_{2}^{1}}& \frac{\partial \overline{{\Sigma}_{1}}}{\partial {M}_{2}^{1}}& \frac{\partial \overline{{\Sigma}_{1}}}{\partial {h}^{{m}_{2}^{1}}}& \frac{\partial \overline{{\Sigma}_{1}}}{\partial {m}_{2}^{2}}& \frac{\partial \overline{{\Sigma}_{1}}}{\partial {M}_{2}^{2}}& \frac{\partial \overline{{\Sigma}_{1}}}{\partial {h}^{{m}_{2}^{2}}}& \frac{\partial \overline{{\Sigma}_{1}}}{\partial Q'}& \frac{\partial \overline{{\Sigma}_{1}}}{\partial q}\\ \frac{\partial \overline{{\Sigma}_{2}}}{\partial \sigma '}& \frac{\partial \overline{{\Sigma}_{2}}}{\partial {\sigma}_{p}}& \frac{\partial \overline{{\Sigma}_{2}}}{\partial {m}_{1}^{1}}& \frac{\partial \overline{{\Sigma}_{2}}}{\partial {M}_{1}^{1}}& \frac{\partial \overline{{\Sigma}_{2}}}{\partial {h}^{{m}_{1}^{1}}}& \frac{\partial \overline{{\Sigma}_{2}}}{\partial {m}_{1}^{2}}& \frac{\partial \overline{{\Sigma}_{2}}}{\partial {M}_{1}^{2}}& \frac{\partial \overline{{\Sigma}_{2}}}{\partial {h}^{{m}_{1}^{2}}}& \frac{\partial \overline{{\Sigma}_{2}}}{\partial {m}_{2}^{1}}& \frac{\partial \overline{{\Sigma}_{2}}}{\partial {M}_{2}^{1}}& \frac{\partial \overline{{\Sigma}_{2}}}{\partial {h}^{{m}_{2}^{1}}}& \frac{\partial \overline{{\Sigma}_{2}}}{\partial {m}_{2}^{2}}& \frac{\partial \overline{{\Sigma}_{2}}}{\partial {M}_{2}^{2}}& \frac{\partial \overline{{\Sigma}_{2}}}{\partial {h}^{{m}_{2}^{2}}}& \frac{\partial \overline{{\Sigma}_{2}}}{\partial Q'}& \frac{\partial \overline{{\Sigma}_{2}}}{\partial q}\\ \frac{\partial \overline{{\Sigma}_{3}}}{\partial \sigma '}& \frac{\partial \overline{{\Sigma}_{3}}}{\partial {\sigma}_{p}}& \frac{\partial \overline{{\Sigma}_{3}}}{\partial {m}_{1}^{1}}& \frac{\partial \overline{{\Sigma}_{3}}}{\partial {M}_{1}^{1}}& \frac{\partial \overline{{\Sigma}_{3}}}{\partial {h}^{{m}_{1}^{1}}}& \frac{\partial \overline{{\Sigma}_{3}}}{\partial {m}_{1}^{2}}& \frac{\partial \overline{{\Sigma}_{3}}}{\partial {M}_{1}^{2}}& \frac{\partial \overline{{\Sigma}_{3}}}{\partial {h}^{{m}_{1}^{2}}}& \frac{\partial \overline{{\Sigma}_{3}}}{\partial {m}_{2}^{1}}& \frac{\partial \overline{{\Sigma}_{3}}}{\partial {M}_{2}^{1}}& \frac{\partial \overline{{\Sigma}_{3}}}{\partial {h}^{{m}_{2}^{1}}}& \frac{\partial \overline{{\Sigma}_{3}}}{\partial {m}_{2}^{2}}& \frac{\partial \overline{{\Sigma}_{3}}}{\partial {M}_{2}^{2}}& \frac{\partial \overline{{\Sigma}_{3}}}{\partial {h}^{{m}_{2}^{2}}}& \frac{\partial \overline{{\Sigma}_{3}}}{\partial Q'}& \frac{\partial \overline{{\Sigma}_{3}}}{\partial q}\\ \frac{\partial \overline{{\Sigma}_{4}}}{\partial \sigma '}& \frac{\partial \overline{{\Sigma}_{4}}}{\partial {\sigma}_{p}}& \frac{\partial \overline{{\Sigma}_{4}}}{\partial {m}_{1}^{1}}& \frac{\partial \overline{{\Sigma}_{4}}}{\partial {M}_{1}^{1}}& \frac{\partial \overline{{\Sigma}_{4}}}{\partial {h}^{{m}_{1}^{1}}}& \frac{\partial \overline{{\Sigma}_{4}}}{\partial {m}_{1}^{2}}& \frac{\partial \overline{{\Sigma}_{4}}}{\partial {M}_{1}^{2}}& \frac{\partial \overline{{\Sigma}_{4}}}{\partial {h}^{{m}_{1}^{2}}}& \frac{\partial \overline{{\Sigma}_{4}}}{\partial {m}_{2}^{1}}& \frac{\partial \overline{{\Sigma}_{4}}}{\partial {M}_{2}^{1}}& \frac{\partial \overline{{\Sigma}_{4}}}{\partial {h}^{{m}_{2}^{1}}}& \frac{\partial \overline{{\Sigma}_{4}}}{\partial {m}_{2}^{2}}& \frac{\partial \overline{{\Sigma}_{4}}}{\partial {M}_{2}^{2}}& \frac{\partial \overline{{\Sigma}_{4}}}{\partial {h}^{{m}_{2}^{2}}}& \frac{\partial \overline{{\Sigma}_{4}}}{\partial Q'}& \frac{\partial \overline{{\Sigma}_{4}}}{\partial q}\\ \frac{\partial \overline{{\Sigma}_{5}}}{\partial \sigma '}& \frac{\partial \overline{{\Sigma}_{5}}}{\partial {\sigma}_{p}}& \frac{\partial \overline{{\Sigma}_{5}}}{\partial {m}_{1}^{1}}& \frac{\partial \overline{{\Sigma}_{5}}}{\partial {M}_{1}^{1}}& \frac{\partial \overline{{\Sigma}_{5}}}{\partial {h}^{{m}_{1}^{1}}}& \frac{\partial \overline{{\Sigma}_{5}}}{\partial {m}_{1}^{2}}& \frac{\partial \overline{{\Sigma}_{5}}}{\partial {M}_{1}^{2}}& \frac{\partial \overline{{\Sigma}_{5}}}{\partial {h}^{{m}_{1}^{2}}}& \frac{\partial \overline{{\Sigma}_{5}}}{\partial {m}_{2}^{1}}& \frac{\partial \overline{{\Sigma}_{5}}}{\partial {M}_{2}^{1}}& \frac{\partial \overline{{\Sigma}_{5}}}{\partial {h}^{{m}_{2}^{1}}}& \frac{\partial \overline{{\Sigma}_{5}}}{\partial {m}_{2}^{2}}& \frac{\partial \overline{{\Sigma}_{5}}}{\partial {M}_{2}^{2}}& \frac{\partial \overline{{\Sigma}_{5}}}{\partial {h}^{{m}_{2}^{2}}}& \frac{\partial \overline{{\Sigma}_{5}}}{\partial Q'}& \frac{\partial \overline{{\Sigma}_{5}}}{\partial q}\\ \frac{\partial \overline{{\Sigma}_{6}}}{\partial \sigma '}& \frac{\partial \overline{{\Sigma}_{6}}}{\partial {\sigma}_{p}}& \frac{\partial \overline{{\Sigma}_{6}}}{\partial {m}_{1}^{1}}& \frac{\partial \overline{{\Sigma}_{6}}}{\partial {M}_{1}^{1}}& \frac{\partial \overline{{\Sigma}_{6}}}{\partial {h}^{{m}_{1}^{1}}}& \frac{\partial \overline{{\Sigma}_{6}}}{\partial {m}_{1}^{2}}& \frac{\partial \overline{{\Sigma}_{6}}}{\partial {M}_{1}^{2}}& \frac{\partial \overline{{\Sigma}_{6}}}{\partial {h}^{{m}_{1}^{2}}}& \frac{\partial \overline{{\Sigma}_{6}}}{\partial {m}_{2}^{1}}& \frac{\partial \overline{{\Sigma}_{6}}}{\partial {M}_{2}^{1}}& \frac{\partial \overline{{\Sigma}_{6}}}{\partial {h}^{{m}_{2}^{1}}}& \frac{\partial \overline{{\Sigma}_{6}}}{\partial {m}_{2}^{2}}& \frac{\partial \overline{{\Sigma}_{6}}}{\partial {M}_{2}^{2}}& \frac{\partial \overline{{\Sigma}_{6}}}{\partial {h}^{{m}_{2}^{2}}}& \frac{\partial \overline{{\Sigma}_{6}}}{\partial Q'}& \frac{\partial \overline{{\Sigma}_{6}}}{\partial q}\\ \frac{\partial \overline{{\Sigma}_{7}}}{\partial \sigma '}& \frac{\partial \overline{{\Sigma}_{7}}}{\partial {\sigma}_{p}}& \frac{\partial \overline{{\Sigma}_{7}}}{\partial {m}_{1}^{1}}& \frac{\partial \overline{{\Sigma}_{7}}}{\partial {M}_{1}^{1}}& \frac{\partial \overline{{\Sigma}_{7}}}{\partial {h}^{{m}_{1}^{1}}}& \frac{\partial \overline{{\Sigma}_{7}}}{\partial {m}_{1}^{2}}& \frac{\partial \overline{{\Sigma}_{7}}}{\partial {M}_{1}^{2}}& \frac{\partial \overline{{\Sigma}_{7}}}{\partial {h}^{{m}_{1}^{2}}}& \frac{\partial \overline{{\Sigma}_{7}}}{\partial {m}_{2}^{1}}& \frac{\partial \overline{{\Sigma}_{7}}}{\partial {M}_{2}^{1}}& \frac{\partial \overline{{\Sigma}_{7}}}{\partial {h}^{{m}_{2}^{1}}}& \frac{\partial \overline{{\Sigma}_{7}}}{\partial {m}_{2}^{2}}& \frac{\partial \overline{{\Sigma}_{7}}}{\partial {M}_{2}^{2}}& \frac{\partial \overline{{\Sigma}_{7}}}{\partial {h}^{{m}_{2}^{2}}}& \frac{\partial \overline{{\Sigma}_{7}}}{\partial Q'}& \frac{\partial \overline{{\Sigma}_{7}}}{\partial q}\\ \frac{\partial \overline{{\Sigma}_{8}}}{\partial \sigma '}& \frac{\partial \overline{{\Sigma}_{8}}}{\partial {\sigma}_{p}}& \frac{\partial \overline{{\Sigma}_{8}}}{\partial {m}_{1}^{1}}& \frac{\partial \overline{{\Sigma}_{8}}}{\partial {M}_{1}^{1}}& \frac{\partial \overline{{\Sigma}_{8}}}{\partial {h}^{{m}_{1}^{1}}}& \frac{\partial \overline{{\Sigma}_{8}}}{\partial {m}_{1}^{2}}& \frac{\partial \overline{{\Sigma}_{8}}}{\partial {M}_{1}^{2}}& \frac{\partial \overline{{\Sigma}_{8}}}{\partial {h}^{{m}_{1}^{2}}}& \frac{\partial \overline{{\Sigma}_{8}}}{\partial {m}_{2}^{1}}& \frac{\partial \overline{{\Sigma}_{8}}}{\partial {M}_{2}^{1}}& \frac{\partial \overline{{\Sigma}_{8}}}{\partial {h}^{{m}_{2}^{1}}}& \frac{\partial \overline{{\Sigma}_{8}}}{\partial {m}_{2}^{2}}& \frac{\partial \overline{{\Sigma}_{8}}}{\partial {M}_{2}^{2}}& \frac{\partial \overline{{\Sigma}_{8}}}{\partial {h}^{{m}_{2}^{2}}}& \frac{\partial \overline{{\Sigma}_{8}}}{\partial Q'}& \frac{\partial \overline{{\Sigma}_{8}}}{\partial q}\end{array}\right]\)

Toutes ces quantités n’étant pas forcément calculées, on notera, pour \(i\) de 1 à 8:

\(\begin{array}{cc}\left[D\overline{\Sigma}iD\sigma \right]=\left[\frac{\partial {\overline{\Sigma}}_{i}}{\partial \sigma '},\frac{\partial {\overline{\Sigma}}_{i}}{\partial {\sigma}_{p}}\right]& \left[D\overline{\Sigma}i\text{DP}21\right]=\left[\frac{\partial {\overline{\Sigma}}_{i}}{\partial {m}_{2}^{1}},\frac{\partial {\overline{\Sigma}}_{i}}{\partial {M}_{2}^{1}},\frac{\partial {\overline{\Sigma}}_{i}}{\partial {h}^{{m}_{2}^{1}}}\right]\\ \left[D\overline{\Sigma}i\text{DP}11\right]=\left[\frac{\partial {\overline{\Sigma}}_{i}}{\partial {m}_{1}^{2}},\frac{\partial {\overline{\Sigma}}_{i}}{\partial {M}_{1}^{1}},\frac{\partial {\overline{\Sigma}}_{i}}{\partial {h}^{{m}_{1}^{2}}}\right]& \left[D\overline{\Sigma}i\text{DP}22\right]=\left[\frac{\partial {\overline{\Sigma}}_{i}}{\partial {m}_{2}^{2}},\frac{\partial {\overline{\Sigma}}_{i}}{\partial {M}_{2}^{2}},\frac{\partial {\overline{\Sigma}}_{i}}{\partial {h}^{{m}_{2}^{2}}}\right]\\ \left[D\overline{\Sigma}i\text{DP}12\right]=\left[\frac{\partial {\overline{\Sigma}}_{i}}{\partial {m}_{1}^{2}},\frac{\partial {\overline{\Sigma}}_{i}}{\partial {M}_{1}^{1}},\frac{\partial {\overline{\Sigma}}_{i}}{\partial {h}^{{m}_{1}^{2}}}\right]& \left[D\overline{\Sigma}i\text{DT}\right]=\left[\frac{\partial {\overline{\Sigma}}_{i}}{\partial Q'},\frac{\partial {\overline{\Sigma}}_{i}}{\partial q}\right]\end{array}\)

Il est alors clair que:

\(D\overline{\Sigma}\mathrm{DE}=D\overline{\Sigma}D\Sigma .D\Sigma \mathrm{DE}\)

Et la contribution du point de Gauss à la matrice tangente \(D{F}_{g}^{\text{el}}\) s’obtient par:

\(D{F}_{g}^{\text{el}}={Q}_{g}^{{\text{el}}^{T}}.D\overline{\Sigma}\text{DE}.{Q}_{g}^{\text{el}}\)

Algorithme global

L’algorithme devient alors:

Initialisations:

Calcul de \({L}^{{\text{meca}}^{+}}\) (option CHAR_MECA)

Calcul de \(D{F}^{-}\) (option RIGI_MECA_TANG)

Calcul de \(\Delta {U}_{0}\) par: \(D{F}^{-}.\Delta {U}_{0}={L}^{{\text{meca}}^{+}}-{L}^{{\text{meca}}^{-}}\)

Itérations d’équilibre de Newton

Boucle éléments el

Boucle points de Gauss g

calcul \({Q}_{g}^{\text{el}}\)

calcul \({E}_{g}^{{\text{el}}^{-}}={Q}_{g}^{\text{el}}.{U}^{{\text{el}}^{-}}\) et \({E}_{g}^{{\text{el}}^{+}}={Q}_{g}^{\text{el}}.{U}^{{\text{el}}^{+}}\)

calcul de: \({\Sigma}_{gn}^{{\text{el}}^{+}},{\alpha}_{g}^{{\text{el}}^{+}},\frac{\partial {\Sigma}_{gn}^{{\text{el}}^{+}}}{\partial {E}_{gn}^{{\text{el}}^{+}}}\) (selon option) à partir de \({E}_{g}^{{\text{el}}^{-}},{\Sigma}_{g}^{{\text{el}}^{-}},{\alpha}_{g}^{{\text{el}}^{-}},{E}_{gn}^{{\text{el}}^{+}}\)

calcul de \({\overline{\Sigma}}_{gn}^{{\text{el}}^{+}}\) à partir de \({\Sigma}_{gn}^{{\text{el}}^{+}}\) ; \({R}_{gn}^{{\text{el}}^{+}}={Q}_{g}^{{\text{el}}^{T}}.{\overline{\Sigma}}_{gn}^{{\text{el}}^{+}}\)

calcul de \(\frac{\partial {\overline{\Sigma}}_{gn}^{{\text{el}}^{+}}}{\partial {\Sigma}_{gn}^{{\text{el}}^{+}}}\) à partir de \({\Sigma}_{gn}^{{\text{el}}^{+}}\) ; \(D{F}_{gn}^{{\text{el}}^{+}}={Q}_{g}^{{\text{el}}^{T}}.\frac{\partial {\overline{\Sigma}}_{gn}^{{\text{el}}^{+}}}{\partial {\Sigma}_{gn}^{{\text{el}}^{+}}}.\frac{\partial {\Sigma}_{gn}^{{\text{el}}^{+}}}{\partial {E}_{gn}^{{\text{el}}^{+}}}.{Q}_{g}^{\text{el}}\) (selon option)

Calcul de \(\delta {U}_{n+1}\) par:

\(D{F}_{n}^{+}.\delta {U}_{n+1}=-{R}_{n}^{+}+{L}^{{\text{meca}}^{+}}\)

Actualisation :

\(\Delta {U}_{n+1}=\Delta {U}_{n}+\rho \delta {U}_{n+1}\)

Si test convergence OK

fin Newton: pas de temps suivant

Sinon

\(n=n+1\)

L’option FORC_NODA

Au niveau des équations continues, l’option FORC_NODA correspond au calcul de l’opérateur \(R={Q}^{T}\overline{\Sigma}\) .Au niveau discret, l’option FORC_NODA revient à calculer le vecteur \({R}_{g}^{\text{el}}={Q}_{g}^{{\text{el}}^{T}}.{\overline{\Sigma}}_{g}^{\text{el}}\) .

Comme nous avons déjà noté que \(\overline{\Sigma}\) dépend non seulement de \(\Sigma\) , mais aussi de \(\dot{\Sigma}\) , il ne faut pas s’étonner de voir apparaître le pas de temps \(\Delta t\) et les contraintes à la fois au temps + et au temps 

L’algorithme de Newton-Raphson de la commande STAT_NON_LINE utilise l’option FORC_NODA pour le calcul de la prédiction au début de chaque pas de temps. Il n’est donc pas anodin de calculer correctement tous les termes pour cette option, y compris ceux qui dépendent du pas de temps. Nous illustrons cette question par un exemple simple correspondant à la seule équation de l’hydraulique.

Soit une version simplifiée de l’équation hydraulique:

\(-{\int}_{\Omega}\frac{\text{dm}}{\text{dt}}.{p}^{\text{*}}d\Omega +{\int}_{\Omega}M.\nabla {p}^{\text{*}}d\Omega ={M}_{\text{ext}}.{p}^{\text{*}}\) éq 8-1

Après discrétisation en temps:

\(-{\int}_{\Omega}\Delta \text{m}.{p}^{\text{*}}d\Omega +\Delta t{\int}_{\Omega}(\theta {M}^{+}+(1-\theta ){M}^{-}).\nabla {p}^{\text{*}}d\Omega ={M}_{\text{ext}}^{\theta}.{p}^{\text{*}}\)

Faisant apparaître \(\Delta M={M}^{+}-{M}^{-}\) et \(\Delta p={p}^{+}-{p}^{-}\) , et écrivant une loi de comportement: \(\Delta m =N\Delta p\) , on trouve:

\(-{\int}_{\Omega}N\Delta \text{p}.{p}^{\text{*}}d\Omega +\Delta t{\int}_{\Omega}\Delta M.\nabla {p}^{\text{*}}d\Omega ={M}_{\text{ext}}^{\theta}{p}^{\text{*}}-\Delta t{\int}_{\Omega}{M}^{-}.\nabla {p}^{\text{*}}d\Omega\)

Par définition la phase de prédiction STAT_NON_LINE s’écrit:

\({K}_{0}\Delta {u}^{0}={F}_{\text{ext}}^{1}-{Q}^{T}{\sigma}_{0}\)

Il est alors clair que l’on doit prendre \(\underset{\Omega}{\int}{Q}^{T}{\sigma}_{0}{p}^{\text{*}}d\Omega =-\Delta t \underset{\Omega}{\int}{M}^{-}\nabla {p}^{\text{*}}d\Omega\)

Discrétisation spatiale

Les éléments finis THM sont des éléments mixtes, dans le sens qu’ils ont à la fois des inconnues de déplacements, de pressions et de températures. Un choix de discrétisation où les déplacements, les pressions et les températures sont interpolés avec le même ordre d’approximation conduit à des oscillations, surtout pour des choix de pas de temps trop petits par rapport à la discrétisation en espace On consultera à ce sujet par exemple [bib10]. Ce problème est également en rapport avec la manière de calculer la matrice dite de masse, et on pourra consulter à ce sujet [bib14]. Nous donnons par ailleurs en annexe, pour illustrer notre propos, la solution pour le premier pas de temps d’un problème de consolidation mono dimensionnel avec une interpolation \(\mathit{P1P1}\) . On voit que pour un petit pas de temps, elle est très oscillante.

Pour cette raison, les éléments quadratiques THM sont des éléments en \(\mathit{P2P1}\) , c’est-à-dire que l’interpolation des déplacements est quadratique et celle des températures et pressions est linéaire. Nous avons néanmoins gardé toutes les inconnues sur tous les nœuds, y compris les nœuds milieux, mais nous avons imposé dans le calcul des matrices de rigidité que la pression d’un nœud milieu de segment soit égale à la demi somme des nœuds sommet du segment auquel il appartient.

Par ailleurs, dans la programmation, nous avons tenu compte de la propriété suivante:

Soit \(s\) un nœud sommet, \({w}_{s}^{1}\) sa fonction de forme en tant qu’appartenant à un élément linéaire (par exemple QUAD4), et \({w}_{s}^{2}\) sa fonction de forme en tant qu’appartenant à un élément quadratique (par exemple QUAD8). Soit \(\mathrm{na}\) le nombre d’arêtes ayant \(s\) comme extrémité et \({w}_{\text{ma}}^{2}\) la fonction de forme de l’interpolation quadratique attachée à un nœud milieu d’arête, on a alors la relation:

\({w}_{s}^{1}={w}_{s}^{2}+\sum_{\text{ma}=1}^{\text{na}}\frac{{w}_{\text{ma}}^{2}}{2}\)

Ceci dit, et y compris en interpolation \(\mathit{P2P1}\) , des conditions de non oscillation existent sur le pas de temps. [bib10] donne la relation:

\(\Delta t >\frac{\Delta {x}^{2}}{20{C}_{v}}\)

où \({C}_{v}\) est le coefficient de consolidation: \({C}_{v}=\frac{\text{kE}(1-\nu )}{{\rho}_{\text{lq}}(1+\nu )(1-\mathrm{2\nu })}\) , \(k\) étant la perméabilité mesurée en \(m/s\) .

Bibliographie

1. 15. COUSSY: «Mécanique des milieux poreux». éditions TECHNIP.

2. 16. GERMAIN: «Cours de mécanique des milieux continus».

3. 7. DUVAUT, J. L. LIONS: «Les inéquations en mécanique et en physique».

4. 3. CHAVANT, P. CHARLES, Th. DUFORESTEL, F. VOLDOIRE: «Thermo-hydro-mécanique des milieux poreux non saturés dans le Code_Aster ». Note HI-74/99/011/A.

5. «Algorithme non linéaire quasi statique (opérateur STAT_NON_LINE)» Document Aster [R5.03.01]

6. «Modèles de comportement THHM» Document Aster [R7.01.11] Indice A

7. 1. GIRAUD: «Adaptation au modèle poroélastique non linéaire de Lassabatère-Coussy à la modélisation milieu poreux non saturé», (ENSG).

8. 10. WABINSKI, F. VOLDOIRE: Thermohydromécanique en milieu saturé. Note EDF/DER HI‑74/96/010, de septembre 1996.

9. 20. LASSABATERE: «Couplages hydromécaniques en milieu poreux non saturé avec changement de phase: application au retrait de dessiccation du béton». Thèse ENPC.

10. Ph. MESTAT, M. PRAT: «Ouvrages en interaction». Hermes

11. 10. 6. et al THILUS: «Poro-mechanics», Biot Conférence 1998.

12. CROUZEIX et MIGNOT, MASSON: «Analyse numérique des équations différentielles» 1992

13. «Algorithme de thermique transitoire» Document Aster [R5.02.01].

14. «Diagonalisation de la matrice de masse thermique» Document Aster [R3.06.07].

15. O.C.ZIENKIEWICZ, C.T. CHANG, P. BETTESS: Drained, undrained, consolidating dynamic behaviour assumptions in soils. Géotechnique, 30, pp. 385-395, 1980.

Problème mono dimensionnel P1P1

On considère un problème de consolidation unidimensionnel dont les inconnues ne varient qu’en fonction de la seule variable d’espace \(x\) .

Un domaine rectangulaire de longueur \(L\) , est rempli d’un matériau poreux de coefficients de Lamé \(\lambda\) et \(\mu\) , de coefficient de biot \(b\) et de module de biot \(N=\frac{1}{M}\) et de conductivité hydraulique \({\lambda}_{h}\) . La masse volumique du fluide est notée \(\rho\) . On note \(\sigma\) la contrainte \({\sigma}_{xx}\) , \(u\) le déplacement dans la direction \(x\) , \(\varepsilon =\frac{\partial u}{\partial x}\) la déformation, \(p\) la pression, \(m\) l’apport massique en fluide, \(M\) le flux de fluide.

Les conditions aux limites sont:

en \(x=0\) : \(\sigma =0;p=0\)

en \(x=L\) : \(u=0;p=1\) pour \(t>0\)

Les conditions initiales en \(t=0\) sont \(\sigma =u=p=0\) .

La mise en équation donne:

Equilibre mécanique et élasticité linéaire: \(\sigma =(\lambda +2\mu )\varepsilon -\text{bp}=0\)

Conservation de la masse: \(\frac{\partial m}{\partial t}+\frac{\partial M}{\partial x}=0\)

Loi de Darcy: \(M=-{\lambda}_{h}\rho \frac{\partial p}{\partial x}\)

Couplages: \(\frac{m}{\rho}=\text{Np}+b\varepsilon\)

Nous faisons une discrétisation en temps implicite et nous nous intéressons au calcul du premier pas de temps :

On obtient le système de deux équations:

\(\lbrace \begin{array}{}\sigma =(\lambda +2\mu )\varepsilon (u)-\text{bp}=0\\ \text{Np}+b\varepsilon -{\lambda}_{h}\Delta t \Delta p =0\end{array}\)

La formulation variationnelle mixte de ce problème est:

\(\lbrace \begin{array}{}{\int}_{\Omega}\left[(\lambda +2\mu )\varepsilon (u)\varepsilon ({u}^{\text{*}})-\text{bp}\varepsilon ({u}^{\text{*}})\right]=0\forall {u}^{\text{*}}\\ {\int}_{\Omega}\left[{\text{Npp}}^{\text{*}}+b\varepsilon (u){p}^{\text{*}}+{\lambda}_{h}\Delta t \nabla p\nabla {p}^{\text{*}}\right]=0\forall {p}^{\text{*}}\end{array}\) éq An 1-1

Concernant la discrétisation spatiale, nous découpons le milieux en \(n\) éléments finis. Les nœuds de l’élément \(i\) sont \(i\) et \(i+1\) . On note \({u}_{e}^{1}\) et \({p}_{e}^{1}\) le déplacement et la pression du premier nœud de l’élément \(e\) , \({u}_{e}^{2}\) et \({p}_{e}^{2}\) le déplacement et la pression de son deuxième nœud. On suppose que l’on utilise des éléments finis P1P1, c’est à dire que les déplacements comme les pressions sont interpolés linéairement. La discrétisation de la première équation donne alors:

\(\sum_{e}\frac{{u}_{e}^{\text{*}2}-{u}_{e}^{\text{*}1}}{\Delta x }\left[(\lambda +2\mu )\frac{{u}_{e}^{2}-{u}_{e}^{1}}{\Delta x }-b\frac{{p}_{e}^{1}+{p}_{e}^{2}}{2}\right]=0\)

D’où l’on déduit:

\(({u}_{e}^{2}-{u}_{e}^{1})=\frac{b\Delta x }{\lambda +2\mu }\frac{{p}_{e}^{1}+{p}_{e}^{2}}{2}\text{}\forall e\) éq An 1-2

La discrétisation de la deuxième équation donne:

\(\sum_{e}\frac{{p}_{e}^{\text{*}2}+{p}_{e}^{\text{*}1}}{2}\left[b\frac{{u}_{e}^{2}-{u}_{e}^{1}}{\Delta x }+N\frac{{p}_{e}^{1}+{p}_{e}^{2}}{2}\right]+{\lambda}_{h}\Delta t \sum_{e}\frac{{p}_{e}^{\text{*}2}-{p}_{e}^{{\text{*}}^{1}}}{\Delta x }\frac{{p}_{e}^{2}-{p}_{e}^{1}}{\Delta x }=0\)

En y portant [éq An 1-2], on trouve:

\(\sum_{e}\frac{{p}_{e}^{\text{*}2}+{p}_{e}^{\text{*}1}}{2}(N+\frac{{b}^{2}}{\lambda +2\mu })\frac{{p}_{e}^{1}+{p}_{e}^{2}}{2}+{\lambda}_{h}\frac{\Delta t }{\Delta {x}^{2}}\sum_{e}({p}_{e}^{\text{*}2}-{p}_{e}^{\text{*}1})({p}_{e}^{2}-{p}_{e}^{1})=0\) éq An 1-3

Si maintenant, nous faisons tendre le pas de temps vers zéro à pas d’espace constant, \({\lambda}_{h}\frac{\Delta t }{\Delta {x}^{2}}\text{<<}(N+\frac{{b}^{2}}{\lambda +2\mu })\) , et [éq An 1-3] se réduit à:

\(\sum_{e}({p}_{e}^{\text{*}2}+{p}_{e}^{\text{*}1})({p}_{e}^{1}+{p}_{e}^{2})=0\)

On introduit une numérotation globale des nœuds et inconnues de pression:

\({p}_{j}^{1}={p}_{j};{p}_{j}^{2}={p}_{j+1}\)

On peut voir que cet ensemble de relations donne finalement:

\(\lbrace \begin{array}{}{p}_{1}=0\\ {p}_{i}+2{p}_{i+1}+{p}_{i+2}=0\forall i\in \left[1,n-1\right]\\ {p}_{n+1}=1\end{array}\)

La solution de cette suite est :

\(\lbrace \begin{array}{}{p}_{2j-1}=\frac{j-1}{n}\\ {p}_{2j}=-\frac{2j-1}{2n}\end{array}\)

Ce qui donne la répartition de pression suivante:

A titre indicatif, nous donnons ci dessous une comparaison de résultats numériques obtenus avec des éléments \(\mathit{P1P1}\) , \(\mathit{P2P2}\) et \(\mathit{P2P1}\) .

On voit que l’élément \(\mathit{P2P1}\) n’enlève pas l’oscillation, mais l’atténue sensiblement.