diff --git a/MaterialLib/MPL/Properties/PermeabilityMohrCoulombFailureIndexModel.cpp b/MaterialLib/MPL/Properties/PermeabilityMohrCoulombFailureIndexModel.cpp
index 48ccccaf6b9421bfea13a808565d76f29e80feba..0aedac2079f36651decbf1ac65b8056e61e06348 100644
--- a/MaterialLib/MPL/Properties/PermeabilityMohrCoulombFailureIndexModel.cpp
+++ b/MaterialLib/MPL/Properties/PermeabilityMohrCoulombFailureIndexModel.cpp
@@ -76,7 +76,7 @@ PermeabilityMohrCoulombFailureIndexModel<DisplacementDim>::value(
double const /*dt*/) const
{
auto const& stress_vector = std::get<SymmetricTensor<DisplacementDim>>(
- variable_array[static_cast<int>(Variable::stress)]);
+ variable_array[static_cast<int>(Variable::total_stress)]);
auto const& stress_tensor =
formEigenTensor<3>(static_cast<PropertyDataType>(stress_vector));
diff --git a/MaterialLib/MPL/VariableType.h b/MaterialLib/MPL/VariableType.h
index 3871530f3efeabdfc4a46682fb550f34cc8fa925..27638b595d698aad20c6d70f8eedcf5ebed12ec6 100644
--- a/MaterialLib/MPL/VariableType.h
+++ b/MaterialLib/MPL/VariableType.h
@@ -56,6 +56,7 @@ enum class Variable : int
strain,
stress,
temperature,
+ total_stress,
transport_porosity,
volumetric_strain,
number_of_variables
diff --git a/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h b/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h
index d3b664d7f2837185f9043cf1d20df7ad116d54b4..17be2b67eb6a9dfa57b24f3d23cd131cd1b81060 100644
--- a/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h
+++ b/ProcessLib/HydroMechanics/HydroMechanicsFEM-impl.h
@@ -234,7 +234,7 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
.template value<double>(vars, x_position, t, dt);
// For stress dependent permeability.
- vars[static_cast<int>(MPL::Variable::stress)].emplace<SymmetricTensor>(
+ vars[static_cast<int>(MPL::Variable::total_stress)].emplace<SymmetricTensor>(
MathLib::KelvinVector::kelvinVectorToSymmetricTensor(
(_ip_data[ip].sigma_eff - alpha * identity2 * p_int_pt)
.eval()));
@@ -399,7 +399,7 @@ HydroMechanicsLocalAssembler<ShapeFunctionDisplacement, ShapeFunctionPressure,
.template value<double>(vars, x_position, t, dt);
// For stress dependent permeability.
- vars[static_cast<int>(MPL::Variable::stress)].emplace<SymmetricTensor>(
+ vars[static_cast<int>(MPL::Variable::total_stress)].emplace<SymmetricTensor>(
MathLib::KelvinVector::kelvinVectorToSymmetricTensor(
(_ip_data[ip].sigma_eff - alpha * identity2 * p_int_pt)
.eval()));
@@ -508,7 +508,7 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
.template value<double>(vars, x_position, t, dt);
// For stress dependent permeability.
- vars[static_cast<int>(MPL::Variable::stress)].emplace<SymmetricTensor>(
+ vars[static_cast<int>(MPL::Variable::total_stress)].emplace<SymmetricTensor>(
MathLib::KelvinVector::kelvinVectorToSymmetricTensor(
(_ip_data[ip].sigma_eff - alpha_b * identity2 * p_int_pt)
.eval()));
@@ -872,6 +872,7 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
_integration_method.getNumberOfPoints();
auto const& medium = _process_data.media_map->getMedium(elem_id);
+ auto const& solid = medium->phase("Solid");
MPL::VariableArray vars;
SymmetricTensor k_sum = SymmetricTensor::Zero(KelvinVectorSize);
@@ -886,12 +887,22 @@ void HydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
auto const& sigma_eff = _ip_data[ip].sigma_eff;
sigma_eff_sum += sigma_eff;
+ auto const alpha_b =
+ solid.property(MPL::PropertyType::biot_coefficient)
+ .template value<double>(vars, x_position, t, dt);
+ auto const& identity2 = MathLib::KelvinVector::Invariants<
+ MathLib::KelvinVector::KelvinVectorDimensions<
+ DisplacementDim>::value>::identity2;
+
+ double const p_int_pt = _ip_data[ip].N_p.dot(p);
+ vars[static_cast<int>(MPL::Variable::phase_pressure)] = p_int_pt;
+
vars[static_cast<int>(MPL::Variable::strain)].emplace<SymmetricTensor>(
MathLib::KelvinVector::kelvinVectorToSymmetricTensor(eps));
- vars[static_cast<int>(MPL::Variable::stress)].emplace<SymmetricTensor>(
- MathLib::KelvinVector::kelvinVectorToSymmetricTensor(sigma_eff));
- vars[static_cast<int>(MPL::Variable::phase_pressure)] =
- _ip_data[ip].N_p.dot(p);
+ vars[static_cast<int>(MPL::Variable::total_stress)]
+ .emplace<SymmetricTensor>(
+ MathLib::KelvinVector::kelvinVectorToSymmetricTensor(
+ (sigma_eff - alpha_b * identity2 * p_int_pt).eval()));
k_sum += MPL::getSymmetricTensor<DisplacementDim>(
medium->property(MPL::PropertyType::permeability)
.value(vars, x_position, t, dt));
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
index 249e2f8bd99390384a6824a274503d729df56152..415ebe9e6348fbc97bbd851720d2d41de4cddcf5 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
@@ -472,13 +472,14 @@ void RichardsMechanicsLocalAssembler<
liquid_phase.property(MPL::PropertyType::viscosity)
.template value<double>(variables, x_position, t, dt);
+ auto const sigma_total = (_ip_data[ip].sigma_eff + sigma_sw +
+ alpha * chi_S_L * identity2 * p_cap_ip)
+ .eval();
// For stress dependent permeability.
- variables[static_cast<int>(MPL::Variable::stress)]
+ variables[static_cast<int>(MPL::Variable::total_stress)]
.emplace<SymmetricTensor>(
MathLib::KelvinVector::kelvinVectorToSymmetricTensor(
- (_ip_data[ip].sigma_eff +
- alpha * chi_S_L * identity2 * p_cap_ip)
- .eval()));
+ sigma_total));
auto const K_intrinsic = MPL::formEigenTensor<DisplacementDim>(
solid_phase.property(MPL::PropertyType::permeability)
@@ -809,13 +810,14 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
liquid_phase.property(MPL::PropertyType::viscosity)
.template value<double>(variables, x_position, t, dt);
+ auto const sigma_total = (_ip_data[ip].sigma_eff + sigma_sw +
+ alpha * chi_S_L * identity2 * p_cap_ip)
+ .eval();
// For stress dependent permeability.
- variables[static_cast<int>(MPL::Variable::stress)]
+ variables[static_cast<int>(MPL::Variable::total_stress)]
.emplace<SymmetricTensor>(
MathLib::KelvinVector::kelvinVectorToSymmetricTensor(
- (_ip_data[ip].sigma_eff +
- alpha * chi_S_L * identity2 * p_cap_ip)
- .eval()));
+ sigma_total));
auto const K_intrinsic = MPL::formEigenTensor<DisplacementDim>(
solid_phase.property(MPL::PropertyType::permeability)
.value(variables, x_position, t, dt));
@@ -1540,13 +1542,14 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
liquid_phase.property(MPL::PropertyType::density)
.template value<double>(variables, x_position, t, dt);
+ auto const sigma_total = (_ip_data[ip].sigma_eff + sigma_sw +
+ alpha * chi_S_L * identity2 * p_cap_ip)
+ .eval();
// For stress dependent permeability.
- variables[static_cast<int>(MPL::Variable::stress)]
+ variables[static_cast<int>(MPL::Variable::total_stress)]
.emplace<SymmetricTensor>(
MathLib::KelvinVector::kelvinVectorToSymmetricTensor(
- (_ip_data[ip].sigma_eff +
- alpha * chi_S_L * identity2 * p_cap_ip)
- .eval()));
+ sigma_total));
auto const K_intrinsic = MPL::formEigenTensor<DisplacementDim>(
solid_phase.property(MPL::PropertyType::permeability)
.value(variables, x_position, t, dt));
diff --git a/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsFEM-impl.h b/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsFEM-impl.h
index ef09e40b13e61740549bff4c4b3a432ca2d7fb83..bfc81f6f3011ca6a5d1462d1e43705cb4bc7cd39 100644
--- a/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsFEM-impl.h
+++ b/ProcessLib/ThermoHydroMechanics/ThermoHydroMechanicsFEM-impl.h
@@ -263,7 +263,7 @@ void ThermoHydroMechanicsLocalAssembler<ShapeFunctionDisplacement,
auto const beta_SR = (1 - alpha) * solid_skeleton_compressibility;
// For stress dependent permeability.
- vars[static_cast<int>(MaterialPropertyLib::Variable::stress)]
+ vars[static_cast<int>(MaterialPropertyLib::Variable::total_stress)]
.emplace<SymmetricTensor>(
MathLib::KelvinVector::kelvinVectorToSymmetricTensor(
(_ip_data[ip].sigma_eff - alpha * identity2 * p_int_pt)
@@ -529,7 +529,7 @@ std::vector<double> const& ThermoHydroMechanicsLocalAssembler<
.property(MaterialPropertyLib::PropertyType::biot_coefficient)
.template value<double>(vars, x_position, t, dt);
// For stress dependent permeability.
- vars[static_cast<int>(MaterialPropertyLib::Variable::stress)]
+ vars[static_cast<int>(MaterialPropertyLib::Variable::total_stress)]
.emplace<SymmetricTensor>(
MathLib::KelvinVector::kelvinVectorToSymmetricTensor(
(_ip_data[ip].sigma_eff - alpha * identity2 * p_int_pt)
diff --git a/Tests/MaterialLib/TestMPLPermeabilityMohrCoulombFailureIndexModel.cpp b/Tests/MaterialLib/TestMPLPermeabilityMohrCoulombFailureIndexModel.cpp
index 778812ef0d06baea3e095dbc8e82b2550e83bb91..dc70a7304e2c172f34d3260bac642748b06abaca 100644
--- a/Tests/MaterialLib/TestMPLPermeabilityMohrCoulombFailureIndexModel.cpp
+++ b/Tests/MaterialLib/TestMPLPermeabilityMohrCoulombFailureIndexModel.cpp
@@ -49,7 +49,7 @@ TEST(MaterialPropertyLib, PermeabilityMohrCoulombFailureIndexModel)
stress[5] = -1.39944e+7;
MPL::VariableArray vars;
- vars[static_cast<int>(MPL::Variable::stress)].emplace<SymmetricTensor>(
+ vars[static_cast<int>(MPL::Variable::total_stress)].emplace<SymmetricTensor>(
stress);
ParameterLib::SpatialPosition const pos;
@@ -71,7 +71,7 @@ TEST(MaterialPropertyLib, PermeabilityMohrCoulombFailureIndexModel)
stress[3] = -1e+5;
stress[4] = -2200.2;
stress[5] = -8e+5;
- vars[static_cast<int>(MPL::Variable::stress)] = stress;
+ vars[static_cast<int>(MPL::Variable::total_stress)] = stress;
auto const k_non_f =
MPL::formEigenTensor<3>(k_model.value(vars, pos, t, dt));
@@ -91,7 +91,7 @@ TEST(MaterialPropertyLib, PermeabilityMohrCoulombFailureIndexModel)
stress[3] = 0.0;
stress[4] = 0.0;
stress[5] = 0.0;
- vars[static_cast<int>(MPL::Variable::stress)] = stress;
+ vars[static_cast<int>(MPL::Variable::total_stress)] = stress;
auto const k_apex_f =
MPL::formEigenTensor<3>(k_model.value(vars, pos, t, dt));