KratosMultiphysics / Kratos

Kratos Multiphysics (A.K.A Kratos) is a framework for building parallel multi-disciplinary simulation software. Modularity, extensibility and HPC are the main objectives. Kratos has BSD license and is written in C++ with extensive Python interface.
https://kratosmultiphysics.github.io/Kratos/
Other
1.05k stars 246 forks source link

{GeoMechanicsApplication] INVESTIGATION - Partially saturated flow does not work as expected #12842

Open mnabideltares opened 2 weeks ago

mnabideltares commented 2 weeks ago

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 kr 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.

For some reason we are still getting results from this configuration., but the results are incorrect.

Image

Image

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.

Image

Image

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,