Skip to content

Commit

Permalink
[GeoMechanicsApplication] Partially saturated flow does not work as e…
Browse files Browse the repository at this point in the history 
…xpected (#12914)

Take saturation into account in fluid body flow, avoid zero divisions, streamline some code
Co-authored-by: Wijtze Pieter Kikstra <wijtzepieter.kikstra@deltares.nl>
Co-authored-by: Anne van de Graaf <anne.vandegraaf@deltares.nl>
  • Loading branch information
@mnabideltares
mnabideltares authored Dec 21, 2024
1 parent 9735f25 commit 5e77366
Show file tree
Hide file tree
Showing 45 changed files with 4,234 additions and 980 deletions.
11 changes: 4 additions & 7 deletions applications/GeoMechanicsApplication/custom_elements/U_Pw_base_element.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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(
r_properties, r_geometry, row(r_geometry.ShapeFunctionsValues(mThisIntegrationMethod), i));
for (auto& r_retention_law : mRetentionLawVector) {
r_retention_law = RetentionLawFactory::Clone(r_properties);
}

if (mStressVector.size() != number_of_integration_points) {
Expand All @@ -151,9 +149,8 @@ void UPwBaseElement::Initialize(const ProcessInfo& rCurrentProcessInfo)

mStateVariablesFinalized.resize(number_of_integration_points);
for (unsigned int i = 0; i < mConstitutiveLawVector.size(); ++i) {
int nStateVariables = 0;
nStateVariables = mConstitutiveLawVector[i]->GetValue(NUMBER_OF_UMAT_STATE_VARIABLES, nStateVariables);
if (nStateVariables > 0) {
if (auto nStateVariables = 0;
mConstitutiveLawVector[i]->GetValue(NUMBER_OF_UMAT_STATE_VARIABLES, nStateVariables) > 0) {
mConstitutiveLawVector[i]->SetValue(STATE_VARIABLES, mStateVariablesFinalized[i], rCurrentProcessInfo);
}
}
Expand Down
75 changes: 30 additions & 45 deletions applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -258,8 +258,6 @@ void UPwSmallStrainElement<TDim, TNumNodes>::InitializeSolutionStep(const Proces
ElementVariables Variables;
this->InitializeElementVariables(Variables, rCurrentProcessInfo);

RetentionLaw::Parameters RetentionParameters(this->GetProperties());

const auto b_matrices = CalculateBMatrices(Variables.DN_DXContainer, Variables.NContainer);
const auto deformation_gradients = CalculateDeformationGradients();
const auto determinants_of_deformation_gradients =
Expand All @@ -283,8 +281,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);
}

// Reset hydraulic discharge
Expand All @@ -299,9 +295,8 @@ void UPwSmallStrainElement<TDim, TNumNodes>::ResetHydraulicDischarge()
KRATOS_TRY

// Reset hydraulic discharge
GeometryType& rGeom = this->GetGeometry();
for (unsigned int i = 0; i < TNumNodes; ++i) {
ThreadSafeNodeWrite(rGeom[i], HYDRAULIC_DISCHARGE, 0.0);
for (auto& r_node : this->GetGeometry()) {
ThreadSafeNodeWrite(r_node, HYDRAULIC_DISCHARGE, 0.0);
}

KRATOS_CATCH("")
Expand Down Expand Up @@ -357,8 +352,8 @@ void UPwSmallStrainElement<TDim, TNumNodes>::InitializeNonLinearIteration(const
{
KRATOS_TRY

const GeometryType& rGeom = this->GetGeometry();
ConstitutiveLaw::Parameters ConstitutiveParameters(rGeom, this->GetProperties(), rCurrentProcessInfo);
ConstitutiveLaw::Parameters ConstitutiveParameters(this->GetGeometry(), this->GetProperties(),
rCurrentProcessInfo);
ConstitutiveParameters.Set(ConstitutiveLaw::COMPUTE_STRESS);
ConstitutiveParameters.Set(ConstitutiveLaw::USE_ELEMENT_PROVIDED_STRAIN);
ConstitutiveParameters.Set(ConstitutiveLaw::INITIALIZE_MATERIAL_RESPONSE); // Note: this is for nonlocal damage
Expand Down Expand Up @@ -403,8 +398,6 @@ void UPwSmallStrainElement<TDim, TNumNodes>::FinalizeSolutionStep(const ProcessI
ElementVariables Variables;
this->InitializeElementVariables(Variables, rCurrentProcessInfo);

RetentionLaw::Parameters RetentionParameters(this->GetProperties());

const auto b_matrices = CalculateBMatrices(Variables.DN_DXContainer, Variables.NContainer);
const auto deformation_gradients = CalculateDeformationGradients();
const auto determinants_of_deformation_gradients =
Expand Down Expand Up @@ -434,8 +427,6 @@ void UPwSmallStrainElement<TDim, TNumNodes>::FinalizeSolutionStep(const ProcessI
mStateVariablesFinalized[GPoint] =
mConstitutiveLawVector[GPoint]->GetValue(STATE_VARIABLES, mStateVariablesFinalized[GPoint]);

mRetentionLawVector[GPoint]->FinalizeSolutionStep(RetentionParameters);

if (rCurrentProcessInfo[NODAL_SMOOTHING])
this->SaveGPStress(StressContainer, mStressVector[GPoint], GPoint);
}
Expand Down Expand Up @@ -622,19 +613,8 @@ void UPwSmallStrainElement<TDim, TNumNodes>::CalculateOnIntegrationPoints(const

RetentionParameters.SetFluidPressure(GeoTransportEquationUtilities::CalculateFluidPressure(
Variables.Np, Variables.PressureVector));

if (rVariable == DEGREE_OF_SATURATION)
rOutput[GPoint] = mRetentionLawVector[GPoint]->CalculateSaturation(RetentionParameters);
else if (rVariable == EFFECTIVE_SATURATION)
rOutput[GPoint] = mRetentionLawVector[GPoint]->CalculateEffectiveSaturation(RetentionParameters);
else if (rVariable == BISHOP_COEFFICIENT)
rOutput[GPoint] = mRetentionLawVector[GPoint]->CalculateBishopCoefficient(RetentionParameters);
else if (rVariable == DERIVATIVE_OF_SATURATION)
rOutput[GPoint] =
mRetentionLawVector[GPoint]->CalculateDerivativeOfSaturation(RetentionParameters);
else if (rVariable == RELATIVE_PERMEABILITY)
rOutput[GPoint] =
mRetentionLawVector[GPoint]->CalculateRelativePermeability(RetentionParameters);
rOutput[GPoint] = mRetentionLawVector[GPoint]->CalculateValue(
RetentionParameters, rVariable, rOutput[GPoint]);
}
} else if (rVariable == HYDRAULIC_HEAD) {
const PropertiesType& rProp = this->GetProperties();
Expand Down Expand Up @@ -1246,16 +1226,17 @@ void UPwSmallStrainElement<TDim, TNumNodes>::CalculateAndAddLHS(MatrixType& rLef
KRATOS_TRY

this->CalculateAndAddStiffnessMatrix(rLeftHandSideMatrix, rVariables);
this->CalculateAndAddCompressibilityMatrix(rLeftHandSideMatrix, rVariables);

if (!rVariables.IgnoreUndrained) {
this->CalculateAndAddCouplingMatrix(rLeftHandSideMatrix, rVariables);
this->CalculateAndAddCouplingMatrix(rLeftHandSideMatrix, rVariables);

if (!rVariables.IgnoreUndrained) {
const auto permeability_matrix =
GeoTransportEquationUtilities::CalculatePermeabilityMatrix<TDim, TNumNodes>(
rVariables.GradNpT, rVariables.DynamicViscosityInverse, rVariables.PermeabilityMatrix,
rVariables.RelativePermeability, rVariables.IntegrationCoefficient);
GeoElementUtilities::AssemblePPBlockMatrix(rLeftHandSideMatrix, permeability_matrix);

this->CalculateAndAddCompressibilityMatrix(rLeftHandSideMatrix, rVariables);
}

KRATOS_CATCH("")
Expand All @@ -1280,17 +1261,18 @@ void UPwSmallStrainElement<TDim, TNumNodes>::CalculateAndAddCouplingMatrix(Matri
{
KRATOS_TRY

const auto coupling_matrix = GeoTransportEquationUtilities::CalculateCouplingMatrix(
rVariables.B, this->GetStressStatePolicy().GetVoigtVector(), rVariables.Np,
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);

if (!rVariables.IgnoreUndrained) {
const double SaturationCoefficient = rVariables.DegreeOfSaturation / rVariables.BishopCoefficient;
const BoundedMatrix<double, TNumNodes, TNumNodes * TDim> transposed_coupling_matrix =
PORE_PRESSURE_SIGN_FACTOR * SaturationCoefficient * rVariables.VelocityCoefficient *
trans(coupling_matrix);
GeoElementUtilities::AssemblePUBlockMatrix(rLeftHandSideMatrix, transposed_coupling_matrix);
const Matrix coupling_matrix_pu = GeoTransportEquationUtilities::CalculateCouplingMatrix(
rVariables.B, GetStressStatePolicy().GetVoigtVector(), rVariables.Np,
rVariables.BiotCoefficient, rVariables.DegreeOfSaturation, rVariables.IntegrationCoefficient);
GeoElementUtilities::AssemblePUBlockMatrix(
rLeftHandSideMatrix,
PORE_PRESSURE_SIGN_FACTOR * rVariables.VelocityCoefficient * trans(coupling_matrix_pu));
}

KRATOS_CATCH("")
Expand Down Expand Up @@ -1389,15 +1371,18 @@ void UPwSmallStrainElement<TDim, TNumNodes>::CalculateAndAddCouplingTerms(Vector
rVariables.B, this->GetStressStatePolicy().GetVoigtVector(), rVariables.Np,
rVariables.BiotCoefficient, rVariables.BishopCoefficient, rVariables.IntegrationCoefficient);

const array_1d<double, TNumNodes * TDim> coupling_terms = prod(coupling_matrix, rVariables.PressureVector);
GeoElementUtilities::AssembleUBlockVector(rRightHandSideVector, coupling_terms);
const array_1d<double, TNumNodes * TDim> coupling_force = prod(coupling_matrix, rVariables.PressureVector);
GeoElementUtilities::AssembleUBlockVector(rRightHandSideVector, coupling_force);

if (!rVariables.IgnoreUndrained) {
const double SaturationCoefficient = rVariables.DegreeOfSaturation / rVariables.BishopCoefficient;
const array_1d<double, TNumNodes> transposed_coupling_matrix =
PORE_PRESSURE_SIGN_FACTOR * SaturationCoefficient *
prod(trans(coupling_matrix), rVariables.VelocityVector);
GeoElementUtilities::AssemblePBlockVector(rRightHandSideVector, transposed_coupling_matrix);
const Matrix p_coupling_matrix =
(-1.0) * GeoTransportEquationUtilities::CalculateCouplingMatrix(
rVariables.B, GetStressStatePolicy().GetVoigtVector(), rVariables.Np,
rVariables.BiotCoefficient, rVariables.DegreeOfSaturation,
rVariables.IntegrationCoefficient);
const array_1d<double, TNumNodes> coupling_flow =
PORE_PRESSURE_SIGN_FACTOR * prod(trans(p_coupling_matrix), rVariables.VelocityVector);
GeoElementUtilities::AssemblePBlockVector(rRightHandSideVector, coupling_flow);
}

KRATOS_CATCH("")
Expand Down Expand Up @@ -1493,7 +1478,7 @@ array_1d<double, TNumNodes> UPwSmallStrainElement<TDim, TNumNodes>::CalculateFlu
const BoundedMatrix<double, TNumNodes, TDim> temp_matrix =
prod(rVariables.GradNpT, rVariables.PermeabilityMatrix) * rVariables.IntegrationCoefficient;

return rVariables.DynamicViscosityInverse * this->GetProperties()[DENSITY_WATER] *
return rVariables.DynamicViscosityInverse * this->GetProperties()[DENSITY_WATER] * rVariables.BishopCoefficient *
rVariables.RelativePermeability * prod(temp_matrix, rVariables.BodyAcceleration);

KRATOS_CATCH("")
Expand Down
10 changes: 0 additions & 10 deletions applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_interface_element.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -192,8 +192,6 @@ void UPwSmallStrainInterfaceElement<TDim, TNumNodes>::InitializeSolutionStep(con
unsigned int NumGPoints = mConstitutiveLawVector.size();
std::vector<double> JointWidthContainer(NumGPoints);

RetentionLaw::Parameters RetentionParameters(this->GetProperties());

// Loop over integration points
for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) {
InterfaceElementUtilities::CalculateNuMatrix(Nu, NContainer, GPoint);
Expand All @@ -206,9 +204,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);
}

KRATOS_CATCH("")
Expand Down Expand Up @@ -253,8 +248,6 @@ void UPwSmallStrainInterfaceElement<TDim, TNumNodes>::FinalizeSolutionStep(const
unsigned int NumGPoints = mConstitutiveLawVector.size();
std::vector<double> JointWidthContainer(NumGPoints);

RetentionLaw::Parameters RetentionParameters(this->GetProperties());

// Loop over integration points
for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) {
InterfaceElementUtilities::CalculateNuMatrix(Nu, NContainer, GPoint);
Expand All @@ -279,9 +272,6 @@ void UPwSmallStrainInterfaceElement<TDim, TNumNodes>::FinalizeSolutionStep(const

mStateVariablesFinalized[GPoint] =
mConstitutiveLawVector[GPoint]->GetValue(STATE_VARIABLES, mStateVariablesFinalized[GPoint]);

// retention law
mRetentionLawVector[GPoint]->FinalizeSolutionStep(RetentionParameters);
}

if (rCurrentProcessInfo[NODAL_SMOOTHING]) this->ExtrapolateGPValues(JointWidthContainer);
Expand Down
Loading

0 comments on commit 5e77366

Please sign in to comment.