2. Benchmark solution#
2.1. Calculation method used for the reference solution#
The equation for the conservation of fluid mass is given by the following expression:
\(\frac{{\mathrm{dm}}_{i}}{\mathrm{dt}}+{\mathrm{DivM}}_{i}=0\) \(i\) varying from \(1\) to the number of components (1)
In our example, the model consists of two fluids: liquid water \((e)\) and dry air \((a)\). Equation (1) is therefore divided in two:
\(\mathrm{\{}\begin{array}{c}\frac{{\mathit{dm}}_{e}}{\mathit{dt}}+{\mathit{DivM}}_{e}\mathrm{=}0\\ \frac{{\mathit{dm}}_{a}}{\mathit{dt}}+{\mathit{DivM}}_{a}\mathrm{=}0\end{array}\) (2)
Fluid flows are expressed as:
\(\mathrm{\{}\begin{array}{c}{M}_{e}\mathrm{=}{\rho }_{e}{\lambda }_{e}(\mathrm{-}\mathrm{\nabla }{p}_{e}+{\rho }_{e}g)\\ {M}_{a}\mathrm{=}{\rho }_{a}{\lambda }_{a}(\mathrm{-}\mathrm{\nabla }{p}_{a}+{\rho }_{a}g)\end{array}\) (3)
However, the mass supply of fluid is defined by equations (4) where \(N\mathrm{=}\left[\begin{array}{cc}{N}_{\mathit{ee}}& {N}_{\mathit{ea}}\\ {N}_{\mathit{ae}}& {N}_{\mathit{aa}}\end{array}\right]\) is a symmetric matrix whose terms (equations (5)) depend on the degree of saturation \(S\), on the porosity \(\phi\), on the coefficient of \(\mathrm{Biot}b\), on the permeability of the liquid and \({K}_{e}\) on the elasticity of the solid matrix \({K}_{s}\).
\(\mathrm{\{}\begin{array}{c}\frac{{\mathit{dm}}_{e}}{\mathit{dt}}\mathrm{=}{\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}}\mathrm{=}{\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)
\(\mathrm{\{}\begin{array}{c}{N}_{\mathit{ee}}\mathrm{=}\mathrm{-}\varphi \frac{\mathrm{\partial }S}{\mathrm{\partial }{p}_{c}}+S(\frac{\varphi }{{K}_{e}}+\frac{b\mathrm{-}\varphi }{{K}_{s}}S)\\ {N}_{\mathit{aa}}\mathrm{=}\mathrm{-}\varphi \frac{\mathrm{\partial }S}{\mathrm{\partial }{p}_{c}}+(1\mathrm{-}S)(\frac{\varphi }{{p}_{a}}+\frac{b\mathrm{-}\varphi }{{K}_{s}}(1\mathrm{-}S))\\ {N}_{\mathit{ea}}\mathrm{=}{N}_{\mathit{ae}}\mathrm{=}\varphi \frac{\mathrm{\partial }S}{\mathrm{\partial }{p}_{c}}+(1\mathrm{-}S)(\frac{b\mathrm{-}\varphi }{{K}_{s}}S)\end{array}\) (5)
The variational formulation of equations (2), taking into account (3) and (4), is:
\(\{\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}\mathrm{.}\nabla {p}_{e}^{\text{*}}={\int }_{\Omega }{\lambda }_{e}{\rho }_{e}g\mathrm{.}\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}\mathrm{.}\nabla {p}_{a}^{\text{*}}={\int }_{\Omega }{\lambda }_{a}{\rho }_{a}g\mathrm{.}\nabla {p}_{a}^{\text{*}}-{\int }_{\partial \Omega }\frac{{M}_{a}^{\mathrm{ext}}}{{\rho }_{a}}{p}_{a}^{\text{*}}\end{array}\) (6)
Discretization
To calculate the analytical solution, we use a one-dimensional case with a single element of degree \(1\).
Note:
In \(\mathrm{THM}\), all meshes must be quadratic, but in the hydraulic case the integration is always linear (middle nodes are ignored).
In both cases, we assume that gravity is oriented according to negative \(z\).
On the other hand, it is assumed that the nonlinearities are low and that the coefficients \(N,\lambda ,\rho\) are constant. The pressure variations must therefore be small enough for \(N\) and \(\rho\) to be assumed to be constant.
In hydraulics, discretization will always be linear.
Linear discretization:
We will write:
\(p(z,t)=\sum _{i=1}^{2}{p}^{i}(t){\lambda }_{i}(z)\) (7)
With:
\(\{\begin{array}{c}{\lambda }_{1}=\frac{1}{2}-z\\ {\lambda }_{2}=\frac{1}{2}+z\end{array}\) (8)
By then introducing the matrices and vectors:
\(\{\begin{array}{c}[A]=[{A}_{\mathrm{ij}}]\text{;}{A}_{\mathrm{ij}}=\underset{-1/2}{\overset{1/2}{\int }}{\lambda }_{i}{\lambda }_{j}\mathrm{dz}\\ [B]=[{B}_{\mathrm{ij}}]\text{;}{B}_{\mathrm{ij}}=\underset{-1/2}{\overset{1/2}{\int }}\frac{d{\lambda }_{i}}{\mathrm{dz}}\frac{d{\lambda }_{j}}{\mathrm{dz}}\mathrm{dz}\\ \{{F}_{g}\}=\{{F}_{\mathrm{gi}}\}\text{;}{F}_{\mathrm{gj}}=\underset{-1/2}{\overset{1/2}{\int }}\frac{d{\lambda }_{i}}{\mathrm{dz}}\mathrm{dz}\end{array}\) (9)
And noting:
\(\{{p}_{e}\}=\left\{\begin{array}{c}{P}_{e}^{1}\\ {p}_{e}^{2}\end{array}\right\}\text{;}\{{p}_{a}\}=\left\{\begin{array}{c}{P}_{a}^{1}\\ {p}_{a}^{2}\end{array}\right\}\) (10)
\(\{{M}_{e}^{\mathrm{ext}}\}=\left\{\begin{array}{c}{{M}_{e}^{\mathrm{ext}}}_{1}\\ {{M}_{e}^{\mathrm{ext}}}_{2}\end{array}\right\}\text{;}\{{M}_{a}^{\mathrm{ext}}\}=\left\{\begin{array}{c}{{M}_{a}^{\mathrm{ext}}}_{1}\\ {{M}_{a}^{\mathrm{ext}}}_{2}\end{array}\right\}\) (11)
Equations (6) become:
\(\{\begin{array}{c}\frac{{N}_{\mathrm{ee}}}{{\lambda }_{e}}[A]\left\{\frac{{\mathrm{dp}}_{e}}{\mathrm{dt}}\right\}+\frac{{N}_{\mathrm{ea}}}{{\lambda }_{e}}[A]\left\{\frac{{\mathrm{dp}}_{a}}{\mathrm{dt}}\right\}+[B]\{{p}_{e}\}={\rho }_{e}\{{F}_{g}\}-\frac{1}{{\lambda }_{e}{\rho }_{e}}\{{M}_{e}^{\mathrm{ext}}\}\\ \frac{{N}_{\mathrm{ae}}}{{\lambda }_{a}}[A]\left\{\frac{{\mathrm{dp}}_{e}}{\mathrm{dt}}\right\}+\frac{{N}_{\mathrm{aa}}}{{\lambda }_{a}}[A]\left\{\frac{{\mathrm{dp}}_{a}}{\mathrm{dt}}\right\}+[B]\{{p}_{a}\}={\rho }_{a}\{{F}_{g}\}-\frac{1}{{\lambda }_{a}{\rho }_{a}}\{{M}_{a}^{\mathrm{ext}}\}\end{array}\) (12)
The calculation of the matrices \([A]\) and \([B]\) and the vector \(\{F\}\) gives:
\([A]=\frac{1}{3}\left[\begin{array}{cc}1& 1/2\\ 1/2& 1\end{array}\right]\text{;}[B]=\left[\begin{array}{cc}1& -1\\ -1& 1\end{array}\right]\text{;}\{{F}_{g}\}=\left\{\begin{array}{c}-1\\ 1\end{array}\right\}\) (13)
We then define \(\{{v}_{1}\},\{{v}_{2}\}\) the eigenvectors of \({[A]}^{-1}[B]\).
We have the properties:
\({\{{v}_{i}\}}^{T}[A]\{{v}_{j}\}={\{{v}_{i}\}}^{T}[B]\{{v}_{j}\}=0\text{}\mathrm{si}\text{}i\ne j\) (14)
And we ask:
\({a}_{i}={\{{v}_{i}\}}^{T}[A]\{{v}_{i}\}\text{,}{b}_{i}={\{{v}_{i}\}}^{T}[B]\{{v}_{i}\}\text{,}{f}_{i}={\{{v}_{i}\}}^{T}\{{F}_{g}\}\mathrm{et}{M}^{i}={\{{v}_{i}\}}^{T}\{{M}^{\mathrm{ext}}\}\) (15)
We find:
\(\{{v}_{1}\}=\{\begin{array}{c}1\\ 1\end{array}\}\text{;}\{{v}_{2}\}=\{\begin{array}{c}-1\\ 1\end{array}\}\) (16)
\(\{\begin{array}{c}{a}_{1}=1\text{;}{b}_{1}=0\text{;}{f}_{1}=0\\ {a}_{2}=\frac{1}{3}\text{;}{b}_{2}=4\text{;}{f}_{2}=-\mathrm{2g}\end{array}\) (17)
We then break down \(\{{p}_{e}\}\) and \(\{{p}_{a}\}\) on the basis of \(\{{v}_{i}\}\)
\(\{{p}_{e}\}=\sum _{i=1}^{2}{\alpha }_{e}^{i}\{{v}_{i}\}\text{;}\{{p}_{a}\}=\sum _{i=1}^{2}{\alpha }_{a}^{i}\{{v}_{i}\}\) (18)
Taking into account the properties of orthogonality (14), the system of equations (12) is written:
\(\{\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)
Posing:
\(\{{\alpha }^{i}\}=\{\begin{array}{c}{\alpha }_{e}^{i}\\ {\alpha }_{a}^{i}\end{array}\}\text{;}[N]=\left[\begin{array}{cc}{N}_{\mathrm{ee}}& {N}_{\mathrm{ea}}\\ {N}_{\mathrm{ae}}& {N}_{\mathrm{aa}}\end{array}\right]\text{;}[L]=\left[\begin{array}{cc}{\lambda }_{e}& 0\\ 0& {\lambda }_{a}\end{array}\right]\) (20)
Equation (19) is written as:
\([N]=\left[\frac{d{\alpha }^{i}}{\mathrm{dt}}\right]+\frac{{b}_{i}}{{a}_{i}}[L]\{{\alpha }^{i}\}=\frac{{f}_{i}}{{a}_{i}}\left\{\begin{array}{c}{\rho }_{e}{\lambda }_{e}\\ {\rho }_{a}{\lambda }_{a}\end{array}\right\}-\left\{\begin{array}{c}{M}_{e}^{i}/{\rho }_{e}{a}_{i}\\ {M}_{a}^{i}/{\rho }_{a}{a}_{i}\end{array}\right\}\) (21)
Initial conditions
It is assumed that:
\(\begin{array}{c}{p}_{a}(x,t=0)={p}_{a}^{0}\\ {p}_{e}(x,t=0)={p}_{a}^{0}-{p}_{c}^{0}\end{array}\) uniforms in space;
Given the values of vectors \(\{{v}_{1}\},\{{v}_{2}\}\) (equations (16)), it is easy to see that:
\(\begin{array}{c}{\alpha }_{a}^{1}(t=0)={p}_{a}^{0}\text{;}{\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)
We are in a case where the hydraulic equations are decoupled \(({N}_{\mathrm{ea}}={N}_{\mathrm{ae}}=0)\) and in which the fluid flows are zero \((\{{M}_{e}^{\mathrm{ext}}\}=\{{M}_{a}^{\mathrm{ext}}\}=0)\).
Taking into account (21), \({f}_{1}={f}_{3}=0\) (equations (17)), the system of equations (21) has the following solution:
\(\left\{\begin{array}{c}{\alpha }_{e}^{1}={P}_{a}^{0}-{p}_{c}^{0}\\ {\alpha }_{e}^{2}=\frac{{f}_{2}}{{b}_{2}}{\rho }_{e}(1-\mathrm{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-\mathrm{exp}(-\frac{{b}_{2}}{{a}_{2}}\frac{{\lambda }_{a}}{{N}_{\mathrm{aa}}}t))\end{array}\right\}\) (23)
Returning to the nodal variables, we find:
\(\left\{\begin{array}{c}{P}_{1}\\ {p}_{2}\end{array}\right\}=\left\{\begin{array}{c}{\alpha }_{1}-{\alpha }_{2}\\ {\alpha }_{1}+{\alpha }_{2}\end{array}\right\}\)
\({\left\{\begin{array}{c}{P}_{1}\\ {p}_{2}\end{array}\right\}}_{\mathrm{eau}}=\left\{\begin{array}{c}{P}_{a}^{0}-{p}_{c}^{0}+\frac{{\rho }_{e}g}{2}(1-\mathrm{exp}(-12\frac{{\lambda }_{e}}{{N}_{\mathrm{ee}}}t))\\ {P}_{a}^{0}-{p}_{c}^{0}-\frac{{\rho }_{e}g}{2}(1-\mathrm{exp}(-12\frac{{\lambda }_{e}}{{N}_{\mathrm{ee}}}t))\end{array}\right\}\) (24)
\({\left\{\begin{array}{c}{P}_{1}\\ {p}_{2}\end{array}\right\}}_{\mathrm{air}}=\left\{\begin{array}{c}{P}_{a}^{0}+\frac{{\rho }_{a}g}{2}\mathrm{exp}(-12\frac{{\lambda }_{a}}{{N}_{\mathrm{aa}}}t)\\ {P}_{a}^{0}-\frac{{\rho }_{a}g}{2}\mathrm{exp}(-12\frac{{\lambda }_{a}}{{N}_{\mathrm{aa}}}t)\end{array}\right\}\) (25)
and the capillary pressure, defined as the difference between air pressure and water pressure, has the value:
\({\left\{\begin{array}{c}{P}_{1}\\ {p}_{2}\end{array}\right\}}_{\mathrm{capillaire}}={\left\{\begin{array}{c}{P}_{1}\\ {p}_{2}\end{array}\right\}}_{\mathrm{air}}-{\left\{\begin{array}{c}{P}_{1}\\ {p}_{2}\end{array}\right\}}_{\mathrm{eau}}\)
We considered the following calculation case:
\(\begin{array}{c}S\ne 1\text{;}\frac{\partial S}{\partial {p}_{c}}=0\text{;}{K}_{s}=\infty \\ {N}_{\mathrm{ee}}=S\frac{\varphi }{{K}_{e}}\text{;}{N}_{\mathrm{aa}}=(1-S)\frac{\varphi }{{p}_{a}}\end{array}\)
2.2. Reference quantity#
Evolution of capillary pressure and dry air pressure as a function of time at the points
\(C,D(z=h)\)
\(A,B(z=0)\)
For quadratic discretization, check the constant value of the pressure at nodes \(E,F(z=\frac{h}{2})\).
2.3. Uncertainties#
Analytical solution on hydraulic equations therefore negligible uncertainties for models A, B, C.
Attention these analytical solutions do not apply to selective or lumped models (D and E). Indeed, in the latter case, the integrations are made at the nodes and no longer at the Gauss points. In fact, the integration by Gauss point is exact in \(\mathrm{1D}\) for polynomials of degrees less than or equal to 3 and therefore for all the integrals presented in equation (9). On the other hand, the vertex integration method is only accurate for polynomials of degree 1. So we can see that the terms in matrix \([A]\) will be under-integrated. It is therefore logical that on a unit mesh such as here the results obtained here are not exact. However, these tests are retained but with a « non-regression » result.
2.4. Bibliographical references#
Thermo-hydro-mechanics of porous media in Code_Aster — Note EDF, HI-74/99/011/A