[GeoMechanicsApplication] Implement well constitutive behaviour for line pressure element

applications/GeoMechanicsApplication/custom_elements/filter_compressibility_calculator.cpp

@@ -0,0 +1,107 @@
+// KRATOS___
+//     //   ) )
+//    //         ___      ___
+//   //  ____  //___) ) //   ) )
+//  //    / / //       //   / /
+// ((____/ / ((____   ((___/ /  MECHANICS
+//
+//  License:         geo_mechanics_application/license.txt
+//
+//  Main authors:    Mohamed Nabi
+//                   Richard Faasse
+//
+#include "filter_compressibility_calculator.h"
+#include "custom_utilities/transport_equation_utilities.hpp"
+#include "geo_mechanics_application_variables.h"
+
+namespace Kratos
+{
+
+FilterCompressibilityCalculator::FilterCompressibilityCalculator(InputProvider rInputProvider)
+    : mInputProvider(std::move(rInputProvider))
+{
+}
+
+Matrix FilterCompressibilityCalculator::LHSContribution()
+{
+    return LHSContribution(CalculateCompressibilityMatrix());
+}
+
+Vector FilterCompressibilityCalculator::RHSContribution()
+{
+    return RHSContribution(CalculateCompressibilityMatrix());
+}
+
+Vector FilterCompressibilityCalculator::RHSContribution(const Matrix& rCompressibilityMatrix) const
+{
+    return -prod(rCompressibilityMatrix, mInputProvider.GetNodalValues(DT_WATER_PRESSURE));
+}
+
+Matrix FilterCompressibilityCalculator::LHSContribution(const Matrix& rCompressibilityMatrix) const
+{
+    return mInputProvider.GetMatrixScalarFactor() * rCompressibilityMatrix;
+}
+
+std::pair<Matrix, Vector> FilterCompressibilityCalculator::LocalSystemContribution()
+{
+    const auto compressibility_matrix = CalculateCompressibilityMatrix();
+    return {LHSContribution(compressibility_matrix), RHSContribution(compressibility_matrix)};
+}
+
+Matrix FilterCompressibilityCalculator::CalculateCompressibilityMatrix() const
+{
+    const auto& r_N_container            = mInputProvider.GetNContainer();
+    const auto& integration_coefficients = mInputProvider.GetIntegrationCoefficients();
+    const auto& projected_gravity_on_integration_points =
+        mInputProvider.GetProjectedGravityForIntegrationPoints();
+    auto result = Matrix{ZeroMatrix{r_N_container.size2(), r_N_container.size2()}};
+    for (unsigned int integration_point_index = 0;
+         integration_point_index < integration_coefficients.size(); ++integration_point_index) {
+        const auto N = Vector{row(r_N_container, integration_point_index)};
+        result += GeoTransportEquationUtilities::CalculateCompressibilityMatrix(
+            N, CalculateBiotModulusInverse(projected_gravity_on_integration_points[integration_point_index]),
+            integration_coefficients[integration_point_index]);
+    }
+    return result;
+}
+
+double FilterCompressibilityCalculator::CalculateBiotModulusInverse(double ProjectedGravity) const
+{
+    const auto&  r_properties = mInputProvider.GetElementProperties();
+    const double result =
+        1.0 / (r_properties[DENSITY_WATER] * ProjectedGravity * r_properties[PIPE_ELEMENT_LENGTH]) +
+        1.0 / r_properties[BULK_MODULUS_FLUID];
+    return result;
+}
+
+const Properties& FilterCompressibilityCalculator::InputProvider::GetElementProperties() const
+{
+    return mGetElementProperties();
+}
+
+const Matrix& FilterCompressibilityCalculator::InputProvider::GetNContainer() const
+{
+    return mGetNContainer();
+}
+
+Vector FilterCompressibilityCalculator::InputProvider::GetIntegrationCoefficients() const
+{
+    return mGetIntegrationCoefficients();
+}
+
+Vector FilterCompressibilityCalculator::InputProvider::GetProjectedGravityForIntegrationPoints() const
+{
+    return mGetProjectedGravityForIntegrationPoints();
+}
+
+double FilterCompressibilityCalculator::InputProvider::GetMatrixScalarFactor() const
+{
+    return mGetMatrixScalarFactor();
+}
+
+Vector FilterCompressibilityCalculator::InputProvider::GetNodalValues(const Variable<double>& rVariable) const
+{
+    return mGetNodalValues(rVariable);
+}
+
+} // namespace Kratos

