Add a linear elastic constitutive law for line interfaces (#12654)

Added a linear elastic constitutive law that uses an incremental formulation for line interface elements. This new constitutive law will eventually replace the existing classes `LinearElastic2DInterfaceLaw` and `LinearElastic3DInterfaceLaw`. So far it aims at 2D models only. In addition to the `README` file, the unit tests of class `GeoIncrementalLinearElasticInterfaceLaw` document its intended behavior.

Other changes in this PR include adding two new variables to the GeoMechanicsApplication: `INTERFACE_NORMAL_STIFFNESS` and `INTERFACE_SHEAR_STIFFNESS`.

applications/GeoMechanicsApplication/custom_constitutive/README.md

@@ -0,0 +1,49 @@
+# Constitutive laws
+
+
+## Incremental linear elastic interface law
+
+This constitutive law for an interface element linearly relates increments of tractions $\Delta \tau$ to increments of relative displacement $\Delta \Delta u$.
+Relative displacement for interface element is the differential motion between the two sides of the interface. As a
+consequence the relative displacement has unit of length $[\mathrm{L}]$ and the stiffness has unit of force over cubic length $[\mathrm{F/L^3}]$.
+
+### Relative displacement and traction
+In 2D plane strain $y$ is the opening/closing direction of the interface, while differential motion in the tangential direction
+gives shear. Similar to the continuum, the normal behavior is placed first in the relative displacement and traction vectors. The shear behavior
+follows after the normal motion or traction in these vectors.
+
+$$ \Delta u = \begin{bmatrix} \Delta u_y \\ \Delta u_x \end{bmatrix} $$
+
+$$ \tau = \begin{bmatrix} \tau_{yy} \\ \tau_{xy} \end{bmatrix} $$
+
+Where:
+* $\Delta u_y$: Relative displacement in the $y$-direction (normal to the interface).
+* $\Delta u_x$: Relative displacement in the $x$-direction (tangential to the interface).
+* $\tau_{yy}$: Traction in the $y$-direction (normal to the interface).
+* $\tau_{xy}$: Shear traction in the $x$-direction (tangential to the interface).
+
+### 2D Elastic constitutive tensor
+
+The elastic behavior of the interface is characterized by the 2D elastic constitutive tensor $C$, which relates the traction and relative displacement vectors. The constitutive tensor in 2D is expressed as:
+
+$$ C = \begin{bmatrix} C_{yy} & 0     \\
+                       0     & C_{xy} \end{bmatrix}$$
+
+Where:
+* $C$: Represents the 2D elastic constitutive tensor.
+* $C_{yy}$: Represents the `INTERFACE_NORMAL_STIFFNESS`, which characterizes the stiffness in the normal direction.
+* $C_{xy}$: Represents the `INTERFACE_SHEAR_STIFFNESS`, which characterizes the stiffness in the tangential direction.
+
+Both stiffness values have dimension $\mathrm{F/L^3}$.
+
+### Incremental relation
+
+The relation between the traction and the relative displacement for a given time increment is given by the following incremental equation:
+
+$$ \tau_{t + \Delta t} = \tau_t + C \cdot \Delta \Delta u $$
+
+Where:
+* $\tau_{t + \Delta t}$: The traction at the updated time $t + \Delta t$.
+* $\tau_t$: The traction at the current time $t$.
+* $C$: The elastic constitutive tensor.
+* $\Delta \Delta u$: The incremental relative displacement vector.

applications/GeoMechanicsApplication/custom_constitutive/incremental_linear_elastic_interface_law.cpp

@@ -0,0 +1,131 @@
+// KRATOS___
+//     //   ) )
+//    //         ___      ___
+//   //  ____  //___) ) //   ) )
+//  //    / / //       //   / /
+// ((____/ / ((____   ((___/ /  MECHANICS
+//
+//  License:         geo_mechanics_application/license.txt
+//
+//  Main authors:    Wijtze Pieter Kikstra
+//                   Anne van de Graaf
+//
+
+#include "incremental_linear_elastic_interface_law.h"
+#include "custom_geometries/line_interface_geometry.h"
+#include "geo_mechanics_application_constants.h"
+#include "geo_mechanics_application_variables.h"
+
+namespace Kratos
+{
+
+ConstitutiveLaw::Pointer GeoIncrementalLinearElasticInterfaceLaw::Clone() const
+{
+    return std::make_shared<GeoIncrementalLinearElasticInterfaceLaw>(*this);
+}
+
+ConstitutiveLaw::SizeType GeoIncrementalLinearElasticInterfaceLaw::WorkingSpaceDimension()
+{
+    return 2;
+}
+
+ConstitutiveLaw::SizeType GeoIncrementalLinearElasticInterfaceLaw::GetStrainSize() const
+{
+    return VOIGT_SIZE_2D_INTERFACE;
+}
+
+Vector& GeoIncrementalLinearElasticInterfaceLaw::GetValue(const Variable<Vector>& rThisVariable, Vector& rValue)
+{
+    if (rThisVariable == STRAIN) {
+        rValue = mPreviousRelativeDisplacement;
+    } else if (rThisVariable == CAUCHY_STRESS_VECTOR) {
+        rValue = mPreviousTraction;
+    }
+
+    return rValue;
+}
+
+ConstitutiveLaw::StressMeasure GeoIncrementalLinearElasticInterfaceLaw::GetStressMeasure()
+{
+    return ConstitutiveLaw::StressMeasure_Cauchy;
+}
+
+bool GeoIncrementalLinearElasticInterfaceLaw::IsIncremental() { return true; }
+
+void GeoIncrementalLinearElasticInterfaceLaw::InitializeMaterial(const Properties&,
+                                                                 const ConstitutiveLaw::GeometryType&,
+                                                                 const Vector&)
+{
+    mPreviousRelativeDisplacement =
+        HasInitialState() ? GetInitialState().GetInitialStrainVector() : ZeroVector{GetStrainSize()};
+    mPreviousTraction =
+        HasInitialState() ? GetInitialState().GetInitialStressVector() : ZeroVector{GetStrainSize()};
+}
+
+void GeoIncrementalLinearElasticInterfaceLaw::CalculateMaterialResponseCauchy(ConstitutiveLaw::Parameters& rValues)
+{
+    rValues.GetStressVector() =
+        mPreviousTraction +
+        prod(MakeConstitutiveMatrix(rValues.GetMaterialProperties()[INTERFACE_NORMAL_STIFFNESS],
+                                    rValues.GetMaterialProperties()[INTERFACE_SHEAR_STIFFNESS]),
+             rValues.GetStrainVector() - mPreviousRelativeDisplacement);
+}
+
+bool GeoIncrementalLinearElasticInterfaceLaw::RequiresInitializeMaterialResponse() { return false; }
+
+void GeoIncrementalLinearElasticInterfaceLaw::FinalizeMaterialResponseCauchy(ConstitutiveLaw::Parameters& rValues)
+{
+    mPreviousRelativeDisplacement = rValues.GetStrainVector();
+    mPreviousTraction             = rValues.GetStressVector();
+}
+
+int GeoIncrementalLinearElasticInterfaceLaw::Check(const Properties& rMaterialProperties,
+                                                   const ConstitutiveLaw::GeometryType& rElementGeometry,
+                                                   const ProcessInfo& rCurrentProcessInfo) const
+{
+    const auto result = BaseType::Check(rMaterialProperties, rElementGeometry, rCurrentProcessInfo);
+
+    KRATOS_ERROR_IF_NOT(rMaterialProperties.Has(INTERFACE_NORMAL_STIFFNESS))
+        << "No interface normal stiffness is defined" << std::endl;
+
+    KRATOS_ERROR_IF_NOT(rMaterialProperties[INTERFACE_NORMAL_STIFFNESS] > 0.0)
+        << "Interface normal stiffness must be positive, but got "
+        << rMaterialProperties[INTERFACE_NORMAL_STIFFNESS] << std::endl;
+
+    KRATOS_ERROR_IF_NOT(rMaterialProperties.Has(INTERFACE_SHEAR_STIFFNESS))
+        << "No interface shear stiffness is defined" << std::endl;
+
+    KRATOS_ERROR_IF_NOT(rMaterialProperties[INTERFACE_SHEAR_STIFFNESS] > 0.0)
+        << "Interface shear stiffness must be positive, but got "
+        << rMaterialProperties[INTERFACE_SHEAR_STIFFNESS] << std::endl;
+
+    KRATOS_ERROR_IF_NOT(dynamic_cast<const LineInterfaceGeometry*>(&rElementGeometry))
+        << "Expected a line interface geometry, but got " << rElementGeometry.Info() << std::endl;
+
+    return result;
+}
+
+void GeoIncrementalLinearElasticInterfaceLaw::save(Serializer& rSerializer) const
+{
+    KRATOS_SERIALIZE_SAVE_BASE_CLASS(rSerializer, BaseType)
+    rSerializer.save("PreviousRelativeDisplacement", mPreviousRelativeDisplacement);
+    rSerializer.save("PreviousTraction", mPreviousTraction);
+}
+
+void GeoIncrementalLinearElasticInterfaceLaw::load(Serializer& rSerializer)
+{
+    KRATOS_SERIALIZE_LOAD_BASE_CLASS(rSerializer, BaseType)
+    rSerializer.load("PreviousRelativeDisplacement", mPreviousRelativeDisplacement);
+    rSerializer.load("PreviousTraction", mPreviousTraction);
+}
+
+Matrix GeoIncrementalLinearElasticInterfaceLaw::MakeConstitutiveMatrix(double NormalStiffness,
+                                                                       double ShearStiffness) const
+{
+    auto result  = Matrix{ZeroMatrix{GetStrainSize(), GetStrainSize()}};
+    result(0, 0) = NormalStiffness;
+    result(1, 1) = ShearStiffness;
+    return result;
+}
+
+} // namespace Kratos

applications/GeoMechanicsApplication/custom_constitutive/incremental_linear_elastic_interface_law.h

@@ -0,0 +1,53 @@
+// KRATOS___
+//     //   ) )
+//    //         ___      ___
+//   //  ____  //___) ) //   ) )
+//  //    / / //       //   / /
+// ((____/ / ((____   ((___/ /  MECHANICS
+//
+//  License:         geo_mechanics_application/license.txt
+//
+//  Main authors:    Wijtze Pieter Kikstra
+//                   Anne van de Graaf
+//
+
+#pragma once
+
+#include "includes/constitutive_law.h"
+#include "includes/kratos_export_api.h"
+
+namespace Kratos
+{
+
+class KRATOS_API(GEO_MECHANICS_APPLICATION) GeoIncrementalLinearElasticInterfaceLaw : public ConstitutiveLaw
+{
+public:
+    using BaseType = ConstitutiveLaw;
+
+    [[nodiscard]] Pointer  Clone() const override;
+    SizeType               WorkingSpaceDimension() override;
+    [[nodiscard]] SizeType GetStrainSize() const override;
+    Vector&                GetValue(const Variable<Vector>& rThisVariable, Vector& rValue) override;
+    using BaseType::GetValue;
+    StressMeasure GetStressMeasure() override;
+    bool          IsIncremental() override;
+    void InitializeMaterial(const Properties&, const GeometryType&, const Vector&) override;
+    void CalculateMaterialResponseCauchy(Parameters& rValues) override;
+    bool RequiresInitializeMaterialResponse() override;
+    void FinalizeMaterialResponseCauchy(Parameters& rValues) override;
+    int  Check(const Properties&   rMaterialProperties,
+               const GeometryType& rElementGeometry,
+               const ProcessInfo&  rCurrentProcessInfo) const override;
+
+private:
+    friend class Serializer;
+    void save(Serializer& rSerializer) const override;
+    void load(Serializer& rSerializer) override;
+
+    [[nodiscard]] Matrix MakeConstitutiveMatrix(double NormalStiffness, double ShearStiffness) const;
+
+    Vector mPreviousRelativeDisplacement;
+    Vector mPreviousTraction;
+};
+
+} // namespace Kratos

applications/GeoMechanicsApplication/geo_mechanics_application.cpp

@@ -591,5 +591,8 @@ void KratosGeoMechanicsApplication::Register()
 
     KRATOS_REGISTER_VARIABLE(STRAINS_OF_PIECEWISE_LINEAR_LAW)
     KRATOS_REGISTER_VARIABLE(STRESSES_OF_PIECEWISE_LINEAR_LAW)
+
+    KRATOS_REGISTER_VARIABLE(INTERFACE_NORMAL_STIFFNESS)
+    KRATOS_REGISTER_VARIABLE(INTERFACE_SHEAR_STIFFNESS)
 }
 } // namespace Kratos.

