Merge branch 'THM_M_Freezing' into 'master'

Extend THM to TM-freezing See merge request ogs/ogs!4642
ufz · Jul 19, 2023 · 56ca45e · 56ca45e
2 parents 4bfda59 + 52d55eb
commit 56ca45e
Show file tree

Hide file tree

Showing 35 changed files with 1,529 additions and 1 deletion.
diff --git a/Documentation/.vale/Vocab/ogs/ignore.txt b/Documentation/.vale/Vocab/ogs/ignore.txt
@@ -51,5 +51,6 @@ subproblems
 thermohydromechanics
 thermomechanical
 thermomechanics
+undeformed
 voxel
 wellbore
diff --git a/Documentation/ProjectFile/material/solid/t_ice_constitutive_relation.md b/Documentation/ProjectFile/material/solid/t_ice_constitutive_relation.md
@@ -0,0 +1 @@
+The constitutive relation for ice used in some processes.
diff --git a/MaterialLib/MPL/PropertyType.h b/MaterialLib/MPL/PropertyType.h
@@ -66,6 +66,7 @@ enum PropertyType : int
     molecular_diffusion,
     name,
     permeability,
+    phase_change_expansivity,
     phase_velocity,
     /// ion diffusivity in the porous medium with account of the effect of
     /// tortuosity and connectivity.
@@ -149,6 +150,7 @@ static const std::array<std::string, PropertyType::number_of_properties>
                              "molecular_diffusion",
                              "name",
                              "permeability",
+                             "phase_change_expansivity",
                              "phase_velocity",
                              "pore_diffusion",
                              "poissons_ratio",

diff --git a/MaterialLib/SolidModels/CreateConstitutiveRelation.cpp b/MaterialLib/SolidModels/CreateConstitutiveRelation.cpp
@@ -82,6 +82,20 @@ createConstitutiveRelation(
               type);
 }
 
+template std::unique_ptr<MaterialLib::Solids::MechanicsBase<2>>
+createConstitutiveRelation(
+    std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
+    std::optional<ParameterLib::CoordinateSystem> const&
+        local_coordinate_system,
+    BaseLib::ConfigTree const& config);
+
+template std::unique_ptr<MaterialLib::Solids::MechanicsBase<3>>
+createConstitutiveRelation(
+    std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
+    std::optional<ParameterLib::CoordinateSystem> const&
+        local_coordinate_system,
+    BaseLib::ConfigTree const& config);
+
 template <int DisplacementDim>
 std::map<int,
          std::unique_ptr<MaterialLib::Solids::MechanicsBase<DisplacementDim>>>

diff --git a/MaterialLib/SolidModels/CreateConstitutiveRelation.h b/MaterialLib/SolidModels/CreateConstitutiveRelation.h
@@ -34,6 +34,28 @@ namespace Solids
 template <int DisplacementDim>
 struct MechanicsBase;
 
+template <int DisplacementDim>
+std::unique_ptr<MaterialLib::Solids::MechanicsBase<DisplacementDim>>
+createConstitutiveRelation(
+    std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
+    std::optional<ParameterLib::CoordinateSystem> const&
+        local_coordinate_system,
+    BaseLib::ConfigTree const& config);
+
+extern template std::unique_ptr<MaterialLib::Solids::MechanicsBase<2>>
+createConstitutiveRelation(
+    std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
+    std::optional<ParameterLib::CoordinateSystem> const&
+        local_coordinate_system,
+    BaseLib::ConfigTree const& config);
+
+extern template std::unique_ptr<MaterialLib::Solids::MechanicsBase<3>>
+createConstitutiveRelation(
+    std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
+    std::optional<ParameterLib::CoordinateSystem> const&
+        local_coordinate_system,
+    BaseLib::ConfigTree const& config);
+
 template <int DisplacementDim>
 std::map<int,
          std::unique_ptr<MaterialLib::Solids::MechanicsBase<DisplacementDim>>>

