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


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

The reference solution is one-dimensional because it only depends on the vertical coordinate. We will therefore reason on a linear axis along z.

If we neglect gravity and consider that there is no external heat source, the thermal equation is given by the following expression:


.. math::
: label: EQ-None

    {h} _ {l} ^ {m}\ frac {d {m} _ {l}} {\ mathrm {dt}} +\ frac {\ delta Q\ text {'}} {\ mathrm {dt}} {\ mathrm {dt}}}} D {m} _ {l} +\ mathrm {divq} =0

Where :math:`{h}_{l}^{m}` is the mass enthalpy of water, :math:`{m}_{l}` its mass, its mass, :math:`{M}_{l}` its flow, q the heat flow, and :math:`Q\text{'}` the non-convected heat.

The mass conservation equation is as follows:


.. math::
: label: EQ-None

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

In the absence of water flow, equation (1) is simplified to become:


.. math::
: label: EQ-None

    \ frac {\ delta Q\ text {'}}} {\ mathrm {dt}} +\ mathrm {divq} =0

In this equation, quantity :math:`Q\text{'}` represents the heat received by the system in a transformation for which there is no heat input per fluid inlet. A low pressure variation and parameters are selected so that the temperature varies little and therefore that:

:math:`\frac{dQ\text{'}}{T}=\frac{dQ\text{'}}{{T}_{0}}`

So :math:`\mathrm{dQ}\text{'}` has this expression:

:math:`\mathrm{dQ}\text{'}=3{\alpha }_{0}{K}_{0}{T}_{0}d{\varepsilon }_{v}+{C}_{0}^{\varepsilon }\mathrm{dT}-3{\alpha }_{l}^{m}{T}_{0}{\mathrm{dp}}_{l}`

with:


* :math:`{\alpha }_{0}` the homogenized thermal expansion coefficient equivalent to that of solid :math:`{\alpha }_{s}`.


* :math:`{K}_{0}` the drained elasticity coefficient.


* :math:`{C}_{0}^{\varepsilon }` heat has constant deformation which has the expression :math:`{C}_{0}^{\varepsilon }={C}_{0}^{\sigma }-9{T}_{0}{K}_{0}{\alpha }_{0}^{2}` and :math:`{C}_{0}^{\sigma }` the heat at constant stress.


* :math:`{\alpha }_{l}^{m}` The coefficient of relative thermal expansion of the liquid with respect to the skeleton, it has for expression: :math:`{\alpha }_{l}^{m}=(b-\phi ){\alpha }_{s}+\phi {\alpha }_{l}` with b the Biot coefficient, :math:`\phi` the porosity and :math:`{\alpha }_{l}` the coefficient of expansion of the liquid.


* :math:`d{\varepsilon }_{v}` The variation in volume deformation.


* :math:`{\mathrm{dp}}_{l}` The change in liquid pressure.

Moreover, the heat flow has the following expression: :math:`q=-{\lambda }_{T}\partial \frac{T}{\partial z}` where :math:`{\lambda }_{T}` is the thermal conductivity coefficient.

By replacing :math:`\mathrm{dQ}\text{'}` and q by their value in equation (3) and limiting ourselves to only time step :math:`\Delta t`, we obtain:

:math:`3{\alpha }_{0}{K}_{0}{T}_{0}\Delta {\varepsilon }_{v}-3{\alpha }_{l}^{m}{T}_{0}\Delta {p}_{l}+{C}_{0}^{\varepsilon }\Delta T=-\Delta tdq={\lambda }_{l}\Delta t\frac{{\partial }^{2}T}{\partial {z}^{2}}`

Considering temperatures and displacements that are initially zero, we can write:

:math:`\Delta T=T(t)-T(0)=T` and :math:`\Delta \varepsilon =\varepsilon (t)-\varepsilon (0)=\varepsilon` and :math:`\Delta {p}_{l}={p}_{l}(t)-{p}_{l}(0)={p}_{l}`

So at the first step of time, we will have:

:math:`T-a\frac{{\partial }^{2}T}{\partial {z}^{2}}=b`

with :math:`a=\frac{{\lambda }_{l}\Delta T}{{C}_{0}^{\varepsilon }}` and :math:`b=-\frac{3{\alpha }_{0}{K}_{0}{T}_{0}\varepsilon -3{\alpha }_{l}^{m}{T}_{0}{p}_{l}}{{C}_{0}^{\varepsilon }}`

The variational formulation of this expression (in the one-dimensional case) is then as follows:


.. math::
: label: EQ-None

    {\ int} _ {\ omega} T\ mathrm {.} \ stackrel {} {T}\ mathrm {dz} +a {\ int} +a {\ int} _ {\ omega}\ frac {\ partial\ stackrel {} {T} {T}} {T}} {\ partial z}} {\ partial z}} {\ partial z}\ mathrm {dz} -a {\ int} _ {\ partial\ Omega}\ frac {\ partial T} {\ partial T} {\ partial T} {\ partial T} {\ partial T} {\ partial z}}\ stackrel {} {T}\ mathrm {dz} = {\ int} _ {\ omega} b\ mathrm {.} \ stackrel {} {T}\ mathrm {dz}

To establish the analytical solution, we consider a single element of degree 1 since in modeling THM the thermo-hydraulic part is treated by linear elements.

Let it be a linear element:


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


.. _RefImage_Object_23.svg:

We recall the boundary conditions: :math:`\frac{\partial T}{\partial z}=0` at both ends (z=0.5 and z = -0.5)

The temperature based on form functions is written as follows:

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

with

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

and

: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 z}\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:`\{T\}=\left\{\begin{array}{c}{T}^{1}\\ {T}^{2}\end{array}\right\}`

and :math:`\{\frac{\partial T}{\partial t}\}=\left\{\begin{array}{c}\frac{\partial {T}^{1}}{\partial t}\\ \frac{\partial {T}^{2}}{\partial t}\end{array}\right\}`

Equation (4) then becomes

:math:`[A]\{T\}+a[B]\{T\}-a[A]\{\frac{\partial T}{\partial t}\}=[A]\{b\}`

with the imposed boundary conditions (zero displacement at the bottom and zero temperature flow at both ends), we have:

:math:`\{b\}=\left\{\begin{array}{c}0\\ b\end{array}\right\}\mathrm{et}\{\frac{\partial T}{\partial t}\}=\left\{\begin{array}{c}0\\ 0\end{array}\right\}`

Which in the end gives us

:math:`\left\{\begin{array}{c}{T}^{1}\\ {T}^{2}\end{array}\right\}+a{[A]}^{-1}[B]\left\{\begin{array}{c}{T}^{1}\\ {T}^{2}\end{array}\right\}=\left\{\begin{array}{c}0\\ b\end{array}\right\}`

In the end we get:

:math:`\left[\begin{array}{cc}1+\mathrm{6a}& -\mathrm{6a}\\ -\mathrm{6a}& 1+\mathrm{6a}\end{array}\right]\{T\}=\left\{\begin{array}{c}{T}^{1}\\ {T}^{2}\end{array}\right\}`


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

For a short time of 100s, we will have :math:`\left\{\begin{array}{c}{T}^{1}\\ {T}^{2}\end{array}\right\}=\left\{\begin{array}{c}{2.10}^{-14}\\ \mathrm{1,045}{.10}^{-7}\end{array}\right\}`

We will consider :math:`{T}^{1}` to be null.


Uncertainty about the solution
---------------------------

None