diff --git a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_FIC_element.cpp b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_FIC_element.cpp index d1006505f1b8..0db0f2ace1eb 100644 --- a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_FIC_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_FIC_element.cpp @@ -14,6 +14,7 @@ // Application includes #include "custom_elements/U_Pw_small_strain_FIC_element.hpp" #include "custom_utilities/math_utilities.h" +#include "custom_utilities/transport_equation_utilities.hpp" namespace Kratos { @@ -443,8 +444,6 @@ void UPwSmallStrainFICElement::CalculateAll(MatrixType& rLeftHa // create general parameters of retention law RetentionLaw::Parameters RetentionParameters(this->GetProperties(), CurrentProcessInfo); - const bool hasBiotCoefficient = Prop.Has(BIOT_COEFFICIENT); - const auto b_matrices = this->CalculateBMatrices(Variables.DN_DXContainer, Variables.NContainer); const auto integration_coefficients = this->CalculateIntegrationCoefficients(IntegrationPoints, Variables.detJContainer); @@ -455,10 +454,10 @@ void UPwSmallStrainFICElement::CalculateAll(MatrixType& rLeftHa this->CalculateAnyOfMaterialResponse(deformation_gradients, ConstitutiveParameters, Variables.NContainer, Variables.DN_DXContainer, strain_vectors, mStressVector, constitutive_matrices); + const auto biot_coefficients = + GeoTransportEquationUtilities::CalculateBiotCoefficients(constitutive_matrices, Prop); - // Loop over integration points for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) { - // Compute Np, GradNpT, B and StrainVector this->CalculateKinematics(Variables, GPoint); Variables.B = b_matrices[GPoint]; Variables.F = deformation_gradients[GPoint]; @@ -469,16 +468,14 @@ void UPwSmallStrainFICElement::CalculateAll(MatrixType& rLeftHa GeoElementUtilities::InterpolateVariableWithComponents( Variables.BodyAcceleration, Variables.NContainer, Variables.VolumeAcceleration, GPoint); - // Compute ShapeFunctionsSecondOrderGradients this->CalculateShapeFunctionsSecondOrderGradients(FICVariables, Variables); - this->CalculateRetentionResponse(Variables, RetentionParameters, GPoint); - // set shear modulus from stiffness matrix FICVariables.ShearModulus = CalculateShearModulus(Variables.ConstitutiveMatrix); - // calculate Bulk modulus from stiffness matrix - this->InitializeBiotCoefficients(Variables, hasBiotCoefficient); + Variables.BiotCoefficient = biot_coefficients[GPoint]; + Variables.BiotModulusInverse = GeoTransportEquationUtilities::CalculateBiotModulusInverse( + Variables.BiotCoefficient, Variables.DegreeOfSaturation, Variables.DerivativeOfSaturation, Prop); Variables.IntegrationCoefficient = integration_coefficients[GPoint]; diff --git a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_FIC_element.hpp b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_FIC_element.hpp index 9699e2fcc011..69822626f834 100644 --- a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_FIC_element.hpp +++ b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_FIC_element.hpp @@ -47,7 +47,6 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwSmallStrainFICElement using UPwBaseElement::mStateVariablesFinalized; using UPwBaseElement::mThisIntegrationMethod; - using UPwSmallStrainElement::CalculateBulkModulus; using UPwSmallStrainElement::VoigtSize; using ElementVariables = typename UPwSmallStrainElement::ElementVariables; diff --git a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp index edcd25595a76..7fb109bde8ff 100644 --- a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.cpp @@ -697,8 +697,6 @@ void UPwSmallStrainElement::CalculateOnIntegrationPoints(const for (unsigned int i = 0; i < StressTensorSize; ++i) VoigtVector[i] = 1.0; - const bool hasBiotCoefficient = rProp.Has(BIOT_COEFFICIENT); - const auto b_matrices = CalculateBMatrices(Variables.DN_DXContainer, Variables.NContainer); const auto deformation_gradients = CalculateDeformationGradients(); auto strain_vectors = StressStrainUtilities::CalculateStrains( @@ -707,23 +705,18 @@ void UPwSmallStrainElement::CalculateOnIntegrationPoints(const this->CalculateAnyOfMaterialResponse(deformation_gradients, ConstitutiveParameters, Variables.NContainer, Variables.DN_DXContainer, strain_vectors, mStressVector, constitutive_matrices); + const auto biot_coefficients = GeoTransportEquationUtilities::CalculateBiotCoefficients( + constitutive_matrices, this->GetProperties()); - // Loop over integration points for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) { this->CalculateKinematics(Variables, GPoint); - Variables.B = b_matrices[GPoint]; - Variables.F = deformation_gradients[GPoint]; - Variables.StrainVector = strain_vectors[GPoint]; - Variables.ConstitutiveMatrix = constitutive_matrices[GPoint]; - - Variables.BiotCoefficient = CalculateBiotCoefficient(Variables, hasBiotCoefficient); - - const auto fluid_pressure = - GeoTransportEquationUtilities::CalculateFluidPressure(Variables.Np, Variables.PressureVector); + const auto fluid_pressure = GeoTransportEquationUtilities::CalculateFluidPressure( + Variables.Np, Variables.PressureVector); RetentionParameters.SetFluidPressure(fluid_pressure); - const auto bishop_coefficient = mRetentionLawVector[GPoint]->CalculateBishopCoefficient(RetentionParameters); + const auto bishop_coefficient = + mRetentionLawVector[GPoint]->CalculateBishopCoefficient(RetentionParameters); - rOutput[GPoint] = mStressVector[GPoint] + PORE_PRESSURE_SIGN_FACTOR * Variables.BiotCoefficient * + rOutput[GPoint] = mStressVector[GPoint] + PORE_PRESSURE_SIGN_FACTOR * biot_coefficients[GPoint] * bishop_coefficient * fluid_pressure * VoigtVector; } } else if (rVariable == ENGINEERING_STRAIN_VECTOR) { @@ -973,8 +966,6 @@ void UPwSmallStrainElement::CalculateAll(MatrixType& rLe // Create general parameters of retention law RetentionLaw::Parameters RetentionParameters(this->GetProperties(), rCurrentProcessInfo); - const bool hasBiotCoefficient = rProp.Has(BIOT_COEFFICIENT); - const auto b_matrices = CalculateBMatrices(Variables.DN_DXContainer, Variables.NContainer); const auto integration_coefficients = this->CalculateIntegrationCoefficients(IntegrationPoints, Variables.detJContainer); @@ -985,6 +976,8 @@ void UPwSmallStrainElement::CalculateAll(MatrixType& rLe this->CalculateAnyOfMaterialResponse(deformation_gradients, ConstitutiveParameters, Variables.NContainer, Variables.DN_DXContainer, strain_vectors, mStressVector, constitutive_matrices); + const auto biot_coefficients = GeoTransportEquationUtilities::CalculateBiotCoefficients( + constitutive_matrices, this->GetProperties()); for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) { this->CalculateKinematics(Variables, GPoint); @@ -1000,7 +993,10 @@ void UPwSmallStrainElement::CalculateAll(MatrixType& rLe CalculateRetentionResponse(Variables, RetentionParameters, GPoint); - this->InitializeBiotCoefficients(Variables, hasBiotCoefficient); + Variables.BiotCoefficient = biot_coefficients[GPoint]; + Variables.BiotModulusInverse = GeoTransportEquationUtilities::CalculateBiotModulusInverse( + Variables.BiotCoefficient, Variables.DegreeOfSaturation, + Variables.DerivativeOfSaturation, this->GetProperties()); Variables.PermeabilityUpdateFactor = this->CalculatePermeabilityUpdateFactor(Variables.StrainVector); Variables.IntegrationCoefficient = integration_coefficients[GPoint]; @@ -1019,52 +1015,6 @@ void UPwSmallStrainElement::CalculateAll(MatrixType& rLe KRATOS_CATCH("") } -template -double UPwSmallStrainElement::CalculateBiotCoefficient(const ElementVariables& rVariables, - bool hasBiotCoefficient) const -{ - KRATOS_TRY - - const PropertiesType& rProp = this->GetProperties(); - - // Properties variables - if (hasBiotCoefficient) { - return rProp[BIOT_COEFFICIENT]; - } else { - // Calculate Bulk modulus from stiffness matrix - const double BulkModulus = CalculateBulkModulus(rVariables.ConstitutiveMatrix); - return 1.0 - BulkModulus / rProp[BULK_MODULUS_SOLID]; - } - - KRATOS_CATCH("") -} - -template -void UPwSmallStrainElement::InitializeBiotCoefficients(ElementVariables& rVariables, - bool hasBiotCoefficient) -{ - KRATOS_TRY - - const PropertiesType& rProp = this->GetProperties(); - - // Properties variables - rVariables.BiotCoefficient = CalculateBiotCoefficient(rVariables, hasBiotCoefficient); - - if (!rProp[IGNORE_UNDRAINED]) { - rVariables.BiotModulusInverse = - (rVariables.BiotCoefficient - rProp[POROSITY]) / rProp[BULK_MODULUS_SOLID] + - rProp[POROSITY] / rProp[BULK_MODULUS_FLUID]; - } else { - rVariables.BiotModulusInverse = - (rVariables.BiotCoefficient - rProp[POROSITY]) / rProp[BULK_MODULUS_SOLID] + rProp[POROSITY] / TINY; - } - - rVariables.BiotModulusInverse *= rVariables.DegreeOfSaturation; - rVariables.BiotModulusInverse -= rVariables.DerivativeOfSaturation * rProp[POROSITY]; - - KRATOS_CATCH("") -} - template std::vector> UPwSmallStrainElement::CalculateFluidFluxes( const std::vector& rPermeabilityUpdateFactors, const ProcessInfo& rCurrentProcessInfo) @@ -1126,17 +1076,6 @@ double UPwSmallStrainElement::CalculatePermeabilityUpdateFactor KRATOS_CATCH("") } -template -double UPwSmallStrainElement::CalculateBulkModulus(const Matrix& ConstitutiveMatrix) const -{ - KRATOS_TRY - - const int IndexG = ConstitutiveMatrix.size1() - 1; - return ConstitutiveMatrix(0, 0) - (4.0 / 3.0) * ConstitutiveMatrix(IndexG, IndexG); - - KRATOS_CATCH("") -} - template void UPwSmallStrainElement::InitializeElementVariables(ElementVariables& rVariables, const ProcessInfo& rCurrentProcessInfo) diff --git a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.hpp b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.hpp index e9397115b85d..7ae1a9a1d6c2 100644 --- a/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.hpp +++ b/applications/GeoMechanicsApplication/custom_elements/U_Pw_small_strain_element.hpp @@ -220,8 +220,6 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwSmallStrainElement : public UPwBa virtual void CalculateKinematics(ElementVariables& rVariables, unsigned int PointNumber); - void InitializeBiotCoefficients(ElementVariables& rVariables, bool hasBiotCoefficient = false); - double CalculatePermeabilityUpdateFactor(const Vector& rStrainVector) const; Matrix CalculateBMatrix(const Matrix& rDN_DX, const Vector& rN) const; @@ -266,9 +264,6 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwSmallStrainElement : public UPwBa array_1d& rPVector, const ElementVariables& rVariables) const; - double CalculateBulkModulus(const Matrix& ConstitutiveMatrix) const; - double CalculateBiotCoefficient(const ElementVariables& rVariables, bool hasBiotCoefficient) const; - virtual Vector CalculateGreenLagrangeStrain(const Matrix& rDeformationGradient) const; Matrix CalculateDeformationGradient(unsigned int GPoint) const; diff --git a/applications/GeoMechanicsApplication/custom_elements/U_Pw_updated_lagrangian_FIC_element.cpp b/applications/GeoMechanicsApplication/custom_elements/U_Pw_updated_lagrangian_FIC_element.cpp index 24335ab04e1c..c344305bb2de 100644 --- a/applications/GeoMechanicsApplication/custom_elements/U_Pw_updated_lagrangian_FIC_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/U_Pw_updated_lagrangian_FIC_element.cpp @@ -15,6 +15,7 @@ // Project includes #include "custom_elements/U_Pw_updated_lagrangian_FIC_element.hpp" #include "custom_utilities/math_utilities.h" +#include "custom_utilities/transport_equation_utilities.hpp" #include "utilities/math_utils.h" namespace Kratos @@ -71,8 +72,6 @@ void UPwUpdatedLagrangianFICElement::CalculateAll(MatrixType& r // create general parameters of retention law RetentionLaw::Parameters RetentionParameters(this->GetProperties(), rCurrentProcessInfo); - const bool hasBiotCoefficient = Prop.Has(BIOT_COEFFICIENT); - const auto b_matrices = this->CalculateBMatrices(Variables.DN_DXContainer, Variables.NContainer); const auto integration_coefficients = this->CalculateIntegrationCoefficients(IntegrationPoints, Variables.detJContainer); @@ -84,6 +83,8 @@ void UPwUpdatedLagrangianFICElement::CalculateAll(MatrixType& r this->CalculateAnyOfMaterialResponse(deformation_gradients, ConstitutiveParameters, Variables.NContainer, Variables.DN_DXContainer, strain_vectors, mStressVector, constitutive_matrices); + const auto biot_coefficients = + GeoTransportEquationUtilities::CalculateBiotCoefficients(constitutive_matrices, Prop); // Computing in all integrations points for (IndexType GPoint = 0; GPoint < IntegrationPoints.size(); ++GPoint) { @@ -111,8 +112,10 @@ void UPwUpdatedLagrangianFICElement::CalculateAll(MatrixType& r // set shear modulus from stiffness matrix FICVariables.ShearModulus = CalculateShearModulus(Variables.ConstitutiveMatrix); - // calculate Bulk modulus from stiffness matrix - this->InitializeBiotCoefficients(Variables, hasBiotCoefficient); + Variables.BiotCoefficient = biot_coefficients[GPoint]; + Variables.BiotModulusInverse = GeoTransportEquationUtilities::CalculateBiotModulusInverse( + Variables.BiotCoefficient, Variables.DegreeOfSaturation, + Variables.DerivativeOfSaturation, this->GetProperties()); Variables.IntegrationCoefficient = integration_coefficients[GPoint]; diff --git a/applications/GeoMechanicsApplication/custom_elements/U_Pw_updated_lagrangian_FIC_element.hpp b/applications/GeoMechanicsApplication/custom_elements/U_Pw_updated_lagrangian_FIC_element.hpp index ff6a2371ab8a..0053cefbe994 100644 --- a/applications/GeoMechanicsApplication/custom_elements/U_Pw_updated_lagrangian_FIC_element.hpp +++ b/applications/GeoMechanicsApplication/custom_elements/U_Pw_updated_lagrangian_FIC_element.hpp @@ -77,7 +77,6 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwUpdatedLagrangianFICElement using UPwBaseElement::CalculateDerivativesOnInitialConfiguration; using UPwBaseElement::mThisIntegrationMethod; using UPwSmallStrainFICElement::CalculateShearModulus; - using UPwSmallStrainElement::CalculateBulkModulus; using ElementVariables = typename UPwSmallStrainElement::ElementVariables; using FICElementVariables = typename UPwSmallStrainFICElement::FICElementVariables; diff --git a/applications/GeoMechanicsApplication/custom_elements/U_Pw_updated_lagrangian_element.cpp b/applications/GeoMechanicsApplication/custom_elements/U_Pw_updated_lagrangian_element.cpp index bc02c59201ec..c5dc9997c3f1 100644 --- a/applications/GeoMechanicsApplication/custom_elements/U_Pw_updated_lagrangian_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/U_Pw_updated_lagrangian_element.cpp @@ -68,7 +68,6 @@ void UPwUpdatedLagrangianElement::CalculateAll(MatrixType& rLef // create general parameters of retention law RetentionLaw::Parameters RetentionParameters(this->GetProperties(), rCurrentProcessInfo); - const bool hasBiotCoefficient = this->GetProperties().Has(BIOT_COEFFICIENT); const auto b_matrices = this->CalculateBMatrices(Variables.DN_DXContainer, Variables.NContainer); const auto integration_coefficients = this->CalculateIntegrationCoefficients(IntegrationPoints, Variables.detJContainer); @@ -80,6 +79,8 @@ void UPwUpdatedLagrangianElement::CalculateAll(MatrixType& rLef this->CalculateAnyOfMaterialResponse(deformation_gradients, ConstitutiveParameters, Variables.NContainer, Variables.DN_DXContainer, strain_vectors, mStressVector, constitutive_matrices); + const auto biot_coefficients = GeoTransportEquationUtilities::CalculateBiotCoefficients( + constitutive_matrices, this->GetProperties()); for (IndexType GPoint = 0; GPoint < IntegrationPoints.size(); ++GPoint) { this->CalculateKinematics(Variables, GPoint); @@ -98,8 +99,10 @@ void UPwUpdatedLagrangianElement::CalculateAll(MatrixType& rLef this->CalculateRetentionResponse(Variables, RetentionParameters, GPoint); - // calculate Bulk modulus from stiffness matrix - this->InitializeBiotCoefficients(Variables, hasBiotCoefficient); + Variables.BiotCoefficient = biot_coefficients[GPoint]; + Variables.BiotModulusInverse = GeoTransportEquationUtilities::CalculateBiotModulusInverse( + Variables.BiotCoefficient, Variables.DegreeOfSaturation, + Variables.DerivativeOfSaturation, this->GetProperties()); Variables.IntegrationCoefficient = integration_coefficients[GPoint]; diff --git a/applications/GeoMechanicsApplication/custom_elements/U_Pw_updated_lagrangian_element.hpp b/applications/GeoMechanicsApplication/custom_elements/U_Pw_updated_lagrangian_element.hpp index 7a1a96f1aad8..eeedc3d86ef3 100644 --- a/applications/GeoMechanicsApplication/custom_elements/U_Pw_updated_lagrangian_element.hpp +++ b/applications/GeoMechanicsApplication/custom_elements/U_Pw_updated_lagrangian_element.hpp @@ -79,7 +79,6 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) UPwUpdatedLagrangianElement using UPwBaseElement::mThisIntegrationMethod; using ElementVariables = typename UPwSmallStrainElement::ElementVariables; - using UPwSmallStrainElement::CalculateBulkModulus; /// Counted pointer of UPwUpdatedLagrangianElement KRATOS_CLASS_INTRUSIVE_POINTER_DEFINITION(UPwUpdatedLagrangianElement); diff --git a/applications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.cpp b/applications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.cpp index 4e1447819a95..a2850323edfc 100644 --- a/applications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.cpp @@ -1128,8 +1128,6 @@ void SmallStrainUPwDiffOrderElement::CalculateOnIntegrationPoints(const Variable ConstitutiveParameters.GetOptions().Set(ConstitutiveLaw::COMPUTE_CONSTITUTIVE_TENSOR); ConstitutiveParameters.GetOptions().Set(ConstitutiveLaw::USE_ELEMENT_PROVIDED_STRAIN); - const PropertiesType& rProp = this->GetProperties(); - const SizeType VoigtSize = mStressVector[0].size(); Vector VoigtVector = ZeroVector(VoigtSize); const SizeType StressTensorSize = (Dim == N_DIM_3D ? STRESS_TENSOR_SIZE_3D : STRESS_TENSOR_SIZE_2D); @@ -1139,8 +1137,6 @@ void SmallStrainUPwDiffOrderElement::CalculateOnIntegrationPoints(const Variable // create general parameters of retention law RetentionLaw::Parameters RetentionParameters(GetProperties(), rCurrentProcessInfo); - const bool hasBiotCoefficient = rProp.Has(BIOT_COEFFICIENT); - const auto b_matrices = CalculateBMatrices(Variables.DNu_DXContainer, Variables.NuContainer); const auto deformation_gradients = CalculateDeformationGradients(); auto strain_vectors = StressStrainUtilities::CalculateStrains( @@ -1150,22 +1146,18 @@ void SmallStrainUPwDiffOrderElement::CalculateOnIntegrationPoints(const Variable this->CalculateAnyOfMaterialResponse(deformation_gradients, ConstitutiveParameters, Variables.NuContainer, Variables.DNu_DXContainer, strain_vectors, mStressVector, constitutive_matrices); + const auto biot_coefficients = GeoTransportEquationUtilities::CalculateBiotCoefficients( + constitutive_matrices, this->GetProperties()); - // Loop over integration points for (unsigned int GPoint = 0; GPoint < mConstitutiveLawVector.size(); ++GPoint) { this->CalculateKinematics(Variables, GPoint); - Variables.B = b_matrices[GPoint]; - Variables.F = deformation_gradients[GPoint]; - Variables.StrainVector = strain_vectors[GPoint]; - Variables.ConstitutiveMatrix = constitutive_matrices[GPoint]; - Variables.BiotCoefficient = CalculateBiotCoefficient(Variables, hasBiotCoefficient); - - const auto fluid_pressure = - GeoTransportEquationUtilities::CalculateFluidPressure(Variables.Np, Variables.PressureVector); + const auto fluid_pressure = GeoTransportEquationUtilities::CalculateFluidPressure( + Variables.Np, Variables.PressureVector); RetentionParameters.SetFluidPressure(fluid_pressure); - const auto bishop_coefficient = mRetentionLawVector[GPoint]->CalculateBishopCoefficient(RetentionParameters); + const auto bishop_coefficient = + mRetentionLawVector[GPoint]->CalculateBishopCoefficient(RetentionParameters); - rOutput[GPoint] = mStressVector[GPoint] + PORE_PRESSURE_SIGN_FACTOR * Variables.BiotCoefficient * + rOutput[GPoint] = mStressVector[GPoint] + PORE_PRESSURE_SIGN_FACTOR * biot_coefficients[GPoint] * bishop_coefficient * fluid_pressure * VoigtVector; } } else { @@ -1275,7 +1267,6 @@ void SmallStrainUPwDiffOrderElement::CalculateAll(MatrixType& rLeftHandSi const GeometryType::IntegrationPointsArrayType& IntegrationPoints = rGeom.IntegrationPoints(this->GetIntegrationMethod()); - const bool hasBiotCoefficient = rProp.Has(BIOT_COEFFICIENT); const auto b_matrices = CalculateBMatrices(Variables.DNu_DXContainer, Variables.NuContainer); const auto deformation_gradients = CalculateDeformationGradients(); const auto integration_coefficients = @@ -1287,20 +1278,22 @@ void SmallStrainUPwDiffOrderElement::CalculateAll(MatrixType& rLeftHandSi this->CalculateAnyOfMaterialResponse(deformation_gradients, ConstitutiveParameters, Variables.NuContainer, Variables.DNu_DXContainer, strain_vectors, mStressVector, constitutive_matrices); + const auto biot_coefficients = GeoTransportEquationUtilities::CalculateBiotCoefficients( + constitutive_matrices, this->GetProperties()); for (unsigned int GPoint = 0; GPoint < IntegrationPoints.size(); ++GPoint) { - // compute element kinematics (Np, gradNpT, |J|, B, strains) this->CalculateKinematics(Variables, GPoint); - Variables.B = b_matrices[GPoint]; - - // Compute infinitesimal strain + Variables.B = b_matrices[GPoint]; Variables.F = deformation_gradients[GPoint]; Variables.StrainVector = strain_vectors[GPoint]; Variables.ConstitutiveMatrix = constitutive_matrices[GPoint]; CalculateRetentionResponse(Variables, RetentionParameters, GPoint); - this->InitializeBiotCoefficients(Variables, hasBiotCoefficient); + Variables.BiotCoefficient = biot_coefficients[GPoint]; + Variables.BiotModulusInverse = GeoTransportEquationUtilities::CalculateBiotModulusInverse( + Variables.BiotCoefficient, Variables.DegreeOfSaturation, + Variables.DerivativeOfSaturation, this->GetProperties()); Variables.PermeabilityUpdateFactor = this->CalculatePermeabilityUpdateFactor(Variables.StrainVector); Variables.IntegrationCoefficient = integration_coefficients[GPoint]; @@ -1359,16 +1352,6 @@ void SmallStrainUPwDiffOrderElement::CalculateMaterialStiffnessMatrix(MatrixType KRATOS_CATCH("") } -double SmallStrainUPwDiffOrderElement::CalculateBulkModulus(const Matrix& ConstitutiveMatrix) const -{ - KRATOS_TRY - - const SizeType IndexG = ConstitutiveMatrix.size1() - 1; - return ConstitutiveMatrix(0, 0) - (4.0 / 3.0) * ConstitutiveMatrix(IndexG, IndexG); - - KRATOS_CATCH("") -} - void SmallStrainUPwDiffOrderElement::InitializeElementVariables(ElementVariables& rVariables, const ProcessInfo& rCurrentProcessInfo) { @@ -1517,48 +1500,6 @@ void SmallStrainUPwDiffOrderElement::InitializeNodalVariables(ElementVariables& KRATOS_CATCH("") } -double SmallStrainUPwDiffOrderElement::CalculateBiotCoefficient(const ElementVariables& rVariables, - const bool& hasBiotCoefficient) const -{ - KRATOS_TRY - - const PropertiesType& rProp = this->GetProperties(); - - // Properties variables - if (hasBiotCoefficient) { - return rProp[BIOT_COEFFICIENT]; - } else { - // calculate Bulk modulus from stiffness matrix - return 1.0 - CalculateBulkModulus(rVariables.ConstitutiveMatrix) / rProp[BULK_MODULUS_SOLID]; - } - - KRATOS_CATCH("") -} - -void SmallStrainUPwDiffOrderElement::InitializeBiotCoefficients(ElementVariables& rVariables, - const bool& hasBiotCoefficient) -{ - KRATOS_TRY - - const PropertiesType& rProp = this->GetProperties(); - - rVariables.BiotCoefficient = CalculateBiotCoefficient(rVariables, hasBiotCoefficient); - - if (rProp[IGNORE_UNDRAINED]) { - rVariables.BiotModulusInverse = - (rVariables.BiotCoefficient - rProp[POROSITY]) / rProp[BULK_MODULUS_SOLID] + rProp[POROSITY] / TINY; - } else { - rVariables.BiotModulusInverse = - (rVariables.BiotCoefficient - rProp[POROSITY]) / rProp[BULK_MODULUS_SOLID] + - rProp[POROSITY] / rProp[BULK_MODULUS_FLUID]; - } - - rVariables.BiotModulusInverse *= rVariables.DegreeOfSaturation; - rVariables.BiotModulusInverse -= rVariables.DerivativeOfSaturation * rProp[POROSITY]; - - KRATOS_CATCH("") -} - double SmallStrainUPwDiffOrderElement::CalculatePermeabilityUpdateFactor(const Vector& rStrainVector) const { KRATOS_TRY @@ -1836,9 +1777,8 @@ void SmallStrainUPwDiffOrderElement::CalculateAndAddMixBodyForce(VectorType& rRi for (SizeType i = 0; i < NumUNodes; ++i) { Index = i * Dim; for (SizeType idim = 0; idim < Dim; ++idim) { - rRightHandSideVector[Index + idim] += - rVariables.Nu[i] * soil_density * BodyAcceleration[idim] * - rVariables.IntegrationCoefficientInitialConfiguration; + rRightHandSideVector[Index + idim] += rVariables.Nu[i] * soil_density * BodyAcceleration[idim] * + rVariables.IntegrationCoefficientInitialConfiguration; } } diff --git a/applications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.hpp b/applications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.hpp index c5f22acbf3ae..8b178ba36a15 100644 --- a/applications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.hpp +++ b/applications/GeoMechanicsApplication/custom_elements/small_strain_U_Pw_diff_order_element.hpp @@ -238,8 +238,6 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) SmallStrainUPwDiffOrderElement : pub void InitializeProperties(ElementVariables& rVariables); - void InitializeBiotCoefficients(ElementVariables& rVariables, const bool& hasBiotCoefficient = false); - double CalculatePermeabilityUpdateFactor(const Vector& rStrainVector) const; virtual void CalculateKinematics(ElementVariables& rVariables, unsigned int GPoint); @@ -279,9 +277,6 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) SmallStrainUPwDiffOrderElement : pub void CalculateAndAddFluidBodyFlow(VectorType& rRightHandSideVector, ElementVariables& rVariables); - double CalculateBulkModulus(const Matrix& ConstitutiveMatrix) const; - double CalculateBiotCoefficient(const ElementVariables& rVariables, const bool& hasBiotCoefficient) const; - Matrix CalculateBMatrix(const Matrix& rDN_DX, const Vector& rN) const; std::vector CalculateBMatrices(const GeometryType::ShapeFunctionsGradientsType& rDN_DXContainer, const Matrix& rNContainer) const; diff --git a/applications/GeoMechanicsApplication/custom_elements/transient_Pw_element.cpp b/applications/GeoMechanicsApplication/custom_elements/transient_Pw_element.cpp index 250af15e7405..971d1bca4c21 100644 --- a/applications/GeoMechanicsApplication/custom_elements/transient_Pw_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/transient_Pw_element.cpp @@ -453,10 +453,9 @@ void TransientPwElement::CalculateAll(MatrixType& rLeftH // create general parameters of retention law RetentionLaw::Parameters RetentionParameters(this->GetProperties(), rCurrentProcessInfo); - const bool hasBiotCoefficient = Prop.Has(BIOT_COEFFICIENT); - const auto integration_coefficients = this->CalculateIntegrationCoefficients(IntegrationPoints, Variables.detJContainer); + std::vector biot_coefficients(NumGPoints, Prop[BIOT_COEFFICIENT]); // Loop over integration points for (unsigned int GPoint = 0; GPoint < NumGPoints; ++GPoint) { @@ -470,7 +469,10 @@ void TransientPwElement::CalculateAll(MatrixType& rLeftH this->CalculateRetentionResponse(Variables, RetentionParameters, GPoint); - this->InitializeBiotCoefficients(Variables, hasBiotCoefficient); + Variables.BiotCoefficient = biot_coefficients[GPoint]; + Variables.BiotModulusInverse = GeoTransportEquationUtilities::CalculateBiotModulusInverse( + Variables.BiotCoefficient, Variables.DegreeOfSaturation, + Variables.DerivativeOfSaturation, this->GetProperties()); Variables.IntegrationCoefficient = integration_coefficients[GPoint]; diff --git a/applications/GeoMechanicsApplication/custom_elements/updated_lagrangian_U_Pw_diff_order_element.cpp b/applications/GeoMechanicsApplication/custom_elements/updated_lagrangian_U_Pw_diff_order_element.cpp index d2de2fa02705..02e4dac3fbbb 100644 --- a/applications/GeoMechanicsApplication/custom_elements/updated_lagrangian_U_Pw_diff_order_element.cpp +++ b/applications/GeoMechanicsApplication/custom_elements/updated_lagrangian_U_Pw_diff_order_element.cpp @@ -15,6 +15,7 @@ // Project includes #include "custom_elements/updated_lagrangian_U_Pw_diff_order_element.hpp" #include "custom_utilities/math_utilities.h" +#include "custom_utilities/transport_equation_utilities.hpp" #include "utilities/math_utils.h" namespace Kratos @@ -70,7 +71,6 @@ void UpdatedLagrangianUPwDiffOrderElement::CalculateAll(MatrixType& rLeft const GeometryType::IntegrationPointsArrayType& IntegrationPoints = rGeom.IntegrationPoints(this->GetIntegrationMethod()); - const bool hasBiotCoefficient = rProp.Has(BIOT_COEFFICIENT); const auto b_matrices = this->CalculateBMatrices(Variables.DNu_DXContainer, Variables.NuContainer); const auto deformation_gradients = this->CalculateDeformationGradients(); const auto integration_coefficients = @@ -82,22 +82,23 @@ void UpdatedLagrangianUPwDiffOrderElement::CalculateAll(MatrixType& rLeft this->CalculateAnyOfMaterialResponse(deformation_gradients, ConstitutiveParameters, Variables.NuContainer, Variables.DNu_DXContainer, strain_vectors, mStressVector, constitutive_matrices); + const auto biot_coefficients = GeoTransportEquationUtilities::CalculateBiotCoefficients( + constitutive_matrices, this->GetProperties()); - // Computing in all integrations points for (IndexType GPoint = 0; GPoint < IntegrationPoints.size(); ++GPoint) { - // Compute element kinematics B, F, DNu_DX ... this->CalculateKinematics(Variables, GPoint); Variables.B = b_matrices[GPoint]; - // Compute strain Variables.F = deformation_gradients[GPoint]; Variables.StrainVector = strain_vectors[GPoint]; Variables.ConstitutiveMatrix = constitutive_matrices[GPoint]; CalculateRetentionResponse(Variables, RetentionParameters, GPoint); - // calculate Bulk modulus from stiffness matrix - this->InitializeBiotCoefficients(Variables, hasBiotCoefficient); + Variables.BiotCoefficient = biot_coefficients[GPoint]; + Variables.BiotModulusInverse = GeoTransportEquationUtilities::CalculateBiotModulusInverse( + Variables.BiotCoefficient, Variables.DegreeOfSaturation, + Variables.DerivativeOfSaturation, this->GetProperties()); Variables.IntegrationCoefficient = integration_coefficients[GPoint]; diff --git a/applications/GeoMechanicsApplication/custom_utilities/transport_equation_utilities.hpp b/applications/GeoMechanicsApplication/custom_utilities/transport_equation_utilities.hpp index 69f98888f409..5b55ef8a9f1a 100644 --- a/applications/GeoMechanicsApplication/custom_utilities/transport_equation_utilities.hpp +++ b/applications/GeoMechanicsApplication/custom_utilities/transport_equation_utilities.hpp @@ -118,5 +118,49 @@ class GeoTransportEquationUtilities return inner_prod(rN, rPressureVector); } + [[nodiscard]] static double CalculateBiotModulusInverse(double BiotCoefficient, + double DegreeOfSaturation, + double DerivativeOfSaturation, + const Properties& rProperties) + { + const auto bulk_modulus_fluid = rProperties[IGNORE_UNDRAINED] ? TINY : rProperties[BULK_MODULUS_FLUID]; + double result = (BiotCoefficient - rProperties[POROSITY]) / rProperties[BULK_MODULUS_SOLID] + + rProperties[POROSITY] / bulk_modulus_fluid; + + result *= DegreeOfSaturation; + result -= DerivativeOfSaturation * rProperties[POROSITY]; + + return result; + } + + [[nodiscard]] static double CalculateBulkModulus(const Matrix& rConstitutiveMatrix) + { + KRATOS_ERROR_IF(rConstitutiveMatrix.size1() == 0) + << "Constitutive matrix is empty, aborting bulk modulus calculation.\n"; + const SizeType index_G = rConstitutiveMatrix.size1() - 1; + return rConstitutiveMatrix(0, 0) - (4.0 / 3.0) * rConstitutiveMatrix(index_G, index_G); + } + + [[nodiscard]] static std::vector CalculateBiotCoefficients(const std::vector& rConstitutiveMatrices, + const Properties& rProperties) + { + std::vector result; + std::transform(rConstitutiveMatrices.begin(), rConstitutiveMatrices.end(), + std::back_inserter(result), [&rProperties](const Matrix& rConstitutiveMatrix) { + return CalculateBiotCoefficient(rConstitutiveMatrix, rProperties); + }); + + return result; + } + +private: + [[nodiscard]] static double CalculateBiotCoefficient(const Matrix& rConstitutiveMatrix, + const Properties& rProperties) + { + return rProperties.Has(BIOT_COEFFICIENT) + ? rProperties[BIOT_COEFFICIENT] + : 1.0 - CalculateBulkModulus(rConstitutiveMatrix) / rProperties[BULK_MODULUS_SOLID]; + } + }; /* Class GeoTransportEquationUtilities*/ } /* namespace Kratos.*/ diff --git a/applications/GeoMechanicsApplication/tests/cpp_tests/test_transport_equation_utilities.cpp b/applications/GeoMechanicsApplication/tests/cpp_tests/test_transport_equation_utilities.cpp new file mode 100644 index 000000000000..0944c29bde02 --- /dev/null +++ b/applications/GeoMechanicsApplication/tests/cpp_tests/test_transport_equation_utilities.cpp @@ -0,0 +1,137 @@ +// KRATOS___ +// // ) ) +// // ___ ___ +// // ____ //___) ) // ) ) +// // / / // // / / +// ((____/ / ((____ ((___/ / MECHANICS +// +// License: geo_mechanics_application/license.txt +// +// Main authors: Richard Faasse +// + +#include "custom_utilities/transport_equation_utilities.hpp" +#include "testing/testing.h" +#include + +namespace Kratos::Testing +{ + +KRATOS_TEST_CASE_IN_SUITE(CalculateBiotModulusInverse_GivesExpectedResult, KratosGeoMechanicsFastSuite) +{ + Properties properties; + properties[IGNORE_UNDRAINED] = false; + properties[POROSITY] = 0.5; + properties[BULK_MODULUS_SOLID] = 1.0e9; + properties[BULK_MODULUS_FLUID] = 2.0e6; + const double biot_coefficient = 1.0; + const double degree_of_saturation = 0.3; + const double derivative_of_saturation = 0.2; + + const double expected_value = -0.09999992485; + KRATOS_EXPECT_DOUBLE_EQ(GeoTransportEquationUtilities::CalculateBiotModulusInverse( + biot_coefficient, degree_of_saturation, derivative_of_saturation, properties), + expected_value); +} + +KRATOS_TEST_CASE_IN_SUITE(CalculateBiotModulusInverse_ReturnsLargeNumber_WhenIgnoreUndrained, KratosGeoMechanicsFastSuite) +{ + Properties properties; + properties[IGNORE_UNDRAINED] = true; + properties[POROSITY] = 0.5; + properties[BULK_MODULUS_SOLID] = 1.0e9; + properties[BULK_MODULUS_FLUID] = 2.0e6; + const double biot_coefficient = 1.0; + const double degree_of_saturation = 0.3; + const double derivative_of_saturation = 0.2; + + const auto large_number = 1e10; + KRATOS_EXPECT_GT(GeoTransportEquationUtilities::CalculateBiotModulusInverse( + biot_coefficient, degree_of_saturation, derivative_of_saturation, properties), + large_number); +} + +KRATOS_TEST_CASE_IN_SUITE(CalculateBiotModulusInverse_DoesNotThrow_ForEmptyProperties, KratosGeoMechanicsFastSuite) +{ + const Properties properties; + const double biot_coefficient = 1.0; + const double degree_of_saturation = 0.3; + const double derivative_of_saturation = 0.2; + + KRATOS_EXPECT_TRUE(std::isnan(GeoTransportEquationUtilities::CalculateBiotModulusInverse( + biot_coefficient, degree_of_saturation, derivative_of_saturation, properties))) +} + +KRATOS_TEST_CASE_IN_SUITE(CalculateBulkModulus_ReturnsZeroForZeroConstitutiveMatrix, KratosGeoMechanicsFastSuite) +{ + ZeroMatrix constitutive_matrix(3, 3); + KRATOS_EXPECT_DOUBLE_EQ(GeoTransportEquationUtilities::CalculateBulkModulus(constitutive_matrix), 0.0); +} + +KRATOS_TEST_CASE_IN_SUITE(CalculateBulkModulus_Throws_WhenConstitutiveMatrixIsEmpty, KratosGeoMechanicsFastSuite) +{ + const Matrix constitutive_matrix; + KRATOS_EXPECT_EXCEPTION_IS_THROWN( + [[maybe_unused]] const auto bulk_modulus = + GeoTransportEquationUtilities::CalculateBulkModulus(constitutive_matrix), + "Constitutive matrix is empty, aborting bulk modulus calculation.") +} + +KRATOS_TEST_CASE_IN_SUITE(CalculateBulkModulus_GivesExpectedResult_ForFilledConstitutiveMatrix, KratosGeoMechanicsFastSuite) +{ + Matrix constitutive_matrix(3, 3); + constitutive_matrix <<= 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0; + + KRATOS_EXPECT_DOUBLE_EQ(GeoTransportEquationUtilities::CalculateBulkModulus(constitutive_matrix), -11.0); +} + +KRATOS_TEST_CASE_IN_SUITE(CalculateBiotCoefficients_GivesExpectedResults_WhenPropertyHasBiotCoefficient, + KratosGeoMechanicsFastSuite) +{ + const std::vector constitutive_matrices(2, ZeroMatrix(3, 3)); + + const double biot_coefficient = 1.2; + Properties properties; + properties[BIOT_COEFFICIENT] = biot_coefficient; + + const auto actual_values = + GeoTransportEquationUtilities::CalculateBiotCoefficients(constitutive_matrices, properties); + + const std::vector expected_values(2, biot_coefficient); + KRATOS_CHECK_VECTOR_NEAR(expected_values, actual_values, 1e-12) +} + +KRATOS_TEST_CASE_IN_SUITE(CalculateBiotCoefficients_GivesExpectedResults, KratosGeoMechanicsFastSuite) +{ + std::vector constitutive_matrices; + constitutive_matrices.emplace_back(ScalarMatrix(3, 3, 1.0)); + constitutive_matrices.emplace_back(ScalarMatrix(3, 3, 3.0)); + + Properties properties; + properties[BULK_MODULUS_SOLID] = 1.0; + + const auto actual_values = + GeoTransportEquationUtilities::CalculateBiotCoefficients(constitutive_matrices, properties); + + const std::vector expected_values = {4.0 / 3.0, 2.0}; + KRATOS_CHECK_VECTOR_NEAR(expected_values, actual_values, 1e-12) +} + +KRATOS_TEST_CASE_IN_SUITE(CalculateBiotCoefficients_GivesInfResults_ForZeroBulkModulus, KratosGeoMechanicsFastSuite) +{ + std::vector constitutive_matrices; + constitutive_matrices.emplace_back(ScalarMatrix(3, 3, 1.0)); + constitutive_matrices.emplace_back(ScalarMatrix(3, 3, 3.0)); + + Properties properties; + properties[BULK_MODULUS_SOLID] = 0.0; + + const auto actual_values = + GeoTransportEquationUtilities::CalculateBiotCoefficients(constitutive_matrices, Properties()); + + // Due to a division by zero, we end up with inf + KRATOS_EXPECT_TRUE(std::all_of(actual_values.begin(), actual_values.end(), + [](const double value) { return std::isinf(value); })) +} + +} // namespace Kratos::Testing \ No newline at end of file