[PL/THM] Correct history variable and push-back

ufz · Aug 31, 2023 · 38ac4be · 38ac4be
1 parent 699959a
commit 38ac4be
Show file tree

Hide file tree

Showing 3 changed files with 17 additions and 26 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;
@@ -81,7 +80,7 @@ struct IntegrationPointData final
     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;

diff --git a/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsFEM-impl.h b/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsFEM-impl.h
@@ -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) *

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