v7.31.109 WTNV109 - Chargement hydrique et mécanique d’un milieu poreux saturé#
Résumé :
On considère un problème tridimensionnel de couplage thermo-hydro-mécanique d’un milieu poreux saturé.
Ce test consiste à étudier l’effet de la mécanique et de l’hydraulique sur la thermique. On étire l’élément en lui imposant un déplacement dans la direction \(z\) , on impose une pression hydraulique sur l’ensemble du domaine et on étudie l’effet de ces chargements sur la température du modèle. Il s’agit de regarder le très faible couplage de la poromécanique vers la thermique. On se limite au premier pas de temps.
Les modèles étudiés sont \(\mathrm{2D}\) plans (DPQ8 et DPTR6) et \(\mathrm{3D}\) volumique (HEXA20) avec un comportement linéaire pour l’hydraulique et la thermique.
Solution de référence#
Méthode de calcul#
La solution de référence est unidimensionnelle car elle ne dépend que de la coordonnée verticale. On raisonnera donc sur un axe linéaire selon z.
Si on néglige la gravité et qu’on considère qu’il n’y a pas de source de chaleur externe, l’équation de la thermique est donnée par l’expression suivante :
Où \({h}_{l}^{m}\) est l’enthalpie massique de l’eau, \({m}_{l}\) sa masse, \({M}_{l}\) son flux, q le flux de chaleur et \(Q'\) la chaleur non convectée.
L’équation de conservation de la masse est la suivante :
En l’absence de flux hydrique, l’équation (1) se simplifie pour devenir :
Dans cette équation, la quantité \(Q'\) 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. On choisit une variation de pression faible et des paramètres de manière à ce que la température varie peu et donc que :
\(\frac{dQ'}{T}=\frac{dQ'}{{T}_{0}}\)
\(\mathrm{dQ}'\) a par conséquent cette expression :
\(\mathrm{dQ}'=3{\alpha}_{0}{K}_{0}{T}_{0}d{\varepsilon}_{v}+{C}_{0}^{\varepsilon}\mathrm{dT}-3{\alpha}_{l}^{m}{T}_{0}{\mathrm{dp}}_{l}\)
avec :
\({\alpha}_{0}\) le coefficient de dilatation thermique homogénéisée équivalent à celui du solide \({\alpha}_{s}\) .
\({K}_{0}\) le coefficient drainé d’élasticité.
\({C}_{0}^{\varepsilon}\) la chaleur a déformation constante qui a pour expression \({C}_{0}^{\varepsilon}={C}_{0}^{\sigma}-9{T}_{0}{K}_{0}{\alpha}_{0}^{2}\) et \({C}_{0}^{\sigma}\) la chaleur à contrainte constante.
\({\alpha}_{l}^{m}\) Le coefficient de dilatation thermique relative du liquide par rapport au squelette, il a pour expression : \({\alpha}_{l}^{m}=(b-\phi ){\alpha}_{s}+\phi {\alpha}_{l}\) avec b le coefficient de Biot, \(\phi\) la porosité et \({\alpha}_{l}\) le coefficient de dilatation du liquide.
\(d{\varepsilon}_{v}\) La variation de déformation volumique.
\({\mathrm{dp}}_{l}\) La variation de pression de liquide.
Par ailleurs, le flux de chaleur a l’expression suivante : \(q=-{\lambda}_{T}\partial \frac{T}{\partial z}\) où \({\lambda}_{T}\) est le coefficient de conductivité thermique.
En remplaçant \(\mathrm{dQ}'\) et q par leur valeur dans l’équation (3) et se limitant à seul pas de temps \(\Delta t\) , on obtient :
\(3{\alpha}_{0}{K}_{0}{T}_{0}\Delta {\varepsilon}_{v}-3{\alpha}_{l}^{m}{T}_{0}\Delta {p}_{l}+{C}_{0}^{\varepsilon}\Delta T=-\Delta tdq={\lambda}_{l}\Delta t\frac{{\partial}^{2}T}{\partial {z}^{2}}\)
En considérant des températures et les déplacements initialement nuls, on peut écrire :
\(\Delta T=T(t)-T(0)=T\) et \(\Delta \varepsilon =\varepsilon (t)-\varepsilon (0)=\varepsilon\) et \(\Delta {p}_{l}={p}_{l}(t)-{p}_{l}(0)={p}_{l}\)
Donc au premier pas de temps, on aura :
\(T-a\frac{{\partial}^{2}T}{\partial {z}^{2}}=b\)
avec \(a=\frac{{\lambda}_{l}\Delta T}{{C}_{0}^{\varepsilon}}\) et \(b=-\frac{3{\alpha}_{0}{K}_{0}{T}_{0}\varepsilon -3{\alpha}_{l}^{m}{T}_{0}{p}_{l}}{{C}_{0}^{\varepsilon}}\)
La formulation variationnelle de cette expression (dans le cas unidimensionnel) est alors la suivante :
Pour établir la solution analytique on considère un unique élément de degré 1 puisque en modélisation THM la partie thermo-hydraulique est traitée par des éléments linéaires.
Soit un élément linéaire :
On rappelle les conditions aux limites : \(\frac{\partial T}{\partial z}=0\) aux 2 extrémités (z=0,5 et z = -0,5)
La température sur la base des fonctions de formes s’écrit de la manière suivante :
\(T(z,t)=\sum_{i=1}^{i=2}{T}^{i}(t){\lambda}_{i}(z)\)
avec
\({\lambda}_{1}(z)=0,5(1+\mathrm{2z})\)
et
\({\lambda}_{2}(z)=0,5(1-\mathrm{2z})\)
On introduit alors les matrices suivantes :
\(\left[A\right]=\left[{A}_{ij}\right];{A}_{ij}={\int}_{-0,5}^{0,5}{\lambda}_{i}{\lambda}_{j}\mathrm{dz}\)
\(\left[B\right]=\left[{B}_{ij}\right];{B}_{ij}={\int}_{-0,5}^{0,5}\frac{\partial {\lambda}_{i}}{\partial z}\frac{\partial {\lambda}_{j}}{\partial z}\mathrm{dz}\)
ce qui conduit à :
\(\left[A\right]=\frac{1}{6}\left[\begin{array}{cc}2& 1\\ 1& 2\end{array}\right]\) et \(\left[B\right]=\left[\begin{array}{cc}1& -1\\ -1& 1\end{array}\right]\)
On note ensuite classiquement :
\(\lbrace T\rbrace =\left\lbrace \begin{array}{c}{T}^{1}\\ {T}^{2}\end{array}\right\rbrace\)
et \(\lbrace \frac{\partial T}{\partial t}\rbrace =\left\lbrace \begin{array}{c}\frac{\partial {T}^{1}}{\partial t}\\ \frac{\partial {T}^{2}}{\partial t}\end{array}\right\rbrace\)
L’équation (4) devient alors
\([A]\lbrace T\rbrace +a[B]\lbrace T\rbrace -a[A]\lbrace \frac{\partial T}{\partial t}\rbrace =[A]\lbrace b\rbrace\)
avec les conditions aux limites imposées (déplacement nul en bas et flux de températures nuls aux deux extrémités), on a :
\(\lbrace b\rbrace =\left\lbrace \begin{array}{c}0\\ b\end{array}\right\rbrace \mathrm{et}\lbrace \frac{\partial T}{\partial t}\rbrace =\left\lbrace \begin{array}{c}0\\ 0\end{array}\right\rbrace\)
Ce qui au final nous donne
\(\left\lbrace \begin{array}{c}{T}^{1}\\ {T}^{2}\end{array}\right\rbrace +a{[A]}^{-1}[B]\left\lbrace \begin{array}{c}{T}^{1}\\ {T}^{2}\end{array}\right\rbrace =\left\lbrace \begin{array}{c}0\\ b\end{array}\right\rbrace\)
Finalement on obtient :
\(\left[\begin{array}{cc}1+\mathrm{6a}& -\mathrm{6a}\\ -\mathrm{6a}& 1+\mathrm{6a}\end{array}\right]\lbrace T\rbrace =\left\lbrace \begin{array}{c}{T}^{1}\\ {T}^{2}\end{array}\right\rbrace\)
Grandeurs et résultats de référence#
Pour un temps court de 100s, on aura \(\left\lbrace \begin{array}{c}{T}^{1}\\ {T}^{2}\end{array}\right\rbrace =\left\lbrace \begin{array}{c}{2.10}^{-14}\\ 1,045{.10}^{-7}\end{array}\right\rbrace\)
On considérera \({T}^{1}\) nul.
Incertitude sur la solution#
Aucune
Modélisation A#
Caractéristiques de la modélisation#
Modélisation volumique 3D_THM
Caractéristiques du maillage#
1 maille HEXA20.
Grandeurs testés et résultats#
Discrétisation en temps: un seul pas de temps de \(100s\) . Le schéma en temps est implicite \((\theta =1)\) . Les résultats correspondent parfaitement à la solution analytique.
Nœud |
Type de valeur |
Instant \((s)\) |
Référence (analytique) |
Tolérance \((\text{\%})\) |
\(\mathrm{N1},\mathrm{N3}\) |
\(\mathrm{TEMP}\) |
\({10}^{2}\) |
\(-1.045\times {10}^{-7}\) |
\(0,1\) |
\(\mathrm{N6},\mathrm{N8}\) |
\(\mathrm{TEMP}\) |
\({10}^{2}\) |
\(-1.045\times {10}^{-7}\) |
\(0,1\) |
Modélisation B#
Caractéristiques de la modélisation#
Modélisation plane: D_PLAN_THM
Caractéristiques du maillage#
1 maille QUAD8.
Grandeurs testés et résultats#
Discrétisation en temps : un seul pas de temps de 100s. Le schéma en temps est implicite \((\vartheta =1)\) . Les résultats correspondent parfaitement à la solution analytique.
Nœud |
Type de valeur |
Instant (s) |
Référence (analytique) |
Tolérance \((\text{\%})\) |
\(\mathit{N3}\) |
\(\mathit{TEMP}\) |
\({10}^{2}\) |
\(-1.045\times {10}^{-7}\) |
\(0,1\) |
\(\mathit{N4}\) |
\(\mathit{TEMP}\) |
\({10}^{2}\) |
\(-1.045\times {10}^{-7}\) |
\(0,1\) |
Modélisation C#
Caractéristiques de la modélisation#
Modélisation plane :D_PLAN_THM
9
Caractéristiques du maillage#
2 mailles TRIA6.
Grandeurs testés et résultats#
Discrétisation en temps : un seul pas de temps : \({10}^{2}s\) . Le schéma en temps est implicite \((\vartheta =1)\) . Les résultats correspondent parfaitement à la solution analytique.
Nœud |
Type de valeur |
Instant (s) |
Référence (analytique) |
Tolérance \((\text{\%})\) |
\(\mathrm{N3}\) |
\(\mathrm{TEMP}\) |
\({10}^{2}\) |
\(-1.045\times {10}^{-7}\) |
\(0,1\) |
\(\mathrm{N4}\) |
\(\mathrm{TEMP}\) |
\({10}^{2}\) |
\(-1.045\times {10}^{-7}\) |
\(0,1\) |
Modélisation D#
Caractéristiques de la modélisation#
Modélisation volumique 3D_THM.
Caractéristiques du maillage#
2 maille PENTA15.
Grandeurs testés et résultats#
Discrétisation en temps: un seul pas de temps : \({10}^{2}s\) . Le schéma en temps est implicite \((\theta =1)\) . Les résultats correspondent parfaitement à la solution analytique.
Nœud |
Type de valeur |
Instant (s) |
Référence (analytique) |
Tolérance \((\text{\%})\) |
\(\mathrm{N3}\) |
\(\mathrm{TEMP}\) |
\({10}^{2}\) |
\(-1.045\times {10}^{-7}\) |
\(0,1\) |
\(\mathrm{N4}\) |
\(\mathrm{TEMP}\) |
\({10}^{2}\) |
\(-1.045\times {10}^{-7}\) |
\(0,1\) |
Synthèse des résultats#
Les résultats sont en cohérence avec la solution analytique.