diff --git a/MaterialLib/SolidModels/CreateConstitutiveRelationIce.cpp b/MaterialLib/SolidModels/CreateConstitutiveRelationIce.cpp
@@ -0,0 +1,55 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2023, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ */
+
+#include "CreateConstitutiveRelationIce.h"
+
+#include "BaseLib/ConfigTree.h"
+#include "MaterialLib/SolidModels/CreateConstitutiveRelation.h"
+#include "MechanicsBase.h"
+
+namespace MaterialLib
+{
+namespace Solids
+{
+template <int DisplacementDim>
+std::unique_ptr<MaterialLib::Solids::MechanicsBase<DisplacementDim>>
+createConstitutiveRelationIce(
+    std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
+    std::optional<ParameterLib::CoordinateSystem> const&
+        local_coordinate_system,
+    BaseLib::ConfigTree const& config)
+{
+    auto const ice_constitutive_relation_config =
+        //! \ogs_file_param{material__solid__ice_constitutive_relation}
+        config.getConfigSubtreeOptional("ice_constitutive_relation");
+
+    if (!ice_constitutive_relation_config)
+    {
+        return nullptr;
+    }
+
+    return MaterialLib::Solids::createConstitutiveRelation<DisplacementDim>(
+        parameters, local_coordinate_system, *ice_constitutive_relation_config);
+}
+
+template std::unique_ptr<MaterialLib::Solids::MechanicsBase<2>>
+createConstitutiveRelationIce(
+    std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
+    std::optional<ParameterLib::CoordinateSystem> const&
+        local_coordinate_system,
+    BaseLib::ConfigTree const& config);
+
+template std::unique_ptr<MaterialLib::Solids::MechanicsBase<3>>
+createConstitutiveRelationIce(
+    std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
+    std::optional<ParameterLib::CoordinateSystem> const&
+        local_coordinate_system,
+    BaseLib::ConfigTree const& config);
+}  // namespace Solids
+}  // namespace MaterialLib
diff --git a/MaterialLib/SolidModels/CreateConstitutiveRelationIce.h b/MaterialLib/SolidModels/CreateConstitutiveRelationIce.h
@@ -0,0 +1,56 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2023, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ */
+
+#pragma once
+
+#include <memory>
+#include <optional>
+#include <vector>
+
+namespace BaseLib
+{
+class ConfigTree;
+}
+
+namespace ParameterLib
+{
+struct ParameterBase;
+struct CoordinateSystem;
+}  // namespace ParameterLib
+
+namespace MaterialLib
+{
+namespace Solids
+{
+template <int DisplacementDim>
+struct MechanicsBase;
+
+template <int DisplacementDim>
+std::unique_ptr<MaterialLib::Solids::MechanicsBase<DisplacementDim>>
+createConstitutiveRelationIce(
+    std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
+    std::optional<ParameterLib::CoordinateSystem> const&
+        local_coordinate_system,
+    BaseLib::ConfigTree const& config);
+
+extern template std::unique_ptr<MaterialLib::Solids::MechanicsBase<2>>
+createConstitutiveRelationIce(
+    std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
+    std::optional<ParameterLib::CoordinateSystem> const&
+        local_coordinate_system,
+    BaseLib::ConfigTree const& config);
+
+extern template std::unique_ptr<MaterialLib::Solids::MechanicsBase<3>>
+createConstitutiveRelationIce(
+    std::vector<std::unique_ptr<ParameterLib::ParameterBase>> const& parameters,
+    std::optional<ParameterLib::CoordinateSystem> const&
+        local_coordinate_system,
+    BaseLib::ConfigTree const& config);
+}  // namespace Solids
+}  // namespace MaterialLib
diff --git a/ProcessLib/ThermoHydroMechanics/CreateThermoHydroMechanicsProcess.cpp b/ProcessLib/ThermoHydroMechanics/CreateThermoHydroMechanicsProcess.cpp
@@ -15,6 +15,7 @@
 #include "MaterialLib/MPL/CreateMaterialSpatialDistributionMap.h"
 #include "MaterialLib/MPL/MaterialSpatialDistributionMap.h"
 #include "MaterialLib/SolidModels/CreateConstitutiveRelation.h"
