[GeoMechanicsApplication] Implement well constitutive behaviour for line pressure element

applications/GeoMechanicsApplication/custom_elements/README.md

@@ -175,5 +175,54 @@ where
 The superscripts $^l$ and $^r$ for Robin boundary condition indicate the left hands side (matrix) and right hand side (vector), respectively. The superscripts $^e$ and $^{ep}$ for $\Omega$ and $\Gamma$ indicate values in the element volume and perpendicular to element boundaries, respectively.  
 
 
+# Pressure Filter Line Element
+For laminar flow along the axis of the well and uniformly distributed storage over the length of the well bore (Diersch, 2014)
+
+```math
+\pi R^2 \left( \frac{1}{l_w} + \rho_0 g \beta \right) \frac{\partial h}{\partial t} - \pi R^2 K_w \frac{\partial}{\partial y} \left[ f_{\mu} \left( \frac{\partial h}{\partial y} + \chi e \right)\right] = -Q_w \delta \left( y - y_w \right)
+```
+
+where
+
+```math
+K_w = \frac{R^2 \rho_0 g}{8 \mu_0} \quad \quad \quad \quad f_u = \frac{\mu_0}{\mu} \quad \quad \quad \quad \chi = \frac{\rho - \rho_0}{\rho_0} \quad \quad \quad \quad h = \frac{p}{\rho_0 g} + y
+```
+
+- $Q_w$			= pumping rate sink $\mathrm{\left[ m^3/s \right]}$
+- $t$				= time $\mathrm{\left[ s \right]}$
+- $y$				= vertical coordinate $\mathrm{\left[ m \right]}$
+- $y_w$				= location of the discharge point $\mathrm{\left[ m \right]}$
+- $h$				= hydraulic head in the well $\mathrm{\left[ m \right]}$
+- $l_w$				= total length of liquid filled well bore $\mathrm{\left[ m \right]}$
+- $R$				= radius of the well casing $\mathrm{\left[ m \right]}$
+- $K_w$				= Hagen-Poiseuille permeability $\mathrm{\left[ m \right]}$
+- $\delta$			= Dirichlet delta function $\mathrm{\left[ - \right]}$
+- $\beta$			= compressibility of the liquid $\mathrm{\left[ m^2/N \right]}$
+- $f_{\mu}$			= viscosity relation function of liquid $\mathrm{\left[ - \right]}$
+- $\chi$			= buoyancy coefficient $\mathrm{\left[ - \right]}$
+- $e$				= gravitational unit vector $\mathrm{\left[ - \right]}$
+- $g$				= gravitational acceleration $\mathrm{\left[ m/s^2 \right]}$
+- $\rho_0$			= reference density of the fluid $\mathrm{\left[ kg/m^3 \right]}$
+- $\mu_0$			= reference viscosity of the fluid $\mathrm{\left[ Pas \right]}$
+- $p$				= liquid pressure $\mathrm{\left[ Pa \right]}$
+
+reformulated in pressure and preserving mass
+
+```math
+\rho^w \left( \frac{1}{\rho^w g l_w} + \beta \right) \frac{\partial p}{\partial t} - \rho^w \frac{\partial}{\partial y} \left[ \frac{R^2}{8 \mu} \left( \frac{\partial p}{\partial y} - \rho^w g \right)\right] = -\frac{\rho^w Q_w}{\pi R^2}
+```
+
+The one dimensionalpressure equation with $\alpha = 1$ and $n = 1$ captures the well flow equation.
+
+```math
+\rho^w \beta^* \frac{\partial p}{\partial t} - \frac{\partial}{\partial y} \left[ \frac{\rho^w K^*}{\mu} \left( \frac{\partial p}{\partial y} - \rho^w g \right) \right] = - \rho^w Q^*_w
+```
+
+where the intrinsic permeability and boundary condition are given by
+
+```math
+\beta^* = \frac{1}{\rho^w g l_w} + \beta \quad \quad \quad \quad K^* = \frac{R^2}{8} \quad \quad \quad \quad Q^*_w = \frac{Q_w}{\pi R^2}
+```
+
 ## Bibliography
 Diersch, H.-J. G., 2014. FEFLOW; Finite Element Modeling of Flow, Mass and Heat Transport in Porous and Fractured Media. Springer.

applications/GeoMechanicsApplication/custom_elements/filter_compressibility_calculator.cpp

@@ -0,0 +1,105 @@
+// 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 AnInputProvider)
+    : mInputProvider(std::move(AnInputProvider))
+{
+}
+
+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, CalculateElasticCapacity(projected_gravity_on_integration_points[integration_point_index]),
+            integration_coefficients[integration_point_index]);
+    }
+    return result;
+}
+
+double FilterCompressibilityCalculator::CalculateElasticCapacity(double ProjectedGravity) const
+{
+    const auto&  r_properties = mInputProvider.GetElementProperties();
+    return 1.0 / (r_properties[DENSITY_WATER] * ProjectedGravity * r_properties[FILTER_LENGTH]) +
+           1.0 / r_properties[BULK_MODULUS_FLUID];
+}
+
+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 AnInputProvider);
+
+    Matrix                    LHSContribution() override;
+    Vector                    RHSContribution() override;
+    std::pair<Matrix, Vector> LocalSystemContribution() override;
+
+private:
+    [[nodiscard]] Matrix CalculateCompressibilityMatrix() const;
+    [[nodiscard]] double CalculateElasticCapacity(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,57 +253,69 @@ 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 +=
                 r_properties[DENSITY_WATER] * RelativePermeability *
-                prod(prod(shape_function_gradients[integration_point_index], constitutive_matrix), projected_gravity) *
+                prod(prod(shape_function_gradients[integration_point_index], constitutive_matrix),
+                     ScalarVector(1, projected_gravity[integration_point_index])) *
                 integration_coefficients[integration_point_index] / r_properties[DYNAMIC_VISCOSITY];
         }
         return fluid_body_vector;
     }
 
+    Vector CalculateProjectedGravityAtIntegrationPoints(const Matrix& rNContainer) const
+    {
+        const auto number_integration_points = GetGeometry().IntegrationPointsNumber(GetIntegrationMethod());
+        GeometryType::JacobiansType J_container{number_integration_points};
+        for (auto& container : J_container) {
+            container.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<>::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 +330,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 +371,14 @@ 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_python/geo_mechanics_python_application.cpp

@@ -145,6 +145,8 @@ PYBIND11_MODULE(KratosGeoMechanicsApplication, m)
     KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, UMAT_PARAMETERS)
 
     KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, DT_TEMPERATURE_COEFFICIENT)
+
+    KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, FILTER_LENGTH)
 }
 
 } // namespace Kratos::Python.