[GeoMechanicsApplication] Add micro climate condition

applications/GeoMechanicsApplication/custom_conditions/T_microclimate_flux_condition.h

@@ -0,0 +1,150 @@
+// KRATOS___
+//     //   ) )
+//    //         ___      ___
+//   //  ____  //___) ) //   ) )
+//  //    / / //       //   / /
+// ((____/ / ((____   ((___/ /  MECHANICS
+//
+//  License:         geo_mechanics_application/license.txt
+//
+//
+//  Main authors:    Mohamed Nabi
+//
+
+#pragma once
+
+// Project includes
+#include "includes/serializer.h"
+
+// Application includes
+#include "custom_conditions/T_condition.h"
+#include "custom_utilities/condition_utilities.hpp"
+#include "custom_utilities/element_utilities.hpp"
+#include "geo_mechanics_application_variables.h"
+
+namespace Kratos
+{
+
+struct WaterFluxes
+{
+    double precipitation = 0.0;
+    double evaporation = 0.0;
+};
+
+template <unsigned int TDim, unsigned int TNumNodes>
+class KRATOS_API(GEO_MECHANICS_APPLICATION) TMicroClimateFluxCondition
+    : public GeoTCondition<TDim, TNumNodes>
+{
+public:
+    KRATOS_CLASS_INTRUSIVE_POINTER_DEFINITION(TMicroClimateFluxCondition);
+
+    using IndexType = std::size_t;
+    using GeometryType = Geometry<Node>;
+    using NodesArrayType = GeometryType::PointsArrayType;
+
+    TMicroClimateFluxCondition() : GeoTCondition<TDim, TNumNodes>() {}
+
+    TMicroClimateFluxCondition(IndexType NewId, GeometryType::Pointer pGeometry)
+        : GeoTCondition<TDim, TNumNodes>(NewId, pGeometry)
+    {
+    }
+
+    TMicroClimateFluxCondition(IndexType NewId,
+                               GeometryType::Pointer pGeometry,
+                               Properties::Pointer pProperties)
+        : GeoTCondition<TDim, TNumNodes>(NewId, pGeometry, pProperties)
+    {
+    }
+
+    Condition::Pointer Create(IndexType NewId,
+                              const NodesArrayType& rNodes,
+                              Properties::Pointer pProperties) const override;
+
+    void Initialize(const ProcessInfo& rCurrentProcessInfo) override;
+
+    void InitializeSolutionStep(const ProcessInfo& rCurrentProcessInfo) override;
+
+    void CalculateLocalSystem(Matrix& rLeftHandSideMatrix,
+                              Vector& rRightHandSideVector,
+                              const ProcessInfo& rCurrentProcessInfo) override;
+
+private:
+    void InitializeProperties();
+
+    void CalculateAndAddLHS(Matrix& rLeftHandSideMatrix,
+                            const array_1d<double, TNumNodes>& rN,
+                            double IntegrationCoefficient,
+                            const array_1d<double, TNumNodes>& rLeftHandSideFluxes);
+
+    void CalculateAndAddRHS(Vector& rRightHandSideVector,
+                            const array_1d<double, TNumNodes>& rN,
+                            double IntegrationCoefficient,
+                            const Vector& rNodalTemperatures,
+                            const array_1d<double, TNumNodes>& rLeftHandSideFluxes,
+                            const array_1d<double, TNumNodes>& rRightHandSideFluxes);
+
+    array_1d<double, TNumNodes> CalculateLeftHandSideFluxes() const;
+    array_1d<double, TNumNodes> CalculateRightHandSideFluxes(double TimeStepSize,
+                                                             double PreviousStorage,
+                                                             double PreviousRadiation) const;
+    double CalculateRightHandSideFlux(double NetRadiation,
+                                      double SurfaceHeatStorage,
+                                      double ActualEvaporation) const;
+
+    double CalculateCurrentWaterStorage(double TimeStepSize,
+                                        double PreviousStorage,
+                                        double PreviousRadiation) const;
+
+    double CalculateCurrentNetRadiation() const;
+    double CalculateNetRadiation(unsigned int NodeIndex) const;
+
+    double CalculateSurfaceHeatStorage(double TimeStepSize,
+                                       double PreviousRadiation,
+                                       double NetRadiation) const;
+
+    WaterFluxes CalculateWaterFluxes(unsigned int NodeIndex,
+                                     double TimeStepSize,
+                                     double PreviousStorage,
+                                     double NetRadiation,
+                                     double SurfaceHeatStorage) const;
+
+    double CalculatePotentialEvaporation(unsigned int NodeIndex,
+                                         double NetRadiation,
+                                         double SurfaceHeatStorage) const;
+
+    void CalculateRoughness(const ProcessInfo& rCurrentProcessInfo);
+
+    double CalculateSurfaceRoughnessFactor(double CurrentAirTemperature,
+                                           double PreviousRoughnessTemperature,
+                                           double RichardsonBulkModulus,
+                                           double FrictionDragCoefficient) const;
+
+    // Serialization
+    friend class Serializer;
+
+    void save(Serializer& rSerializer) const override
+    {
+        KRATOS_SERIALIZE_SAVE_BASE_CLASS(rSerializer, Condition)
+    }
+
+    void load(Serializer& rSerializer) override
+    {
+        KRATOS_SERIALIZE_LOAD_BASE_CLASS(rSerializer, Condition)
+    }
+
+    bool mIsInitialized = false;
+    double mAlbedoCoefficient = 0.0;
+    double mFirstCoverStorageCoefficient = 0.0;
+    double mSecondCoverStorageCoefficient = 0.0;
+    double mThirdCoverStorageCoefficient = 0.0;
+    double mBuildEnvironmentRadiation = 0.0;
+    double mMinimalStorage = 0.0;
+    double mMaximalStorage = 0.0;
+    double mRoughnessTemperature = 0.0;
+    double mNetRadiation = 0.0;
+    double mWaterStorage = 0.0;
+};
+
+// class TMicroClimateFluxCondition.
+
+} // namespace Kratos.