+#include "MaterialLib/SolidModels/CreateConstitutiveRelationIce.h"
 #include "MaterialLib/SolidModels/MechanicsBase.h"
 #include "NumLib/NumericalStability/CreateNumericalStabilization.h"
 #include "ParameterLib/Utils.h"
@@ -145,6 +146,10 @@ std::unique_ptr<Process> createThermoHydroMechanicsProcess(
         MaterialLib::Solids::createConstitutiveRelations<DisplacementDim>(
             parameters, local_coordinate_system, config);
 
+    auto ice_constitutive_relation =
+        MaterialLib::Solids::createConstitutiveRelationIce<DisplacementDim>(
+            parameters, local_coordinate_system, config);
+
     // Specific body force
     Eigen::Matrix<double, DisplacementDim, 1> specific_body_force;
     {
@@ -181,6 +186,7 @@ std::unique_ptr<Process> createThermoHydroMechanicsProcess(
         materialIDs(mesh),
         std::move(media_map),
         std::move(solid_constitutive_relations),
+        std::move(ice_constitutive_relation),
         initial_stress,
         specific_body_force,
         std::move(stabilizer)};

diff --git a/ProcessLib/ThermoHydroMechanics/IntegrationPointData.h b/ProcessLib/ThermoHydroMechanics/IntegrationPointData.h
@@ -37,19 +37,29 @@ struct IntegrationPointData final
         static const int kelvin_vector_size =
             MathLib::KelvinVector::kelvin_vector_dimensions(DisplacementDim);
         sigma_eff.setZero(kelvin_vector_size);
+        sigma_eff_ice.setZero(kelvin_vector_size);
         eps.setZero(kelvin_vector_size);
         eps_prev.setZero(kelvin_vector_size);
+        eps0_prev.setZero(kelvin_vector_size);
+        eps0_prev2.setZero(kelvin_vector_size);
         eps_m.setZero(kelvin_vector_size);
         eps_m_prev.resize(kelvin_vector_size);
+        eps_m_ice.setZero(kelvin_vector_size);
+        eps_m_ice_prev.resize(kelvin_vector_size);
 
         // Previous time step values are not initialized and are set later.
         sigma_eff_prev.resize(kelvin_vector_size);
+        sigma_eff_ice_prev.resize(kelvin_vector_size);
     }
 
     typename BMatricesType::KelvinVectorType sigma_eff, sigma_eff_prev;
-    typename BMatricesType::KelvinVectorType eps, eps_prev;
+    typename BMatricesType::KelvinVectorType eps, eps_prev, eps0_prev,
+        eps0_prev2;
     typename BMatricesType::KelvinVectorType eps_m, eps_m_prev;
 
+    typename BMatricesType::KelvinVectorType sigma_eff_ice, sigma_eff_ice_prev;
+    typename BMatricesType::KelvinVectorType eps_m_ice, eps_m_ice_prev;
+
     typename ShapeMatrixTypeDisplacement::NodalRowVectorType N_u;
     typename ShapeMatrixTypeDisplacement::GlobalDimNodalMatrixType dNdx_u;
 
@@ -61,13 +71,20 @@ struct IntegrationPointData final
     std::unique_ptr<typename MaterialLib::Solids::MechanicsBase<
         DisplacementDim>::MaterialStateVariables>
         material_state_variables;
+
+    double phi_fr = 0;
+    double phi_fr_prev;
     double integration_weight;
 
     void pushBackState()
     {
+        phi_fr_prev = phi_fr;
+        eps0_prev2 = eps0_prev;
         eps_prev = eps;
         eps_m_prev = eps_m;
         sigma_eff_prev = sigma_eff;
+        sigma_eff_ice_prev = sigma_eff_ice;
+        eps_m_ice_prev = eps_m_ice;
         material_state_variables->pushBackState();
     }
 
@@ -141,6 +158,41 @@ struct IntegrationPointData final
         return C;
     }
 
