[GeoMechanicsApplication] Partially saturated flow does not work as expected #12914

mnabideltares · 2024-12-09T09:58:55Z

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,

…est for the errors.

…mall_strain_upw.

…try loops.

…istent.

…ltiphysics/Kratos into geo/partial_saturation2

… test. New results.

…ltiphysics/Kratos into geo/partial_saturation2

markelov208

Hi Mohamed and Wijtze-Pieter, the code looks very professional! Sorry I missed what was done to fix the problem. Based on the recent discussion, a higher order shall be used; however, I do not see it here. Does it mean that this PR's goal is to make the code professional? Many thanks.

markelov208 · 2024-12-13T13:56:22Z

applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp

@@ -283,8 +283,6 @@ void UPwSmallStrainElement<TDim, TNumNodes>::InitializeSolutionStep(const Proces
        noalias(Variables.StressVector) = mStressVector[GPoint];
        ConstitutiveParameters.SetStressVector(Variables.StressVector);
        mConstitutiveLawVector[GPoint]->InitializeMaterialResponseCauchy(ConstitutiveParameters);
-
-        mRetentionLawVector[GPoint]->InitializeSolutionStep(RetentionParameters);


RetentionParameters is not used any longer.

markelov208 · 2024-12-13T13:59:22Z

applications/GeoMechanicsApplication/custom_elements/U_Pw_base_element.cpp

@@ -135,10 +135,8 @@ void UPwBaseElement::Initialize(const ProcessInfo& rCurrentProcessInfo)
    }

    mRetentionLawVector.resize(number_of_integration_points);