applications/GeoMechanicsApplication/geo_mechanics_application_variables.cpp

@@ -243,4 +243,7 @@ KRATOS_CREATE_VARIABLE(double, STATE_VARIABLE_50)
 KRATOS_CREATE_VARIABLE(Vector, STRAINS_OF_PIECEWISE_LINEAR_LAW)
 KRATOS_CREATE_VARIABLE(Vector, STRESSES_OF_PIECEWISE_LINEAR_LAW)
 
+KRATOS_CREATE_VARIABLE(double, INTERFACE_NORMAL_STIFFNESS)
+KRATOS_CREATE_VARIABLE(double, INTERFACE_SHEAR_STIFFNESS)
+
 } // namespace Kratos

applications/GeoMechanicsApplication/geo_mechanics_application_variables.h

@@ -259,6 +259,10 @@ KRATOS_DEFINE_APPLICATION_VARIABLE(GEO_MECHANICS_APPLICATION, double, STATE_VARI
 
 KRATOS_DEFINE_APPLICATION_VARIABLE(GEO_MECHANICS_APPLICATION, Vector, STRAINS_OF_PIECEWISE_LINEAR_LAW)
 KRATOS_DEFINE_APPLICATION_VARIABLE(GEO_MECHANICS_APPLICATION, Vector, STRESSES_OF_PIECEWISE_LINEAR_LAW)
+
+KRATOS_DEFINE_APPLICATION_VARIABLE(GEO_MECHANICS_APPLICATION, double, INTERFACE_NORMAL_STIFFNESS)
+KRATOS_DEFINE_APPLICATION_VARIABLE(GEO_MECHANICS_APPLICATION, double, INTERFACE_SHEAR_STIFFNESS)
+
 } // namespace Kratos
 
 #endif /* KRATOS_GEO_MECHANICS_APPLICATION_VARIABLES_H_INCLUDED */