applications/GeoMechanicsApplication/custom_conditions/T_normal_flux_condition.cpp

@@ -14,6 +14,7 @@
 
 // Application includes
 #include "custom_conditions/T_normal_flux_condition.h"
+#include "custom_utilities/condition_utilities.hpp"
 #include "utilities/math_utils.h"
 
 namespace Kratos {
@@ -73,8 +74,10 @@ void GeoTNormalFluxCondition<TDim, TNumNodes>::CalculateRHS(Vector& rRightHandSi
         // Interpolation of nodal normal flux to integration point normal flux.
         auto normal_flux_on_integration_point = MathUtils<>::Dot(N, normal_flux_vector);
 
-        auto weighting_integration_coefficient = CalculateIntegrationCoefficient(j_container[integration_point],
-                                                                                        r_integration_points[integration_point].Weight());
+        auto weighting_integration_coefficient =
+            ConditionUtilities::CalculateIntegrationCoefficient<TDim, TNumNodes>(
+                j_container[integration_point],
+                r_integration_points[integration_point].Weight());
 
         // Contributions to the right hand side
         auto normal_flux_on_DOF = normal_flux_on_integration_point * N * weighting_integration_coefficient;
@@ -83,22 +86,6 @@ void GeoTNormalFluxCondition<TDim, TNumNodes>::CalculateRHS(Vector& rRightHandSi
     }
 }
 
-template <unsigned int TDim, unsigned int TNumNodes>
-double GeoTNormalFluxCondition<TDim, TNumNodes>::CalculateIntegrationCoefficient(const Matrix& rJacobian,
-                                                                                 double Weight)
-{
-    Vector normal_vector = ZeroVector(TDim);
-
-    if constexpr (TDim == 2) {
-        normal_vector = column(rJacobian, 0);
-    }
-    else if constexpr (TDim == 3) {
-        MathUtils<double>::CrossProduct(normal_vector, column(rJacobian, 0),
-                                        column(rJacobian, 1));
-    }
-    return Weight * MathUtils<double>::Norm(normal_vector);
-}
-
 template class GeoTNormalFluxCondition<2, 2>;
 template class GeoTNormalFluxCondition<2, 3>;
 template class GeoTNormalFluxCondition<2, 4>;

applications/GeoMechanicsApplication/custom_conditions/T_normal_flux_condition.h

@@ -53,9 +53,6 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) GeoTNormalFluxCondition
     void CalculateRHS(Vector& rRightHandSideVector,
                       const ProcessInfo& CurrentProcessInfo) override;
 