-    for (unsigned int i = 0; i < mRetentionLawVector.size(); ++i) {
-        mRetentionLawVector[i] = RetentionLawFactory::Clone(r_properties);
-        mRetentionLawVector[i]->InitializeMaterial(


As far as I understand, InitializeMaterial and other retention law functions were doing nothing after certain changes. Nice to remove them finally.

markelov208 · 2024-12-13T21:50:03Z

applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_interface_element.cpp

@@ -206,9 +206,6 @@ void UPwSmallStrainInterfaceElement<TDim, TNumNodes>::InitializeSolutionStep(con
        noalias(StressVector) = mStressVector[GPoint];
        ConstitutiveParameters.SetStressVector(StressVector);
        mConstitutiveLawVector[GPoint]->InitializeMaterialResponseCauchy(ConstitutiveParameters);
-
-        // Initialize retention law
-        mRetentionLawVector[GPoint]->InitializeSolutionStep(RetentionParameters);


RetentionParameters should be removed as well.

markelov208 · 2024-12-13T21:50:49Z

applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_interface_element.cpp

@@ -279,9 +276,6 @@ void UPwSmallStrainInterfaceElement<TDim, TNumNodes>::FinalizeSolutionStep(const

        mStateVariablesFinalized[GPoint] =
            mConstitutiveLawVector[GPoint]->GetValue(STATE_VARIABLES, mStateVariablesFinalized[GPoint]);
-
-        // retention law
-        mRetentionLawVector[GPoint]->FinalizeSolutionStep(RetentionParameters);


Sorry, again about RetentionParameters

markelov208 · 2024-12-13T22:32:43Z

applications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.cpp

+        const double p0 = r_geom[0].FastGetSolutionStepValue(WATER_PRESSURE);
+        const double p1 = r_geom[1].FastGetSolutionStepValue(WATER_PRESSURE);
+        const double p2 = r_geom[2].FastGetSolutionStepValue(WATER_PRESSURE);


There are many assignments for pressure values. What do you think on using
Vector pressure(3);
std::transform(r_geom.begin(), r_geom.begin()+3, pressure.begin(), [](const auto& node) {
return node.FastGetSolutionStepValue(WATER_PRESSURE);
});
It can be a separate function that returns pressure vector of given size.

Tried, see the version on github now.
Eventually, I would like to replace the whole summation below with shapefunctions of the pressure geometry multiplied by the p gotten from this pressure geometry. That is what it should be, there is no real need to write out this interpolation in the way it is done here. ( This is only very error prone )

markelov208 · 2024-12-13T23:22:20Z

...ly_saturated/test_saturated_below_phreatic_level_pw_triangle3N/ProjectParameters_stage2.json

+    "domain_size"                      : 2,
+    "model_import_settings": {
+      "input_type"                     : "mdpa",
+      "input_filename"                 : "kolom"


Should the file name be column?

markelov208 · 2024-12-13T23:23:35Z

...y_saturated/test_saturated_below_phreatic_level_pw_triangle6N/MaterialParameters_stage1.json

+          "K0_MAIN_DIRECTION": 1,
+          "K0_NC": 0.6,


markelov208 · 2024-12-13T23:23:53Z

...y_saturated/test_saturated_below_phreatic_level_pw_triangle6N/MaterialParameters_stage2.json

+          "K0_MAIN_DIRECTION": 1,
+          "K0_NC": 0.6,


markelov208 · 2024-12-13T23:24:13Z

...ly_saturated/test_saturated_below_phreatic_level_pw_triangle6N/ProjectParameters_stage1.json

+    "domain_size"                      : 2,
+    "model_import_settings": {
+      "input_type"                     : "mdpa",
+      "input_filename"                 : "kolom"


markelov208 · 2024-12-13T23:24:44Z

...ly_saturated/test_saturated_below_phreatic_level_pw_triangle6N/ProjectParameters_stage2.json

+    "domain_size"                      : 2,
+    "model_import_settings": {
+      "input_type"                     : "mdpa",
+      "input_filename"                 : "kolom"


…ions

…ion2' into geo/partial_saturation2

WPK4FEM and others added 6 commits December 6, 2024 10:16

Avoid zero divisions + a lot of printing

93e1263

Removed an unused variable and some commented out code

6fcc7fa

Fix effective saturation. Commented print statements.

37758a7

start of unit tests for SaturatedBelowPhreaticLevelLaw

409851f

Improved error messages for saturated below phreatic level and unit t…

bf5bd8d

…est for the errors.

removed unwanted include

188099f

mnabideltares requested a review from a team as a code owner December 9, 2024 09:58

mnabideltares linked an issue Dec 9, 2024 that may be closed by this pull request

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

Closed

mnabideltares assigned mnabideltares and WPK4FEM Dec 9, 2024

mnabideltares added the Bug label Dec 9, 2024

WPK4FEM and others added 9 commits December 9, 2024 11:17

ternary operator

7e10e65

Unit test for saturated law

85cb890

ternary operator

d31ed60

Checks and unit tests for Van Genuchten added.

c4412d9

Fixes for partially saturated flow

60eb250

Attempt to only place includes there where they are used.

e203f4c

Cleanup attempts

124346e

Trying to get rid of empty functions.

62f388d

Remove added print statements from core

28ff169

WPK4FEM removed the request for review from a team December 9, 2024 16:07

WPK4FEM added 4 commits December 10, 2024 11:46

sonar cloud remarks and naming style corrections

6574b78

Coupling matrix computation without potential zero division also in s…

c5871a2

…mall_strain_upw.

Also avoid zero division in RHS setup of coupling terms.

ef05087

van Genuchten gl >= 0, clean unused parameters from tests, some geome…

2a98e5c

…try loops.

mnabideltares linked an issue Dec 11, 2024 that may be closed by this pull request

[GeoMechanicsApplication] Fix Partially saturated flow for Pw element #12929

Open

WPK4FEM and others added 4 commits December 11, 2024 17:09

IgnoreUndrained setup of LHS and RHS made symmetric.

83debf8

No water solution intended in these tests, made IGNORE_UNDRAINED cons…

04ce103

…istent.

Added test cases for Pw 2D3N and 2D6N, with documentation

dd719eb

Merge branch 'geo/partial_saturation2' of https://github.com/KratosMu…

021cc95

…ltiphysics/Kratos into geo/partial_saturation2

fix for Readme

0a1cdba

mnabideltares linked an issue Dec 12, 2024 that may be closed by this pull request

[GeoMechanicsApplication] Fix Partially saturated flow for triangular U-Pw diff order element #12931

Open

A test xase for U-Pw diff order 2D6N is added

03c1ba7

mnabideltares linked an issue Dec 12, 2024 that may be closed by this pull request

[GeoMechanicsApplication] Fix Partially saturated flow for U-Pw small strain element #12930

Open

WPK4FEM and others added 5 commits December 12, 2024 15:52

fix of revision 60eb250 changed results of transient groundwater flow…

fb3b09e

… test. New results.

2 tests are added for U-Pw small strain element

eb14f80

Merge branch 'geo/partial_saturation2' of https://github.com/KratosMu…

5825e04

…ltiphysics/Kratos into geo/partial_saturation2

fixes

437df1f

fix

bd43cc4

WPK4FEM requested review from avdg81, rfaasse and markelov208 December 13, 2024 09:49

markelov208 requested changes Dec 13, 2024

View reviewed changes

WPK4FEM and others added 9 commits December 16, 2024 09:24

Sonar Cloud unused stuff complaint resolved

162c3f4

Removed unused retention parameters ( review Gennady )

f7c0261

std::transform for pressure getting ( review comment of Gennady )

be608a3

Removed unused K0 things ( review comment Gennady )

83d65fe

Renamed kolom to the english column ( review comment Gennady )

35b9872

avoid code repetition

27b69bc

kp format on request of Anne

856256e

Set material parameters in a common directory, to avoid file duplicat…

5e78ddd

…ions

Merge remote-tracking branch 'refs/remotes/origin/geo/partial_saturat…

eec6c84

…ion2' into geo/partial_saturation2

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

[GeoMechanicsApplication] Partially saturated flow does not work as expected #12914

[GeoMechanicsApplication] Partially saturated flow does not work as expected #12914

mnabideltares commented Dec 9, 2024

markelov208 left a comment

markelov208 Dec 13, 2024

WPK4FEM Dec 17, 2024

markelov208 Dec 13, 2024

WPK4FEM Dec 17, 2024

markelov208 Dec 13, 2024

WPK4FEM Dec 17, 2024

markelov208 Dec 13, 2024

WPK4FEM Dec 17, 2024

markelov208 Dec 13, 2024

WPK4FEM Dec 17, 2024

markelov208 Dec 13, 2024

WPK4FEM Dec 17, 2024

markelov208 Dec 13, 2024

WPK4FEM Dec 17, 2024

markelov208 Dec 13, 2024

WPK4FEM Dec 17, 2024

markelov208 Dec 13, 2024

WPK4FEM Dec 17, 2024

markelov208 Dec 13, 2024

WPK4FEM Dec 17, 2024

[GeoMechanicsApplication] Partially saturated flow does not work as expected #12914

Are you sure you want to change the base?

[GeoMechanicsApplication] Partially saturated flow does not work as expected #12914

Conversation

mnabideltares commented Dec 9, 2024

This issue is time-boxed to 10 developer days.

Description

markelov208 left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment