Benchmark solution
=====================


Calculation method
-----------------

The reference solution is one-dimensional because it only depends on the vertical coordinate.

The water mass conservation equation :math:`{m}_{l}` is given by the following expression:

.. math::
: label: EQ-None

    \ frac {d {m} _ {l}}} {\ mathrm {dt}} +D {M} _ {l} =0

with :math:`{M}_{l}` the flow of water such that, if we neglect gravity:

:math:`{M}_{l}={\rho }_{l}{\lambda }_{l}(-\nabla {p}_{l})` with :math:`{p}_{l}` the liquid pressure, :math:`{\rho }_{l}` its density, and :math:`{\lambda }_{l}` its hydraulic conductivity such as :math:`{\lambda }_{l}=\frac{{K}_{i}}{\mu }`. :math:`{K}_{i}` is the intrinsic permeability and :math:`\mu` the viscosity of the liquid.

Considering an undeformable solid, we can write that

:math:`\frac{d{m}_{l}}{\mathrm{dt}}={\rho }_{l}\frac{\phi }{{K}_{l}}\frac{d{p}_{l}}{\mathrm{dt}}` with :math:`\phi` the porosity and :math:`{K}_{l}` the compressibility of water. We note :math:`N=\frac{\phi }{{K}_{l}}`.

Equation (1) then becomes:

:math:`{\rho }_{l}N\frac{d{p}_{l}}{\mathrm{dt}}+D({\rho }_{l}{\lambda }_{l}(-\nabla {p}_{l}))=0`

whose variational formulation is as follows:

.. math::
: label: EQ-None

    {\ int} _ {\ omega} N\ frac {d {p} _ {l}} {\ mathrm {dt}}\ mathrm {.} \ stackrel {} {{p} _ {l}} + {\ int}} + {\ int}} _ {\ omega} {\ lambda} _ {l}\nabla\ stackrel {}} {{p}} {{p}}} {{p}} _ {l}} =- {\ int}} =- {\ int}} _ {\ partial\ omega}}\ frac {{M}} ^ {\ mathrm {ext}}} {{\ rho} _ {l}}\ stackrel {} {{p} _ {l}}

with :math:`{M}^{\mathrm{ext}}` the exterior load.

To establish the solution, we use the one-dimensional case and we adopt a discretization corresponding to a single element of degree 1 since in modeling THM the hydraulic part is treated on linear elements. It is also assumed that the nonlinearities are low in this case and that the coefficients :math:`N` and :math:`{\rho }_{l}` are constant, which implies a relatively small variation in pressure.

Let it be a linear element:


.. image:: images/Object_103.svg
    :width: 489
    :height: 206


.. _RefImage_Object_103.svg:

We then write the pressure based on form functions such as:

:math:`{p}_{l}(z,t)=\sum _{i=1}^{i=2}{p}^{i}(t){\lambda }_{i}(z)`

with

:math:`{\lambda }_{1}(z)=\mathrm{0,5}(1+\mathrm{2z})`

:math:`{\lambda }_{2}(z)=\mathrm{0,5}(1-\mathrm{2z})`

The following matrices are then introduced:

:math:`\left[A\right]=\left[{A}_{\mathrm{ij}}\right];{A}_{\mathrm{ij}}={\int }_{-\mathrm{0,5}}^{\mathrm{0,5}}{\lambda }_{i}{\lambda }_{j}\mathrm{dz}`

:math:`\left[B\right]=\left[{B}_{\mathrm{ij}}\right];{B}_{\mathrm{ij}}={\int }_{-\mathrm{0,5}}^{\mathrm{0,5}}\frac{\partial {\lambda }_{i}}{\partial z}\frac{\partial {\lambda }_{j}}{\partial \mathrm{dz}}\mathrm{dz}`

Which leads to

:math:`\left[A\right]=\frac{1}{6}\left[\begin{array}{cc}2& 1\\ 1& 2\end{array}\right]` and :math:`\left[B\right]=\left[\begin{array}{cc}1& -1\\ -1& 1\end{array}\right]`

We then note classically:

:math:`\{{p}_{l}\}=\left\{\begin{array}{c}{P}_{l}^{1}\\ {P}_{l}^{2}\end{array}\right\}` :math:`\{{M}^{\mathrm{ext}}\}=\left\{\begin{array}{c}{{M}^{\mathrm{ext}}}^{1}\\ {{M}^{\mathrm{ext}}}^{2}\end{array}\right\}`

Knowing that we are injecting a flow of water Q on the upper face, we therefore have

:math:`\{{M}^{\mathrm{ext}}\}=\left\{\begin{array}{c}0\\ -Q\end{array}\right\}`

Equation (2) becomes :math:`\frac{N}{{\lambda }_{l}}\left[A\right]\{\frac{{\mathrm{dp}}_{l}}{\mathrm{dt}}\}+\left[B\right]\{{p}_{l}\}=\frac{-1}{{\lambda }_{l}\rho }\{{M}^{\mathrm{ext}}\}`

Considering that for short time variations Dt, the evolution of p is almost linear, we can write that

:math:`\{\frac{{\mathrm{dp}}_{l}}{\mathrm{dt}}\}=\frac{1}{\Delta t}\left\{\begin{array}{c}{P}_{l}^{1}-{p}_{0}^{1}\\ {P}_{l}^{2}-{p}_{0}^{2}\end{array}\right\}`

and since we start here from an initial state of zero pressure, we will have

:math:`\{\frac{{\mathrm{dp}}_{l}}{\mathrm{dt}}\}=\frac{1}{\Delta t}\left\{\begin{array}{c}{P}_{l}^{1}\\ {P}_{l}^{2}\end{array}\right\}`

Finally, we obtain the system of two equations with two unknowns:

:math:`{p}_{l}^{1}(\frac{N}{{\lambda }_{l}\Delta t}+6)-{\mathrm{6p}}_{l}^{2}=\frac{-\mathrm{2Q}}{{\lambda }_{l}{\rho }_{l}}`

:math:`{p}_{l}^{2}(\frac{N}{{\lambda }_{l}\Delta t}+6)-{\mathrm{6p}}_{l}^{1}=\frac{-\mathrm{4Q}}{{\lambda }_{l}{\rho }_{l}}`

We then have the following result:

:math:`\left\{\begin{array}{c}{P}_{l}^{1}\\ {P}_{l}^{2}\end{array}\right\}={\left[\begin{array}{cc}\frac{N}{{\lambda }_{l}\Delta t}+6& -6\\ -6& \frac{N}{{\lambda }_{l}\Delta t}+6\end{array}\right]}^{-1}\mathrm{.}\left\{\begin{array}{c}\frac{-\mathrm{2Q}}{{\lambda }_{l}{\rho }_{l}}\\ \frac{-\mathrm{4Q}}{{\lambda }_{l}{\rho }_{l}}\end{array}\right\}`


Reference quantities and results
-----------------------------------

The value of the pore pressure on the lower and upper faces of the mass is tested.


Uncertainties about the solution
----------------------------

No uncertainty about the solution. The solution is analytical.


Bibliographical references
---------------------------

1) *Mechanics of porous environments*. O. Coussy, Editions Technip, 2000.