Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Merged
merged 53 commits into from
Dec 21, 2024

Conversation

mnabideltares
Copy link
Contributor

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,

Copy link
Contributor

@markelov208 markelov208 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Contributor

@avdg81 avdg81 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Gentlemen, you two have done a lot of work to clean up the code and to fix a few issues. I'm very satisfied with the result. It's yet another step towards code that is easier to understand and more maintainable.

I do have quite a few suggestions and some questions, but none of them is blocking in my opinion.

else if (rVariable == RELATIVE_PERMEABILITY)
rOutput[GPoint] =
mRetentionLawVector[GPoint]->CalculateRelativePermeability(RetentionParameters);
rOutput[GPoint] = mRetentionLawVector[GPoint]->CalculateValue(RetentionParameters, rVariable, rOutput[GPoint]);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I have noticed that the switch that has been taken out here, is repeated for each retention law (in their overrides of CalculateValue). We shouldn't try to address that in this PR, but I think it would be nice to pick that up soon after merging this one. I would propose to implement the switch in class RetentionLaw and to not make CalculateValue virtual. In that way, the switch is implemented only once, in a place where you would expect it to be.

Comment on lines +1264 to +1267
const Matrix coupling_matrix_up = GeoTransportEquationUtilities::CalculateCouplingMatrix(
rVariables.B, GetStressStatePolicy().GetVoigtVector(), rVariables.Np,
rVariables.BiotCoefficient, rVariables.BishopCoefficient, rVariables.IntegrationCoefficient);
GeoElementUtilities::AssembleUPBlockMatrix(rLeftHandSideMatrix, coupling_matrix);
GeoElementUtilities::AssembleUPBlockMatrix(rLeftHandSideMatrix, coupling_matrix_up);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am confused here: I had expected a minus sign in front of the coupling matrix, but I haven't found it yet... Perhaps @WPK4FEM can shed some light here, please?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

As we discussed, I have the same doubt and have been looking for the same. Trying to change it leads at least to some tests failing. We need to sort this out building up a really validated test suite.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Perhaps we can try to apply the Rule of zero here?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Tried with some help.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Perhaps we can try to apply the Rule of zero here?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Tried.

Comment on lines 19 to 21
SaturatedBelowPhreaticLevelLaw::SaturatedBelowPhreaticLevelLaw() : RetentionLaw() {}

SaturatedBelowPhreaticLevelLaw::SaturatedBelowPhreaticLevelLaw(const SaturatedBelowPhreaticLevelLaw& rOther)
: RetentionLaw(rOther)
{
}
SaturatedBelowPhreaticLevelLaw::SaturatedBelowPhreaticLevelLaw(const SaturatedBelowPhreaticLevelLaw& rOther) = default;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We should be able to remove these when we apply the Rule of zero.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

tried.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

tried.

Comment on lines 256 to 262
Begin SubModelPart DISPLACEMENT_Bottom // Group Bottom // Subtree DISPLACEMENT
Begin SubModelPartNodes
1
2
43
End SubModelPartNodes
End SubModelPart
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is an unused sub-model part.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Done

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Removed

End SubModelPartNodes
End SubModelPart

Begin SubModelPart DISPLACEMENT_Sides // Group Sides // Subtree DISPLACEMENT
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

And also this sub-model part is not being used.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Removed

End SubModelPartNodes
End SubModelPart

Begin SubModelPart DISPLACEMENT_Sides // Group Sides // Subtree DISPLACEMENT
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Also this sub-model part is not being used...

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Removed

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am missing displacement boundary conditions here. I had expected to see that the bottom edge would be fixed, and sliders along the sides of the column. Shouldn't we better add these?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Removed the bottom boundary, as it is not necessary

End SubModelPartNodes
End SubModelPart

Begin SubModelPart DISPLACEMENT_Bottom // Group Bottom // Subtree DISPLACEMENT
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This sub-model part should be referenced from the ProjectParameters.json files to fix the displacements at the bottom, but that is currently not being done.

UPDATE: Now I see that you fix all displacement DOFs, by prescribing zeroes for them. In that case, we don't need this sub-model part.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Removed

@markelov208 markelov208 self-requested a review December 20, 2024 10:19
@WPK4FEM WPK4FEM force-pushed the geo/partial_saturation2 branch from 6ace9cb to a8cd4ca Compare December 20, 2024 10:46
avdg81
avdg81 previously approved these changes Dec 20, 2024
Copy link
Contributor

@avdg81 avdg81 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks for all the hard work. I believe this PR is ready to be merged.

markelov208
markelov208 previously approved these changes Dec 20, 2024
@WPK4FEM WPK4FEM dismissed stale reviews from markelov208 and avdg81 via a741cd1 December 20, 2024 15:56
Copy link
Contributor

@avdg81 avdg81 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Let me approve it once more :-)

@WPK4FEM WPK4FEM enabled auto-merge (squash) December 20, 2024 16:30
@WPK4FEM WPK4FEM merged commit 5e77366 into master Dec 21, 2024
11 checks passed
@WPK4FEM WPK4FEM deleted the geo/partial_saturation2 branch December 21, 2024 15:13
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment