Merge branch 'THMF_CorrectHistoryVariableUpdate' into 'master'

THM/F Correct history variable update See merge request ogs/ogs!4723
ufz · Aug 31, 2023 · 4a5ce09 · 4a5ce09
2 parents 073d0b9 + 0ab7091
commit 4a5ce09
Show file tree

Hide file tree

Showing 33 changed files with 1,273 additions and 325 deletions.
diff --git a/ProcessLib/ThermoHydroMechanics/IntegrationPointData.h b/ProcessLib/ThermoHydroMechanics/IntegrationPointData.h
@@ -40,8 +40,8 @@ struct IntegrationPointData final
  sigma_eff_ice.setZero(kelvin_vector_size);
  eps.setZero(kelvin_vector_size);
  eps_prev.setZero(kelvin_vector_size);
+ eps0.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);
@@ -53,8 +53,7 @@ struct IntegrationPointData final
  }
 
  typename BMatricesType::KelvinVectorType sigma_eff, sigma_eff_prev;
- typename BMatricesType::KelvinVectorType eps, eps_prev, eps0_prev,
- eps0_prev2;
+ typename BMatricesType::KelvinVectorType eps, eps_prev, eps0, eps0_prev;
  typename BMatricesType::KelvinVectorType eps_m, eps_m_prev;
 
  typename BMatricesType::KelvinVectorType sigma_eff_ice, sigma_eff_ice_prev;
@@ -76,10 +75,12 @@ struct IntegrationPointData final
  double phi_fr_prev;
  double integration_weight;
 
+ double porosity;
+
  void pushBackState()
  {
  phi_fr_prev = phi_fr;
- eps0_prev2 = eps0_prev;
+ eps0_prev = eps0;
  eps_prev = eps;
  eps_m_prev = eps_m;
  sigma_eff_prev = sigma_eff;
@@ -180,6 +181,8 @@ struct IntegrationPointData final
  // Extend for non-linear ice materials if necessary.
  auto const null_state =
  ice_constitutive_relation.createMaterialStateVariables();
+ ice_constitutive_relation.initializeInternalStateVariables(
+ t, x_position, *null_state);
  auto&& solution = ice_constitutive_relation.integrateStress(
  variable_array_prev, variable_array, t, x_position, dt,
  *null_state);
@@ -227,7 +230,6 @@ struct ConstitutiveRelationsValues
  double c_f;
  double effective_volumetric_heat_capacity;
  double fluid_compressibility;
- double porosity;
  double rho;
 
  // Freezing related values.

diff --git a/ProcessLib/ThermoHydroMechanics/Tests.cmake b/ProcessLib/ThermoHydroMechanics/Tests.cmake
@@ -6,6 +6,7 @@ if (NOT OGS_USE_MPI)
  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)
+ OgsTest(PROJECTFILE ThermoHydroMechanics/9percentWaterFreezingExpansion/UnitSquare.prj RUNTIME 1)
 endif()
 
 AddTest(

diff --git a/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsFEM-impl.h b/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsFEM-impl.h
@@ -213,7 +213,7 @@ ConstitutiveRelationsValues<DisplacementDim> ThermoHydroMechanicsLocalAssembler<
  medium->property(MaterialPropertyLib::PropertyType::porosity)
  .template value<double>(vars, x_position, t, dt);
  vars.porosity = porosity;
- crv.porosity = porosity;
+ ip_data.porosity = porosity;
 
  crv.alpha_biot =
  medium->property(MaterialPropertyLib::PropertyType::biot_coefficient)
@@ -413,29 +413,21 @@ ConstitutiveRelationsValues<DisplacementDim> ThermoHydroMechanicsLocalAssembler<
  phase_change_expansivity)
  .value(vars, x_position, t, dt));
 
- // Heaviside(S_I - 0.5) = Heaviside(phi_I / phi - 0.5)
- auto heaviside_I = [porosity](double const phi_fr)
- { return phi_fr >= 0.5 * porosity ? 1. : 0.; };
-
  MathLib::KelvinVector::KelvinVectorType<DisplacementDim> const
- dphase_change_strain =
- phase_change_expansion_coefficient *
- (heaviside_I(phi_fr) - heaviside_I(phi_fr_prev));
+ dphase_change_strain = phase_change_expansion_coefficient *
+ (phi_fr - phi_fr_prev) / porosity;
 
  // eps0 ia a 'history variable' -- a solid matrix strain accrued
  // prior to the onset of ice forming
