diff --git a/MaterialLib/MPL/Properties/PorosityFromMassBalance.cpp b/MaterialLib/MPL/Properties/PorosityFromMassBalance.cpp
index 5ef873bdb6d5b208bf11c8e81c60a7cc3470f1b7..8109d2a3517b5c2ea01baa3199bc2ebe17059144 100644
--- a/MaterialLib/MPL/Properties/PorosityFromMassBalance.cpp
+++ b/MaterialLib/MPL/Properties/PorosityFromMassBalance.cpp
@@ -34,10 +34,21 @@ void PorosityFromMassBalance::checkScale() const
}
}
+PropertyDataType PorosityFromMassBalance::value(
+ VariableArray const& /*variable_array*/,
+ ParameterLib::SpatialPosition const& /*pos*/, double const /*t*/,
+ double const /*dt*/) const
+{
+ OGS_FATAL(
+ "PorosityFromMassBalance value call requires previous time step "
+ "values.");
+}
+
PropertyDataType PorosityFromMassBalance::value(
VariableArray const& variable_array,
- ParameterLib::SpatialPosition const& pos,
- double const t, double const dt) const
+ VariableArray const& variable_array_prev,
+ ParameterLib::SpatialPosition const& pos, double const t,
+ double const dt) const
{
double const beta_SR = std::get<double>(
variable_array[static_cast<int>(Variable::grain_compressibility)]);
@@ -46,17 +57,24 @@ PropertyDataType PorosityFromMassBalance::value(
->property(PropertyType::biot_coefficient)
.template value<double>(variable_array, pos, t, dt);
- double const e_dot = std::get<double>(
- variable_array[static_cast<int>(Variable::volumetric_strain_rate)]);
+ double const e = std::get<double>(
+ variable_array[static_cast<int>(Variable::volumetric_strain)]);
+ double const e_prev = std::get<double>(
+ variable_array_prev[static_cast<int>(Variable::volumetric_strain)]);
+ double const delta_e = e - e_prev;
- double const p_eff_dot = std::get<double>(variable_array[static_cast<int>(
- Variable::effective_pore_pressure_rate)]);
+ double const p_eff = std::get<double>(
+ variable_array[static_cast<int>(Variable::effective_pore_pressure)]);
+ double const p_eff_prev =
+ std::get<double>(variable_array_prev[static_cast<int>(
+ Variable::effective_pore_pressure)]);
+ double const delta_p_eff = p_eff - p_eff_prev;
- double const phi = std::get<double>(
- variable_array[static_cast<int>(Variable::porosity)]);
+ double const phi_prev = std::get<double>(
+ variable_array_prev[static_cast<int>(Variable::porosity)]);
- double const w = dt * (e_dot + p_eff_dot * beta_SR);
- return std::clamp((phi + alpha_b * w) / (1 + w), phi_min_, phi_max_);
+ double const w = delta_e + delta_p_eff * beta_SR;
+ return std::clamp((phi_prev + alpha_b * w) / (1 + w), phi_min_, phi_max_);
}
PropertyDataType PorosityFromMassBalance::dValue(
diff --git a/MaterialLib/MPL/Properties/PorosityFromMassBalance.h b/MaterialLib/MPL/Properties/PorosityFromMassBalance.h
index 20a91e0ea27e08491dd678e6273ddda277cf71a0..1bf2e0a95897a37090745fdb3c8cf9400ce453aa 100644
--- a/MaterialLib/MPL/Properties/PorosityFromMassBalance.h
+++ b/MaterialLib/MPL/Properties/PorosityFromMassBalance.h
@@ -51,6 +51,10 @@ public:
PropertyDataType value(VariableArray const& variable_array,
ParameterLib::SpatialPosition const& pos,
double const t, double const dt) const override;
+ PropertyDataType value(VariableArray const& variable_array,
+ VariableArray const& variable_array_prev,
+ ParameterLib::SpatialPosition const& pos,
+ double const t, double const dt) const override;
PropertyDataType dValue(VariableArray const& variable_array,
Variable const variable,
ParameterLib::SpatialPosition const& pos,
diff --git a/MaterialLib/MPL/Properties/SaturationDependentSwelling.cpp b/MaterialLib/MPL/Properties/SaturationDependentSwelling.cpp
index 41137b68d3c6636a1e7a80e375572d520da3c441..f8e662ba0f8b78392591728ac793542b24ba420a 100644
--- a/MaterialLib/MPL/Properties/SaturationDependentSwelling.cpp
+++ b/MaterialLib/MPL/Properties/SaturationDependentSwelling.cpp
@@ -55,13 +55,14 @@ void SaturationDependentSwelling::checkScale() const
PropertyDataType SaturationDependentSwelling::value(
VariableArray const& variable_array,
+ VariableArray const& variable_array_prev,
ParameterLib::SpatialPosition const& pos, double const /*t*/,
double const dt) const
{
auto const S_L = std::get<double>(
variable_array[static_cast<int>(Variable::liquid_saturation)]);
- auto const S_L_dot = std::get<double>(
- variable_array[static_cast<int>(Variable::liquid_saturation_rate)]);
+ auto const S_L_prev = std::get<double>(
+ variable_array_prev[static_cast<int>(Variable::liquid_saturation)]);
Eigen::Matrix<double, 3, 3> const e =
local_coordinate_system_ == nullptr
@@ -76,8 +77,6 @@ PropertyDataType SaturationDependentSwelling::value(
return delta_sigma_sw; // still being zero.
}
- double const S_L_prev = S_L - S_L_dot * dt;
-
double const S_eff = std::clamp((S_L - S_min_) / (S_max_ - S_min_), 0., 1.);
double const S_eff_prev =
std::clamp((S_L_prev - S_min_) / (S_max_ - S_min_), 0., 1.);
@@ -100,9 +99,10 @@ PropertyDataType SaturationDependentSwelling::value(
}
PropertyDataType SaturationDependentSwelling::dValue(
- VariableArray const& variable_array, Variable const primary_variable,
+ VariableArray const& variable_array,
+ VariableArray const& variable_array_prev, Variable const primary_variable,
ParameterLib::SpatialPosition const& pos, double const /*t*/,
- double const dt) const
+ double const /*dt*/) const
{
(void)primary_variable;
assert((primary_variable == Variable::liquid_saturation) &&
@@ -111,8 +111,8 @@ PropertyDataType SaturationDependentSwelling::dValue(
auto const S_L = std::get<double>(
variable_array[static_cast<int>(Variable::liquid_saturation)]);
- auto const S_L_dot = std::get<double>(
- variable_array[static_cast<int>(Variable::liquid_saturation_rate)]);
+ auto const S_L_prev = std::get<double>(
+ variable_array_prev[static_cast<int>(Variable::liquid_saturation)]);
Eigen::Matrix<double, 3, 3> const e =
local_coordinate_system_ == nullptr
@@ -127,8 +127,6 @@ PropertyDataType SaturationDependentSwelling::dValue(
return delta_sigma_sw; // still being zero.
}
- double const S_L_prev = S_L - S_L_dot * dt;
-
double const S_eff = std::clamp((S_L - S_min_) / (S_max_ - S_min_), 0., 1.);
double const S_eff_prev =
std::clamp((S_L_prev - S_min_) / (S_max_ - S_min_), 0., 1.);
diff --git a/MaterialLib/MPL/Properties/SaturationDependentSwelling.h b/MaterialLib/MPL/Properties/SaturationDependentSwelling.h
index eeaf4d0158f1e35d27db75a61ea26687fe9ea71e..254a3bda9a83d0b94886e524c6134d73f16a197f 100644
--- a/MaterialLib/MPL/Properties/SaturationDependentSwelling.h
+++ b/MaterialLib/MPL/Properties/SaturationDependentSwelling.h
@@ -56,9 +56,11 @@ public:
void checkScale() const override;
PropertyDataType value(VariableArray const& variable_array,
+ VariableArray const& variable_array_prev,
ParameterLib::SpatialPosition const& pos,
double const t, double const dt) const override;
PropertyDataType dValue(VariableArray const& variable_array,
+ VariableArray const& variable_array_prev,
Variable const variable,
ParameterLib::SpatialPosition const& pos,
double const t, double const dt) const override;
diff --git a/MaterialLib/MPL/Properties/SwellingStress/LinearSaturationSwellingStress.cpp b/MaterialLib/MPL/Properties/SwellingStress/LinearSaturationSwellingStress.cpp
index 7897b12697dd315d38df91689b8072740b5ed997..3b1ddbca7c1c652130578f262c4375b3b0727698 100644
--- a/MaterialLib/MPL/Properties/SwellingStress/LinearSaturationSwellingStress.cpp
+++ b/MaterialLib/MPL/Properties/SwellingStress/LinearSaturationSwellingStress.cpp
@@ -24,10 +24,21 @@ LinearSaturationSwellingStress::LinearSaturationSwellingStress(
name_ = std::move(name);
}
+PropertyDataType LinearSaturationSwellingStress::value(
+ const VariableArray& /*variable_array*/,
+ const ParameterLib::SpatialPosition& /*pos*/, const double /*t*/,
+ const double /*dt*/) const
+{
+ OGS_FATAL(
+ "LinearSaturationSwellingStress value call requires previous time step "
+ "values.");
+}
+
PropertyDataType LinearSaturationSwellingStress::value(
const VariableArray& variable_array,
+ const VariableArray& variable_array_prev,
const ParameterLib::SpatialPosition& /*pos*/, const double /*t*/,
- const double dt) const
+ const double /*dt*/) const
{
// Sl <= S_max is guaranteed by the saturation property or
// the saturation calculation.
@@ -39,12 +50,10 @@ PropertyDataType LinearSaturationSwellingStress::value(
return 0.0;
}
- const double dS =
- dt *
- std::get<double>(
- variable_array[static_cast<int>(Variable::liquid_saturation_rate)]);
+ const double Sl_prev = std::get<double>(
+ variable_array_prev[static_cast<int>(Variable::liquid_saturation)]);
- return coefficient_ * dS;
+ return coefficient_ * (Sl - Sl_prev);
}
PropertyDataType LinearSaturationSwellingStress::dValue(
diff --git a/MaterialLib/MPL/Properties/SwellingStress/LinearSaturationSwellingStress.h b/MaterialLib/MPL/Properties/SwellingStress/LinearSaturationSwellingStress.h
index 2837f5e807ed0a4a2e6119c7940f966c2d92e108..dbb4784583907f3931900ba8f170553c0c0af849 100644
--- a/MaterialLib/MPL/Properties/SwellingStress/LinearSaturationSwellingStress.h
+++ b/MaterialLib/MPL/Properties/SwellingStress/LinearSaturationSwellingStress.h
@@ -77,8 +77,14 @@ public:
}
}
+ PropertyDataType value(VariableArray const& variable_array,
+ ParameterLib::SpatialPosition const& pos,
+ double const t,
+ double const dt) const override;
+
/// \return \f$\Delta {{\sigma}}^{\text{sw}} \f$.
PropertyDataType value(VariableArray const& variable_array,
+ VariableArray const& variable_array_prev,
ParameterLib::SpatialPosition const& pos,
double const t,
double const dt) const override;
diff --git a/MaterialLib/MPL/Properties/TransportPorosityFromMassBalance.cpp b/MaterialLib/MPL/Properties/TransportPorosityFromMassBalance.cpp
index 90d618132e75452b86804378ec330d6c2aceb5f0..26ed86ba10eded696407444e6d286b09d79672df 100644
--- a/MaterialLib/MPL/Properties/TransportPorosityFromMassBalance.cpp
+++ b/MaterialLib/MPL/Properties/TransportPorosityFromMassBalance.cpp
@@ -36,6 +36,7 @@ void TransportPorosityFromMassBalance::checkScale() const
PropertyDataType TransportPorosityFromMassBalance::value(
VariableArray const& variable_array,
+ VariableArray const& variable_array_prev,
ParameterLib::SpatialPosition const& pos, double const t,
double const dt) const
{
@@ -46,22 +47,39 @@ PropertyDataType TransportPorosityFromMassBalance::value(
->property(PropertyType::biot_coefficient)
.template value<double>(variable_array, pos, t, dt);
- double const e_dot = std::get<double>(
- variable_array[static_cast<int>(Variable::volumetric_strain_rate)]);
+ double const e = std::get<double>(
+ variable_array[static_cast<int>(Variable::volumetric_strain)]);
+ double const e_prev = std::get<double>(
+ variable_array_prev[static_cast<int>(Variable::volumetric_strain)]);
+ double const delta_e = e - e_prev;
- double const p_eff_dot = std::get<double>(variable_array[static_cast<int>(
- Variable::effective_pore_pressure_rate)]);
+ double const p_eff = std::get<double>(
+ variable_array[static_cast<int>(Variable::effective_pore_pressure)]);
+ double const p_eff_prev =
+ std::get<double>(variable_array_prev[static_cast<int>(
+ Variable::effective_pore_pressure)]);
+ double const delta_p_eff = p_eff - p_eff_prev;
double const phi =
std::get<double>(variable_array[static_cast<int>(Variable::porosity)]);
double const phi_tr_prev = std::get<double>(
- variable_array[static_cast<int>(Variable::transport_porosity)]);
+ variable_array_prev[static_cast<int>(Variable::transport_porosity)]);
- double const w = dt * (e_dot + p_eff_dot * beta_SR);
+ double const w = delta_e + delta_p_eff * beta_SR;
return std::clamp(phi_tr_prev + (alpha_b - phi) * w, phi_min_, phi_max_);
}
+PropertyDataType TransportPorosityFromMassBalance::value(
+ VariableArray const& /*variable_array*/,
+ ParameterLib::SpatialPosition const& /*pos*/, double const /*t*/,
+ double const /*dt*/) const
+{
+ OGS_FATAL(
+ "TransportPorosityFromMassBalance value call requires previous time "
+ "step values.");
+}
+
PropertyDataType TransportPorosityFromMassBalance::dValue(
VariableArray const& /*variable_array*/,
Variable const /*primary_variable*/,
diff --git a/MaterialLib/MPL/Properties/TransportPorosityFromMassBalance.h b/MaterialLib/MPL/Properties/TransportPorosityFromMassBalance.h
index 824a6d3ef0fee7807e7068203b1cbf7dc9dbd714..64f89b9ff2e62871ab05b4360e3ad8388e475dbf 100644
--- a/MaterialLib/MPL/Properties/TransportPorosityFromMassBalance.h
+++ b/MaterialLib/MPL/Properties/TransportPorosityFromMassBalance.h
@@ -48,6 +48,11 @@ public:
}
PropertyDataType value(VariableArray const& variable_array,
+ ParameterLib::SpatialPosition const& pos,
+ double const t,
+ double const dt) const override;
+ PropertyDataType value(VariableArray const& variable_array,
+ VariableArray const& variable_array_prev,
ParameterLib::SpatialPosition const& pos,
double const t, double const dt) const override;
PropertyDataType dValue(VariableArray const& variable_array,
diff --git a/MaterialLib/MPL/Property.cpp b/MaterialLib/MPL/Property.cpp
index 6c500292442c7dde865df4a56c6a76081bbbc031..677d2adfb30a5a194bded17099f33e2846e23b97 100644
--- a/MaterialLib/MPL/Property.cpp
+++ b/MaterialLib/MPL/Property.cpp
@@ -74,18 +74,23 @@ PropertyDataType Property::value() const
return value_;
}
-/// The default implementation of this method only returns the property value
-/// without altering it.
PropertyDataType Property::value(VariableArray const& /*variable_array*/,
+ VariableArray const& /*variable_array_prev*/,
ParameterLib::SpatialPosition const& /*pos*/,
double const /*t*/, double const /*dt*/) const
{
return value_;
}
-/// The default implementation of this method only returns the
-/// property value derivative without altering it.
+PropertyDataType Property::value(VariableArray const& variable_array,
+ ParameterLib::SpatialPosition const& pos,
+ double const t, double const dt) const
+{
+ return value(variable_array, VariableArray{}, pos, t, dt);
+}
+
PropertyDataType Property::dValue(VariableArray const& /*variable_array*/,
+ VariableArray const& /*variable_array_prev*/,
Variable const /*variable*/,
ParameterLib::SpatialPosition const& /*pos*/,
double const /*t*/, double const /*dt*/) const
@@ -93,6 +98,16 @@ PropertyDataType Property::dValue(VariableArray const& /*variable_array*/,
return dvalue_;
}
+/// The default implementation of this method only returns the
+/// property value derivative without altering it.
+PropertyDataType Property::dValue(VariableArray const& variable_array,
+ Variable const variable,
+ ParameterLib::SpatialPosition const& pos,
+ double const t, double const dt) const
+{
+ return dValue(variable_array, VariableArray{}, variable, pos, t, dt);
+}
+
/// Default implementation: 2nd derivative of any constant property is zero.
PropertyDataType Property::d2Value(VariableArray const& /*variable_array*/,
Variable const /*variable*/,
diff --git a/MaterialLib/MPL/Property.h b/MaterialLib/MPL/Property.h
index 84c2b75b9823e67434353606c431a3f97988dae4..ff1799bd454da121ce9268093a99a5a55191778a 100644
--- a/MaterialLib/MPL/Property.h
+++ b/MaterialLib/MPL/Property.h
@@ -54,13 +54,30 @@ public:
/// This virtual method simply returns the private value_ attribute without
/// changing it.
virtual PropertyDataType value() const;
- /// This virtual method will compute the property value based on the primary
- /// variables that are passed as arguments.
+ /// This virtual method will compute the property value based on the
+ /// variables that are passed as arguments and the variables from the
+ /// previous time step.
virtual PropertyDataType value(VariableArray const& variable_array,
+ VariableArray const& variable_array_prev,
ParameterLib::SpatialPosition const& pos,
double const t, double const dt) const;
- /// This virtual method will compute the derivative of a property
- /// with respect to the given variable pv.
+ /// This virtual method will compute the property value based on the
+ /// variables that are passed as arguments with the default implementation
+ /// using empty variables array for the previous time step.
+ virtual PropertyDataType value(VariableArray const& variable_array,
+ ParameterLib::SpatialPosition const& pos,
+ double const t, double const dt) const;
+ /// This virtual method will compute the property derivative value based on
+ /// the variables that are passed as arguments and the variables from the
+ /// previous time step.
+ virtual PropertyDataType dValue(VariableArray const& variable_array,
+ VariableArray const& variable_array_prev,
+ Variable const variable,
+ ParameterLib::SpatialPosition const& pos,
+ double const t, double const dt) const;
+ /// This virtual method will compute the property derivative value based on
+ /// the variables that are passed as arguments with the default
+ /// implementation using empty variables array for the previous time step.
virtual PropertyDataType dValue(VariableArray const& variable_array,
Variable const variable,
ParameterLib::SpatialPosition const& pos,
@@ -116,6 +133,29 @@ public:
}
}
+ template <typename T>
+ T value(VariableArray const& variable_array,
+ VariableArray const& variable_array_prev,
+ ParameterLib::SpatialPosition const& pos, double const t,
+ double const dt) const
+ {
+ try
+ {
+ return std::get<T>(
+ value(variable_array, variable_array_prev, pos, t, dt));
+ }
+ catch (std::bad_variant_access const&)
+ {
+ OGS_FATAL(
+ "The value of {:s} is not of the requested type '{:s}' but a "
+ "{:s}.",
+ description(),
+ typeid(T).name(),
+ property_data_type_names_[value(variable_array,
+ variable_array_prev, pos, t, dt)
+ .index()]);
+ }
+ }
template <typename T>
T value(VariableArray const& variable_array,
ParameterLib::SpatialPosition const& pos, double const t,
@@ -137,6 +177,28 @@ public:
}
}
+ template <typename T>
+ T dValue(VariableArray const& variable_array,
+ VariableArray const& variable_array_prev, Variable const variable,
+ ParameterLib::SpatialPosition const& pos, double const t,
+ double const dt) const
+ {
+ try
+ {
+ return std::get<T>(dValue(variable_array, variable_array_prev,
+ variable, pos, t, dt));
+ }
+ catch (std::bad_variant_access const&)
+ {
+ OGS_FATAL(
+ "The first derivative value of {:s} is not of the requested "
+ "type '{:s}' but a {:s}.",
+ description(),
+ typeid(T).name(),
+ property_data_type_names_
+ [dValue(variable_array, variable, pos, t, dt).index()]);
+ }
+ }
template <typename T>
T dValue(VariableArray const& variable_array, Variable const variable,
ParameterLib::SpatialPosition const& pos, double const t,
diff --git a/MaterialLib/MPL/VariableType.h b/MaterialLib/MPL/VariableType.h
index 85102143c87772f1201b49ed443f39b61484b495..b5c16c996fa324088b25e0c9a4de9ef24126d554 100644
--- a/MaterialLib/MPL/VariableType.h
+++ b/MaterialLib/MPL/VariableType.h
@@ -46,10 +46,9 @@ enum class Variable : int
concentration,
density,
displacement,
- effective_pore_pressure_rate,
+ effective_pore_pressure,
grain_compressibility,
liquid_saturation,
- liquid_saturation_rate,
phase_pressure,
porosity,
solid_grain_pressure,
@@ -57,7 +56,7 @@ enum class Variable : int
stress,
temperature,
transport_porosity,
- volumetric_strain_rate,
+ volumetric_strain,
number_of_variables
};
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
index 4ec1da3122e0ff00354d2ac1f0d208d58b841140..4f1b15bff1e2ba1eb675b1c594d6d91508edc664 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
@@ -283,6 +283,7 @@ void RichardsMechanicsLocalAssembler<
auto const& liquid_phase = medium->phase("AqueousLiquid");
auto const& solid_phase = medium->phase("Solid");
MPL::VariableArray variables;
+ MPL::VariableArray variables_prev;
ParameterLib::SpatialPosition x_position;
x_position.setElementID(_element.getID());
@@ -361,8 +362,8 @@ void RichardsMechanicsLocalAssembler<
S_L = medium->property(MPL::PropertyType::saturation)
.template value<double>(variables, x_position, t, dt);
variables[static_cast<int>(MPL::Variable::liquid_saturation)] = S_L;
- variables[static_cast<int>(MPL::Variable::liquid_saturation_rate)] =
- (S_L - S_L_prev) / dt;
+ variables_prev[static_cast<int>(MPL::Variable::liquid_saturation)] =
+ S_L_prev;
double const dS_L_dp_cap =
medium->property(MPL::PropertyType::saturation)
@@ -379,29 +380,27 @@ void RichardsMechanicsLocalAssembler<
double const chi_S_L = chi(S_L);
double const chi_S_L_prev = chi(S_L_prev);
- variables[static_cast<int>(
- MPL::Variable::effective_pore_pressure_rate)] =
- (chi_S_L * (-p_cap_ip) -
- chi_S_L_prev * (-p_cap_ip + p_cap_dot_ip * dt)) /
- dt;
+ variables[static_cast<int>(MPL::Variable::effective_pore_pressure)] =
+ -chi_S_L * p_cap_ip;
+ variables_prev[static_cast<int>(
+ MPL::Variable::effective_pore_pressure)] =
+ -chi_S_L_prev * (p_cap_ip - p_cap_dot_ip * dt);
// Set volumetric strain rate for the general case without swelling.
- double const div_u_dot = identity2.transpose() * B * u_dot;
- variables[static_cast<int>(MPL::Variable::volumetric_strain_rate)]
- .emplace<double>(div_u_dot);
+ variables[static_cast<int>(MPL::Variable::volumetric_strain)]
+ .emplace<double>(identity2.transpose() * B * u);
+ variables_prev[static_cast<int>(MPL::Variable::volumetric_strain)]
+ .emplace<double>(identity2.transpose() * B * (u - u_dot * dt));
- auto& porosity = _ip_data[ip].porosity;
+ auto& phi = _ip_data[ip].porosity;
{ // Porosity update
-
- // Use previous time step porosity for porosity update, ...
- variables[static_cast<int>(MPL::Variable::porosity)] =
+ variables_prev[static_cast<int>(MPL::Variable::porosity)] =
_ip_data[ip].porosity_prev;
- porosity =
- solid_phase.property(MPL::PropertyType::porosity)
- .template value<double>(variables, x_position, t, dt);
- // ... then use new porosity.
- variables[static_cast<int>(MPL::Variable::porosity)] = porosity;
+ phi = solid_phase.property(MPL::PropertyType::porosity)
+ .template value<double>(variables, variables_prev,
+ x_position, t, dt);
}
+ variables[static_cast<int>(MPL::Variable::porosity)] = phi;
// Swelling and possibly volumetric strain rate update.
auto& sigma_sw = _ip_data[ip].sigma_sw;
@@ -419,39 +418,40 @@ void RichardsMechanicsLocalAssembler<
MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
solid_phase
.property(MPL::PropertyType::swelling_stress_rate)
- .template value<DimMatrix>(variables, x_position, t,
- dt));
+ .template value<DimMatrix>(
+ variables, variables_prev, x_position, t, dt));
sigma_sw += sigma_sw_dot * dt;
auto const C_el = _ip_data[ip].computeElasticTangentStiffness(
t, x_position, dt, temperature);
- variables[static_cast<int>(
- MPL::Variable::volumetric_strain_rate)]
- .emplace<double>(div_u_dot + identity2.transpose() *
- C_el.inverse() *
- sigma_sw_dot);
+ // !!! Misusing volumetric strain for mechanical volumetric
+ // strain just to update the transport porosity !!!
+ std::get<double>(variables[static_cast<int>(
+ MPL::Variable::volumetric_strain)]) +=
+ identity2.transpose() * C_el.inverse() * sigma_sw;
+ std::get<double>(variables_prev[static_cast<int>(
+ MPL::Variable::volumetric_strain)]) +=
+ identity2.transpose() * C_el.inverse() * sigma_sw_prev;
}
if (solid_phase.hasProperty(MPL::PropertyType::transport_porosity))
{
- double& phi_tr = _ip_data[ip].transport_porosity;
- // Use previous time step transport_porosity for
- // transport_porosity update, ...
- variables[static_cast<int>(MPL::Variable::transport_porosity)] =
+ variables_prev[static_cast<int>(MPL::Variable::transport_porosity)] =
_ip_data[ip].transport_porosity_prev;
- // ... then use new transport_porosity.
- phi_tr =
+
+ _ip_data[ip].transport_porosity =
solid_phase.property(MPL::PropertyType::transport_porosity)
- .template value<double>(variables, x_position, t, dt);
+ .template value<double>(variables, variables_prev,
+ x_position, t, dt);
variables[static_cast<int>(MPL::Variable::transport_porosity)] =
- phi_tr;
+ _ip_data[ip].transport_porosity;
}
else
{
variables[static_cast<int>(MPL::Variable::transport_porosity)] =
- porosity;
+ phi;
}
}
@@ -492,18 +492,17 @@ void RichardsMechanicsLocalAssembler<
displacement_index, displacement_index)
.noalias() += B.transpose() * C * B * w;
- double const rho = rho_SR * (1 - porosity) + S_L * porosity * rho_LR;
+ double const rho = rho_SR * (1 - phi) + S_L * phi * rho_LR;
rhs.template segment<displacement_size>(displacement_index).noalias() +=
N_u_op.transpose() * rho * b * w - B.transpose() * sigma_sw * w;
//
// pressure equation, pressure part.
//
- double const a0 = S_L * (alpha - porosity) * beta_SR;
+ double const a0 = S_L * (alpha - phi) * beta_SR;
// Volumetric average specific storage of the solid and fluid phases.
double const specific_storage =
- dS_L_dp_cap * (p_cap_ip * a0 - porosity) +
- S_L * (porosity / K_LR + a0);
+ dS_L_dp_cap * (p_cap_ip * a0 - phi) + S_L * (phi / K_LR + a0);
M.template block<pressure_size, pressure_size>(pressure_index,
pressure_index)
.noalias() += N_p.transpose() * rho_LR * specific_storage * N_p * w;
@@ -622,6 +621,7 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
auto const& liquid_phase = medium->phase("AqueousLiquid");
auto const& solid_phase = medium->phase("Solid");
MPL::VariableArray variables;
+ MPL::VariableArray variables_prev;
ParameterLib::SpatialPosition x_position;
x_position.setElementID(_element.getID());
@@ -694,8 +694,8 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
S_L = medium->property(MPL::PropertyType::saturation)
.template value<double>(variables, x_position, t, dt);
variables[static_cast<int>(MPL::Variable::liquid_saturation)] = S_L;
- variables[static_cast<int>(MPL::Variable::liquid_saturation_rate)] =
- (S_L - S_L_prev) / dt;
+ variables_prev[static_cast<int>(MPL::Variable::liquid_saturation)] =
+ S_L_prev;
double const dS_L_dp_cap =
medium->property(MPL::PropertyType::saturation)
@@ -712,28 +712,27 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
double const chi_S_L = chi(S_L);
double const chi_S_L_prev = chi(S_L_prev);
- variables[static_cast<int>(
- MPL::Variable::effective_pore_pressure_rate)] =
- (chi_S_L * (-p_cap_ip) -
- chi_S_L_prev * (-p_cap_ip + p_cap_dot_ip * dt)) /
- dt;
+ variables[static_cast<int>(MPL::Variable::effective_pore_pressure)] =
+ -chi_S_L * p_cap_ip;
+ variables_prev[static_cast<int>(
+ MPL::Variable::effective_pore_pressure)] =
+ -chi_S_L_prev * (p_cap_ip - p_cap_dot_ip * dt);
// Set volumetric strain rate for the general case without swelling.
- double const div_u_dot = identity2.transpose() * B * u_dot;
- variables[static_cast<int>(MPL::Variable::volumetric_strain_rate)]
- .emplace<double>(div_u_dot);
+ variables[static_cast<int>(MPL::Variable::volumetric_strain)]
+ .emplace<double>(identity2.transpose() * B * u);
+ variables_prev[static_cast<int>(MPL::Variable::volumetric_strain)]
+ .emplace<double>(identity2.transpose() * B * (u - u_dot * dt));
- auto& porosity = _ip_data[ip].porosity;
+ auto& phi = _ip_data[ip].porosity;
{ // Porosity update
- // Use previous time step porosity for porosity update, ...
- variables[static_cast<int>(MPL::Variable::porosity)] =
+ variables_prev[static_cast<int>(MPL::Variable::porosity)] =
_ip_data[ip].porosity_prev;
- porosity =
- solid_phase.property(MPL::PropertyType::porosity)
- .template value<double>(variables, x_position, t, dt);
- // ... then use new porosity.
- variables[static_cast<int>(MPL::Variable::porosity)] = porosity;
+ phi = solid_phase.property(MPL::PropertyType::porosity)
+ .template value<double>(variables, variables_prev,
+ x_position, t, dt);
+ variables[static_cast<int>(MPL::Variable::porosity)] = phi;
}
// Swelling and possibly volumetric strain rate update.
@@ -752,36 +751,36 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
solid_phase
.property(MPL::PropertyType::swelling_stress_rate)
- .template value<DimMatrix>(variables, x_position, t,
- dt));
+ .template value<DimMatrix>(
+ variables, variables_prev, x_position, t, dt));
sigma_sw += sigma_sw_dot * dt;
- variables[static_cast<int>(
- MPL::Variable::volumetric_strain_rate)]
- .emplace<double>(div_u_dot + identity2.transpose() *
- C_el.inverse() *
- sigma_sw_dot);
+ // !!! Misusing volumetric strain for mechanical volumetric
+ // strain just to update the transport porosity !!!
+ std::get<double>(variables[static_cast<int>(
+ MPL::Variable::volumetric_strain)]) +=
+ identity2.transpose() * C_el.inverse() * sigma_sw;
+ std::get<double>(variables_prev[static_cast<int>(
+ MPL::Variable::volumetric_strain)]) +=
+ identity2.transpose() * C_el.inverse() * sigma_sw_prev;
}
if (solid_phase.hasProperty(MPL::PropertyType::transport_porosity))
{
- double& phi_tr = _ip_data[ip].transport_porosity;
-
- // Use previous time step transport_porosity for
- // transport_porosity update, ...
- variables[static_cast<int>(MPL::Variable::transport_porosity)] =
+ variables_prev[static_cast<int>(MPL::Variable::transport_porosity)] =
_ip_data[ip].transport_porosity_prev;
- // ... then use new transport_porosity.
- phi_tr =
+
+ _ip_data[ip].transport_porosity =
solid_phase.property(MPL::PropertyType::transport_porosity)
- .template value<double>(variables, x_position, t, dt);
+ .template value<double>(variables, variables_prev,
+ x_position, t, dt);
variables[static_cast<int>(MPL::Variable::transport_porosity)] =
- phi_tr;
+ _ip_data[ip].transport_porosity;
}
else
{
variables[static_cast<int>(MPL::Variable::transport_porosity)] =
- porosity;
+ phi;
}
}
@@ -818,7 +817,6 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
displacement_index, displacement_index)
.noalias() += B.transpose() * C * B * w;
- auto const phi = _ip_data[ip].porosity;
double const p_FR = -chi_S_L * p_cap_ip;
// p_SR
variables[static_cast<int>(MPL::Variable::solid_grain_pressure)] =
@@ -827,7 +825,7 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
solid_phase.property(MPL::PropertyType::density)
.template value<double>(variables, x_position, t, dt);
- double const rho = rho_SR * (1 - porosity) + S_L * porosity * rho_LR;
+ double const rho = rho_SR * (1 - phi) + S_L * phi * rho_LR;
local_rhs.template segment<displacement_size>(displacement_index)
.noalias() -= (B.transpose() * (sigma_eff + sigma_sw) -
N_u_op.transpose() * rho * b) *
@@ -854,7 +852,7 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
.template block<displacement_size, pressure_size>(
displacement_index, pressure_index)
.noalias() +=
- N_u_op.transpose() * porosity * rho_LR * dS_L_dp_cap * b * N_p * w;
+ N_u_op.transpose() * phi * rho_LR * dS_L_dp_cap * b * N_p * w;
if (solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate))
{
@@ -864,8 +862,9 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
solid_phase
.property(MPL::PropertyType::swelling_stress_rate)
.template dValue<DimMatrix>(
- variables, MPL::Variable::liquid_saturation,
- x_position, t, dt));
+ variables, variables_prev,
+ MPL::Variable::liquid_saturation, x_position, t,
+ dt));
local_Jac
.template block<displacement_size, pressure_size>(
displacement_index, pressure_index)
@@ -884,12 +883,12 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
laplace_p.noalias() +=
dNdx_p.transpose() * k_rel * rho_Ki_over_mu * dNdx_p * w;
- double const a0 = (alpha - porosity) * beta_SR;
- double const specific_storage_a_p = S_L * (porosity / K_LR + S_L * a0);
- double const specific_storage_a_S = porosity - p_cap_ip * S_L * a0;
+ double const a0 = (alpha - phi) * beta_SR;
+ double const specific_storage_a_p = S_L * (phi / K_LR + S_L * a0);
+ double const specific_storage_a_S = phi - p_cap_ip * S_L * a0;
double const dspecific_storage_a_p_dp_cap =
- dS_L_dp_cap * (porosity / K_LR + 2 * S_L * a0);
+ dS_L_dp_cap * (phi / K_LR + 2 * S_L * a0);
double const dspecific_storage_a_S_dp_cap =
-a0 * (S_L + p_cap_ip * dS_L_dp_cap);
@@ -1365,6 +1364,7 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
auto const& liquid_phase = medium->phase("AqueousLiquid");
auto const& solid_phase = medium->phase("Solid");
MPL::VariableArray variables;
+ MPL::VariableArray variables_prev;
ParameterLib::SpatialPosition x_position;
x_position.setElementID(_element.getID());
@@ -1416,8 +1416,8 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
S_L = medium->property(MPL::PropertyType::saturation)
.template value<double>(variables, x_position, t, dt);
variables[static_cast<int>(MPL::Variable::liquid_saturation)] = S_L;
- variables[static_cast<int>(MPL::Variable::liquid_saturation_rate)] =
- (S_L - S_L_prev) / dt;
+ variables_prev[static_cast<int>(MPL::Variable::liquid_saturation)] =
+ S_L_prev;
auto const chi = [medium, x_position, t, dt](double const S_L) {
MPL::VariableArray variables;
@@ -1442,33 +1442,30 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
variables[static_cast<int>(MPL::Variable::grain_compressibility)] =
beta_SR;
- variables[static_cast<int>(
- MPL::Variable::effective_pore_pressure_rate)] =
- (chi_S_L * (-p_cap_ip) -
- chi_S_L_prev * (-p_cap_ip + p_cap_dot_ip * dt)) /
- dt;
+ variables[static_cast<int>(MPL::Variable::effective_pore_pressure)] =
+ -chi_S_L * p_cap_ip;
+ variables_prev[static_cast<int>(
+ MPL::Variable::effective_pore_pressure)] =
+ -chi_S_L_prev * (p_cap_ip - p_cap_dot_ip * dt);
// Set volumetric strain rate for the general case without swelling.
- double const div_u_dot = identity2.transpose() * B * u_dot;
- variables[static_cast<int>(MPL::Variable::volumetric_strain_rate)]
- .emplace<double>(div_u_dot);
+ variables[static_cast<int>(MPL::Variable::volumetric_strain)]
+ .emplace<double>(identity2.transpose() * B * u);
+ variables_prev[static_cast<int>(MPL::Variable::volumetric_strain)]
+ .emplace<double>(identity2.transpose() * B * (u - u_dot * dt));
- auto& porosity = _ip_data[ip].porosity;
+ auto& phi = _ip_data[ip].porosity;
{ // Porosity update
-
- // Use previous time step porosity for porosity update, ...
- variables[static_cast<int>(MPL::Variable::porosity)] =
+ variables_prev[static_cast<int>(MPL::Variable::porosity)] =
_ip_data[ip].porosity_prev;
- porosity =
- solid_phase.property(MPL::PropertyType::porosity)
- .template value<double>(variables, x_position, t, dt);
- // ... then use new porosity.
- variables[static_cast<int>(MPL::Variable::porosity)] = porosity;
+ phi = solid_phase.property(MPL::PropertyType::porosity)
+ .template value<double>(variables, variables_prev,
+ x_position, t, dt);
+ variables[static_cast<int>(MPL::Variable::porosity)] = phi;
}
// Swelling and possibly volumetric strain rate update.
auto& sigma_sw = _ip_data[ip].sigma_sw;
- double& phi_tr = _ip_data[ip].transport_porosity;
{
auto const& sigma_sw_prev = _ip_data[ip].sigma_sw_prev;
@@ -1483,34 +1480,37 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
MathLib::KelvinVector::tensorToKelvin<DisplacementDim>(
solid_phase
.property(MPL::PropertyType::swelling_stress_rate)
- .template value<DimMatrix>(variables, x_position, t,
- dt));
+ .template value<DimMatrix>(
+ variables, variables_prev, x_position, t, dt));
sigma_sw += sigma_sw_dot * dt;
- variables[static_cast<int>(
- MPL::Variable::volumetric_strain_rate)]
- .emplace<double>(div_u_dot + identity2.transpose() *
- C_el.inverse() *
- sigma_sw_dot);
+ // !!! Misusing volumetric strain for mechanical volumetric
+ // strain just to update the transport porosity !!!
+ std::get<double>(variables[static_cast<int>(
+ MPL::Variable::volumetric_strain)]) +=
+ identity2.transpose() * C_el.inverse() * sigma_sw;
+ std::get<double>(variables_prev[static_cast<int>(
+ MPL::Variable::volumetric_strain)]) +=
+ identity2.transpose() * C_el.inverse() * sigma_sw_prev;
}
if (solid_phase.hasProperty(MPL::PropertyType::transport_porosity))
{
- // Use previous time step transport_porosity for
- // transport_porosity update, ...
- variables[static_cast<int>(MPL::Variable::transport_porosity)] =
+ variables_prev[static_cast<int>(
+ MPL::Variable::transport_porosity)] =
_ip_data[ip].transport_porosity_prev;
- // ... then use new transport_porosity.
- phi_tr =
+
+ _ip_data[ip].transport_porosity =
solid_phase.property(MPL::PropertyType::transport_porosity)
- .template value<double>(variables, x_position, t, dt);
+ .template value<double>(variables, variables_prev,
+ x_position, t, dt);
variables[static_cast<int>(MPL::Variable::transport_porosity)] =
- phi_tr;
+ _ip_data[ip].transport_porosity;
}
else
{
variables[static_cast<int>(MPL::Variable::transport_porosity)] =
- porosity;
+ phi;
}
}
auto const mu =
@@ -1537,7 +1537,6 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
GlobalDimMatrixType const K_over_mu = k_rel * K_intrinsic / mu;
- auto const phi = _ip_data[ip].porosity;
auto const& sigma_eff = _ip_data[ip].sigma_eff;
double const p_FR = -chi_S_L * p_cap_ip;
// p_SR
@@ -1561,9 +1560,7 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
-K_over_mu * dNdx_p * p_L + rho_LR * K_over_mu * b;
saturation_avg += S_L;
- porosity_avg +=
- _ip_data[ip].porosity; // Note, this is not updated, because needs
- // xdot and dt to be passed.
+ porosity_avg += phi;
sigma_avg += sigma_eff;
}
saturation_avg /= n_integration_points;
diff --git a/Tests/MaterialLib/TestMPLLinearSaturationSwellingStress.cpp b/Tests/MaterialLib/TestMPLLinearSaturationSwellingStress.cpp
index eb2b060c32a473410ed9a5e34ad62831061a501b..2640befd917de5cbcf34015b95b8c8d0e6e3494e 100644
--- a/Tests/MaterialLib/TestMPLLinearSaturationSwellingStress.cpp
+++ b/Tests/MaterialLib/TestMPLLinearSaturationSwellingStress.cpp
@@ -37,7 +37,7 @@ TEST(MaterialPropertyLib, LinearSaturationSwellingStress)
" <reference_saturation> 0.65 </reference_saturation>"
"</property>";
- auto const swelling_stress_rate =
+ auto const swelling_stress_increment =
createLinearSaturationSwellingStressModel(xml);
double const coefficient = 1.0e+6;
@@ -49,20 +49,21 @@ TEST(MaterialPropertyLib, LinearSaturationSwellingStress)
const double dS = S1 - S0;
MaterialPropertyLib::VariableArray variable_array;
- variable_array[static_cast<int>(
- MaterialPropertyLib::Variable::liquid_saturation_rate)] = dS / dt;
+ MaterialPropertyLib::VariableArray variable_array_prev;
variable_array[static_cast<int>(
MaterialPropertyLib::Variable::liquid_saturation)] = S1;
+ variable_array_prev[static_cast<int>(
+ MaterialPropertyLib::Variable::liquid_saturation)] = S0;
ParameterLib::SpatialPosition const pos;
double const time = std::numeric_limits<double>::quiet_NaN();
double d_sw_expected = coefficient * dS;
- double d_sw = std::get<double>(
- swelling_stress_rate->value(variable_array, pos, time, dt));
+ double d_sw = std::get<double>(swelling_stress_increment->value(
+ variable_array, variable_array_prev, pos, time, dt));
ASSERT_LE(std::fabs(d_sw_expected - d_sw), 1e-19)
- << "for expected swelling stress rate" << d_sw_expected
- << " for computed swelling stress rate." << d_sw_expected;
+ << "for expected swelling stress increment" << d_sw_expected
+ << " for computed swelling stress increment." << d_sw_expected;
- double dsw_dS = std::get<double>(swelling_stress_rate->dValue(
+ double dsw_dS = std::get<double>(swelling_stress_increment->dValue(
variable_array, MaterialPropertyLib::Variable::liquid_saturation, pos,
time, dt));
double dsw_dS_expected = coefficient;
@@ -74,15 +75,17 @@ TEST(MaterialPropertyLib, LinearSaturationSwellingStress)
double const S2 = 0.3;
variable_array[static_cast<int>(
MaterialPropertyLib::Variable::liquid_saturation)] = S2;
+ variable_array_prev[static_cast<int>(
+ MaterialPropertyLib::Variable::liquid_saturation)] = S1;
- d_sw = std::get<double>(
- swelling_stress_rate->value(variable_array, pos, time, dt));
+ d_sw = std::get<double>(swelling_stress_increment->value(
+ variable_array, variable_array_prev, pos, time, dt));
d_sw_expected = 0.0;
ASSERT_LE(std::fabs(d_sw_expected - d_sw), 1e-19)
- << "for expected swelling stress rate" << d_sw_expected
- << " for computed swelling stress rate." << d_sw_expected;
+ << "for expected swelling stress increment" << d_sw_expected
+ << " for computed swelling stress increment." << d_sw_expected;
- dsw_dS = std::get<double>(swelling_stress_rate->dValue(
+ dsw_dS = std::get<double>(swelling_stress_increment->dValue(
variable_array, MaterialPropertyLib::Variable::liquid_saturation, pos,
time, dt));
dsw_dS_expected = 0.0;