From beb1ec175ab41eb88d173f00d5cd4e966c3ac37c Mon Sep 17 00:00:00 2001 From: mnabideltares Date: Tue, 14 May 2024 13:28:37 +0200 Subject: [PATCH 1/2] Added pw_filter_law --- .../custom_constitutive/pw_filter_law.cpp | 45 ++++++++++++++ .../custom_constitutive/pw_filter_law.h | 62 +++++++++++++++++++ 2 files changed, 107 insertions(+) create mode 100644 applications/GeoMechanicsApplication/custom_constitutive/pw_filter_law.cpp create mode 100644 applications/GeoMechanicsApplication/custom_constitutive/pw_filter_law.h diff --git a/applications/GeoMechanicsApplication/custom_constitutive/pw_filter_law.cpp b/applications/GeoMechanicsApplication/custom_constitutive/pw_filter_law.cpp new file mode 100644 index 000000000000..0f6a66755d02 --- /dev/null +++ b/applications/GeoMechanicsApplication/custom_constitutive/pw_filter_law.cpp @@ -0,0 +1,45 @@ +// KRATOS___ +// // ) ) +// // ___ ___ +// // ____ //___) ) // ) ) +// // / / // // / / +// ((____/ / ((____ ((___/ / MECHANICS +// +// License: geo_mechanics_application/license.txt +// +// Main authors: Mohamed Nabi +// John van Esch +// Gennady Markelov +// + +#include "custom_constitutive/pw_filter_law.h" +#include "geo_mechanics_application_variables.h" +#include "includes/global_variables.h" + +namespace Kratos +{ + +GeoPwFilterLaw::GeoPwFilterLaw() : mNumberOfDimensions{2} {} + +GeoPwFilterLaw::GeoPwFilterLaw(std::size_t NumberOfDimensions) + : mNumberOfDimensions{NumberOfDimensions} +{ + KRATOS_ERROR_IF(mNumberOfDimensions != 1) + << "Got invalid number of dimensions. The dimension has to be 1, but got: " << mNumberOfDimensions + << std::endl; +} + +ConstitutiveLaw::Pointer GeoPwFilterLaw::Clone() const +{ + return Kratos::make_shared(*this); +} + +SizeType GeoPwFilterLaw::WorkingSpaceDimension() { return mNumberOfDimensions; } + +Matrix GeoPwFilterLaw::CalculatePwFilterMatrix(const Properties& rProp, const ProcessInfo& rProcessInfo) const +{ + const double equivalent_radius_square = rProp[CROSS_AREA] / Globals::Pi; + return ScalarMatrix(1, 1, equivalent_radius_square * 0.125); +} + +} // Namespace Kratos diff --git a/applications/GeoMechanicsApplication/custom_constitutive/pw_filter_law.h b/applications/GeoMechanicsApplication/custom_constitutive/pw_filter_law.h new file mode 100644 index 000000000000..3bf8075e5b25 --- /dev/null +++ b/applications/GeoMechanicsApplication/custom_constitutive/pw_filter_law.h @@ -0,0 +1,62 @@ +// KRATOS___ +// // ) ) +// // ___ ___ +// // ____ //___) ) // ) ) +// // / / // // / / +// ((____/ / ((____ ((___/ / MECHANICS +// +// License: geo_mechanics_application/license.txt +// +// Main authors: Mohamed Nabi +// John van Esch +// Gennady Markelov +// + +#pragma once + +#include "includes/constitutive_law.h" + +namespace Kratos +{ + +class RetentionLaw; + +/** + * @class GeoPwFilterLaw + * @ingroup GeoMechanicsApplication + * @brief This class defines the Pw filter for water pressure cases + */ +class KRATOS_API(GEO_MECHANICS_APPLICATION) GeoPwFilterLaw : public ConstitutiveLaw +{ +public: + /// Counted pointer of GeoPwFilterLaw + KRATOS_CLASS_POINTER_DEFINITION(GeoPwFilterLaw); + + GeoPwFilterLaw(); + + explicit GeoPwFilterLaw(SizeType NumberOfDimensions); + + ConstitutiveLaw::Pointer Clone() const override; + + SizeType WorkingSpaceDimension() override; + + Matrix CalculatePwFilterMatrix(const Properties& rProp, const ProcessInfo& rProcessInfo) const; + +private: + SizeType mNumberOfDimensions; + + friend class Serializer; + + void save(Serializer& rSerializer) const override + { + KRATOS_SERIALIZE_SAVE_BASE_CLASS(rSerializer, ConstitutiveLaw) + rSerializer.save("NumberOfDimensions", mNumberOfDimensions); + } + + void load(Serializer& rSerializer) override + { + KRATOS_SERIALIZE_LOAD_BASE_CLASS(rSerializer, ConstitutiveLaw) + rSerializer.load("NumberOfDimensions", mNumberOfDimensions); + } +}; +} // namespace Kratos. From 222b88c98bc9fc392b71860540cd5c0ae7f9ef59 Mon Sep 17 00:00:00 2001 From: mnabideltares Date: Fri, 7 Jun 2024 17:25:42 +0200 Subject: [PATCH 2/2] Added filter for pw --- .../custom_constitutive/pw_filter_law.cpp | 45 ----------- .../custom_constitutive/pw_filter_law.h | 62 --------------- .../transient_Pw_line_element.h | 77 ++++++++++++------- .../custom_retention/retention_law_factory.h | 3 + .../custom_utilities/element_utilities.hpp | 8 +- 5 files changed, 59 insertions(+), 136 deletions(-) delete mode 100644 applications/GeoMechanicsApplication/custom_constitutive/pw_filter_law.cpp delete mode 100644 applications/GeoMechanicsApplication/custom_constitutive/pw_filter_law.h diff --git a/applications/GeoMechanicsApplication/custom_constitutive/pw_filter_law.cpp b/applications/GeoMechanicsApplication/custom_constitutive/pw_filter_law.cpp deleted file mode 100644 index 0f6a66755d02..000000000000 --- a/applications/GeoMechanicsApplication/custom_constitutive/pw_filter_law.cpp +++ /dev/null @@ -1,45 +0,0 @@ -// KRATOS___ -// // ) ) -// // ___ ___ -// // ____ //___) ) // ) ) -// // / / // // / / -// ((____/ / ((____ ((___/ / MECHANICS -// -// License: geo_mechanics_application/license.txt -// -// Main authors: Mohamed Nabi -// John van Esch -// Gennady Markelov -// - -#include "custom_constitutive/pw_filter_law.h" -#include "geo_mechanics_application_variables.h" -#include "includes/global_variables.h" - -namespace Kratos -{ - -GeoPwFilterLaw::GeoPwFilterLaw() : mNumberOfDimensions{2} {} - -GeoPwFilterLaw::GeoPwFilterLaw(std::size_t NumberOfDimensions) - : mNumberOfDimensions{NumberOfDimensions} -{ - KRATOS_ERROR_IF(mNumberOfDimensions != 1) - << "Got invalid number of dimensions. The dimension has to be 1, but got: " << mNumberOfDimensions - << std::endl; -} - -ConstitutiveLaw::Pointer GeoPwFilterLaw::Clone() const -{ - return Kratos::make_shared(*this); -} - -SizeType GeoPwFilterLaw::WorkingSpaceDimension() { return mNumberOfDimensions; } - -Matrix GeoPwFilterLaw::CalculatePwFilterMatrix(const Properties& rProp, const ProcessInfo& rProcessInfo) const -{ - const double equivalent_radius_square = rProp[CROSS_AREA] / Globals::Pi; - return ScalarMatrix(1, 1, equivalent_radius_square * 0.125); -} - -} // Namespace Kratos diff --git a/applications/GeoMechanicsApplication/custom_constitutive/pw_filter_law.h b/applications/GeoMechanicsApplication/custom_constitutive/pw_filter_law.h deleted file mode 100644 index 3bf8075e5b25..000000000000 --- a/applications/GeoMechanicsApplication/custom_constitutive/pw_filter_law.h +++ /dev/null @@ -1,62 +0,0 @@ -// KRATOS___ -// // ) ) -// // ___ ___ -// // ____ //___) ) // ) ) -// // / / // // / / -// ((____/ / ((____ ((___/ / MECHANICS -// -// License: geo_mechanics_application/license.txt -// -// Main authors: Mohamed Nabi -// John van Esch -// Gennady Markelov -// - -#pragma once - -#include "includes/constitutive_law.h" - -namespace Kratos -{ - -class RetentionLaw; - -/** - * @class GeoPwFilterLaw - * @ingroup GeoMechanicsApplication - * @brief This class defines the Pw filter for water pressure cases - */ -class KRATOS_API(GEO_MECHANICS_APPLICATION) GeoPwFilterLaw : public ConstitutiveLaw -{ -public: - /// Counted pointer of GeoPwFilterLaw - KRATOS_CLASS_POINTER_DEFINITION(GeoPwFilterLaw); - - GeoPwFilterLaw(); - - explicit GeoPwFilterLaw(SizeType NumberOfDimensions); - - ConstitutiveLaw::Pointer Clone() const override; - - SizeType WorkingSpaceDimension() override; - - Matrix CalculatePwFilterMatrix(const Properties& rProp, const ProcessInfo& rProcessInfo) const; - -private: - SizeType mNumberOfDimensions; - - friend class Serializer; - - void save(Serializer& rSerializer) const override - { - KRATOS_SERIALIZE_SAVE_BASE_CLASS(rSerializer, ConstitutiveLaw) - rSerializer.save("NumberOfDimensions", mNumberOfDimensions); - } - - void load(Serializer& rSerializer) override - { - KRATOS_SERIALIZE_LOAD_BASE_CLASS(rSerializer, ConstitutiveLaw) - rSerializer.load("NumberOfDimensions", mNumberOfDimensions); - } -}; -} // namespace Kratos. diff --git a/applications/GeoMechanicsApplication/custom_elements/transient_Pw_line_element.h b/applications/GeoMechanicsApplication/custom_elements/transient_Pw_line_element.h index 156a67a0d0e3..333e6fe12a7a 100644 --- a/applications/GeoMechanicsApplication/custom_elements/transient_Pw_line_element.h +++ b/applications/GeoMechanicsApplication/custom_elements/transient_Pw_line_element.h @@ -74,11 +74,13 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) TransientPwLineElement : public Elem std::transform(dN_dX_container.begin(), dN_dX_container.end(), det_J_container.begin(), dN_dX_container.begin(), std::divides<>()); const Matrix& r_N_container = GetGeometry().ShapeFunctionsValues(GetIntegrationMethod()); + Vector projected_gravity = this->CalculateProjectedGravityAtIntegrationPoints(r_N_container); + const auto integration_coefficients = CalculateIntegrationCoefficients(det_J_container); const auto permeability_matrix = CalculatePermeabilityMatrix(dN_dX_container, integration_coefficients); - const auto compressibility_matrix = CalculateCompressibilityMatrix(r_N_container, integration_coefficients); + const auto compressibility_matrix = CalculateCompressibilityMatrix(r_N_container, integration_coefficients, projected_gravity); - const auto fluid_body_vector = CalculateFluidBodyVector(r_N_container, dN_dX_container, integration_coefficients); + const auto fluid_body_vector = CalculateFluidBodyVector(r_N_container, dN_dX_container, integration_coefficients, projected_gravity); AddContributionsToLhsMatrix(rLeftHandSideMatrix, permeability_matrix, compressibility_matrix, rCurrentProcessInfo[DT_PRESSURE_COEFFICIENT]); AddContributionsToRhsVector(rRightHandSideVector, permeability_matrix, compressibility_matrix, fluid_body_vector); @@ -241,7 +243,8 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) TransientPwLineElement : public Elem } BoundedMatrix CalculateCompressibilityMatrix(const Matrix& rNContainer, - const Vector& rIntegrationCoefficients) const + const Vector& rIntegrationCoefficients, + const Vector& rProjectedGravity) const { const auto& r_properties = GetProperties(); RetentionLaw::Parameters parameters(r_properties); @@ -252,7 +255,7 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) TransientPwLineElement : public Elem integration_point_index < GetGeometry().IntegrationPointsNumber(GetIntegrationMethod()); ++integration_point_index) { const auto N = Vector{row(rNContainer, integration_point_index)}; - const double BiotModulusInverse = CalculateBiotModulusInverse(integration_point_index); + const double BiotModulusInverse = CalculateBiotModulusInverse(integration_point_index, rProjectedGravity); result += GeoTransportEquationUtilities::CalculateCompressibilityMatrix( N, BiotModulusInverse, rIntegrationCoefficients[integration_point_index]); } @@ -284,17 +287,26 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) TransientPwLineElement : public Elem } - double CalculateBiotModulusInverse(const unsigned int integrationPointIndex) const + double CalculateBiotModulusInverse(const unsigned int integrationPointIndex, const Vector& rProjectedGravity) const { const auto& r_properties = GetProperties(); + + double result = 0.0; + if (r_properties[RETENTION_LAW] == "PressureFilterLaw") { + result = 1.0 / (r_properties[DENSITY_WATER] * rProjectedGravity[integrationPointIndex] * + r_properties[PIPE_ELEMENT_LENGTH]) + + 1.0 / r_properties[BULK_MODULUS_FLUID]; + return result; + } + const double biot_coefficient = r_properties[BIOT_COEFFICIENT]; double bulk_fluid = TINY; if (!r_properties[IGNORE_UNDRAINED]) { bulk_fluid = r_properties[BULK_MODULUS_FLUID]; } - double result = (biot_coefficient - r_properties[POROSITY]) / r_properties[BULK_MODULUS_SOLID] + - r_properties[POROSITY] / bulk_fluid; + result = (biot_coefficient - r_properties[POROSITY]) / r_properties[BULK_MODULUS_SOLID] + + r_properties[POROSITY] / bulk_fluid; RetentionLaw::Parameters RetentionParameters(GetProperties()); const double degree_of_saturation = mRetentionLawVector[integrationPointIndex]->CalculateSaturation(RetentionParameters); @@ -327,38 +339,21 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) TransientPwLineElement : public Elem array_1d CalculateFluidBodyVector( const Matrix& rNContainer, const GeometryType::ShapeFunctionsGradientsType& rShapeFunctionGradients, - const Vector& rIntegrationCoefficients) const + const Vector& rIntegrationCoefficients, + const Vector& rProjectedGravity) const { - const std::size_t number_integration_points = GetGeometry().IntegrationPointsNumber(GetIntegrationMethod()); - GeometryType::JacobiansType J_container; - J_container.resize(number_integration_points, false); - for (std::size_t i = 0; i < number_integration_points; ++i) { - J_container[i].resize(GetGeometry().WorkingSpaceDimension(), GetGeometry().LocalSpaceDimension(), false); - } - GetGeometry().Jacobian(J_container, this->GetIntegrationMethod()); - const auto& r_properties = GetProperties(); BoundedMatrix constitutive_matrix; GeoElementUtilities::FillPermeabilityMatrix(constitutive_matrix, r_properties); RetentionLaw::Parameters RetentionParameters(GetProperties()); - - array_1d volume_acceleration; - GeoElementUtilities::GetNodalVariableVector(volume_acceleration, GetGeometry(), VOLUME_ACCELERATION); - array_1d body_acceleration; + array_1d projected_gravity = ZeroVector(1); array_1d fluid_body_vector = ZeroVector(TNumNodes); for (unsigned int integration_point_index = 0; integration_point_index < GetGeometry().IntegrationPointsNumber(GetIntegrationMethod()); ++integration_point_index) { - GeoElementUtilities::InterpolateVariableWithComponents( - body_acceleration, rNContainer, volume_acceleration, integration_point_index); - - array_1d tangent_vector = column(J_container[integration_point_index], 0); - tangent_vector /= norm_2(tangent_vector); - - array_1d projected_gravity = ZeroVector(1); - projected_gravity(0) = MathUtils::Dot(tangent_vector, body_acceleration); + projected_gravity(0) = rProjectedGravity[integration_point_index]; const auto N = Vector{row(rNContainer, integration_point_index)}; double RelativePermeability = mRetentionLawVector[integration_point_index]->CalculateRelativePermeability(RetentionParameters); fluid_body_vector += r_properties[DENSITY_WATER] * RelativePermeability @@ -369,6 +364,32 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) TransientPwLineElement : public Elem } + Vector CalculateProjectedGravityAtIntegrationPoints(const Matrix& rNContainer) + { + const std::size_t number_integration_points = GetGeometry().IntegrationPointsNumber(GetIntegrationMethod()); + GeometryType::JacobiansType J_container; + J_container.resize(number_integration_points, false); + for (std::size_t i = 0; i < number_integration_points; ++i) { + J_container[i].resize(GetGeometry().WorkingSpaceDimension(), GetGeometry().LocalSpaceDimension(), false); + } + GetGeometry().Jacobian(J_container, this->GetIntegrationMethod()); + + array_1d volume_acceleration; + GeoElementUtilities::GetNodalVariableVector(volume_acceleration, GetGeometry(), VOLUME_ACCELERATION); + array_1d body_acceleration; + + Vector projected_gravity = ZeroVector(number_integration_points); + + for (unsigned int integration_point_index = 0; integration_point_index < number_integration_points; ++integration_point_index) { + GeoElementUtilities::InterpolateVariableWithComponents(body_acceleration, rNContainer, volume_acceleration, integration_point_index); + array_1d tangent_vector = column(J_container[integration_point_index], 0); + tangent_vector /= norm_2(tangent_vector); + projected_gravity(integration_point_index) = MathUtils::Dot(tangent_vector, body_acceleration); + } + return projected_gravity; + } + + [[nodiscard]] DofsVectorType GetDofs() const { return Geo::DofUtilities::ExtractDofsFromNodes(GetGeometry(), WATER_PRESSURE); diff --git a/applications/GeoMechanicsApplication/custom_retention/retention_law_factory.h b/applications/GeoMechanicsApplication/custom_retention/retention_law_factory.h index 158cee8aacd1..7422a43ae3dd 100644 --- a/applications/GeoMechanicsApplication/custom_retention/retention_law_factory.h +++ b/applications/GeoMechanicsApplication/custom_retention/retention_law_factory.h @@ -56,6 +56,9 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) RetentionLawFactory if (RetentionLawName == "SaturatedBelowPhreaticLevelLaw") return make_unique(); + if (RetentionLawName == "PressureFilterLaw") + return make_unique(); + KRATOS_ERROR << "Undefined RETENTION_LAW! " << RetentionLawName << std::endl; diff --git a/applications/GeoMechanicsApplication/custom_utilities/element_utilities.hpp b/applications/GeoMechanicsApplication/custom_utilities/element_utilities.hpp index c37b10db6816..27a9195f44b9 100644 --- a/applications/GeoMechanicsApplication/custom_utilities/element_utilities.hpp +++ b/applications/GeoMechanicsApplication/custom_utilities/element_utilities.hpp @@ -19,6 +19,7 @@ // Application includes #include "geo_mechanics_application_variables.h" +#include "includes/global_variables.h" namespace Kratos { @@ -142,7 +143,12 @@ class GeoElementUtilities const Element::PropertiesType& Prop) { // 1D - rPermeabilityMatrix(0, 0) = Prop[PERMEABILITY_XX]; + if (Prop[RETENTION_LAW] == "PressureFilterLaw") { + const double equivalent_radius_square = Prop[CROSS_AREA] / Globals::Pi; + rPermeabilityMatrix(0, 0) = equivalent_radius_square * 0.125; + } else { + rPermeabilityMatrix(0, 0) = Prop[PERMEABILITY_XX]; + } } static inline void FillPermeabilityMatrix(BoundedMatrix& rPermeabilityMatrix,