- auto& eps0_prev = ip_data.eps0_prev;
- auto const& eps0_prev2 = ip_data.eps0_prev2;
- // definition of eps_m_ice
- eps0_prev = eps0_prev2 +
- (1 - heaviside_I(phi_fr_prev)) * (eps_prev - eps0_prev2);
+ auto& eps0 = ip_data.eps0;
+ auto const& eps0_prev = ip_data.eps0_prev;
 
- // ice
+ // definition of eps_m_ice
  auto& eps_m_ice = ip_data.eps_m_ice;
  auto const& eps_m_ice_prev = ip_data.eps_m_ice_prev;
 
  eps_m_ice.noalias() = eps_m_ice_prev + eps - eps_prev -
- (eps0_prev - eps0_prev2) - dthermal_strain_ice -
+ (eps0 - eps0_prev) - dthermal_strain_ice -
  dphase_change_strain;
 
  vars_ice.mechanical_strain
@@ -456,16 +448,12 @@ ConstitutiveRelationsValues<DisplacementDim> ThermoHydroMechanicsLocalAssembler<
  MaterialPropertyLib::Variable::temperature, x_position, t,
  dt);
 
- crv.J_uu_fr = heaviside_I(phi_fr) * porosity * C_IR;
+ crv.J_uu_fr = phi_fr * C_IR;
 
  auto const& sigma_eff_ice = ip_data.sigma_eff_ice;
- crv.r_u_fr = heaviside_I(phi_fr) * porosity * sigma_eff_ice;
- crv.J_uT_fr =
- dphi_fr_dT * sigma_eff_ice + // here, there must be a Dirac-delta
- heaviside_I(phi_fr) * porosity * C_IR *
- (ice_linear_thermal_expansion_coefficient +
- phase_change_expansion_coefficient * dphi_fr_dT /
- porosity); // here, there must be a Dirac-delta
+ crv.r_u_fr = phi_fr * sigma_eff_ice;
+
+ crv.J_uT_fr = phi_fr * C_IR * ice_linear_thermal_expansion_coefficient;
 
  crv.J_TT_fr = ((rho_fr * c_fr - fluid_density * crv.c_f) * dphi_fr_dT +
  l_fr * rho_fr * d2phi_fr_dT2) *
@@ -633,10 +621,11 @@ void ThermoHydroMechanicsLocalAssembler<
  //
  laplace_p.noalias() += dNdx.transpose() * crv.K_over_mu * dNdx * w;
 
- storage_p.noalias() += N.transpose() *
- (crv.porosity * crv.fluid_compressibility +
- (crv.alpha_biot - crv.porosity) * crv.beta_SR) *
- N * w;
+ storage_p.noalias() +=
+ N.transpose() *
+ (_ip_data[ip].porosity * crv.fluid_compressibility +
+ (crv.alpha_biot - _ip_data[ip].porosity) * crv.beta_SR) *
+ N * w;
 
  laplace_T.noalias() +=
  dNdx.transpose() * crv.K_pT_thermal_osmosis * dNdx * w;

diff --git a/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsFEM.h b/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsFEM.h
@@ -174,6 +174,10 @@ class ThermoHydroMechanicsLocalAssembler : public LocalAssemblerInterface
 
  for (unsigned ip = 0; ip < n_integration_points; ip++)
  {
+ _ip_data[ip].eps0 =
+ _ip_data[ip].eps0_prev +
+ (1 - _ip_data[ip].phi_fr_prev / _ip_data[ip].porosity) *
+ (_ip_data[ip].eps_prev - _ip_data[ip].eps0_prev);
  _ip_data[ip].pushBackState();
  }
  }

