Merge branch 'THM_M_Freezing' into 'master'

Extend THM  to TM-freezing

See merge request ogs/ogs!4642

Documentation/.vale/Vocab/ogs/ignore.txt

@@ -51,5 +51,6 @@ subproblems
 thermohydromechanics
 thermomechanical
 thermomechanics
+undeformed
 voxel
 wellbore

Documentation/ProjectFile/material/solid/t_ice_constitutive_relation.md

@@ -0,0 +1 @@
+The constitutive relation for ice used in some processes.

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",

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>>>

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>>>

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

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

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)};

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

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,

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()