-    double CalculateIntegrationCoefficient(const Matrix& rJacobian,
-                                           double Weight);
-
 private:
     friend class Serializer;
 

applications/GeoMechanicsApplication/custom_conditions/micro_climate_constants.h

@@ -0,0 +1,34 @@
+// KRATOS___
+//     //   ) )
+//    //         ___      ___
+//   //  ____  //___) ) //   ) )
+//  //    / / //       //   / /
+// ((____/ / ((____   ((___/ /  MECHANICS
+//
+//  License:         geo_mechanics_application/license.txt
+//
+//  Main authors:    Richard Faasse
+//
+#pragma once
+
+namespace Kratos::MicroClimateConstants
+{
+
+constexpr static auto AirDensity = 1.18;
+constexpr static auto AirHeatCapacity = 1004.67;
+constexpr static auto RoughnessLayerResistance = 30.0;
+constexpr static auto LatentEvaporationHeat = 2.45e6;
+constexpr static auto WaterDensity = 1e3;
+constexpr static auto PsychometricConstant = 0.63;
+constexpr static auto SurfaceResistance = 30.0;
+
+constexpr static auto RoughnessLayerHeight = 10.0;
+constexpr static auto VonNeumannCoefficient = 0.4;
+constexpr static auto MeasurementHeight = 10.0;
+constexpr static auto RoughnessHeight = 1.0;
+constexpr static auto GravitationalAcceleration = 9.81;
+
+constexpr static auto EffectiveEmissivity = 0.95;
+constexpr static auto BoltzmannCoefficient = 5.67e-8;
+
+} // namespace Kratos::MicroClimateConstants

applications/GeoMechanicsApplication/custom_python/geo_mechanics_python_application.cpp

@@ -62,6 +62,12 @@ PYBIND11_MODULE(KratosGeoMechanicsApplication,m)
 
     KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, DT_TEMPERATURE)
     KRATOS_REGISTER_IN_PYTHON_VARIABLE(m, NORMAL_HEAT_FLUX)
+	KRATOS_REGISTER_IN_PYTHON_VARIABLE( m, AIR_TEMPERATURE )
+    KRATOS_REGISTER_IN_PYTHON_VARIABLE( m, SOLAR_RADIATION )
+    KRATOS_REGISTER_IN_PYTHON_VARIABLE( m, AIR_HUMIDITY )
+    KRATOS_REGISTER_IN_PYTHON_VARIABLE( m, PRECIPITATION )
+    KRATOS_REGISTER_IN_PYTHON_VARIABLE( m, WIND_SPEED )
+
 
     KRATOS_REGISTER_IN_PYTHON_3D_VARIABLE_WITH_COMPONENTS( m, FLUID_FLUX_VECTOR )
     KRATOS_REGISTER_IN_PYTHON_3D_VARIABLE_WITH_COMPONENTS( m, LOCAL_FLUID_FLUX_VECTOR )

applications/GeoMechanicsApplication/custom_utilities/condition_utilities.hpp

@@ -198,6 +198,23 @@ class ConditionUtilities
         }
     }
 
+    template <unsigned int TDim, unsigned int TNumNodes>
+    static double CalculateIntegrationCoefficient(const Matrix& rJacobian, double Weight)
+    {
+        auto normal_vector = Vector{TDim, 0.0};
+
+        if constexpr (TDim == 2)
+        {
+            normal_vector = column(rJacobian, 0);
+        }
+        else if constexpr (TDim == 3)
+        {
+            MathUtils<>::CrossProduct(normal_vector, column(rJacobian, 0),
+                                      column(rJacobian, 1));
+        }
+        return Weight * MathUtils<>::Norm(normal_vector);
+    }
+
 }; /* Class ConditionUtilities*/
 } /* namespace Kratos.*/
 

