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{{\mathit{dm}}_{i}}{\mathit{dt}}+\text{Div}{M}_{i}\mathrm{=}0\) \(i\) varying from \(1\) to the number of components (1)
In our example, the model consists of a fluid: liquid water. Equation (1) therefore applies to this constituent:
\(\frac{{\mathrm{dm}}_{e}}{\mathrm{dt}}+\text{Div}{M}_{e}=0\) (2)
The expression of fluid flow is:
\({M}_{e}={\rho }_{e}{\lambda }_{e}(-\nabla {p}_{e}+{\rho }_{e}g)\) (3)
However, the mass supply of fluid is defined by equation (4) where the terms \({N}_{\mathrm{ee}}\) and \({N}_{\mathrm{ea}}\) (equation (5)) depend on the degree of saturation \(S\), on the porosity \(\phi\), on the Biot coefficient \(b\), on the permeability of the liquid \({K}_{e}\), on the permeability of the liquid and on the elasticity of the solid matrix \({K}_{s}\).
\(\frac{{\mathrm{dm}}_{e}}{\mathrm{dt}}={\rho }_{e}{N}_{\mathrm{ee}}\frac{{\mathrm{dp}}_{e}}{\mathrm{dt}}+{\rho }_{e}{N}_{\mathrm{ea}}\frac{{\mathrm{dp}}_{a}}{\mathrm{dt}}\) (4)
\(\{\begin{array}{c}{N}_{\mathrm{ee}}=-\varphi \frac{\partial S}{\partial {p}_{c}}+S(\frac{\varphi }{{K}_{e}}+\frac{b-\varphi }{{K}_{s}}S)\\ {N}_{\mathrm{ea}}={N}_{\mathrm{ae}}=\varphi \frac{\partial S}{\partial {p}_{c}}+(1-S)(\frac{b-\varphi }{{K}_{s}}S)\end{array}\) (5)
The material is saturated, \(S=1\) and \(\frac{\partial S}{\partial {P}_{c}}=0\). \(\Rightarrow {N}_{\mathrm{ee}}=S(\frac{\varphi }{{K}_{e}}+\frac{b-\varphi }{{K}_{s}}S)\) and \({N}_{\mathrm{ea}}=0\).
The variational formulation of equation (2), taking into account (3) and (4), is:
\(\forall {P}_{e}^{\text{*}}\) verifying the pressure boundary conditions:
\({\int }_{\Omega }{N}_{\mathrm{ee}}\frac{{\mathrm{dp}}_{e}}{\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{*}}\) (6)
Discretization
To calculate the analytical solution, we use a one-dimensional case and consider a discretization with a single element of degree 2 (HEXA20, DPQ8). It should be noted that since all modeling \(\mathrm{HM}\) is of type \(\mathrm{P2P1}\), even if the mesh is quadratic, hydraulic modeling is linear.
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.
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}\mathrm{[}A\mathrm{]}\mathrm{=}\mathrm{[}{A}_{\mathit{ij}}\mathrm{]}\text{;}{A}_{\mathit{ij}}\mathrm{=}\underset{\mathrm{-}1\mathrm{/}2}{\overset{1\mathrm{/}2}{\mathrm{\int }}}{\lambda }_{i}{\lambda }_{j}\mathit{dz}\\ \mathrm{[}B\mathrm{]}\mathrm{=}\mathrm{[}{B}_{\mathit{ij}}\mathrm{]}\text{;}{B}_{\mathit{ij}}\mathrm{=}\underset{\mathrm{-}1\mathrm{/}2}{\overset{1\mathrm{/}2}{\mathrm{\int }}}\frac{d{\lambda }_{i}}{\mathit{dz}}\frac{d{\lambda }_{j}}{\mathit{dz}}\mathit{dz}\\ \mathrm{\{}{F}_{g}\mathrm{\}}\mathrm{=}\mathrm{\{}{F}_{\mathit{gi}}\mathrm{\}}\text{;}{F}_{\mathit{gi}}\mathrm{=}\underset{\mathrm{-}1\mathrm{/}2}{\overset{1\mathrm{/}2}{\mathrm{\int }}}\frac{d{\lambda }_{i}}{\mathit{dz}}\mathit{dz}\end{array}\)
And noting:
\(\{{p}_{e}\}=\left\{\begin{array}{c}{P}_{e}^{1}\\ {p}_{e}^{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\}\) (11)
Equation (6) becomes:
\(\frac{{N}_{\mathit{ee}}}{{\lambda }_{e}}\mathrm{[}A\mathrm{]}\left\{\frac{{\mathit{dp}}_{e}}{\mathit{dt}}\right\}+\mathrm{[}B\mathrm{]}\mathrm{\{}{p}_{e}\mathrm{\}}\mathrm{=}{\rho }_{e}\mathrm{\{}{F}_{g}\mathrm{\}}\mathrm{-}\frac{1}{{\lambda }_{e}{\rho }_{e}}\mathrm{\{}{M}_{e}^{\mathit{ext}}\mathrm{\}}\) (12)
The calculation of the matrices \(\mathrm{[}A\mathrm{]}\) and \(\mathrm{[}B\mathrm{]}\) and the vector \(\mathrm{\{}{f}_{g}\mathrm{\}}\) gives:
\(\mathrm{[}A\mathrm{]}\mathrm{=}\frac{1}{3}\left[\begin{array}{cc}1& 1\mathrm{/}2\\ 1\mathrm{/}2& 1\end{array}\right]\text{;}\mathrm{[}B\mathrm{]}\mathrm{=}\left[\begin{array}{cc}1& \mathrm{-}1\\ \mathrm{-}1& 1\end{array}\right]\text{;}\mathrm{\{}F\mathrm{\}}\mathrm{=}\left\{\begin{array}{c}\mathrm{-}1\\ 1\end{array}\right\}\) (13)
We then define the eigenvectors \({\mathrm{[}A\mathrm{]}}^{\mathrm{-}1}\mathrm{[}B\mathrm{]}\) of: \(\{{v}_{1}\},\{{v}_{2}\}\) which have the following orthogonality 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}\}\) on the basis of \(\{{v}_{i}\}\):
\(\mathrm{\{}{p}_{e}\mathrm{\}}\mathrm{=}\mathrm{\sum }_{i\mathrm{=}1}^{2}{\alpha }_{e}^{i}\mathrm{\{}{v}_{i}\mathrm{\}}\) (18)
Taking into account the orthogonality properties (14), equation (12) is written:
\(\frac{{N}_{\mathit{ee}}}{{\lambda }_{e}}{a}_{i}\frac{d{\alpha }_{e}^{i}}{\mathit{dt}}+{b}_{i}{\alpha }_{e}^{i}\mathrm{=}{\rho }_{e}{f}_{i}\mathrm{-}\frac{1}{{\lambda }_{e}{\rho }_{e}}{M}_{e}^{i}\) (19)
Initial conditions
It is assumed that:
\({p}_{e}(x,t\mathrm{=}0)\mathrm{=}{p}_{a}^{0}\mathrm{-}{p}_{c}^{0}\) uniforms in space;
Given the values of vectors \(\{{v}_{1}\},\{{v}_{2}\}\) (equations (16)), it is easy to see that:
\(\mathrm{\{}\begin{array}{c}{\alpha }_{e}^{1}(t\mathrm{=}0)\mathrm{=}{P}_{a}^{0}\mathrm{-}{p}_{c}^{0}\\ {\alpha }_{e}^{2}(t\mathrm{=}0)\mathrm{=}0\end{array}\) (20)
This is a case where the fluid flow is zero (\(\mathrm{\{}{M}_{e}^{\mathit{ext}}\mathrm{\}}\mathrm{=}0\)).
Taking into account (20), \({f}_{1}=0\) (equations (17)), the solution of the system of equations (19) is:
\(\mathrm{\{}\begin{array}{c}{\alpha }_{e}^{1}\mathrm{=}{P}_{a}^{0}\mathrm{-}{p}_{c}^{0}\\ {\alpha }_{e}^{2}\mathrm{=}\frac{{f}_{2}}{{b}_{2}}{\rho }_{e}(1\mathrm{-}\mathrm{exp}(\mathrm{-}\frac{{b}_{2}}{{a}_{2}}\frac{{\lambda }_{e}}{{N}_{\mathit{ee}}}t))\end{array}\) (21)
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\}}_{\mathit{eau}}\mathrm{=}\left\{\begin{array}{c}{P}_{a}^{0}\mathrm{-}{p}_{c}^{0}+\frac{{\rho }_{e}g}{2}(1\mathrm{-}\mathrm{exp}(\mathrm{-}12\frac{{\lambda }_{e}}{{N}_{\mathit{ee}}}t))\\ {P}_{a}^{0}\mathrm{-}{p}_{c}^{0}\mathrm{-}\frac{{\rho }_{e}g}{2}(1\mathrm{-}\mathrm{exp}(\mathrm{-}12\frac{{\lambda }_{e}}{{N}_{\mathit{ee}}}t))\end{array}\right\}\) (22)
2.2. Reference quantity#
Evolution of capillary pressure as a function of time at the points:
\(C,D\) \((z\mathrm{=}h)\)
\(A,B\) \((z\mathrm{=}0)\)
For quadratic discretization: Verification of the constant value of the pressure at nodes \(E,F(z\mathrm{=}\frac{h}{2})\).
2.3. Uncertainties#
Analytical solution on the hydraulic equation, the uncertainties are therefore negligible.