applications/GeoMechanicsApplication/custom_elements/filter_compressibility_calculator.h

@@ -0,0 +1,78 @@
+// KRATOS___
+//     //   ) )
+//    //         ___      ___
+//   //  ____  //___) ) //   ) )
+//  //    / / //       //   / /
+// ((____/ / ((____   ((___/ /  MECHANICS
+//
+//  License:         geo_mechanics_application/license.txt
+//
+//  Main authors:    Mohamed Nabi
+//                   Richard Faasse
+
+#pragma once
+
+#include "contribution_calculator.h"
+#include "custom_retention/retention_law.h"
+#include "includes/properties.h"
+#include "includes/ublas_interface.h"
+
+#include <utility>
+#include <vector>
+
+namespace Kratos
+{
+
+class FilterCompressibilityCalculator : public ContributionCalculator
+{
+public:
+    class InputProvider
+    {
+    public:
+        InputProvider(std::function<const Properties&()> GetElementProperties,
+                      std::function<const Matrix&()>                             GetNContainer,
+                      std::function<Vector()> GetIntegrationCoefficients,
+                      std::function<Vector()> GetProjectedGravityForIntegrationPoints,
+                      std::function<double()> GetMatrixScalarFactor,
+                      std::function<Vector(const Variable<double>&)> GetNodalValuesOf)
+            : mGetElementProperties(std::move(GetElementProperties)),
+              mGetNContainer(std::move(GetNContainer)),
+              mGetIntegrationCoefficients(std::move(GetIntegrationCoefficients)),
+              mGetProjectedGravityForIntegrationPoints(std::move(GetProjectedGravityForIntegrationPoints)),
+              mGetMatrixScalarFactor(std::move(GetMatrixScalarFactor)),
+              mGetNodalValues(std::move(GetNodalValuesOf))
+        {
+        }
+
+        [[nodiscard]] const Properties&                         GetElementProperties() const;
+        [[nodiscard]] const Matrix&                             GetNContainer() const;
+        [[nodiscard]] Vector                                    GetIntegrationCoefficients() const;
+        [[nodiscard]] Vector GetProjectedGravityForIntegrationPoints() const;
+        [[nodiscard]] double GetMatrixScalarFactor() const;
+        [[nodiscard]] Vector GetNodalValues(const Variable<double>& rVariable) const;
+
+    private:
+        std::function<const Properties&()>                         mGetElementProperties;
+        std::function<const Matrix&()>                             mGetNContainer;
+        std::function<Vector()>                                    mGetIntegrationCoefficients;
+        std::function<Vector()>                        mGetProjectedGravityForIntegrationPoints;
+        std::function<double()>                        mGetMatrixScalarFactor;
+        std::function<Vector(const Variable<double>&)> mGetNodalValues;
+    };
+
+    explicit FilterCompressibilityCalculator(InputProvider rInputProvider);
+
+    Matrix                    LHSContribution() override;
+    Vector                    RHSContribution() override;
+    std::pair<Matrix, Vector> LocalSystemContribution() override;
+
+private:
+    [[nodiscard]] Matrix CalculateCompressibilityMatrix() const;
+    [[nodiscard]] double CalculateBiotModulusInverse(double ProjectedGravity) const;
+    [[nodiscard]] Vector RHSContribution(const Matrix& rCompressibilityMatrix) const;
+    [[nodiscard]] Matrix LHSContribution(const Matrix& rCompressibilityMatrix) const;
+
+    InputProvider mInputProvider;
+};
+
+} // namespace Kratos

applications/GeoMechanicsApplication/custom_elements/transient_Pw_line_element.h

@@ -18,6 +18,7 @@
 #include "custom_retention/retention_law_factory.h"
 #include "custom_utilities/dof_utilities.h"
 #include "custom_utilities/element_utilities.hpp"
+#include "filter_compressibility_calculator.h"
 #include "geo_mechanics_application_variables.h"
 #include "includes/cfd_variables.h"
 #include "includes/element.h"
@@ -252,40 +253,19 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) TransientPwLineElement : public Elem
 
         const auto integration_coefficients = CalculateIntegrationCoefficients(det_J_container);
 
