2. Solution de référence#
2.1. Méthode de calcul utilisée pour la solution de référence#
L’équation de conservation de la masse fluide est donnée par l’expression suivante :
\(\frac{{\mathrm{dm}}_{i}}{\mathrm{dt}}+{\mathrm{DivM}}_{i}=0\) \(i\) variant de \(1\) au nombre de constituants (1)
Dans notre exemple, le modèle est constitué de deux fluides : l’eau liquide \((e)\) et l’air sec \((a)\) . L’équation (1) se divise donc en deux :
\(\lbrace \begin{array}{c}\frac{{\mathit{dm}}_{e}}{\mathit{dt}}+{\mathit{DivM}}_{e}=0\\ \frac{{\mathit{dm}}_{a}}{\mathit{dt}}+{\mathit{DivM}}_{a}=0\end{array}\) (2)
Les flux de fluide ont pour expression :
\(\lbrace \begin{array}{c}{M}_{e}={\rho}_{e}{\lambda}_{e}(-\nabla {p}_{e}+{\rho}_{e}g)\\ {M}_{a}={\rho}_{a}{\lambda}_{a}(-\nabla {p}_{a}+{\rho}_{a}g)\end{array}\) (3)
Or l’apport massique de fluide est défini par les équations (4) où \(N=\left[\begin{array}{cc}{N}_{\mathit{ee}}& {N}_{\mathit{ea}}\\ {N}_{\mathit{ae}}& {N}_{\mathit{aa}}\end{array}\right]\) est une matrice symétrique dont les termes (équations (5)) dépendent du degré de saturation \(S\) , de la porosité \(\phi\) , du coefficient de \(\mathrm{Biot}b\) , de la perméabilité du liquide et \({K}_{e}\) de l’élasticité de la matrice solide \({K}_{s}\) .
\(\lbrace \begin{array}{c}\frac{{\mathit{dm}}_{e}}{\mathit{dt}}={\rho}_{e}{N}_{\mathit{ee}}\frac{{\mathit{dp}}_{e}}{\mathit{dt}}+{\rho}_{e}{N}_{\mathit{ea}}\frac{{\mathit{dp}}_{a}}{\mathit{dt}}\\ \frac{{\mathit{dm}}_{a}}{\mathit{dt}}={\rho}_{a}{N}_{\mathit{ae}}\frac{{\mathit{dp}}_{e}}{\mathit{dt}}+{\rho}_{a}{N}_{\mathit{aa}}\frac{{\mathit{dp}}_{a}}{\mathit{dt}}\end{array}\) (4)
\(\lbrace \begin{array}{c}{N}_{\mathit{ee}}=-\varphi \frac{\partial S}{\partial {p}_{c}}+S(\frac{\varphi}{{K}_{e}}+\frac{b-\varphi }{{K}_{s}}S)\\ {N}_{\mathit{aa}}=-\varphi \frac{\partial S}{\partial {p}_{c}}+(1-S)(\frac{\varphi}{{p}_{a}}+\frac{b-\varphi }{{K}_{s}}(1-S))\\ {N}_{\mathit{ea}}={N}_{\mathit{ae}}=\varphi \frac{\partial S}{\partial {p}_{c}}+(1-S)(\frac{b-\varphi }{{K}_{s}}S)\end{array}\) (5)
La formulation variationnelle des équations (2), en tenant compte de (3) et (4) est :
\(\lbrace \begin{array}{c}{\int}_{\Omega}{N}_{\mathrm{ee}}\frac{{\mathrm{dp}}_{e}}{\mathrm{dt}}{p}_{e}^{\text{*}}+{\int}_{\Omega}{N}_{\mathrm{ea}}\frac{{\mathrm{dp}}_{a}}{\mathrm{dt}}{p}_{e}^{\text{*}}+{\int}_{\Omega}{\lambda}_{e}\nabla {p}_{e}.\nabla {p}_{e}^{\text{*}}={\int}_{\Omega}{\lambda}_{e}{\rho}_{e}g.\nabla {p}_{e}^{\text{*}}-{\int}_{\partial \Omega }\frac{{M}_{e}^{\mathrm{ext}}}{{\rho}_{e}}{p}_{e}^{\text{*}}\\ {\int}_{\Omega}{N}_{\mathrm{ea}}\frac{{\mathrm{dp}}_{e}}{\mathrm{dt}}{p}_{a}^{\text{*}}+{\int}_{\Omega}{N}_{\mathrm{aa}}\frac{{\mathrm{dp}}_{a}}{\mathrm{dt}}{p}_{a}^{\text{*}}+{\int}_{\Omega}{\lambda}_{a}\nabla {p}_{a}.\nabla {p}_{a}^{\text{*}}={\int}_{\Omega}{\lambda}_{a}{\rho}_{a}g.\nabla {p}_{a}^{\text{*}}-{\int}_{\partial \Omega }\frac{{M}_{a}^{\mathrm{ext}}}{{\rho}_{a}}{p}_{a}^{\text{*}}\end{array}\) (6)
Discrétisation
Pour le calcul de la solution analytique, on se place dans un cas unidimensionnel avec un seul élément de degré \(1\) .
Remarque :
En \(\mathrm{THM}\) , tous les maillages doivent être quadratiques, mais dans le cas hydraulique l’intégration est toujours linéaire (les nœuds milieux sont ignorés).
On suppose dans les deux cas que la pesanteur est orientée selon les \(z\) négatifs.
On suppose d’autre part que les non-linéarités sont faibles et que les coefficients \(N,\lambda ,\rho\) sont constants. Il faut donc que les variations de pression soient suffisamment faibles pour que \(N\) et \(\rho\) puissent être supposés constants.
En hydraulique, la discrétisation sera toujours linéaire.
Discrétisation linéaire :
On écrira :
\(p(z,t)=\sum_{i=1}^{2}{p}^{i}(t){\lambda}_{i}(z)\) (7)
Avec :
\(\lbrace \begin{array}{c}{\lambda}_{1}=\frac{1}{2}-z\\ {\lambda}_{2}=\frac{1}{2}+z\end{array}\) (8)
En introduisant alors les matrices et vecteurs :
\(\lbrace \begin{array}{c}[A]=[{A}_{ij}];{A}_{ij}=\underset{-1/2}{\overset{1/2}{\int}}{\lambda}_{i}{\lambda}_{j}\mathrm{dz}\\ [B]=[{B}_{ij}];{B}_{ij}=\underset{-1/2}{\overset{1/2}{\int}}\frac{d{\lambda}_{i}}{\mathrm{dz}}\frac{d{\lambda}_{j}}{\mathrm{dz}}\mathrm{dz}\\ \lbrace {F}_{g}\rbrace =\lbrace {F}_{\mathrm{gi}}\rbrace ;{F}_{\mathrm{gj}}=\underset{-1/2}{\overset{1/2}{\int}}\frac{d{\lambda}_{i}}{\mathrm{dz}}\mathrm{dz}\end{array}\) (9)
Et en notant :
\(\lbrace {p}_{e}\rbrace =\left\lbrace \begin{array}{c}{P}_{e}^{1}\\ {p}_{e}^{2}\end{array}\right\rbrace ;\lbrace {p}_{a}\rbrace =\left\lbrace \begin{array}{c}{P}_{a}^{1}\\ {p}_{a}^{2}\end{array}\right\rbrace\) (10)
\(\lbrace {M}_{e}^{\mathrm{ext}}\rbrace =\left\lbrace \begin{array}{c}{{M}_{e}^{\mathrm{ext}}}_{1}\\ {{M}_{e}^{\mathrm{ext}}}_{2}\end{array}\right\rbrace ;\lbrace {M}_{a}^{\mathrm{ext}}\rbrace =\left\lbrace \begin{array}{c}{{M}_{a}^{\mathrm{ext}}}_{1}\\ {{M}_{a}^{\mathrm{ext}}}_{2}\end{array}\right\rbrace\) (11)
Les équations (6) deviennent :
\(\lbrace \begin{array}{c}\frac{{N}_{\mathrm{ee}}}{{\lambda}_{e}}[A]\left\lbrace \frac{{\mathrm{dp}}_{e}}{\mathrm{dt}}\right\rbrace +\frac{{N}_{\mathrm{ea}}}{{\lambda}_{e}}[A]\left\lbrace \frac{{\mathrm{dp}}_{a}}{\mathrm{dt}}\right\rbrace +[B]\lbrace {p}_{e}\rbrace ={\rho}_{e}\lbrace {F}_{g}\rbrace -\frac{1}{{\lambda}_{e}{\rho}_{e}}\lbrace {M}_{e}^{\mathrm{ext}}\rbrace \\ \frac{{N}_{\mathrm{ae}}}{{\lambda}_{a}}[A]\left\lbrace \frac{{\mathrm{dp}}_{e}}{\mathrm{dt}}\right\rbrace +\frac{{N}_{\mathrm{aa}}}{{\lambda}_{a}}[A]\left\lbrace \frac{{\mathrm{dp}}_{a}}{\mathrm{dt}}\right\rbrace +[B]\lbrace {p}_{a}\rbrace ={\rho}_{a}\lbrace {F}_{g}\rbrace -\frac{1}{{\lambda}_{a}{\rho}_{a}}\lbrace {M}_{a}^{\mathrm{ext}}\rbrace \end{array}\) (12)
Le calcul des matrices \([A]\) et \([B]\) et du vecteur \(\lbrace F\rbrace\) donne :
\([A]=\frac{1}{3}\left[\begin{array}{cc}1& 1/2\\ 1/2& 1\end{array}\right];[B]=\left[\begin{array}{cc}1& -1\\ -1& 1\end{array}\right];\lbrace {F}_{g}\rbrace =\left\lbrace \begin{array}{c}-1\\ 1\end{array}\right\rbrace\) (13)
On définit alors \(\lbrace {v}_{1}\rbrace ,\lbrace {v}_{2}\rbrace\) les vecteurs propres de \({[A]}^{-1}[B]\) .
On a les propriétés :
\({\lbrace {v}_{i}\rbrace }^{T}[A]\lbrace {v}_{j}\rbrace ={\lbrace {v}_{i}\rbrace }^{T}[B]\lbrace {v}_{j}\rbrace =0\text{}\mathrm{si}\text{}i\ne j\) (14)
Et on pose :
\({a}_{i}={\lbrace {v}_{i}\rbrace }^{T}[A]\lbrace {v}_{i}\rbrace ,{b}_{i}={\lbrace {v}_{i}\rbrace }^{T}[B]\lbrace {v}_{i}\rbrace ,{f}_{i}={\lbrace {v}_{i}\rbrace }^{T}\lbrace {F}_{g}\rbrace \mathrm{et}{M}^{i}={\lbrace {v}_{i}\rbrace }^{T}\lbrace {M}^{\mathrm{ext}}\rbrace\) (15)
On trouve :
\(\lbrace {v}_{1}\rbrace =\lbrace \begin{array}{c}1\\ 1\end{array}\rbrace ;\lbrace {v}_{2}\rbrace =\lbrace \begin{array}{c}-1\\ 1\end{array}\rbrace\) (16)
\(\lbrace \begin{array}{c}{a}_{1}=1;{b}_{1}=0;{f}_{1}=0\\ {a}_{2}=\frac{1}{3};{b}_{2}=4;{f}_{2}=-\mathrm{2g}\end{array}\) (17)
On décompose alors \(\lbrace {p}_{e}\rbrace\) et \(\lbrace {p}_{a}\rbrace\) sur la base des \(\lbrace {v}_{i}\rbrace\)
\(\lbrace {p}_{e}\rbrace =\sum_{i=1}^{2}{\alpha}_{e}^{i}\lbrace {v}_{i}\rbrace ;\lbrace {p}_{a}\rbrace =\sum_{i=1}^{2}{\alpha}_{a}^{i}\lbrace {v}_{i}\rbrace\) (18)
Compte tenu des propriétés d’orthogonalité (14), le système d’équations (12) s’écrit :
\(\lbrace \begin{array}{c}\frac{{N}_{\mathrm{ee}}}{{\lambda}_{e}}{a}_{i}\frac{d{\alpha}_{e}^{i}}{\mathrm{dt}}+\frac{{N}_{\mathrm{ea}}}{{\lambda}_{e}}{a}_{i}\frac{d{\alpha}_{a}^{i}}{\mathrm{dt}}+{b}_{i}{\alpha}_{e}^{i}={\rho}_{e}{f}_{i}-\frac{1}{{\lambda}_{e}{\rho}_{e}}{M}_{e}^{i}\\ \frac{{N}_{\mathrm{ae}}}{{\lambda}_{a}}{a}_{i}\frac{d{\alpha}_{e}^{i}}{\mathrm{dt}}+\frac{{N}_{\mathrm{aa}}}{{\lambda}_{a}}{a}_{i}\frac{d{\alpha}_{a}^{i}}{\mathrm{dt}}+{b}_{i}{\alpha}_{a}^{i}={\rho}_{a}{f}_{i}-\frac{1}{{\lambda}_{a}{\rho}_{a}}{M}_{a}^{i}\end{array}\) (19)
Posant :
\(\lbrace {\alpha}^{i}\rbrace =\lbrace \begin{array}{c}{\alpha}_{e}^{i}\\ {\alpha}_{a}^{i}\end{array}\rbrace ;[N]=\left[\begin{array}{cc}{N}_{\mathrm{ee}}& {N}_{\mathrm{ea}}\\ {N}_{\mathrm{ae}}& {N}_{\mathrm{aa}}\end{array}\right];[L]=\left[\begin{array}{cc}{\lambda}_{e}& 0\\ 0& {\lambda}_{a}\end{array}\right]\) (20)
L’équation (19) s’écrit :
\([N]=\left[\frac{d{\alpha}^{i}}{\mathrm{dt}}\right]+\frac{{b}_{i}}{{a}_{i}}[L]\lbrace {\alpha}^{i}\rbrace =\frac{{f}_{i}}{{a}_{i}}\left\lbrace \begin{array}{c}{\rho}_{e}{\lambda}_{e}\\ {\rho}_{a}{\lambda}_{a}\end{array}\right\rbrace -\left\lbrace \begin{array}{c}{M}_{e}^{i}/{\rho}_{e}{a}_{i}\\ {M}_{a}^{i}/{\rho}_{a}{a}_{i}\end{array}\right\rbrace\) (21)
Conditions initiales
On suppose que :
\(\begin{array}{c}{p}_{a}(x,t=0)={p}_{a}^{0}\\ {p}_{e}(x,t=0)={p}_{a}^{0}-{p}_{c}^{0}\end{array}\) uniformes en espace ;
Compte tenu des valeurs des vecteurs \(\lbrace {v}_{1}\rbrace ,\lbrace {v}_{2}\rbrace\) (équations (16)), on voit facilement que :
\(\begin{array}{c}{\alpha}_{a}^{1}(t=0)={p}_{a}^{0};{\alpha}_{e}^{1}(t=0)={p}_{a}^{0}-{p}_{c}^{0}\\ {\alpha}_{a}^{2}(t=0)={\alpha}_{e}^{2}(t=0)=0\end{array}\) (22)
On se place dans un cas où les équations de l’hydraulique sont découplées \(({N}_{\mathrm{ea}}={N}_{\mathrm{ae}}=0)\) et dans lequel les flux de fluide sont nuls \((\lbrace {M}_{e}^{\mathrm{ext}}\rbrace =\lbrace {M}_{a}^{\mathrm{ext}}\rbrace =0)\) .
Compte tenu de (21) , de \({f}_{1}={f}_{3}=0\) (équations (17)), le système d’équations (21) a pour solution :
\(\left\lbrace \begin{array}{c}{\alpha}_{e}^{1}={P}_{a}^{0}-{p}_{c}^{0}\\ {\alpha}_{e}^{2}=\frac{{f}_{2}}{{b}_{2}}{\rho}_{e}(1-\exp(-\frac{{b}_{2}}{{a}_{2}}\frac{{\lambda}_{e}}{{N}_{\mathrm{ee}}}t))\\ {\alpha}_{a}^{1}={P}_{a}^{0}\\ {\alpha}_{a}^{2}=\frac{{f}_{2}}{{b}_{2}}{\rho}_{a}(1-\exp(-\frac{{b}_{2}}{{a}_{2}}\frac{{\lambda}_{a}}{{N}_{\mathrm{aa}}}t))\end{array}\right\rbrace\) (23)
On trouve en revenant aux variables nodales :
\(\left\lbrace \begin{array}{c}{P}_{1}\\ {p}_{2}\end{array}\right\rbrace =\left\lbrace \begin{array}{c}{\alpha}_{1}-{\alpha}_{2}\\ {\alpha}_{1}+{\alpha}_{2}\end{array}\right\rbrace\)
\({\left\lbrace \begin{array}{c}{P}_{1}\\ {p}_{2}\end{array}\right\rbrace }_{\mathrm{eau}}=\left\lbrace \begin{array}{c}{P}_{a}^{0}-{p}_{c}^{0}+\frac{{\rho}_{e}g}{2}(1-\exp(-12\frac{{\lambda}_{e}}{{N}_{\mathrm{ee}}}t))\\ {P}_{a}^{0}-{p}_{c}^{0}-\frac{{\rho}_{e}g}{2}(1-\exp(-12\frac{{\lambda}_{e}}{{N}_{\mathrm{ee}}}t))\end{array}\right\rbrace\) (24)
\({\left\lbrace \begin{array}{c}{P}_{1}\\ {p}_{2}\end{array}\right\rbrace }_{\mathrm{air}}=\left\lbrace \begin{array}{c}{P}_{a}^{0}+\frac{{\rho}_{a}g}{2}\exp(-12\frac{{\lambda}_{a}}{{N}_{\mathrm{aa}}}t)\\ {P}_{a}^{0}-\frac{{\rho}_{a}g}{2}\exp(-12\frac{{\lambda}_{a}}{{N}_{\mathrm{aa}}}t)\end{array}\right\rbrace\) (25)
et la pression capillaire définie comme la différence entre la pression d’air et la pression d’eau a pour valeur :
\({\left\lbrace \begin{array}{c}{P}_{1}\\ {p}_{2}\end{array}\right\rbrace }_{\mathrm{capillaire}}={\left\lbrace \begin{array}{c}{P}_{1}\\ {p}_{2}\end{array}\right\rbrace }_{\mathrm{air}}-{\left\lbrace \begin{array}{c}{P}_{1}\\ {p}_{2}\end{array}\right\rbrace }_{\mathrm{eau}}\)
Nous avons considéré le cas de calcul suivant :
\(\begin{array}{c}S\ne 1;\frac{\partial S}{\partial {p}_{c}}=0;{K}_{s}=\infty \\ {N}_{\mathrm{ee}}=S\frac{\varphi}{{K}_{e}};{N}_{\mathrm{aa}}=(1-S)\frac{\varphi}{{p}_{a}}\end{array}\)
2.2. Grandeur de référence#
Evolution de la pression capillaire et de la pression d’air sec en fonction du temps aux points
\(C,D(z=h)\)
\(A,B(z=0)\)
Pour la discrétisation quadratique, vérification de la valeur constante de la pression aux nœuds \(E,F(z=\frac{h}{2})\) .
2.3. Incertitudes#
Solution analytique sur les équations d’hydrauliques donc incertitudes négligeables pour les modélisations A, B, C.
Attention ces solutions analytiques ne s’appliquent pas aux modélisations sélectives ou lumpées (D et E). En effet, dans ce dernier cas, les intégrations sont faites aux nœuds et non plus aux points de Gauss. En effet l’intégration par point de Gauss est exacte en \(\mathrm{1D}\) pour des polynôme de degré inférieur ou égal à 3 et donc pour toutes les intégrales présentées dans l’équation (9). En revanche la méthode d’intégration au sommet n’est exacte que pour les polynômes de degré 1. On voit donc que les termes de la matrice \([A]\) seront sous intégrés. Il est donc logique que sur un maillage unitaire comme ici les résultats obtenus ici ne soient pas exacts. On conserve cependant ces tests mais avec un résultat en «non régression».
2.4. Références bibliographiques#
Thermo-hydro-mécanique des milieux poreux dans Code_Aster – Note EDF, HI-74/99/011/A