9. Spatial discretization

Finite elements THM are mixed elements, in the sense that they have both unknowns of movements, pressures, and temperatures. A discretization choice where displacements, pressures and temperatures are interpolated with the same order of approximation leads to oscillations, especially for choices of time steps that are too small compared to discretization in space We will consult on this subject for example [bib10]. This problem is also related to the way to calculate the so-called mass matrix, and we can consult [bib14] on this subject. In addition, to illustrate our point, we give in the appendix the solution for the first time step of a one-dimensional consolidation problem with interpolation \(\mathit{P1P1}\). We can see that for a small step of time, it is very oscillating.

For this reason, quadratic elements THM are elements in \(\mathit{P2P1}\), i.e. the interpolation of the displacements is quadratic and that of the temperatures and pressures is linear. However, we have kept all the unknowns about all the nodes, including the middle nodes, but we have imposed in the calculation of the stiffness matrices that the pressure of a middle segment node be equal to the half sum of the vertex nodes of the segment to which it belongs.

In addition, in programming, we took into account the following property:

Let \(s\) be a vertex node, \({w}_{s}^{1}\) its shape function as belonging to a linear element (for example QUAD4), and \({w}_{s}^{2}\) its shape function as belonging to a quadratic element (for example QUAD8). Let \(\mathrm{na}\) be the number of edges having \(s\) as an end and \({w}_{\text{ma}}^{2}\) the shape function of the quadratic interpolation attached to a mid-edge node, we then have the relationship:

\({w}_{s}^{1}={w}_{s}^{2}+\sum _{\text{ma}=1}^{\text{na}}\frac{{w}_{\text{ma}}^{2}}{2}\)

That said, and including in interpolation \(\mathit{P2P1}\), non-oscillation conditions exist on the time step. [bib10] give the relationship:

\(\Delta t >\frac{\Delta {x}^{2}}{20{C}_{v}}\)

where \({C}_{v}\) is the consolidation coefficient: \({C}_{v}=\frac{\text{kE}(1-\mathrm{\nu })}{{\mathrm{\rho }}_{\text{lq}}(1+\mathrm{\nu })(1-\mathrm{2\nu })}\), \(k\) being the permeability measured in \(m/s\).