diff --git a/Tests/Data/ThermoHydroMechanics/9percentWaterFreezingExpansion/Results_ts_0_t_0_sec.vtu b/Tests/Data/ThermoHydroMechanics/9percentWaterFreezingExpansion/Results_ts_0_t_0_sec.vtu
@@ -0,0 +1,46 @@
+<?xml version="1.0"?>
+<VTKFile type="UnstructuredGrid" version="1.0" byte_order="LittleEndian" header_type="UInt64" compressor="vtkZLibDataCompressor">
+ <UnstructuredGrid>
+ <FieldData>
+ <DataArray type="Int8" Name="IntegrationPointMetaData" NumberOfTuples="166" format="appended" RangeMin="34" RangeMax="125" offset="0" />
+ <DataArray type="Int8" Name="OGS_VERSION" NumberOfTuples="26" format="appended" RangeMin="45" RangeMax="121" offset="176" />
+ <DataArray type="Float64" Name="epsilon_ip" NumberOfComponents="4" NumberOfTuples="4" format="appended" RangeMin="0" RangeMax="0" offset="268" />
+ <DataArray type="Float64" Name="sigma_ip" NumberOfComponents="4" NumberOfTuples="4" format="appended" RangeMin="0" RangeMax="0" offset="328" />
+ </FieldData>
+ <Piece NumberOfPoints="4" NumberOfCells="1" >
+ <PointData>
+ <DataArray type="Float64" Name="HeatFlowRate" format="appended" RangeMin="0" RangeMax="0" offset="388" />
+ <DataArray type="Float64" Name="MassFlowRate" format="appended" RangeMin="0" RangeMax="0" offset="448" />
+ <DataArray type="Float64" Name="NodalForces" NumberOfComponents="2" format="appended" RangeMin="0" RangeMax="0" offset="508" />
+ <DataArray type="UInt64" Name="bulk_node_ids" format="appended" RangeMin="0" RangeMax="3" offset="568" />
+ <DataArray type="Float64" Name="displacement" NumberOfComponents="2" format="appended" RangeMin="0" RangeMax="0" offset="640" />
+ <DataArray type="Float64" Name="epsilon" NumberOfComponents="4" format="appended" RangeMin="0" RangeMax="0" offset="700" />
+ <DataArray type="Float64" Name="ice_volume_fraction" format="appended" RangeMin="6.3169798575e-36" RangeMax="6.3169798575e-36" offset="760" />
+ <DataArray type="Float64" Name="pressure" format="appended" RangeMin="0" RangeMax="0" offset="832" />
+ <DataArray type="Float64" Name="pressure_interpolated" format="appended" RangeMin="0" RangeMax="0" offset="892" />
+ <DataArray type="Float64" Name="sigma" NumberOfComponents="4" format="appended" RangeMin="0" RangeMax="0" offset="952" />
+ <DataArray type="Float64" Name="sigma_ice" NumberOfComponents="4" format="appended" RangeMin="1.5630471519e-26" RangeMax="1.5630471519e-26" offset="1012" />
+ <DataArray type="Float64" Name="temperature" format="appended" RangeMin="277.15" RangeMax="277.15" offset="1104" />
+ <DataArray type="Float64" Name="temperature_interpolated" format="appended" RangeMin="277.15" RangeMax="277.15" offset="1172" />
+ </PointData>
+ <CellData>
+ <DataArray type="Int32" Name="MaterialIDs" format="appended" RangeMin="0" RangeMax="0" offset="1240" />
+ <DataArray type="UInt64" Name="bulk_element_ids" format="appended" RangeMin="0" RangeMax="0" offset="1300" />
+ <DataArray type="Float64" Name="fluid_density_avg" format="appended" RangeMin="1000" RangeMax="1000" offset="1360" />
+ <DataArray type="Float64" Name="stress_avg" NumberOfComponents="4" format="appended" RangeMin="0" RangeMax="0" offset="1424" />
+ <DataArray type="Float64" Name="viscosity_avg" format="appended" RangeMin="0.001" RangeMax="0.001" offset="1484" />
+ </CellData>
+ <Points>
+ <DataArray type="Float64" Name="Points" NumberOfComponents="3" format="appended" RangeMin="0" RangeMax="0.14142135624" offset="1552" />
+ </Points>
+ <Cells>
+ <DataArray type="Int64" Name="connectivity" format="appended" RangeMin="" RangeMax="" offset="1628" />
+ <DataArray type="Int64" Name="offsets" format="appended" RangeMin="" RangeMax="" offset="1700" />
+ <DataArray type="UInt8" Name="types" format="appended" RangeMin="" RangeMax="" offset="1760" />
+ </Cells>
+ </Piece>
+ </UnstructuredGrid>
+ <AppendedData encoding="base64">
+ _AQAAAAAAAAAAgAAAAAAAAKYAAAAAAAAAYwAAAAAAAAA=eF6FzDEKgDAMheG7ZO4kTr2KSIgaJWCTktZBxLvb1UXH9374LhCtvDlVMcVsbSG501kgDtcrmi/sELsASokhQpEtEUqG9hxpYkdbcbaUTVlrA/o7/BCci+wtfCLj/QAlSj0gAQAAAAAAAAAAgAAAAAAAABoAAAAAAAAAIgAAAAAAAAA=eF4z0zPRM9E1NDQ21003NU9KSTFJMtVLySwqqQQAUwMHQw==AQAAAAAAAAAAgAAAAAAAAIAAAAAAAAAADAAAAAAAAAA=eF5jYBhYAAAAgAABAQAAAAAAAAAAgAAAAAAAAIAAAAAAAAAADAAAAAAAAAA=eF5jYBhYAAAAgAABAQAAAAAAAAAAgAAAAAAAACAAAAAAAAAACwAAAAAAAAA=eF5jYMAPAAAgAAE=AQAAAAAAAAAAgAAAAAAAACAAAAAAAAAACwAAAAAAAAA=eF5jYMAPAAAgAAE=AQAAAAAAAAAAgAAAAAAAAEAAAAAAAAAADAAAAAAAAAA=eF5jYKAMAAAAQAABAQAAAAAAAAAAgAAAAAAAACAAAAAAAAAAEwAAAAAAAAA=eF5jYIAARijNBKWZoTQAAHAABw==AQAAAAAAAAAAgAAAAAAAAEAAAAAAAAAADAAAAAAAAAA=eF5jYKAMAAAAQAABAQAAAAAAAAAAgAAAAAAAAIAAAAAAAAAADAAAAAAAAAA=eF5jYBhYAAAAgAABAQAAAAAAAAAAgAAAAAAAACAAAAAAAAAAFQAAAAAAAAA=eF7TaiqWkT29wEIDSuug8QHAtAvrAQAAAAAAAAAAgAAAAAAAACAAAAAAAAAACwAAAAAAAAA=eF5jYMAPAAAgAAE=AQAAAAAAAAAAgAAAAAAAACAAAAAAAAAACwAAAAAAAAA=eF5jYMAPAAAgAAE=AQAAAAAAAAAAgAAAAAAAAIAAAAAAAAAADAAAAAAAAAA=eF5jYBhYAAAAgAABAQAAAAAAAAAAgAAAAAAAAIAAAAAAAAAAIwAAAAAAAAA=eF4TCF6Wfja8bZcAGs0ABTxQPjoNkxeG8tFpYvUDACnwLeI=AQAAAAAAAAAAgAAAAAAAACAAAAAAAAAAEQAAAAAAAAA=eF5LSwOCoEKHNBw0AMiADAU=AQAAAAAAAAAAgAAAAAAAACAAAAAAAAAAEQAAAAAAAAA=eF5LSwOCoEKHNBw0AMiADAU=AQAAAAAAAAAAgAAAAAAAAAQAAAAAAAAADAAAAAAAAAA=eF5jYGBgAAAABAABAQAAAAAAAAAAgAAAAAAAAAgAAAAAAAAACwAAAAAAAAA=eF5jYIAAAAAIAAE=AQAAAAAAAAAAgAAAAAAAAAgAAAAAAAAADgAAAAAAAAA=eF5jYAACh34HAAImARA=AQAAAAAAAAAAgAAAAAAAACAAAAAAAAAACwAAAAAAAAA=eF5jYMAPAAAgAAE=AQAAAAAAAAAAgAAAAAAAAAgAAAAAAAAAEAAAAAAAAAA=eF77s/LjJd+kAHsAGYAEpw==AQAAAAAAAAAAgAAAAAAAAGAAAAAAAAAAFwAAAAAAAAA=eF5jYMAOZs0EgZ32uMQJycP4AOa2Ej0=AQAAAAAAAAAAgAAAAAAAACAAAAAAAAAAEwAAAAAAAAA=eF5jYIAARijNBKWZoTQAAHAABw==AQAAAAAAAAAAgAAAAAAAAAgAAAAAAAAACwAAAAAAAAA=eF5jYYAAAAAoAAU=AQAAAAAAAAAAgAAAAAAAAAEAAAAAAAAACQAAAAAAAAA=eF7jBAAACgAK
+ </AppendedData>
+</VTKFile>