[GeoMechanicsApplication] INVESTIGATION - Partially saturated flow does not work as expected

[mnabideltares]

This issue is time-boxed to 10 developer days.

22/10/2024 - JDN - Comment
I need more information to able to prioritize this for this sprint or beyond. i.e. whats the urgency ? who is it for ? also needs refinement ? I wasn't aware we were using partial saturation anywhere and we are currently awaiting for 2025 to do so, but do I interpret correctly that the mechanism also fails when the zone above the phreatic line is fully unsaturated

Description

The Partially saturated flow is implemented but it has been never tested. Moreover, instability arrised when we tried to apply it.

The equation for partially saturated flow is:

$$\left[ (\alpha + n \beta) \rho^w S + n \rho^w \frac{d S}{dp} \right] \frac{\partial p}{\partial t} = \frac{\partial}{\partial x_i} \left[ \frac{\rho^w k_r K_{ij}}{\mu} \left( \frac{\partial p}{\partial x_j} - \rho^w g_j \right) \right]$$

where

$$k_r = S_e^{g_l} \left[ 1 - \left( 1 -S_e^{1/g_m} \right)^{g_m} \right]$$ $$S_e = \frac{S - S_r}{S_s - S_r}$$

Here we start from a simple test case, with saturated flow below the phreatic line:
In the case of $S_r = 0$, the pressure above the phreatic line needs to be zero. However, from this equation, the left hand side includes the saturation $S$ which is zero. It leads to a zero $C$ matrix. Morover, $k_r = 0$ too, which leads to a zero $K$ matrix. Therefore, a zero matrix is added which makes the diagonal of the generic matrix partly to be zero, and it lead to a singular matrix.
This can be probably avoided to adding "BIOT_COEFFICIENT" to the material parameters to set the $C$ matrix manually as a user defined parameter.

In the case of $S_r > 0$, then $K_r =0$`but the left hand side gets a nonzero value. it means $\frac{dp}{dt} = 0$. This is proposed to lead to unchanged pressure above the phreatic line.

However, here we see a change in pressure. It indicates there is driving force, which leads the water pressure goes up. Which is an unexpected behaviour,

[mnabideltares]

Case 1: case01.zip

The solution of this test case is simplified to include only Pw element. At the first stage, steady state Pw element is forced to reach to an equilibrium state, which a specified pheriatic line at the level of $y=-2$. At the second stage, the pressure is solved. However, flux appears at the pheriatic line level and the water flows to the top. The grid is kept to be quadraterial.

[mnabideltares]

Case 2: case02.zip

The solution of this test case is simplified to include only Pw element. At the first stage, steady state Pw element is forced to reach to an equilibrium state, which a specified pheriatic line at the level of $y=-2$. At the second stage, the pressure is solved. The grid is kept to be triangular 3 noded element.

The pressure works correctly here, and the flux is zero at the pheriatic line.

[mnabideltares]

Case 3: case03.zip

The solution of this test case is simplified to include only Pw element. At the first stage, steady state Pw element is forced to reach to an equilibrium state, which a specified pheriatic line at the level of $y=-2$. At the second stage, the pressure is solved. The grid is set to be triangular 6 noded element.

On high order element the solution goes wrong and flux appears toward to the top of the column.

[mnabideltares]

Case 4: case04.zip

The element here is selected to be SmallStrainUPwDiffOrderElement. This test case is solved on a mesh with high order triangular element, namely 6 noded. The deformation is solved on a 6 noded element while the pressure on 3 noded triangle.

As the pressure is solved on a linear element, the solution behaves correctly.

[mnabideltares]

Update:
A minor change in the code is implemented. Minor in programming, but it may have big effect on the whole results. The number of ghost points are increased for some pressure elements.

In "U_Pw_base_element.cpp" line 206 the intergration metod is changed from 2 to 3.

default:
    GI_GAUSS = GeometryData::IntegrationMethod::GI_GAUSS_3;
    break;

This fixed issues in the quadirateral elements.

[mnabideltares]

Case 5:
For file input, refer to "Case 1"

This case is exactly as case 1, with the new update of mentioned above

The pressure works correctly here, and the flux is zero at the pheriatic line. It is achieved by adding the number of gauss points for this element (quad 4 noded).

[mnabideltares]

case 6: case06.zip

Similar to case 4, but with quadraterial elements with 8 nodes.
The update mentioned above is implemented, meaning the number of integration points are increased.

[mnabideltares]

A final case need to be investigated. It is Pw element on 6 noded triangle.
I added added a modification in "U_Pw_base_element.cpp" in line 197:

case 6:
    GI_GAUSS = GeometryData::IntegrationMethod::GI_GAUSS_3;
    break;

previousely was GI_GAUSS_2.
However, this modification did not lead to correct results. The results look as below:

The pressure field does not even look symmetric.