-        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<double, 1, 1> constitutive_matrix;
         GeoElementUtilities::FillPermeabilityMatrix(constitutive_matrix, r_properties);
 
         RetentionLaw::Parameters RetentionParameters(GetProperties());
 
-        array_1d<double, TNumNodes * TDim> volume_acceleration;
-        GeoElementUtilities::GetNodalVariableVector<TDim, TNumNodes>(
-            volume_acceleration, GetGeometry(), VOLUME_ACCELERATION);
+        const auto projected_gravity = CalculateProjectedGravityAtIntegrationPoints(N_container);
 
         array_1d<double, TNumNodes> fluid_body_vector = ZeroVector(TNumNodes);
         for (unsigned int integration_point_index = 0;
              integration_point_index < GetGeometry().IntegrationPointsNumber(GetIntegrationMethod());
              ++integration_point_index) {
-            array_1d<double, TDim> body_acceleration;
-            GeoElementUtilities::InterpolateVariableWithComponents<TDim, TNumNodes>(
-                body_acceleration, N_container, volume_acceleration, integration_point_index);
-
-            array_1d<double, TDim> tangent_vector = column(J_container[integration_point_index], 0);
-            tangent_vector /= norm_2(tangent_vector);
-
-            array_1d<double, 1> projected_gravity = ZeroVector(1);
-            projected_gravity(0) = MathUtils<double>::Dot(tangent_vector, body_acceleration);
-            const auto N         = Vector{row(N_container, integration_point_index)};
+            const auto N = Vector{row(N_container, integration_point_index)};
             double     RelativePermeability =
                 mRetentionLawVector[integration_point_index]->CalculateRelativePermeability(RetentionParameters);
             fluid_body_vector +=
@@ -296,13 +276,47 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) TransientPwLineElement : public Elem
         return fluid_body_vector;
     }
 
+    Vector CalculateProjectedGravityAtIntegrationPoints(const Matrix& rNContainer) const
+    {
+        const auto 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<double, TNumNodes * TDim> volume_acceleration;
+        GeoElementUtilities::GetNodalVariableVector<TDim, TNumNodes>(
+            volume_acceleration, GetGeometry(), VOLUME_ACCELERATION);
+        array_1d<double, TDim> 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<TDim, TNumNodes>(
+                body_acceleration, rNContainer, volume_acceleration, integration_point_index);
+            array_1d<double, TDim> tangent_vector = column(J_container[integration_point_index], 0);
+            tangent_vector /= norm_2(tangent_vector);
+            projected_gravity(integration_point_index) =
+                MathUtils<double>::Dot(tangent_vector, body_acceleration);
+        }
+        return projected_gravity;
+    }
+
     std::unique_ptr<ContributionCalculator> CreateCalculator(const CalculationContribution& rContribution,
                                                              const ProcessInfo& rCurrentProcessInfo)
     {
         switch (rContribution) {
         case CalculationContribution::Permeability:
             return std::make_unique<PermeabilityCalculator>(CreatePermeabilityInputProvider());
         case CalculationContribution::Compressibility:
+            if (GetProperties()[RETENTION_LAW] == "PressureFilterLaw") {
+                return std::make_unique<FilterCompressibilityCalculator>(CreateFilterCompressibilityInputProvider(rCurrentProcessInfo));
+            }
             return std::make_unique<CompressibilityCalculator>(CreateCompressibilityInputProvider(rCurrentProcessInfo));
         default:
             KRATOS_ERROR << "Unknown contribution" << std::endl;
@@ -317,6 +331,14 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) TransientPwLineElement : public Elem
             MakeNodalVariableGetter());
     }
 
+    FilterCompressibilityCalculator::InputProvider CreateFilterCompressibilityInputProvider(const ProcessInfo& rCurrentProcessInfo)
+    {
+        return FilterCompressibilityCalculator::InputProvider(
+            MakePropertiesGetter(), MakeNContainerGetter(),
+            MakeIntegrationCoefficientsGetter(), MakeProjectedGravityForIntegrationPointsGetter(), MakeMatrixScalarFactorGetter(rCurrentProcessInfo),
+            MakeNodalVariableGetter());
+    }
+
     PermeabilityCalculator::InputProvider CreatePermeabilityInputProvider()
     {
         return PermeabilityCalculator::InputProvider(
@@ -350,6 +372,11 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) TransientPwLineElement : public Elem
         };
     }
 
