v7.31.114 WTNV114 - Flux hydrique sur un milieu poreux saturé#
Résumé :
On étudie le comportement hydraulique d’un milieu poreux saturé. Neuf modélisations sont effectuées: deux modélisations bidimensionnelles (modélisations A et B) et sept modélisations tridimensionnelles (modélisations C, D, E, F, G, H, I).
La distinction entre les modélisations réside dans la loi de comportement du fluide (LIQU_SATU pour les modélisations A, C, G, H, I et LIQU_GAZ_ATM pour les modélisations B, D, E, F).
Ce test consiste à appliquer un flux hydrique sur la face supérieure du modèle et à étudier l’effet de ce flux sur la distribution de la pression du fluide dans le milieu saturé. Il s’agit d’un problème évolutif.
Les huit premiers modèles étudiés sont 2D plans (HM_DPQ8) et 3D volumiques (HM_HEXA20 ou HM_PYRAM13) avec un comportement linéaire. Le dernier modèle est HM-XFEM.
La solution de référence est unidimensionnelle car elle ne dépend que de la coordonnée verticale.
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.
L’équation de conservation de la masse d’eau \({m}_{l}\) est donnée par l’expression suivante :
avec \({M}_{l}\) le flux d’eau tel que, si on néglige la gravité :
\({M}_{l}={\rho}_{l}{\lambda}_{l}(-\nabla {p}_{l})\) avec \({p}_{l}\) la pression de liquide, \({\rho}_{l}\) sa masse volumique et \({\lambda}_{l}\) sa conductivité hydraulique telle que \({\lambda}_{l}=\frac{{K}_{i}}{\mu}\) . \({K}_{i}\) est la perméabilité intrinsèque et \(\mu\) la viscosité du liquide.
En considérant un solide indéformable, on peut écrire que
\(\frac{d{m}_{l}}{\mathrm{dt}}={\rho}_{l}\frac{\phi }{{K}_{l}}\frac{d{p}_{l}}{\mathrm{dt}}\) avec \(\phi\) la porosité et \({K}_{l}\) la compressibilité de l’eau. On note \(N=\frac{\phi }{{K}_{l}}\) .
L’équation (1) devient alors :
\({\rho}_{l}N\frac{d{p}_{l}}{\mathrm{dt}}+D({\rho}_{l}{\lambda}_{l}(-\nabla {p}_{l}))=0\)
dont la formulation variationnelle est la suivante :
avec \({M}^{\mathrm{ext}}\) le chargement extérieur.
Pour établir la solution on se place dans le cas unidimensionnel et on adopte une discrétisation correspondant à un unique élément de degré 1 puisque en modélisation THM la partie hydraulique est traitée sur des éléments linéaires. On suppose également que les non linéarités sont faibles dans ce cas et que les coefficients \(N\) et \({\rho}_{l}\) sont constants ce qui suppose une variation relativement faibles de la pression.
Soit un élément linéaire :
On écrit alors la pression sur la base des fonctions de formes tel que :
\({p}_{l}(z,t)=\sum_{i=1}^{i=2}{p}^{i}(t){\lambda}_{i}(z)\)
avec
\({\lambda}_{1}(z)=0,5(1+\mathrm{2z})\)
\({\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 \mathrm{dz}}\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 {p}_{l}\rbrace =\left\lbrace \begin{array}{c}{P}_{l}^{1}\\ {P}_{l}^{2}\end{array}\right\rbrace\) \(\lbrace {M}^{\mathrm{ext}}\rbrace =\left\lbrace \begin{array}{c}{{M}^{\mathrm{ext}}}^{1}\\ {{M}^{\mathrm{ext}}}^{2}\end{array}\right\rbrace\)
Sachant qu’on injecte un flux d’eau Q sur la face supérieure, on a donc
\(\lbrace {M}^{\mathrm{ext}}\rbrace =\left\lbrace \begin{array}{c}0\\ -Q\end{array}\right\rbrace\)
L’équation (2) devient \(\frac{N}{{\lambda}_{l}}\left[A\right]\lbrace \frac{{\mathrm{dp}}_{l}}{\mathrm{dt}}\rbrace +\left[B\right]\lbrace {p}_{l}\rbrace =\frac{-1}{{\lambda}_{l}\rho }\lbrace {M}^{\mathrm{ext}}\rbrace\)
Considérant que pour des variations de temps courts Dt, l’évolution de p est quasiment linéaire, on peut écrire que
\(\lbrace \frac{{\mathrm{dp}}_{l}}{\mathrm{dt}}\rbrace =\frac{1}{\Delta t}\left\lbrace \begin{array}{c}{P}_{l}^{1}-{p}_{0}^{1}\\ {P}_{l}^{2}-{p}_{0}^{2}\end{array}\right\rbrace\)
et comme on part ici d’un état initial de pression nulle, on aura
\(\lbrace \frac{{\mathrm{dp}}_{l}}{\mathrm{dt}}\rbrace =\frac{1}{\Delta t}\left\lbrace \begin{array}{c}{P}_{l}^{1}\\ {P}_{l}^{2}\end{array}\right\rbrace\)
Finalement on obtient le système de deux équations à deux inconnues :
\({p}_{l}^{1}(\frac{N}{{\lambda}_{l}\Delta t}+6)-{\mathrm{6p}}_{l}^{2}=\frac{-\mathrm{2Q}}{{\lambda}_{l}{\rho}_{l}}\)
\({p}_{l}^{2}(\frac{N}{{\lambda}_{l}\Delta t}+6)-{\mathrm{6p}}_{l}^{1}=\frac{-\mathrm{4Q}}{{\lambda}_{l}{\rho}_{l}}\)
On a alors le résultat suivant:
\(\left\lbrace \begin{array}{c}{P}_{l}^{1}\\ {P}_{l}^{2}\end{array}\right\rbrace ={\left[\begin{array}{cc}\frac{N}{{\lambda}_{l}\Delta t}+6& -6\\ -6& \frac{N}{{\lambda}_{l}\Delta t}+6\end{array}\right]}^{-1}.\left\lbrace \begin{array}{c}\frac{-\mathrm{2Q}}{{\lambda}_{l}{\rho}_{l}}\\ \frac{-\mathrm{4Q}}{{\lambda}_{l}{\rho}_{l}}\end{array}\right\rbrace\)
Grandeurs et résultats de référence#
On teste la valeur de la pression de pore sur les faces inférieures et supérieures du massif.
Incertitudes sur la solution#
Aucune incertitude sur la solution. La solution est analytique.
Références bibliographiques#
Mécanique des milieux poreux . O. Coussy, Editions Technip, 2000.
Modélisation A#
Caractéristiques de la modélisation#
Modélisation en déformations planes D_PLAN_HM.
Comportement hydraulique LIQU_SATU.
Caractéristique du maillage#
1 éléments \(\mathit{Q8}\) .
Grandeurs et résultats testés#
Discrétisation en temps: 7 pas de temps croissant. La liste d’instant en secondes est: \((1,5,10,50,100,500,1000)\) .
Si l’on regarde la géométrie du problème le nœud A (numéro N1) correspond à la pression \({p}^{1}\) et le nœud C à la pression \({p}^{2}\)
Tableau de résultats aux différents instants :
\(N°\) NŒUD |
Numéro d’ordre |
\(\mathrm{PRE1}(\mathrm{Pa})\) |
Tolérance \((\text{\%})\) |
\(\mathrm{N1},A\) |
1 |
\(-6.631\times {10}^{3}\) |
\(5\) |
2 |
\(-3.315\times {10}^{4}\) |
\(5\) |
|
3 |
\(-6,631\times {10}^{4}\) |
\(5\) |
|
4 |
\(-3.314\times {10}^{5}\) |
\(5\) |
|
7 |
\(-6,553\times {10}^{5}\) |
\(5\) |
|
\(\mathrm{N3},C\) |
1 |
\(1.326\times {10}^{4}\) |
\(5\) |
2 |
\(6.631\times {10}^{4}\) |
\(5\) |
|
3 |
\(1.326\times {10}^{5}\) |
\(5\) |
|
4 |
\(6.629\times {10}^{5}\) |
\(5\) |
|
7 |
\(1.318\times {10}^{7}\) |
\(5\) |
Tableau 1 : Résultats
Les résultats sont en parfait accord avec la solution analytique.
Modélisation B#
Il s’agit de la même modélisation mais avec une loi de comportement hydraulique LIQU_GAZ_ATM. On a donc une phase gazeuse constante à 1atm (hypothèse de Richard). Les résultats sont logiquement exactement les mêmes que ci-dessus.
Modélisation C#
Il s’agit de la même modélisation qu’en A mais sur un cube 3D (1 élément HEXA20). Les conditions sur les faces de devant et de derrière sont des flux et déplacements nuls. Il s’agit donc du même cas qu’en 2D et les résultats sont logiquement identiques.
Modélisation D#
Il s’agit de la même modélisation qu’en C mais avec une loi de comportement hydraulique LIQU_GAZ_ATM. On a donc une phase gazeuse constante à 1atm (hypothèse de Richard). Les résultats sont logiquement exactement les mêmes que ci-dessus.
Modélisation E#
Exactement le même cas que la modélisation D mais en testant l’option SYME=”OUI” du solveur linéaire.
Les résultats sont logiquement les mêmes que pour la modélisation D.
Modélisation F#
Cette modélisation est la même que la modélisation D mais en modélisation sélective. L’intégration est donc différente: l’intégration des termes instationnaires se fait aux nœuds et celle des autres termes aux points de Gauss.
Les résultats sur un élément diffèrent donc logiquement.
Tableau de résultats aux différents instants :
\(N°\) NŒUD |
Numéro d’ordre |
\(\mathrm{PRE1}(\mathrm{Pa})\) |
Tolérance \((\text{\%})\) |
\(\mathrm{N20},A\) |
1 |
\(8.79\times {10}^{-6}\) |
\(5\) |
2 |
\(1.84\times {10}^{-4}\) |
\(5\) |
|
3 |
\(6.24\times {10}^{-4}\) |
\(5\) |
|
4 |
\(0.018\) |
\(5\) |
|
7 |
\(6.20\) |
\(5\) |
|
\(\mathrm{N1},C\) |
1 |
\(6.63\times {10}^{3}\) |
\(5\) |
2 |
\(3.31\times {10}^{4}\) |
\(5\) |
|
3 |
\(6.63\times {10}^{4}\) |
\(5\) |
|
4 |
\(3.31\times {10}^{5}\) |
\(5\) |
|
7 |
\(6.68\times {10}^{6}\) |
\(5\) |
Tableau 8-1: Résultats
Modélisation G#
Avec les modélisation G et H, on entreprend de valider la maille PYRAM13pour les modélisations HM. On compare le comportement d’un HEXA20 (modélisation G) par rapport à 6 PYRAM13 (modélisation H).
Caractéristiques de la modélisation#
Modélisation en déformations planes 3D_HM.
Comportement hydraulique LIQU_SATU.
Le régime transitoire d’établissement du flux étant très sensible à la structure, on diminue l’inverse de la compressibilité du fluide qui passe à \(1.0E-18\) et on réalise les tests au pas de temps 8 \((10000)\) de manière à dépasser le régime transitoire.
Caractéristique du maillage#
Le cube est composé d’un HEXA20.#
Grandeurs testées et résultats#
On teste les valeurs minimales et maximales de la pression de pore sur la face inférieure et la face supérieure du cube par rapport à la solution analytique.
Grandeurs testées |
Type de référence |
Valeur de référence |
Tolérance |
PRE1(face supérieure) MIN |
“ANALYTIQUE” |
1,250124996E+16 |
0,05 |
PRE1(face supérieure) MAX |
“ANALYTIQUE” |
1,250124996E+16 |
0,05 |
PRE1(face inférieure) MIN |
“ANALYTIQUE” |
1,249875004E+16 |
0,05 |
PRE1(face inférieure) MAX |
“ANALYTIQUE” |
1,249875004E+16 |
0,05 |
Modélisation H#
Caractéristiques de la modélisation#
Les caractérisiques sont les mêmes que pour la modélisation G.
Caractéristique du maillage#
Le maillage est composé de 6 éléments de type PYRAM13au lieu d’un HEXA20.
Le maillage est représenté sur la figure .
Figure 10-1 : Maillage pyramidal
Grandeurs testées et résultats#
On teste les valeurs minimales et maximales de la pression de pore sur la face inférieure et la face supérieure du cube par rapport à la modélisation précédente.
Ces résultats, très proches de ceux obtenus dans la modélisation G, assurent que le comportement
hydraulique des mailles HM_PYRAM13 est identique à celui des mailles HM_HEXA20.
Grandeurs testées |
Type de référence |
Valeur de référence |
Tolérance |
PRE1(face supérieure) MIN |
“AUTRE ASTER” |
1.20598944073E+16 |
1.0E-05 |
PRE1(face supérieure) MAX |
“AUTRE ASTER” |
1.20598944073E+16 |
1.0E-05 |
PRE1(face inférieure) MIN |
“AUTRE ASTER” |
1.20576822134E+16 |
1.0E-05 |
PRE1(face inférieure) MAX |
“AUTRE ASTER” |
1.20576822135E+16 |
1.0E-05 |
Modélisation I#
Avec la modélisation I, on entreprend de valider l’option CHAR_MECA_FLUX pour les éléments HM-XFEM. La modélisation G sert de référence.
Caractéristiques de la modélisation#
Cette modélisation est une modélisation HM-XFEM.
Les caractéristiques sont les mêmes que pour la modélisation G.
Caractéristiques du maillage#
Le maillage initial est le même que dans la modélisation G mais le cube est coupé en deux par une fissure verticale.
On applique les mêmes conditions aux limites que dans la modélisation G. Dans cette modélisation, les éléments de bords situés sur les faces supérieures et inférieures du cube et sur lesquels sont appliquées les conditions aux limites sont des éléments de bord HM-XFEMnon conformes à la fissure. Ils sont ensuite découpés en sous-éléments de manière à être compatibles avec la fissure. Aucun chargement n’est appliqué au niveau de la fissure.#
Grandeurs et résultats testés#
On teste les valeurs minimales et maximales de la pression de pore sur la face inférieure et la face supérieure du cube par rapport à la modélisation G.
Grandeurs testées |
Type de référence |
Valeur de référence |
Tolérance |
PRE1(face supérieure) MIN |
“AUTRE ASTER” |
1.20598944073E+16 |
1.0E-04 |
PRE1(face supérieure) MAX |
“AUTRE ASTER” |
1.20598944073E+16 |
1.0E-04 |
PRE1(face inférieure) MIN |
“AUTRE ASTER” |
1.20576822134E+16 |
1.0E-04 |
PRE1(face inférieure) MAX |
“AUTRE ASTER” |
1.20576822135E+16 |
1.0E-04 |
Ces résultats, très proches de ceux obtenus dans la modélisation G, assurent que le comportement hydraulique des éléments HM-XFEM est similaire au comportement hydraulique des éléments HM. En particulier, l’option CHAR_MECA_FLUX est validée pour les éléments HM-XFEM.
Synthèse des résultats#
Les résultats sont cohérents physiquement. Les modélisations A à G admettent une solution analytique qui est parfaitement vérifiée et les modélisations H et I sont parfaitement cohérentes avec la modélisations G de référence. Elles valident respectivement le comportement hydraulique des mailles pyramidales HM et l’option CHAR_MECA_FLUX pour les éléments HM-XFEM.