+    typename BMatricesType::KelvinMatrixType updateConstitutiveRelationIce(
+        MaterialLib::Solids::MechanicsBase<DisplacementDim> const&
+            ice_constitutive_relation,
+        MaterialPropertyLib::VariableArray const& variable_array,
+        double const t,
+        ParameterLib::SpatialPosition const& x_position,
+        double const dt,
+        double const temperature_prev)
+    {
+        MaterialPropertyLib::VariableArray variable_array_prev;
+
+        variable_array_prev.stress
+            .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
+                sigma_eff_ice_prev);
+        variable_array_prev.mechanical_strain
+            .emplace<MathLib::KelvinVector::KelvinVectorType<DisplacementDim>>(
+                eps_m_ice_prev);
+        variable_array_prev.temperature = temperature_prev;
+
+        // Extend for non-linear ice materials if necessary.
+        auto const null_state =
+            ice_constitutive_relation.createMaterialStateVariables();
+        auto&& solution = ice_constitutive_relation.integrateStress(
+            variable_array_prev, variable_array, t, x_position, dt,
+            *null_state);
+
+        if (!solution)
+            OGS_FATAL("Computation of local constitutive relation failed.");
+
+        MathLib::KelvinVector::KelvinMatrixType<DisplacementDim> C_IR;
+        std::tie(sigma_eff_ice, material_state_variables, C_IR) =
+            std::move(*solution);
+
+        return C_IR;
+    }
     EIGEN_MAKE_ALIGNED_OPERATOR_NEW;
 };
 
@@ -180,6 +232,9 @@ struct ConstitutiveRelationsValues
 
     // Freezing related values.
     double J_TT_fr;
+    MathLib::KelvinVector::KelvinMatrixType<DisplacementDim> J_uu_fr;
+    MathLib::KelvinVector::KelvinVectorType<DisplacementDim> J_uT_fr;
+    MathLib::KelvinVector::KelvinVectorType<DisplacementDim> r_u_fr;
 };
 
 }  // namespace ThermoHydroMechanics

diff --git a/ProcessLib/ThermoHydroMechanics/LocalAssemblerInterface.h b/ProcessLib/ThermoHydroMechanics/LocalAssemblerInterface.h
@@ -32,6 +32,12 @@ struct LocalAssemblerInterface : public ProcessLib::LocalAssemblerInterface,
         std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
         std::vector<double>& cache) const = 0;
 
+    virtual std::vector<double> const& getIntPtSigmaIce(
+        const double t,
+        std::vector<GlobalVector*> const& x,
+        std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
+        std::vector<double>& cache) const = 0;
+
     virtual std::vector<double> getEpsilon() const = 0;
 
     virtual std::vector<double> const& getIntPtEpsilon(
@@ -40,6 +46,12 @@ struct LocalAssemblerInterface : public ProcessLib::LocalAssemblerInterface,
         std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
         std::vector<double>& cache) const = 0;
 
+    virtual std::vector<double> const& getIntPtIceVolume(
+        const double t,
+        std::vector<GlobalVector*> const& x,
+        std::vector<NumLib::LocalToGlobalIndexMap const*> const& dof_table,
+        std::vector<double>& cache) const = 0;
+
     virtual std::vector<double> const& getIntPtDarcyVelocity(
         const double t,
         std::vector<GlobalVector*> const& x,

diff --git a/ProcessLib/ThermoHydroMechanics/Tests.cmake b/ProcessLib/ThermoHydroMechanics/Tests.cmake
@@ -4,6 +4,7 @@ if (NOT OGS_USE_MPI)
     OgsTest(PROJECTFILE ThermoHydroMechanics/A2/A2.prj RUNTIME 23)
     OgsTest(PROJECTFILE ThermoHydroMechanics/A2/A2_heating.prj RUNTIME 23)
     OgsTest(PROJECTFILE ThermoHydroMechanics/1D_freezing_column_Stefan/Stefan_problem.prj RUNTIME 15)
+    OgsTest(PROJECTFILE ThermoHydroMechanics/ColumnDeformationFreezing/TM.prj RUNTIME 13)
     OgsTest(PROJECTFILE ThermoHydroMechanics/HeatingHomogeneousDomain/hex_THM.prj RUNTIME 30)
 endif()