applications/GeoMechanicsApplication/custom_utilities/variables_utilities.hpp

@@ -0,0 +1,38 @@
+// KRATOS___
+//     //   ) )
+//    //         ___      ___
+//   //  ____  //___) ) //   ) )
+//  //    / / //       //   / /
+// ((____/ / ((____   ((___/ /  MECHANICS
+//
+//  License:         geo_mechanics_application/license.txt
+//
+//
+//  Main authors:    Anne van de Graaf
+//
+
+#pragma once
+
+#include <algorithm>
+
+#include "containers/variable.h"
+
+namespace Kratos
+{
+
+class VariablesUtilities
+{
+public:
+    template <typename GeometryType, typename OutputIt>
+    static OutputIt GetNodalValues(const GeometryType& rGeometry,
+                                   const Variable<double>& rNodalVariable,
+                                   OutputIt FirstOut)
+    {
+        return std::transform(rGeometry.begin(), rGeometry.end(), FirstOut,
+                              [&rNodalVariable](const auto& node) {
+                                  return node.FastGetSolutionStepValue(rNodalVariable);
+                              });
+    }
+};
+
+}

applications/GeoMechanicsApplication/geo_mechanics_application.cpp

@@ -287,6 +287,16 @@ void KratosGeoMechanicsApplication::Register() {
     KRATOS_REGISTER_CONDITION("GeoTNormalFluxCondition3D8N", mGeoTNormalFluxCondition3D8N)
     KRATOS_REGISTER_CONDITION("GeoTNormalFluxCondition3D9N", mGeoTNormalFluxCondition3D9N)
 
+    KRATOS_REGISTER_CONDITION("GeoTMicroClimateFluxCondition2D2N", mGeoTMicroClimateFluxCondition2D2N)
+    KRATOS_REGISTER_CONDITION("GeoTMicroClimateFluxCondition2D3N", mGeoTMicroClimateFluxCondition2D3N)
+    KRATOS_REGISTER_CONDITION("GeoTMicroClimateFluxCondition2D4N", mGeoTMicroClimateFluxCondition2D4N)
+    KRATOS_REGISTER_CONDITION("GeoTMicroClimateFluxCondition2D5N", mGeoTMicroClimateFluxCondition2D5N)
+    KRATOS_REGISTER_CONDITION("GeoTMicroClimateFluxCondition3D3N", mGeoTMicroClimateFluxCondition3D3N)
+    KRATOS_REGISTER_CONDITION("GeoTMicroClimateFluxCondition3D4N", mGeoTMicroClimateFluxCondition3D4N)
+    KRATOS_REGISTER_CONDITION("GeoTMicroClimateFluxCondition3D6N", mGeoTMicroClimateFluxCondition3D6N)
+    KRATOS_REGISTER_CONDITION("GeoTMicroClimateFluxCondition3D8N", mGeoTMicroClimateFluxCondition3D8N)
+    KRATOS_REGISTER_CONDITION("GeoTMicroClimateFluxCondition3D9N", mGeoTMicroClimateFluxCondition3D9N)
+
     //Register Constitutive Laws
     KRATOS_REGISTER_CONSTITUTIVE_LAW("BilinearCohesive3DLaw", mBilinearCohesive3DLaw)
     KRATOS_REGISTER_CONSTITUTIVE_LAW("BilinearCohesive2DLaw", mBilinearCohesive2DLaw)

applications/GeoMechanicsApplication/geo_mechanics_application.h

@@ -29,6 +29,7 @@
 
 // conditions
 #include "custom_conditions/Pw_normal_flux_condition.hpp"
+#include "custom_conditions/T_microclimate_flux_condition.h"
 #include "custom_conditions/T_normal_flux_condition.h"
 #include "custom_conditions/U_Pw_face_load_condition.hpp"
 #include "custom_conditions/U_Pw_face_load_interface_condition.hpp"
@@ -535,6 +536,17 @@ class KRATOS_API(GEO_MECHANICS_APPLICATION) KratosGeoMechanicsApplication : publ
     const GeoTNormalFluxCondition<3, 8> mGeoTNormalFluxCondition3D8N{ 0, Kratos::make_shared< Quadrilateral3D8 <NodeType> >(Condition::GeometryType::PointsArrayType(8)) };
     const GeoTNormalFluxCondition<3, 9> mGeoTNormalFluxCondition3D9N{ 0, Kratos::make_shared< Quadrilateral3D9 <NodeType> >(Condition::GeometryType::PointsArrayType(9)) };
 
+    const TMicroClimateFluxCondition<2, 2> mGeoTMicroClimateFluxCondition2D2N{0, Kratos::make_shared<Line2D2         <NodeType>>(Condition::GeometryType::PointsArrayType(2))};
+    const TMicroClimateFluxCondition<2, 3> mGeoTMicroClimateFluxCondition2D3N{0, Kratos::make_shared<Line2D3         <NodeType>>(Condition::GeometryType::PointsArrayType(3))};
+    const TMicroClimateFluxCondition<2, 4> mGeoTMicroClimateFluxCondition2D4N{0, Kratos::make_shared<Line2D4         <NodeType>>(Condition::GeometryType::PointsArrayType(4))};
+    const TMicroClimateFluxCondition<2, 5> mGeoTMicroClimateFluxCondition2D5N{0, Kratos::make_shared<Line2D5         <NodeType>>(Condition::GeometryType::PointsArrayType(5))};
+    const TMicroClimateFluxCondition<3, 3> mGeoTMicroClimateFluxCondition3D3N{0, Kratos::make_shared<Triangle3D3     <NodeType>>(Condition::GeometryType::PointsArrayType(3))};
+    const TMicroClimateFluxCondition<3, 6> mGeoTMicroClimateFluxCondition3D6N{0, Kratos::make_shared<Triangle3D6     <NodeType>>(Condition::GeometryType::PointsArrayType(6))};
+    const TMicroClimateFluxCondition<3, 4> mGeoTMicroClimateFluxCondition3D4N{0, Kratos::make_shared<Quadrilateral3D4<NodeType>>(Condition::GeometryType::PointsArrayType(4))};
+    const TMicroClimateFluxCondition<3, 8> mGeoTMicroClimateFluxCondition3D8N{0, Kratos::make_shared<Quadrilateral3D8<NodeType>>(Condition::GeometryType::PointsArrayType(8))};
+    const TMicroClimateFluxCondition<3, 9> mGeoTMicroClimateFluxCondition3D9N{0, Kratos::make_shared<Quadrilateral3D9<NodeType>>(Condition::GeometryType::PointsArrayType(9))};
+
+
     // constitutive models
     const BilinearCohesive3DLaw             mBilinearCohesive3DLaw;
     const BilinearCohesive2DLaw             mBilinearCohesive2DLaw;

applications/GeoMechanicsApplication/python_scripts/geomechanics_solver.py

@@ -361,6 +361,11 @@ def _add_temperature_variables(self):
         self.main_model_part.AddNodalSolutionStepVariable(GeoMechanicsApplication.DT_TEMPERATURE)
         # Add variables for the heat conditions
         self.main_model_part.AddNodalSolutionStepVariable(GeoMechanicsApplication.NORMAL_HEAT_FLUX)
+        self.main_model_part.AddNodalSolutionStepVariable(GeoMechanicsApplication.AIR_TEMPERATURE)
+        self.main_model_part.AddNodalSolutionStepVariable(GeoMechanicsApplication.SOLAR_RADIATION)
+        self.main_model_part.AddNodalSolutionStepVariable(GeoMechanicsApplication.AIR_HUMIDITY)
+        self.main_model_part.AddNodalSolutionStepVariable(GeoMechanicsApplication.PRECIPITATION)
+        self.main_model_part.AddNodalSolutionStepVariable(GeoMechanicsApplication.WIND_SPEED)
 
     def _add_smoothing_variables(self):
         self.main_model_part.AddNodalSolutionStepVariable(KratosMultiphysics.NODAL_AREA)