+    auto MakeProjectedGravityForIntegrationPointsGetter()
+    {
+        return [this]() -> Vector {return CalculateProjectedGravityAtIntegrationPoints(GetGeometry().ShapeFunctionsValues(GetIntegrationMethod()));};
+    }
+
     static auto MakeMatrixScalarFactorGetter(const ProcessInfo& rCurrentProcessInfo)
     {
         return [&rCurrentProcessInfo]() { return rCurrentProcessInfo[DT_PRESSURE_COEFFICIENT]; };

applications/GeoMechanicsApplication/custom_retention/retention_law_factory.h

@@ -46,20 +46,22 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) RetentionLawFactory
     {
         if (rMaterialProperties.Has(RETENTION_LAW)) {
             const std::string& RetentionLawName = rMaterialProperties[RETENTION_LAW];
-            if (RetentionLawName == "VanGenuchtenLaw") return make_unique<VanGenuchtenLaw>();
+            if (RetentionLawName == "VanGenuchtenLaw") return std::make_unique<VanGenuchtenLaw>();
 
-            if (RetentionLawName == "SaturatedLaw") return make_unique<SaturatedLaw>();
+            if (RetentionLawName == "SaturatedLaw") return std::make_unique<SaturatedLaw>();
 
             if (RetentionLawName == "SaturatedBelowPhreaticLevelLaw")
-                return make_unique<SaturatedBelowPhreaticLevelLaw>();
+                return std::make_unique<SaturatedBelowPhreaticLevelLaw>();
+
+            if (RetentionLawName == "PressureFilterLaw") return std::make_unique<SaturatedLaw>();
 
             KRATOS_ERROR << "Undefined RETENTION_LAW! " << RetentionLawName << std::endl;
 
             return nullptr;
         }
 
         // default is saturated law
-        return make_unique<SaturatedLaw>();
+        return std::make_unique<SaturatedLaw>();
     }
 
 }; // Class RetentionLawFactory

applications/GeoMechanicsApplication/custom_utilities/element_utilities.cpp

@@ -35,7 +35,12 @@ void GeoElementUtilities::FillPermeabilityMatrix(BoundedMatrix<double, 1, 1>&
                                                  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];
+    }
 }
 
 void GeoElementUtilities::FillPermeabilityMatrix(BoundedMatrix<double, 2, 2>&   rPermeabilityMatrix,

applications/GeoMechanicsApplication/tests/test_pressure_line_element.py

@@ -59,5 +59,11 @@ def test_vertical_line_element3D2N(self):
     def test_vertical_line_element3D3N(self):
         self.check_water_pressure("test_vertical_line_element3D3N", self.etalon_value1)
 
+    def test_oblique_filter_element2D3N(self):
+        self.check_water_pressure("test_oblique_filter_element2D3N", self.etalon_value1)
+
+    def test_oblique_filter_element3D3N(self):
+        self.check_water_pressure("test_oblique_filter_element3D3N", self.etalon_value1)
+
 if __name__ == '__main__':
     KratosUnittest.main()

applications/GeoMechanicsApplication/tests/test_pressure_line_element/README.md

@@ -23,4 +23,13 @@ These results are associated with the final time step after the solution reaches
 
 $P = 10000 y$
 
-In this test case, the result at node number 3 at location $y = -2 \mathrm{[m]}$ is compared with the analytical solution. The value of the pressure at node 3 is -20,000 $\mathrm{[Pa]}$
+In this test case, the result at node number 3 at location $y = -2 \mathrm{[m]}$ is compared with the analytical solution. The value of the pressure at node 3 is -20,000 $\mathrm{[Pa]}$
+
+
+# Test Cases for Water Pore Pressure Line Element with Filter as Retention Law
+Two test cases are repeated with modification for retention law in the material parameters. They are test cases on elements of 2D3N and 3D3N (with slope of 45 degrees). Here, the retention law is replaced by "PressureFilterLaw". The results are then compared with analytical results.
+The results for both test configurations are identical and they are also identical to the results of the pressure line element.
+
+The picture below illustrates the pressure profile resulting from the simulation (as an example the 2D3N test is shown below).
+
+<img src="documentation_data/test_pressure_filter_element_2d3n_result.png" alt="Pressure along depth at the last time step" title="Pressure along the depth at the last time step" width="600">