diff --git a/Documentation/ProjectFile/properties/property/BishopsPowerLaw/c_BishopsPowerLaw.md b/Documentation/ProjectFile/properties/property/BishopsPowerLaw/c_BishopsPowerLaw.md
new file mode 100644
index 0000000000000000000000000000000000000000..0d7d20a7b89eff6ed390ee028b5b6f62f1dfb25b
--- /dev/null
+++ b/Documentation/ProjectFile/properties/property/BishopsPowerLaw/c_BishopsPowerLaw.md
@@ -0,0 +1 @@
+\copydoc MaterialPropertyLib::BishopsPowerLaw
diff --git a/Documentation/ProjectFile/properties/property/BishopsPowerLaw/t_exponent.md b/Documentation/ProjectFile/properties/property/BishopsPowerLaw/t_exponent.md
new file mode 100644
index 0000000000000000000000000000000000000000..6f6002095400ad77e24cbb92d1472d038b8d747f
--- /dev/null
+++ b/Documentation/ProjectFile/properties/property/BishopsPowerLaw/t_exponent.md
@@ -0,0 +1 @@
+\copydoc MaterialPropertyLib::BishopsPowerLaw::_m
diff --git a/Documentation/ProjectFile/properties/property/BishopsSaturationCutoff/c_BishopsSaturationCutoff.md b/Documentation/ProjectFile/properties/property/BishopsSaturationCutoff/c_BishopsSaturationCutoff.md
new file mode 100644
index 0000000000000000000000000000000000000000..d3dc4657adbc8e1274db56790e33a8b44bf72c18
--- /dev/null
+++ b/Documentation/ProjectFile/properties/property/BishopsSaturationCutoff/c_BishopsSaturationCutoff.md
@@ -0,0 +1 @@
+\copydoc MaterialPropertyLib::BishopsSaturationCutoff
diff --git a/Documentation/ProjectFile/properties/property/BishopsSaturationCutoff/t_cutoff_value.md b/Documentation/ProjectFile/properties/property/BishopsSaturationCutoff/t_cutoff_value.md
new file mode 100644
index 0000000000000000000000000000000000000000..d4b89ee3b3f46469e736f2e4007ece743220b47a
--- /dev/null
+++ b/Documentation/ProjectFile/properties/property/BishopsSaturationCutoff/t_cutoff_value.md
@@ -0,0 +1 @@
+\copydoc MaterialPropertyLib::BishopsSaturationCutoff::_S_L_max
diff --git a/MaterialLib/MPL/CreateProperty.cpp b/MaterialLib/MPL/CreateProperty.cpp
index 7b164ac5527f409a3f90f2ec3aa0878640083631..1b17482ce9bcdba31393c1204c27a79b93b3b0db 100644
--- a/MaterialLib/MPL/CreateProperty.cpp
+++ b/MaterialLib/MPL/CreateProperty.cpp
@@ -112,6 +112,16 @@ std::unique_ptr<MaterialPropertyLib::Property> createProperty(
local_coordinate_system);
}
+ if (boost::iequals(property_type, "BishopsPowerLaw"))
+ {
+ return createBishopsPowerLaw(config);
+ }
+
+ if (boost::iequals(property_type, "BishopsSaturationCutoff"))
+ {
+ return createBishopsSaturationCutoff(config);
+ }
+
// If none of the above property types are found, OGS throws an error.
OGS_FATAL("The specified component property type '%s' was not recognized",
property_type.c_str());
diff --git a/MaterialLib/MPL/Properties/BishopsPowerLaw.cpp b/MaterialLib/MPL/Properties/BishopsPowerLaw.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..069d861dfd03fdd8b430f9d6879523550a24ed33
--- /dev/null
+++ b/MaterialLib/MPL/Properties/BishopsPowerLaw.cpp
@@ -0,0 +1,54 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2020, 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 "BishopsPowerLaw.h"
+
+namespace MaterialPropertyLib
+{
+BishopsPowerLaw::BishopsPowerLaw(double const exponent) : _m(exponent) {}
+
+void BishopsPowerLaw::setScale(
+ std::variant<Medium*, Phase*, Component*> scale_pointer)
+{
+ if (!std::holds_alternative<Medium*>(scale_pointer))
+ {
+ OGS_FATAL(
+ "The property 'BishopsPowerLaw' is implemented on the 'media' "
+ "scale only.");
+ }
+ _medium = std::get<Medium*>(scale_pointer);
+}
+
+PropertyDataType BishopsPowerLaw::value(
+ VariableArray const& variable_array,
+ 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)]);
+
+ return std::pow(S_L, _m);
+}
+
+PropertyDataType BishopsPowerLaw::dValue(
+ VariableArray const& variable_array, Variable const variable,
+ ParameterLib::SpatialPosition const& /*pos*/, double const /*t*/,
+ double const /*dt*/) const
+{
+ (void)variable;
+ assert((variable == Variable::liquid_saturation) &&
+ "BishopsPowerLaw::dvalue is implemented for derivatives with "
+ "respect to liquid saturation only.");
+
+ auto const S_L = std::get<double>(
+ variable_array[static_cast<int>(Variable::liquid_saturation)]);
+
+ return _m * std::pow(S_L, _m - 1.);
+}
+} // namespace MaterialPropertyLib
diff --git a/MaterialLib/MPL/Properties/BishopsPowerLaw.h b/MaterialLib/MPL/Properties/BishopsPowerLaw.h
new file mode 100644
index 0000000000000000000000000000000000000000..bee26d4d2693d9ef0e953ca52a71c56df15ee629
--- /dev/null
+++ b/MaterialLib/MPL/Properties/BishopsPowerLaw.h
@@ -0,0 +1,38 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2020, 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 "MaterialLib/MPL/Property.h"
+
+namespace MaterialPropertyLib
+{
+/// Bishop's power law for effective stress.
+class BishopsPowerLaw final : public Property
+{
+public:
+ explicit BishopsPowerLaw(double const exponent);
+
+ void setScale(
+ std::variant<Medium*, Phase*, Component*> scale_pointer) override;
+
+ PropertyDataType value(VariableArray const& variable_array,
+ 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*/,
+ double const /*t*/,
+ double const /*dt*/) const override;
+
+private:
+ Medium* _medium = nullptr;
+ double const _m; //< Exponent.
+};
+} // namespace MaterialPropertyLib
diff --git a/MaterialLib/MPL/Properties/BishopsSaturationCutoff.cpp b/MaterialLib/MPL/Properties/BishopsSaturationCutoff.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..8e25999e41531aa78727dbcf90bce25e3f813b80
--- /dev/null
+++ b/MaterialLib/MPL/Properties/BishopsSaturationCutoff.cpp
@@ -0,0 +1,55 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2020, 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 "BishopsSaturationCutoff.h"
+
+namespace MaterialPropertyLib
+{
+BishopsSaturationCutoff::BishopsSaturationCutoff(double const cutoff_value)
+ : _S_L_max(cutoff_value)
+{
+}
+
+void BishopsSaturationCutoff::setScale(
+ std::variant<Medium*, Phase*, Component*> scale_pointer)
+{
+ if (!std::holds_alternative<Medium*>(scale_pointer))
+ {
+ OGS_FATAL(
+ "The property 'BishopsSaturationCutoff' is implemented on the "
+ "'media' scale only.");
+ }
+ _medium = std::get<Medium*>(scale_pointer);
+}
+
+PropertyDataType BishopsSaturationCutoff::value(
+ VariableArray const& variable_array,
+ 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)]);
+
+ return S_L < _S_L_max ? 0. : 1.;
+}
+
+PropertyDataType BishopsSaturationCutoff::dValue(
+ VariableArray const& /*variable_array*/, Variable const variable,
+ ParameterLib::SpatialPosition const& /*pos*/, double const /*t*/,
+ double const /*dt*/) const
+{
+ (void)variable;
+ assert(
+ (variable == Variable::liquid_saturation) &&
+ "BishopsSaturationCutoff::dvalue is implemented for derivatives with "
+ "respect to liquid saturation only.");
+
+ return 0;
+}
+} // namespace MaterialPropertyLib
diff --git a/MaterialLib/MPL/Properties/BishopsSaturationCutoff.h b/MaterialLib/MPL/Properties/BishopsSaturationCutoff.h
new file mode 100644
index 0000000000000000000000000000000000000000..fff74493adcae40faa0fa88a23ddb7f486532c26
--- /dev/null
+++ b/MaterialLib/MPL/Properties/BishopsSaturationCutoff.h
@@ -0,0 +1,40 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2020, 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 "MaterialLib/MPL/Property.h"
+
+namespace MaterialPropertyLib
+{
+/// Bishop's effective stress model using saturation cutoff. The effective
+/// stress is 1 as long as the saturation does not fall below the saturation
+/// cutoff value.
+class BishopsSaturationCutoff final : public Property
+{
+public:
+ explicit BishopsSaturationCutoff(double const cutoff_value);
+
+ void setScale(
+ std::variant<Medium*, Phase*, Component*> scale_pointer) override;
+
+ PropertyDataType value(VariableArray const& variable_array,
+ 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*/,
+ double const /*t*/,
+ double const /*dt*/) const override;
+
+private:
+ Medium* _medium = nullptr;
+ double const _S_L_max; //< Maximum saturation cutoff value.
+};
+} // namespace MaterialPropertyLib
diff --git a/MaterialLib/MPL/Properties/CreateBishopsPowerLaw.cpp b/MaterialLib/MPL/Properties/CreateBishopsPowerLaw.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..782853ac14acf30547748e50eb536dd67b0b139f
--- /dev/null
+++ b/MaterialLib/MPL/Properties/CreateBishopsPowerLaw.cpp
@@ -0,0 +1,28 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2020, 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 "BaseLib/ConfigTree.h"
+#include "BishopsPowerLaw.h"
+
+namespace MaterialPropertyLib
+{
+std::unique_ptr<BishopsPowerLaw> createBishopsPowerLaw(
+ BaseLib::ConfigTree const& config)
+{
+ //! \ogs_file_param{properties__property__type}
+ config.checkConfigParameter("type", "BishopsPowerLaw");
+ DBUG("Create BishopsPowerLaw property");
+
+ auto const exponent =
+ //! \ogs_file_param{properties__property__BishopsPowerLaw__exponent}
+ config.getConfigParameter<double>("exponent");
+
+ return std::make_unique<MaterialPropertyLib::BishopsPowerLaw>(exponent);
+}
+} // namespace MaterialPropertyLib
diff --git a/MaterialLib/MPL/Properties/CreateBishopsPowerLaw.h b/MaterialLib/MPL/Properties/CreateBishopsPowerLaw.h
new file mode 100644
index 0000000000000000000000000000000000000000..dab374eb831635d61eb343e0f46a064fcf2aa3e5
--- /dev/null
+++ b/MaterialLib/MPL/Properties/CreateBishopsPowerLaw.h
@@ -0,0 +1,28 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2020, 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>
+
+namespace BaseLib
+{
+class ConfigTree;
+}
+
+namespace MaterialPropertyLib
+{
+class BishopsPowerLaw;
+}
+
+namespace MaterialPropertyLib
+{
+std::unique_ptr<BishopsPowerLaw> createBishopsPowerLaw(
+ BaseLib::ConfigTree const& config);
+} // namespace MaterialPropertyLib
diff --git a/MaterialLib/MPL/Properties/CreateBishopsSaturationCutoff.cpp b/MaterialLib/MPL/Properties/CreateBishopsSaturationCutoff.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..68c2728504935a6098616231a031520480d9c029
--- /dev/null
+++ b/MaterialLib/MPL/Properties/CreateBishopsSaturationCutoff.cpp
@@ -0,0 +1,29 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2020, 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 "BaseLib/ConfigTree.h"
+#include "BishopsSaturationCutoff.h"
+
+namespace MaterialPropertyLib
+{
+std::unique_ptr<BishopsSaturationCutoff> createBishopsSaturationCutoff(
+ BaseLib::ConfigTree const& config)
+{
+ //! \ogs_file_param{properties__property__type}
+ config.checkConfigParameter("type", "BishopsSaturationCutoff");
+ DBUG("Create BishopsSaturationCutoff property");
+
+ auto const cutoff_value =
+ //! \ogs_file_param{properties__property__BishopsSaturationCutoff__cutoff_value}
+ config.getConfigParameter<double>("cutoff_value");
+
+ return std::make_unique<MaterialPropertyLib::BishopsSaturationCutoff>(
+ cutoff_value);
+}
+} // namespace MaterialPropertyLib
diff --git a/MaterialLib/MPL/Properties/CreateBishopsSaturationCutoff.h b/MaterialLib/MPL/Properties/CreateBishopsSaturationCutoff.h
new file mode 100644
index 0000000000000000000000000000000000000000..2d728a735a1793ef0d93c65703244747b2aed1da
--- /dev/null
+++ b/MaterialLib/MPL/Properties/CreateBishopsSaturationCutoff.h
@@ -0,0 +1,28 @@
+/**
+ * \file
+ * \copyright
+ * Copyright (c) 2012-2020, 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>
+
+namespace BaseLib
+{
+class ConfigTree;
+}
+
+namespace MaterialPropertyLib
+{
+class BishopsSaturationCutoff;
+}
+
+namespace MaterialPropertyLib
+{
+std::unique_ptr<BishopsSaturationCutoff> createBishopsSaturationCutoff(
+ BaseLib::ConfigTree const& config);
+} // namespace MaterialPropertyLib
diff --git a/MaterialLib/MPL/Properties/CreateProperties.h b/MaterialLib/MPL/Properties/CreateProperties.h
index 970586dbdc73965711eb6c117d0499b8e80e6f67..337092a4fbad0ad183f642db19eb031a91542013 100644
--- a/MaterialLib/MPL/Properties/CreateProperties.h
+++ b/MaterialLib/MPL/Properties/CreateProperties.h
@@ -11,6 +11,8 @@
*/
#pragma once
+#include "CreateBishopsPowerLaw.h"
+#include "CreateBishopsSaturationCutoff.h"
#include "CreateConstant.h"
#include "CreateDupuitPermeability.h"
#include "CreateExponentialProperty.h"
diff --git a/MaterialLib/MPL/Properties/PorosityFromMassBalance.cpp b/MaterialLib/MPL/Properties/PorosityFromMassBalance.cpp
index a674feecfedb7484b40ebaadcafcb4ec2759867b..93fb827453c72ceb7ccbdcd53fa2be13ca30c278 100644
--- a/MaterialLib/MPL/Properties/PorosityFromMassBalance.cpp
+++ b/MaterialLib/MPL/Properties/PorosityFromMassBalance.cpp
@@ -52,13 +52,13 @@ PropertyDataType PorosityFromMassBalance::value(
double const e_dot = std::get<double>(
variable_array[static_cast<int>(Variable::volumetric_strain_rate)]);
- double const p_dot = std::get<double>(
- variable_array[static_cast<int>(Variable::phase_pressure_rate)]);
+ double const p_eff_dot = std::get<double>(variable_array[static_cast<int>(
+ Variable::effective_pore_pressure_rate)]);
double const phi = std::get<double>(
variable_array[static_cast<int>(Variable::porosity)]);
- double const w = dt * (e_dot + p_dot / K_SR);
+ double const w = dt * (e_dot + p_eff_dot / K_SR);
return std::clamp((phi + alpha_b * w) / (1 + w), 0., 1.);
}
diff --git a/MaterialLib/MPL/Properties/Properties.h b/MaterialLib/MPL/Properties/Properties.h
index eede46b6c6e74a1a5eedb352c94011c81c231158..a2924469967b39347f7cbf5f1e8d4165a03abb55 100644
--- a/MaterialLib/MPL/Properties/Properties.h
+++ b/MaterialLib/MPL/Properties/Properties.h
@@ -11,6 +11,8 @@
*/
#pragma once
+#include "BishopsPowerLaw.h"
+#include "BishopsSaturationCutoff.h"
#include "Constant.h"
#include "DupuitPermeability.h"
#include "ExponentialProperty.h"
diff --git a/MaterialLib/MPL/PropertyType.h b/MaterialLib/MPL/PropertyType.h
index 48d0c93e430b71cdfada9dcd837b5a8e90e48728..6c36701d8c8eaa8d60090c5650847a6d136bcc43 100644
--- a/MaterialLib/MPL/PropertyType.h
+++ b/MaterialLib/MPL/PropertyType.h
@@ -36,6 +36,7 @@ enum PropertyType : int
acentric_factor,
binary_interaction_coefficient,
biot_coefficient,
+ bishops_effective_stress,
brooks_corey_exponent,
bulk_modulus,
critical_density,
@@ -101,6 +102,10 @@ inline PropertyType convertStringToProperty(std::string const& inString)
{
return PropertyType::biot_coefficient;
}
+ if (boost::iequals(inString, "bishops_effective_stress"))
+ {
+ return PropertyType::bishops_effective_stress;
+ }
if (boost::iequals(inString, "brooks_corey_exponent"))
{
return PropertyType::brooks_corey_exponent;
@@ -266,6 +271,7 @@ static const std::array<std::string, PropertyType::number_of_properties>
property_enum_to_string{{"acentric_factor",
"binary_interaction_coefficient",
"biot_coefficient",
+ "bishops_effective_stress",
"brooks_corey_exponent",
"bulk_modulus",
"critical_density",
diff --git a/MaterialLib/MPL/VariableType.h b/MaterialLib/MPL/VariableType.h
index f88d1da5299a9b41bb3785285a4a4ac2ef259f41..4a29cf2004614f63b1c2abcb48c41fb6dcd2a27b 100644
--- a/MaterialLib/MPL/VariableType.h
+++ b/MaterialLib/MPL/VariableType.h
@@ -44,7 +44,7 @@ enum class Variable : int
{
concentration,
phase_pressure,
- phase_pressure_rate,
+ effective_pore_pressure_rate,
capillary_pressure,
density,
temperature,
diff --git a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
index 58dc2dfbfa594975ce208a88f5925db964c65aee..c1368c241a0eee7a44f0bc41d4a88730db8d58d5 100644
--- a/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
+++ b/ProcessLib/RichardsMechanics/RichardsMechanicsFEM-impl.h
@@ -251,11 +251,12 @@ void RichardsMechanicsLocalAssembler<
double p_cap_ip;
NumLib::shapeFunctionInterpolate(-p_L, N_p, p_cap_ip);
+ double p_cap_dot_ip;
+ NumLib::shapeFunctionInterpolate(-p_L_dot, N_p, p_cap_dot_ip);
+
variables[static_cast<int>(MPL::Variable::capillary_pressure)] =
p_cap_ip;
variables[static_cast<int>(MPL::Variable::phase_pressure)] = -p_cap_ip;
- variables[static_cast<int>(MPL::Variable::phase_pressure_rate)] =
- N_p.dot(p_L_dot);
auto const temperature =
medium->property(MPL::PropertyType::reference_temperature)
@@ -275,15 +276,6 @@ void RichardsMechanicsLocalAssembler<
liquid_phase.property(MPL::PropertyType::bulk_modulus)
.template value<double>(variables, x_position, t, dt);
- // Use previous time step porosity for porosity update, ...
- variables[static_cast<int>(MPL::Variable::porosity)] =
- _ip_data[ip].porosity_prev;
- auto& porosity = _ip_data[ip].porosity;
- 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;
-
auto const rho_LR =
liquid_phase.property(MPL::PropertyType::density)
.template value<double>(variables, x_position, t, dt);
@@ -302,6 +294,31 @@ void RichardsMechanicsLocalAssembler<
MPL::Variable::capillary_pressure,
x_position, t, dt);
+ auto const chi = [&](double const S_L) {
+ MPL::VariableArray variables;
+ variables[static_cast<int>(MPL::Variable::liquid_saturation)] = S_L;
+ return medium->property(MPL::PropertyType::bishops_effective_stress)
+ .template value<double>(variables, x_position, t, dt);
+ };
+ 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;
+
+ // Use previous time step porosity for porosity update, ...
+ variables[static_cast<int>(MPL::Variable::porosity)] =
+ _ip_data[ip].porosity_prev;
+ auto& porosity = _ip_data[ip].porosity;
+ 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;
+
+
double const k_rel =
medium->property(MPL::PropertyType::relative_permeability)
.template value<double>(variables, x_position, t, dt);
@@ -370,7 +387,8 @@ void RichardsMechanicsLocalAssembler<
//
K.template block<displacement_size, pressure_size>(displacement_index,
pressure_index)
- .noalias() -= B.transpose() * alpha * S_L * identity2 * N_p * w;
+ .noalias() -=
+ B.transpose() * alpha * chi_S_L * identity2 * N_p * w;
//
// pressure equation, displacement part.
@@ -511,8 +529,6 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
variables[static_cast<int>(MPL::Variable::capillary_pressure)] =
p_cap_ip;
variables[static_cast<int>(MPL::Variable::phase_pressure)] = -p_cap_ip;
- variables[static_cast<int>(MPL::Variable::phase_pressure_rate)] =
- N_p.dot(p_L_dot);
variables[static_cast<int>(MPL::Variable::volumetric_strain_rate)]
.emplace<double>(identity2.transpose() * B * u_dot);
auto const temperature =
@@ -539,15 +555,6 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
liquid_phase.property(MPL::PropertyType::bulk_modulus)
.template value<double>(variables, x_position, t, dt);
- // Use previous time step porosity for porosity update, ...
- variables[static_cast<int>(MPL::Variable::porosity)] =
- _ip_data[ip].porosity_prev;
- auto& porosity = _ip_data[ip].porosity;
- 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;
-
auto const rho_LR =
liquid_phase.property(MPL::PropertyType::density)
.template value<double>(variables, x_position, t, dt);
@@ -571,6 +578,15 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
variables, MPL::Variable::capillary_pressure,
MPL::Variable::capillary_pressure, x_position, t, dt);
+ auto const chi = [&](double const S_L) {
+ MPL::VariableArray variables;
+ variables[static_cast<int>(MPL::Variable::liquid_saturation)] = S_L;
+ return medium->property(MPL::PropertyType::bishops_effective_stress)
+ .template value<double>(variables, x_position, t, dt);
+ };
+ double const chi_S_L = chi(S_L);
+ double const chi_S_L_prev = chi(S_L_prev);
+
double const k_rel =
medium->property(MPL::PropertyType::relative_permeability)
.template value<double>(variables, x_position, t, dt);
@@ -582,6 +598,21 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
GlobalDimMatrixType const rho_Ki_over_mu = K_intrinsic * rho_LR / mu;
+ 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;
+
+ // Use previous time step porosity for porosity update, ...
+ variables[static_cast<int>(MPL::Variable::porosity)] =
+ _ip_data[ip].porosity_prev;
+ auto& porosity = _ip_data[ip].porosity;
+ 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;
+
sigma_sw = sigma_sw_prev;
if (solid_phase.hasProperty(MPL::PropertyType::swelling_stress_rate))
{
@@ -617,13 +648,19 @@ void RichardsMechanicsLocalAssembler<ShapeFunctionDisplacement,
//
// displacement equation, pressure part
//
- Kup.noalias() += B.transpose() * alpha * S_L * identity2 * N_p * w;
+ Kup.noalias() += B.transpose() * alpha * chi_S_L * identity2 * N_p * w;
+ auto const dchi_dS_L =
+ medium->property(MPL::PropertyType::bishops_effective_stress)
+ .template dValue<double>(variables,
+ MPL::Variable::liquid_saturation,
+ x_position, t, dt);
local_Jac
.template block<displacement_size, pressure_size>(
displacement_index, pressure_index)
.noalias() -= B.transpose() * alpha *
- (S_L + p_cap_ip * dS_L_dp_cap) * identity2 * N_p * w;
+ (chi_S_L + dchi_dS_L * p_cap_ip * dS_L_dp_cap) *
+ identity2 * N_p * w;
local_Jac
.template block<displacement_size, pressure_size>(
diff --git a/ProcessLib/RichardsMechanics/Tests.cmake b/ProcessLib/RichardsMechanics/Tests.cmake
index f72b6458778898d1ecdb5b396ab23fcde82efeaa..4b55cd167fd482c676aa5054f7d2a1d8b0f4528e 100644
--- a/ProcessLib/RichardsMechanics/Tests.cmake
+++ b/ProcessLib/RichardsMechanics/Tests.cmake
@@ -244,3 +244,39 @@ AddTest(
GLOB orthotropic_swelling_xy_pcs_0_ts_*.vtu porosity porosity 1e-15 1e-15
GLOB orthotropic_swelling_xy_pcs_0_ts_*.vtu porosity_avg porosity_avg 1e-15 1e-15
)
+AddTest(
+ NAME RichardsMechanics_bishops_effective_stress_power_law
+ PATH RichardsMechanics
+ EXECUTABLE ogs
+ EXECUTABLE_ARGS bishops_effective_stress_power_law.prj
+ WRAPPER time
+ TESTER vtkdiff
+ REQUIREMENTS NOT OGS_USE_MPI
+ DIFF_DATA
+ GLOB bishops_effective_stress_power_law_pcs_0_ts_*.vtu pressure pressure 1e-15 0
+ GLOB bishops_effective_stress_power_law_pcs_0_ts_*.vtu saturation saturation 1e-14 1e-15
+ GLOB bishops_effective_stress_power_law_pcs_0_ts_*.vtu displacement displacement 1e-15 0
+ GLOB bishops_effective_stress_power_law_pcs_0_ts_*.vtu sigma sigma 5e-15 0
+ GLOB bishops_effective_stress_power_law_pcs_0_ts_*.vtu epsilon epsilon 2e-15 0
+ GLOB bishops_effective_stress_power_law_pcs_0_ts_*.vtu velocity velocity 1e-15 0
+ GLOB bishops_effective_stress_power_law_pcs_0_ts_*.vtu porosity porosity 5e-15 1e-15
+ GLOB bishops_effective_stress_power_law_pcs_0_ts_*.vtu porosity_avg porosity_avg 5e-15 1e-15
+)
+AddTest(
+ NAME RichardsMechanics_bishops_effective_stress_saturation_cutoff
+ PATH RichardsMechanics
+ EXECUTABLE ogs
+ EXECUTABLE_ARGS bishops_effective_stress_saturation_cutoff.prj
+ WRAPPER time
+ TESTER vtkdiff
+ REQUIREMENTS NOT OGS_USE_MPI
+ DIFF_DATA
+ GLOB bishops_effective_stress_saturation_cutoff_pcs_0_ts_*.vtu pressure pressure 1e-15 0
+ GLOB bishops_effective_stress_saturation_cutoff_pcs_0_ts_*.vtu saturation saturation 1e-14 1e-15
+ GLOB bishops_effective_stress_saturation_cutoff_pcs_0_ts_*.vtu displacement displacement 1e-15 0
+ GLOB bishops_effective_stress_saturation_cutoff_pcs_0_ts_*.vtu sigma sigma 1e-15 0
+ GLOB bishops_effective_stress_saturation_cutoff_pcs_0_ts_*.vtu epsilon epsilon 1e-15 0
+ GLOB bishops_effective_stress_saturation_cutoff_pcs_0_ts_*.vtu velocity velocity 1e-15 0
+ GLOB bishops_effective_stress_saturation_cutoff_pcs_0_ts_*.vtu porosity porosity 5e-15 1e-15
+ GLOB bishops_effective_stress_saturation_cutoff_pcs_0_ts_*.vtu porosity_avg porosity_avg 5e-15 1e-15
+)
diff --git a/Tests/Data/RichardsMechanics/RichardsFlow_2d_quasinewton.prj b/Tests/Data/RichardsMechanics/RichardsFlow_2d_quasinewton.prj
index 07e96694c450527ce2614ed2e0dcbb5eed5db621..fc63969a3217dadaf231055ac1646b0b31cf2b86 100644
--- a/Tests/Data/RichardsMechanics/RichardsFlow_2d_quasinewton.prj
+++ b/Tests/Data/RichardsMechanics/RichardsFlow_2d_quasinewton.prj
@@ -112,6 +112,11 @@
<exponent>0.789029535864979</exponent>
<minimum_relative_permeability_liquid>1e-12</minimum_relative_permeability_liquid>
</property>
+ <property>
+ <name>bishops_effective_stress</name>
+ <type>BishopsPowerLaw</type>
+ <exponent>1</exponent>
+ </property>
</properties>
</medium>
</media>
diff --git a/Tests/Data/RichardsMechanics/RichardsFlow_2d_small.prj b/Tests/Data/RichardsMechanics/RichardsFlow_2d_small.prj
index d630df1ac8b956e8fbf38ab060d86baed05cd2a1..c975f69253a9cbff3f8b3f4dc1cdc2969aafffe0 100644
--- a/Tests/Data/RichardsMechanics/RichardsFlow_2d_small.prj
+++ b/Tests/Data/RichardsMechanics/RichardsFlow_2d_small.prj
@@ -107,6 +107,11 @@
<exponent>0.789029535864979</exponent>
<minimum_relative_permeability_liquid>1e-12</minimum_relative_permeability_liquid>
</property>
+ <property>
+ <name>bishops_effective_stress</name>
+ <type>BishopsPowerLaw</type>
+ <exponent>1</exponent>
+ </property>
</properties>
</medium>
</media>
diff --git a/Tests/Data/RichardsMechanics/RichardsFlow_2d_small_masslumping.prj b/Tests/Data/RichardsMechanics/RichardsFlow_2d_small_masslumping.prj
index 17d82de12b94ca3944cc0b68f30573708428e296..5c0be74da6b75219723dd5ff7a7c6694ebee1612 100644
--- a/Tests/Data/RichardsMechanics/RichardsFlow_2d_small_masslumping.prj
+++ b/Tests/Data/RichardsMechanics/RichardsFlow_2d_small_masslumping.prj
@@ -108,6 +108,11 @@
<exponent>0.789029535864979</exponent>
<minimum_relative_permeability_liquid>1e-12</minimum_relative_permeability_liquid>
</property>
+ <property>
+ <name>bishops_effective_stress</name>
+ <type>BishopsPowerLaw</type>
+ <exponent>1</exponent>
+ </property>
</properties>
</medium>
</media>
diff --git a/Tests/Data/RichardsMechanics/bishops_effective_stress_power_law.prj b/Tests/Data/RichardsMechanics/bishops_effective_stress_power_law.prj
new file mode 100644
index 0000000000000000000000000000000000000000..6e1ec61b27ab8ec280ff96dfc227bb413ec69816
--- /dev/null
+++ b/Tests/Data/RichardsMechanics/bishops_effective_stress_power_law.prj
@@ -0,0 +1,335 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+<OpenGeoSysProject>
+ <mesh>cube_1x1x1_hex20_1e0.vtu</mesh>
+ <geometry>cube_1x1x1.gml</geometry>
+ <processes>
+ <process>
+ <name>RM</name>
+ <type>RICHARDS_MECHANICS</type>
+ <integration_order>3</integration_order>
+ <dimension>3</dimension>
+ <!--
+ <jacobian_assembler>
+ <type>CentralDifferences</type>
+ <component_magnitudes>1e5 1e-5 1e-5 1e-5</component_magnitudes>
+ <relative_epsilons>1e-8 1e-8 1e-8 1e-8</relative_epsilons>
+ </jacobian_assembler>
+ -->
+ <constitutive_relation>
+ <type>LinearElasticIsotropic</type>
+ <youngs_modulus>E</youngs_modulus>
+ <poissons_ratio>nu</poissons_ratio>
+ </constitutive_relation>
+ <process_variables>
+ <displacement>displacement</displacement>
+ <pressure>pressure</pressure>
+ </process_variables>
+ <secondary_variables>
+ <secondary_variable internal_name="sigma" output_name="sigma"/>
+ <secondary_variable internal_name="epsilon" output_name="epsilon"/>
+ <secondary_variable internal_name="velocity" output_name="velocity"/>
+ <secondary_variable internal_name="saturation" output_name="saturation"/>
+ <secondary_variable internal_name="porosity" output_name="porosity"/>
+ </secondary_variables>
+ <specific_body_force>0 0 0</specific_body_force>
+ </process>
+ </processes>
+ <media>
+ <medium>
+ <phases>
+ <phase>
+ <type>AqueousLiquid</type>
+ <properties>
+ <property>
+ <name>bulk_modulus</name>
+ <type>Constant</type>
+ <value>1</value>
+ </property>
+ <property>
+ <name>viscosity</name>
+ <type>Constant</type>
+ <value>1</value>
+ </property>
+ <property>
+ <name>density</name>
+ <type>Constant</type>
+ <value>1</value>
+ </property>
+ </properties>
+ </phase>
+ <phase>
+ <type>Solid</type>
+ <properties>
+ <property>
+ <name>bulk_modulus</name>
+ <type>Constant</type>
+ <value>10</value>
+ </property>
+ <property>
+ <name>biot_coefficient</name>
+ <type>Constant</type>
+ <value>1</value>
+ </property>
+ <property>
+ <name>density</name>
+ <type>Constant</type>
+ <value>1</value>
+ </property>
+ <property>
+ <name>porosity</name>
+ <type>Constant</type>
+ <value>0.5</value>
+ </property>
+ <property>
+ <name>permeability</name>
+ <type>Constant</type>
+ <value>1</value>
+ </property>
+ </properties>
+ </phase>
+ </phases>
+ <properties>
+ <property>
+ <name>reference_temperature</name>
+ <type>Constant</type>
+ <value>293.15</value>
+ </property>
+ <property>
+ <name>saturation</name>
+ <type>SaturationVanGenuchten</type>
+ <residual_liquid_saturation>0.1</residual_liquid_saturation>
+ <residual_gas_saturation>0.05</residual_gas_saturation>
+ <exponent>0.8</exponent>
+ <entry_pressure>0.5</entry_pressure>
+ </property>
+ <property>
+ <name>relative_permeability</name>
+ <type>Constant</type>
+ <value>1</value>
+ </property>
+ <property>
+ <name>bishops_effective_stress</name>
+ <type>BishopsPowerLaw</type>
+ <exponent>3</exponent>
+ </property>
+ </properties>
+ </medium>
+ </media>
+ <time_loop>
+ <processes>
+ <process ref="RM">
+ <nonlinear_solver>nonlinear_solver</nonlinear_solver>
+ <compensate_non_equilibrium_initial_residuum>true</compensate_non_equilibrium_initial_residuum>
+ <convergence_criterion>
+ <type>PerComponentDeltaX</type>
+ <norm_type>NORM2</norm_type>
+ <abstols>1e-13 1e-13 1e-13 1e-13</abstols>
+ </convergence_criterion>
+ <time_discretization>
+ <type>BackwardEuler</type>
+ </time_discretization>
+ <time_stepping>
+ <type>FixedTimeStepping</type>
+ <t_initial>0</t_initial>
+ <t_end>3</t_end>
+ <timesteps>
+ <pair>
+ <repeat>50</repeat>
+ <delta_t>.02</delta_t>
+ </pair>
+
+ <pair>
+ <repeat>4</repeat>
+ <delta_t>.25</delta_t>
+ </pair>
+
+ <pair>
+ <repeat>50</repeat>
+ <delta_t>.02</delta_t>
+ </pair>
+ </timesteps>
+ </time_stepping>
+ </process>
+ </processes>
+ <output>
+ <type>VTK</type>
+ <prefix>bishops_effective_stress_power_law</prefix>
+ <timesteps>
+ <pair>
+ <repeat>1000</repeat>
+ <each_steps>1</each_steps>
+ </pair>
+ </timesteps>
+ <variables>
+ <variable>displacement</variable>
+ <variable>pressure</variable>
+ <variable>sigma</variable>
+ <variable>epsilon</variable>
+ <variable>velocity</variable>
+ <variable>saturation</variable>
+ <variable>porosity</variable>
+ </variables>
+ </output>
+ </time_loop>
+ <local_coordinate_system>
+ <basis_vector_0>e0</basis_vector_0>
+ <basis_vector_1>e1</basis_vector_1>
+ <basis_vector_2>e2</basis_vector_2>
+ </local_coordinate_system>
+ <parameters>
+ <parameter>
+ <name>e0</name>
+ <type>Constant</type>
+ <values>1 0 0</values>
+ </parameter>
+ <parameter>
+ <name>e1</name>
+ <type>Constant</type>
+ <values>0 1 0</values>
+ </parameter>
+ <parameter>
+ <name>e2</name>
+ <type>Constant</type>
+ <values>0 0 1</values>
+ </parameter>
+ <!-- Mechanics -->
+ <parameter>
+ <name>E</name>
+ <type>Constant</type>
+ <value>1</value>
+ </parameter>
+ <parameter>
+ <name>nu</name>
+ <type>Constant</type>
+ <value>.3</value>
+ </parameter>
+ <!-- Model parameters -->
+ <parameter>
+ <name>displacement0</name>
+ <type>Constant</type>
+ <values>0 0 0</values>
+ </parameter>
+ <parameter>
+ <name>p0</name>
+ <type>Constant</type>
+ <value>0.1</value>
+ </parameter>
+ <parameter>
+ <name>dirichlet0</name>
+ <type>Constant</type>
+ <values>0</values>
+ </parameter>
+ <parameter>
+ <name>pressure_bc</name>
+ <type>CurveScaled</type>
+ <curve>pressure_ramp</curve>
+ <parameter>p0</parameter>
+ </parameter>
+ </parameters>
+ <curves>
+ <curve>
+ <name>pressure_ramp</name>
+ <coords>0 1 2 3</coords>
+ <values>1 -10 -10 1</values>
+ </curve>
+ </curves>
+ <process_variables>
+ <process_variable>
+ <name>displacement</name>
+ <components>3</components>
+ <order>2</order>
+ <initial_condition>displacement0</initial_condition>
+ <boundary_conditions>
+ <boundary_condition>
+ <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+ <geometry>left</geometry>
+ <type>Dirichlet</type>
+ <component>0</component>
+ <parameter>dirichlet0</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+ <geometry>front</geometry>
+ <type>Dirichlet</type>
+ <component>1</component>
+ <parameter>dirichlet0</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+ <geometry>bottom</geometry>
+ <type>Dirichlet</type>
+ <component>2</component>
+ <parameter>dirichlet0</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+ <geometry>top</geometry>
+ <type>Dirichlet</type>
+ <component>2</component>
+ <parameter>dirichlet0</parameter>
+ </boundary_condition>
+ </boundary_conditions>
+ </process_variable>
+ <process_variable>
+ <name>pressure</name>
+ <components>1</components>
+ <order>1</order>
+ <initial_condition>p0</initial_condition>
+ <boundary_conditions>
+ <boundary_condition>
+ <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+ <geometry>right</geometry>
+ <type>Dirichlet</type>
+ <parameter>pressure_bc</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+ <geometry>bottom</geometry>
+ <type>Dirichlet</type>
+ <parameter>pressure_bc</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+ <geometry>front</geometry>
+ <type>Dirichlet</type>
+ <parameter>pressure_bc</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+ <geometry>left</geometry>
+ <type>Dirichlet</type>
+ <parameter>pressure_bc</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+ <geometry>top</geometry>
+ <type>Dirichlet</type>
+ <parameter>pressure_bc</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+ <geometry>back</geometry>
+ <type>Dirichlet</type>
+ <parameter>pressure_bc</parameter>
+ </boundary_condition>
+ </boundary_conditions>
+ </process_variable>
+ </process_variables>
+ <nonlinear_solvers>
+ <nonlinear_solver>
+ <name>nonlinear_solver</name>
+ <type>Newton</type>
+ <max_iter>100</max_iter>
+ <linear_solver>general_linear_solver</linear_solver>
+ </nonlinear_solver>
+ </nonlinear_solvers>
+ <linear_solvers>
+ <linear_solver>
+ <name>general_linear_solver</name>
+ <eigen>
+ <solver_type>SparseLU</solver_type>
+ <scaling>true</scaling>
+ </eigen>
+ </linear_solver>
+ </linear_solvers>
+</OpenGeoSysProject>
diff --git a/Tests/Data/RichardsMechanics/bishops_effective_stress_power_law_pcs_0_ts_104_t_3.000000.vtu b/Tests/Data/RichardsMechanics/bishops_effective_stress_power_law_pcs_0_ts_104_t_3.000000.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..6968e497ff126cbd1398dadf11e4548c7f878f0b
--- /dev/null
+++ b/Tests/Data/RichardsMechanics/bishops_effective_stress_power_law_pcs_0_ts_104_t_3.000000.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:f02b79a6560a16d684770632812fe55379674caedfee6a3f5092211c5696f937
+size 8187
diff --git a/Tests/Data/RichardsMechanics/bishops_effective_stress_power_law_pcs_0_ts_25_t_0.500000.vtu b/Tests/Data/RichardsMechanics/bishops_effective_stress_power_law_pcs_0_ts_25_t_0.500000.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..5e149d5aa840af32d5fe887ff4464318cfd3b59e
--- /dev/null
+++ b/Tests/Data/RichardsMechanics/bishops_effective_stress_power_law_pcs_0_ts_25_t_0.500000.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:974584d21ffcdabc7d84ad77b22a4006613ee22b719a552a403b1786f74408b1
+size 7423
diff --git a/Tests/Data/RichardsMechanics/bishops_effective_stress_power_law_pcs_0_ts_50_t_1.000000.vtu b/Tests/Data/RichardsMechanics/bishops_effective_stress_power_law_pcs_0_ts_50_t_1.000000.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..39d85628b4a62dfa38e5bd2a93103991736fe4c8
--- /dev/null
+++ b/Tests/Data/RichardsMechanics/bishops_effective_stress_power_law_pcs_0_ts_50_t_1.000000.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:de0db17f4546c25722993b64087a6b41168be7bbd8c1ab996e68a631ef44470d
+size 7299
diff --git a/Tests/Data/RichardsMechanics/bishops_effective_stress_power_law_pcs_0_ts_52_t_1.500000.vtu b/Tests/Data/RichardsMechanics/bishops_effective_stress_power_law_pcs_0_ts_52_t_1.500000.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..10bd69b5350bf54f6db098a62df0169658042deb
--- /dev/null
+++ b/Tests/Data/RichardsMechanics/bishops_effective_stress_power_law_pcs_0_ts_52_t_1.500000.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:338240d63d18184e0297419a33ecfa2308f866eb705978307334c41dd90fd8e7
+size 7283
diff --git a/Tests/Data/RichardsMechanics/bishops_effective_stress_power_law_pcs_0_ts_54_t_2.000000.vtu b/Tests/Data/RichardsMechanics/bishops_effective_stress_power_law_pcs_0_ts_54_t_2.000000.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..ca8ef43f40e1267b0635d5002f3e852d8ebf6a42
--- /dev/null
+++ b/Tests/Data/RichardsMechanics/bishops_effective_stress_power_law_pcs_0_ts_54_t_2.000000.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:abbb64b4f5da0d1468f73a898b4023150e39e88f748f38b71b512b81c2f29537
+size 7447
diff --git a/Tests/Data/RichardsMechanics/bishops_effective_stress_power_law_pcs_0_ts_79_t_2.500000.vtu b/Tests/Data/RichardsMechanics/bishops_effective_stress_power_law_pcs_0_ts_79_t_2.500000.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..0b5677ebfeb15ce41a3cb485663c130908d4cb23
--- /dev/null
+++ b/Tests/Data/RichardsMechanics/bishops_effective_stress_power_law_pcs_0_ts_79_t_2.500000.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:4d17a6b7751c56d1fda352d13f8211618422b791a6a993025594d4e3f703c5cc
+size 7331
diff --git a/Tests/Data/RichardsMechanics/bishops_effective_stress_saturation_cutoff.prj b/Tests/Data/RichardsMechanics/bishops_effective_stress_saturation_cutoff.prj
new file mode 100644
index 0000000000000000000000000000000000000000..fe65fa784251bf5070761b66459f363cb14d0688
--- /dev/null
+++ b/Tests/Data/RichardsMechanics/bishops_effective_stress_saturation_cutoff.prj
@@ -0,0 +1,335 @@
+<?xml version="1.0" encoding="ISO-8859-1"?>
+<OpenGeoSysProject>
+ <mesh>cube_1x1x1_hex20_1e0.vtu</mesh>
+ <geometry>cube_1x1x1.gml</geometry>
+ <processes>
+ <process>
+ <name>RM</name>
+ <type>RICHARDS_MECHANICS</type>
+ <integration_order>3</integration_order>
+ <dimension>3</dimension>
+ <!--
+ <jacobian_assembler>
+ <type>CentralDifferences</type>
+ <component_magnitudes>1e5 1e-5 1e-5 1e-5</component_magnitudes>
+ <relative_epsilons>1e-8 1e-8 1e-8 1e-8</relative_epsilons>
+ </jacobian_assembler>
+ -->
+ <constitutive_relation>
+ <type>LinearElasticIsotropic</type>
+ <youngs_modulus>E</youngs_modulus>
+ <poissons_ratio>nu</poissons_ratio>
+ </constitutive_relation>
+ <process_variables>
+ <displacement>displacement</displacement>
+ <pressure>pressure</pressure>
+ </process_variables>
+ <secondary_variables>
+ <secondary_variable internal_name="sigma" output_name="sigma"/>
+ <secondary_variable internal_name="epsilon" output_name="epsilon"/>
+ <secondary_variable internal_name="velocity" output_name="velocity"/>
+ <secondary_variable internal_name="saturation" output_name="saturation"/>
+ <secondary_variable internal_name="porosity" output_name="porosity"/>
+ </secondary_variables>
+ <specific_body_force>0 0 0</specific_body_force>
+ </process>
+ </processes>
+ <media>
+ <medium>
+ <phases>
+ <phase>
+ <type>AqueousLiquid</type>
+ <properties>
+ <property>
+ <name>bulk_modulus</name>
+ <type>Constant</type>
+ <value>1</value>
+ </property>
+ <property>
+ <name>viscosity</name>
+ <type>Constant</type>
+ <value>1</value>
+ </property>
+ <property>
+ <name>density</name>
+ <type>Constant</type>
+ <value>1</value>
+ </property>
+ </properties>
+ </phase>
+ <phase>
+ <type>Solid</type>
+ <properties>
+ <property>
+ <name>bulk_modulus</name>
+ <type>Constant</type>
+ <value>10</value>
+ </property>
+ <property>
+ <name>biot_coefficient</name>
+ <type>Constant</type>
+ <value>1</value>
+ </property>
+ <property>
+ <name>density</name>
+ <type>Constant</type>
+ <value>1</value>
+ </property>
+ <property>
+ <name>porosity</name>
+ <type>Constant</type>
+ <value>0.5</value>
+ </property>
+ <property>
+ <name>permeability</name>
+ <type>Constant</type>
+ <value>1</value>
+ </property>
+ </properties>
+ </phase>
+ </phases>
+ <properties>
+ <property>
+ <name>reference_temperature</name>
+ <type>Constant</type>
+ <value>293.15</value>
+ </property>
+ <property>
+ <name>saturation</name>
+ <type>SaturationVanGenuchten</type>
+ <residual_liquid_saturation>0.1</residual_liquid_saturation>
+ <residual_gas_saturation>0.05</residual_gas_saturation>
+ <exponent>0.8</exponent>
+ <entry_pressure>0.5</entry_pressure>
+ </property>
+ <property>
+ <name>relative_permeability</name>
+ <type>Constant</type>
+ <value>1</value>
+ </property>
+ <property>
+ <name>bishops_effective_stress</name>
+ <type>BishopsSaturationCutoff</type>
+ <cutoff_value>0.95</cutoff_value>
+ </property>
+ </properties>
+ </medium>
+ </media>
+ <time_loop>
+ <processes>
+ <process ref="RM">
+ <nonlinear_solver>nonlinear_solver</nonlinear_solver>
+ <compensate_non_equilibrium_initial_residuum>true</compensate_non_equilibrium_initial_residuum>
+ <convergence_criterion>
+ <type>PerComponentDeltaX</type>
+ <norm_type>NORM2</norm_type>
+ <abstols>1e-13 1e-13 1e-13 1e-13</abstols>
+ </convergence_criterion>
+ <time_discretization>
+ <type>BackwardEuler</type>
+ </time_discretization>
+ <time_stepping>
+ <type>FixedTimeStepping</type>
+ <t_initial>0</t_initial>
+ <t_end>3</t_end>
+ <timesteps>
+ <pair>
+ <repeat>50</repeat>
+ <delta_t>.02</delta_t>
+ </pair>
+
+ <pair>
+ <repeat>4</repeat>
+ <delta_t>.25</delta_t>
+ </pair>
+
+ <pair>
+ <repeat>50</repeat>
+ <delta_t>.02</delta_t>
+ </pair>
+ </timesteps>
+ </time_stepping>
+ </process>
+ </processes>
+ <output>
+ <type>VTK</type>
+ <prefix>bishops_effective_stress_saturation_cutoff</prefix>
+ <timesteps>
+ <pair>
+ <repeat>1000</repeat>
+ <each_steps>1</each_steps>
+ </pair>
+ </timesteps>
+ <variables>
+ <variable>displacement</variable>
+ <variable>pressure</variable>
+ <variable>sigma</variable>
+ <variable>epsilon</variable>
+ <variable>velocity</variable>
+ <variable>saturation</variable>
+ <variable>porosity</variable>
+ </variables>
+ </output>
+ </time_loop>
+ <local_coordinate_system>
+ <basis_vector_0>e0</basis_vector_0>
+ <basis_vector_1>e1</basis_vector_1>
+ <basis_vector_2>e2</basis_vector_2>
+ </local_coordinate_system>
+ <parameters>
+ <parameter>
+ <name>e0</name>
+ <type>Constant</type>
+ <values>1 0 0</values>
+ </parameter>
+ <parameter>
+ <name>e1</name>
+ <type>Constant</type>
+ <values>0 1 0</values>
+ </parameter>
+ <parameter>
+ <name>e2</name>
+ <type>Constant</type>
+ <values>0 0 1</values>
+ </parameter>
+ <!-- Mechanics -->
+ <parameter>
+ <name>E</name>
+ <type>Constant</type>
+ <value>1</value>
+ </parameter>
+ <parameter>
+ <name>nu</name>
+ <type>Constant</type>
+ <value>.3</value>
+ </parameter>
+ <!-- Model parameters -->
+ <parameter>
+ <name>displacement0</name>
+ <type>Constant</type>
+ <values>0 0 0</values>
+ </parameter>
+ <parameter>
+ <name>p0</name>
+ <type>Constant</type>
+ <value>0.1</value>
+ </parameter>
+ <parameter>
+ <name>dirichlet0</name>
+ <type>Constant</type>
+ <values>0</values>
+ </parameter>
+ <parameter>
+ <name>pressure_bc</name>
+ <type>CurveScaled</type>
+ <curve>pressure_ramp</curve>
+ <parameter>p0</parameter>
+ </parameter>
+ </parameters>
+ <curves>
+ <curve>
+ <name>pressure_ramp</name>
+ <coords>0 1 2 3</coords>
+ <values>1 -10 -10 1</values>
+ </curve>
+ </curves>
+ <process_variables>
+ <process_variable>
+ <name>displacement</name>
+ <components>3</components>
+ <order>2</order>
+ <initial_condition>displacement0</initial_condition>
+ <boundary_conditions>
+ <boundary_condition>
+ <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+ <geometry>left</geometry>
+ <type>Dirichlet</type>
+ <component>0</component>
+ <parameter>dirichlet0</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+ <geometry>front</geometry>
+ <type>Dirichlet</type>
+ <component>1</component>
+ <parameter>dirichlet0</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+ <geometry>bottom</geometry>
+ <type>Dirichlet</type>
+ <component>2</component>
+ <parameter>dirichlet0</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+ <geometry>top</geometry>
+ <type>Dirichlet</type>
+ <component>2</component>
+ <parameter>dirichlet0</parameter>
+ </boundary_condition>
+ </boundary_conditions>
+ </process_variable>
+ <process_variable>
+ <name>pressure</name>
+ <components>1</components>
+ <order>1</order>
+ <initial_condition>p0</initial_condition>
+ <boundary_conditions>
+ <boundary_condition>
+ <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+ <geometry>right</geometry>
+ <type>Dirichlet</type>
+ <parameter>pressure_bc</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+ <geometry>bottom</geometry>
+ <type>Dirichlet</type>
+ <parameter>pressure_bc</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+ <geometry>front</geometry>
+ <type>Dirichlet</type>
+ <parameter>pressure_bc</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+ <geometry>left</geometry>
+ <type>Dirichlet</type>
+ <parameter>pressure_bc</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+ <geometry>top</geometry>
+ <type>Dirichlet</type>
+ <parameter>pressure_bc</parameter>
+ </boundary_condition>
+ <boundary_condition>
+ <geometrical_set>cube_1x1x1_geometry</geometrical_set>
+ <geometry>back</geometry>
+ <type>Dirichlet</type>
+ <parameter>pressure_bc</parameter>
+ </boundary_condition>
+ </boundary_conditions>
+ </process_variable>
+ </process_variables>
+ <nonlinear_solvers>
+ <nonlinear_solver>
+ <name>nonlinear_solver</name>
+ <type>Newton</type>
+ <max_iter>100</max_iter>
+ <linear_solver>general_linear_solver</linear_solver>
+ </nonlinear_solver>
+ </nonlinear_solvers>
+ <linear_solvers>
+ <linear_solver>
+ <name>general_linear_solver</name>
+ <eigen>
+ <solver_type>SparseLU</solver_type>
+ <scaling>true</scaling>
+ </eigen>
+ </linear_solver>
+ </linear_solvers>
+</OpenGeoSysProject>
diff --git a/Tests/Data/RichardsMechanics/bishops_effective_stress_saturation_cutoff_pcs_0_ts_104_t_3.000000.vtu b/Tests/Data/RichardsMechanics/bishops_effective_stress_saturation_cutoff_pcs_0_ts_104_t_3.000000.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..42658f65522ffdd625f12045287b3b3e1c3a1352
--- /dev/null
+++ b/Tests/Data/RichardsMechanics/bishops_effective_stress_saturation_cutoff_pcs_0_ts_104_t_3.000000.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:4ea27efd0ddc782cf3679849c7bd64f229573166b0a977d00956db5fd1e93a0b
+size 8171
diff --git a/Tests/Data/RichardsMechanics/bishops_effective_stress_saturation_cutoff_pcs_0_ts_25_t_0.500000.vtu b/Tests/Data/RichardsMechanics/bishops_effective_stress_saturation_cutoff_pcs_0_ts_25_t_0.500000.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..06fc2256d7a78c760ded73dcfa799d8866acf107
--- /dev/null
+++ b/Tests/Data/RichardsMechanics/bishops_effective_stress_saturation_cutoff_pcs_0_ts_25_t_0.500000.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:e324c4251c2d820142b1cd9a35b2cf0186123e4c60eeed810f03c5b515519a3e
+size 7379
diff --git a/Tests/Data/RichardsMechanics/bishops_effective_stress_saturation_cutoff_pcs_0_ts_50_t_1.000000.vtu b/Tests/Data/RichardsMechanics/bishops_effective_stress_saturation_cutoff_pcs_0_ts_50_t_1.000000.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..785f5b776b55079e525bd786ad7d297f4a81a2b7
--- /dev/null
+++ b/Tests/Data/RichardsMechanics/bishops_effective_stress_saturation_cutoff_pcs_0_ts_50_t_1.000000.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:e0edbb427dd2793f32ebcebaf8995a20a4e6283da05e4ff02e81f9bedd4499c3
+size 7251
diff --git a/Tests/Data/RichardsMechanics/bishops_effective_stress_saturation_cutoff_pcs_0_ts_52_t_1.500000.vtu b/Tests/Data/RichardsMechanics/bishops_effective_stress_saturation_cutoff_pcs_0_ts_52_t_1.500000.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..bc2577ac6537d64072031cb11ef931a60a5030f8
--- /dev/null
+++ b/Tests/Data/RichardsMechanics/bishops_effective_stress_saturation_cutoff_pcs_0_ts_52_t_1.500000.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:719db305516fec867f2ad7220b2ee75b7d4ac390e173bbea0e7773b294b092e4
+size 7235
diff --git a/Tests/Data/RichardsMechanics/bishops_effective_stress_saturation_cutoff_pcs_0_ts_54_t_2.000000.vtu b/Tests/Data/RichardsMechanics/bishops_effective_stress_saturation_cutoff_pcs_0_ts_54_t_2.000000.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..cea5428a833297926ff15ae58618804ef7bd0da7
--- /dev/null
+++ b/Tests/Data/RichardsMechanics/bishops_effective_stress_saturation_cutoff_pcs_0_ts_54_t_2.000000.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:f1136f1ff8cf7c03b9763863b03c7907e753b4add265600285c1f68e2610af43
+size 7395
diff --git a/Tests/Data/RichardsMechanics/bishops_effective_stress_saturation_cutoff_pcs_0_ts_79_t_2.500000.vtu b/Tests/Data/RichardsMechanics/bishops_effective_stress_saturation_cutoff_pcs_0_ts_79_t_2.500000.vtu
new file mode 100644
index 0000000000000000000000000000000000000000..8f3ab459a4f9e13934d2afae910b66175348824a
--- /dev/null
+++ b/Tests/Data/RichardsMechanics/bishops_effective_stress_saturation_cutoff_pcs_0_ts_79_t_2.500000.vtu
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:ac0aab0bc0b91b2028da9ff6d239cfe0aeaa6a1b48e5125da1004804541f67b4
+size 7291
diff --git a/Tests/Data/RichardsMechanics/confined_compression_fully_saturated.prj b/Tests/Data/RichardsMechanics/confined_compression_fully_saturated.prj
index be056de940752f295eb1fa7890b067c518651e4d..4818ac1ff893633a7866da6597f5d484fa699d26 100644
--- a/Tests/Data/RichardsMechanics/confined_compression_fully_saturated.prj
+++ b/Tests/Data/RichardsMechanics/confined_compression_fully_saturated.prj
@@ -108,6 +108,11 @@
<type>Constant</type>
<value>1</value>
</property>
+ <property>
+ <name>bishops_effective_stress</name>
+ <type>BishopsPowerLaw</type>
+ <exponent>1</exponent>
+ </property>
</properties>
</medium>
</media>
diff --git a/Tests/Data/RichardsMechanics/deformation_dependent_porosity.prj b/Tests/Data/RichardsMechanics/deformation_dependent_porosity.prj
index 6810b7220c6de8d96bc6b06020c0886b1e7ea60b..dc3542198dbcfec4667d4399d2cadd2092fb9219 100644
--- a/Tests/Data/RichardsMechanics/deformation_dependent_porosity.prj
+++ b/Tests/Data/RichardsMechanics/deformation_dependent_porosity.prj
@@ -104,6 +104,11 @@
<type>Constant</type>
<value>1</value>
</property>
+ <property>
+ <name>bishops_effective_stress</name>
+ <type>BishopsPowerLaw</type>
+ <exponent>1</exponent>
+ </property>
</properties>
</medium>
</media>
diff --git a/Tests/Data/RichardsMechanics/flow_fully_saturated.prj b/Tests/Data/RichardsMechanics/flow_fully_saturated.prj
index 8cce040950ae2f5ee13e1d3d5568ca0d8b7192d9..73724958a5b1336d52d0c0e2bf4a5730a25a88ac 100644
--- a/Tests/Data/RichardsMechanics/flow_fully_saturated.prj
+++ b/Tests/Data/RichardsMechanics/flow_fully_saturated.prj
@@ -96,6 +96,11 @@
<type>Constant</type>
<value>1</value>
</property>
+ <property>
+ <name>bishops_effective_stress</name>
+ <type>BishopsPowerLaw</type>
+ <exponent>1</exponent>
+ </property>
<property>
<name>reference_temperature</name>
<type>Constant</type>
diff --git a/Tests/Data/RichardsMechanics/flow_fully_saturated_anisotropic.prj b/Tests/Data/RichardsMechanics/flow_fully_saturated_anisotropic.prj
index 11426a0fd66335691e0076e7b17e5a71e75acbc6..b706337ffdc058e5e37a21f98570d644f550d182 100644
--- a/Tests/Data/RichardsMechanics/flow_fully_saturated_anisotropic.prj
+++ b/Tests/Data/RichardsMechanics/flow_fully_saturated_anisotropic.prj
@@ -96,6 +96,11 @@
<type>Constant</type>
<value>1</value>
</property>
+ <property>
+ <name>bishops_effective_stress</name>
+ <type>BishopsPowerLaw</type>
+ <exponent>1</exponent>
+ </property>
<property>
<name>reference_temperature</name>
<type>Constant</type>
diff --git a/Tests/Data/RichardsMechanics/flow_fully_saturated_coordinate_system.prj b/Tests/Data/RichardsMechanics/flow_fully_saturated_coordinate_system.prj
index 7ce294e3e3927e8bee6026fcece0d4e4e18f5504..bf846106ffa157cbafbec73ddfd5b806ed0c0f58 100644
--- a/Tests/Data/RichardsMechanics/flow_fully_saturated_coordinate_system.prj
+++ b/Tests/Data/RichardsMechanics/flow_fully_saturated_coordinate_system.prj
@@ -96,6 +96,11 @@
<type>Constant</type>
<value>1</value>
</property>
+ <property>
+ <name>bishops_effective_stress</name>
+ <type>BishopsPowerLaw</type>
+ <exponent>1</exponent>
+ </property>
<property>
<name>reference_temperature</name>
<type>Constant</type>
diff --git a/Tests/Data/RichardsMechanics/gravity.prj b/Tests/Data/RichardsMechanics/gravity.prj
index 7394f8044d5cd84ff64c6eb1682f0e38c1dcc499..5b3c0037cf93b9ac656603ec94351bd6a53e782e 100644
--- a/Tests/Data/RichardsMechanics/gravity.prj
+++ b/Tests/Data/RichardsMechanics/gravity.prj
@@ -114,6 +114,11 @@
<exponent>0.789029535864979</exponent>
<minimum_relative_permeability_liquid>1e-12</minimum_relative_permeability_liquid>
</property>
+ <property>
+ <name>bishops_effective_stress</name>
+ <type>BishopsPowerLaw</type>
+ <exponent>1</exponent>
+ </property>
</properties>
</medium>
</media>
diff --git a/Tests/Data/RichardsMechanics/mechanics_linear.prj b/Tests/Data/RichardsMechanics/mechanics_linear.prj
index 3729849f21dbfc76283f2d506ac5651db6e121c6..0199ca2bceadfa8198e1d04ebd432a984802ce2a 100644
--- a/Tests/Data/RichardsMechanics/mechanics_linear.prj
+++ b/Tests/Data/RichardsMechanics/mechanics_linear.prj
@@ -109,6 +109,11 @@
<type>Constant</type>
<value>1</value>
</property>
+ <property>
+ <name>bishops_effective_stress</name>
+ <type>BishopsPowerLaw</type>
+ <exponent>1</exponent>
+ </property>
<property>
<name>reference_temperature</name>
<type>Constant</type>
diff --git a/Tests/Data/RichardsMechanics/orthotropic_power_law_permeability_xyz.prj b/Tests/Data/RichardsMechanics/orthotropic_power_law_permeability_xyz.prj
index 4f87c2eb8652cfa83cbba439bae0fd05ff4340d1..b131f794c34370117e6374dc82fa8f47f7397203 100644
--- a/Tests/Data/RichardsMechanics/orthotropic_power_law_permeability_xyz.prj
+++ b/Tests/Data/RichardsMechanics/orthotropic_power_law_permeability_xyz.prj
@@ -103,6 +103,11 @@
<type>Constant</type>
<value>1</value>
</property>
+ <property>
+ <name>bishops_effective_stress</name>
+ <type>BishopsPowerLaw</type>
+ <exponent>1</exponent>
+ </property>
</properties>
</medium>
</media>
diff --git a/Tests/Data/RichardsMechanics/orthotropic_swelling_xy.prj b/Tests/Data/RichardsMechanics/orthotropic_swelling_xy.prj
index 3f467b63cc863413c80dc6167eb1807418c7ecc1..30dbec1c3bab7d7815eb80530c2f08fba13e7c59 100644
--- a/Tests/Data/RichardsMechanics/orthotropic_swelling_xy.prj
+++ b/Tests/Data/RichardsMechanics/orthotropic_swelling_xy.prj
@@ -116,6 +116,11 @@
<type>Constant</type>
<value>1</value>
</property>
+ <property>
+ <name>bishops_effective_stress</name>
+ <type>BishopsPowerLaw</type>
+ <exponent>1</exponent>
+ </property>
</properties>
</medium>
</media>
diff --git a/Tests/Data/RichardsMechanics/orthotropic_swelling_xyz.prj b/Tests/Data/RichardsMechanics/orthotropic_swelling_xyz.prj
index 74bc46d475d17de10b675108394f6e5f74ae8d8a..04867c4c4dd71725700910ce1ab6eb8b9202714c 100644
--- a/Tests/Data/RichardsMechanics/orthotropic_swelling_xyz.prj
+++ b/Tests/Data/RichardsMechanics/orthotropic_swelling_xyz.prj
@@ -116,6 +116,11 @@
<type>Constant</type>
<value>1</value>
</property>
+ <property>
+ <name>bishops_effective_stress</name>
+ <type>BishopsPowerLaw</type>
+ <exponent>1</exponent>
+ </property>
</properties>
</medium>
</media>
diff --git a/web/config.toml b/web/config.toml
index 6b6d6a507d669b596ca9c6053fd02d6f73aa8a3d..a4d459f2bfe955ee3ee186e75b3c9240583f2b7c 100644
--- a/web/config.toml
+++ b/web/config.toml
@@ -76,35 +76,39 @@ disableKinds = ["taxonomyTerm"]
[[menu.benchmarks]]
name = "Heat Transport BHE"
identifier = "heat-transport-bhe"
- weight = 11
+ weight = 12
[[menu.benchmarks]]
name = "Richards Flow"
identifier = "richards-flow"
weight = 4
+[[menu.benchmarks]]
+ name = "Richards-Mechanics"
+ identifier = "richards-mechanics"
+ weight = 5
[[menu.benchmarks]]
name = "Hydro-Mechanics"
identifier = "hydro-mechanics"
- weight = 5
+ weight = 6
[[menu.benchmarks]]
name = "Thermo-Hydro-Mechanics"
identifier = "thermo-hydro-mechanics"
- weight = 7
+ weight = 8
[[menu.benchmarks]]
name = "Hydro-Thermal"
identifier = "hydro-thermal"
- weight = 6
+ weight = 7
[[menu.benchmarks]]
name = "Liquid Flow"
identifier = "liquid-flow"
- weight = 8
+ weight = 9
[[menu.benchmarks]]
name = "Two-phase Flow"
identifier = "two-phase-flow"
- weight = 9
+ weight = 10
[[menu.benchmarks]]
name = "Python Boundary Conditions"
identifier = "python-bc"
- weight = 10
+ weight = 11
# Quickstart sidebar top-level categories
[[menu.userguide]]
diff --git a/web/content/docs/benchmarks/richards-mechanics/BishopsEffectiveStress.png b/web/content/docs/benchmarks/richards-mechanics/BishopsEffectiveStress.png
new file mode 100644
index 0000000000000000000000000000000000000000..13fe53eef21022c0ea63e399ccc67bdbb8033f01
--- /dev/null
+++ b/web/content/docs/benchmarks/richards-mechanics/BishopsEffectiveStress.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:891fc59993dfba7610a73b4f49f8ce310419a679711c427fa1dc1d5783bc05fe
+size 65602
diff --git a/web/content/docs/benchmarks/richards-mechanics/bishops-effective-stress.pandoc b/web/content/docs/benchmarks/richards-mechanics/bishops-effective-stress.pandoc
new file mode 100644
index 0000000000000000000000000000000000000000..b35e15f13ad572add18e702c6e354f9d02ea696c
--- /dev/null
+++ b/web/content/docs/benchmarks/richards-mechanics/bishops-effective-stress.pandoc
@@ -0,0 +1,33 @@
++++
+project = "RichardsMechanics/bishops_effective_stress_power_law.prj"
+author = "Dmitri Naumov"
+date = "2020-02-27"
+title = "Bishop's effective stress models comparison"
+weight = 153
+
+[menu]
+ [menu.benchmarks]
+ parent = "richards-mechanics"
+
++++
+
+{{< data-link >}}
+
+Two models for the Bishop's effective stress computation are presented; the
+power-law model, and saturation cut-off model. The models are:
+$$
+\chi(S_\mathrm{L}) = S_\mathrm{L}^{m_\chi}
+\qquad \mbox{and}\qquad
+\chi(S_\mathrm{L}) =
+ \chi = \begin{cases}
+ 1 & \mbox{for $S_\text{L} \geq S_\text{cutoff}$}
+ \\
+ 0 & \mbox{for $S_\text{L} < S_\text{cutoff}$.}
+ \end{cases}
+$$
+Simulation result shows different influence of the effective stress on the
+displacement. In the test the medium is desaturated and then saturated again,
+which causes shrinkage and expansion of the domain. Power law with exponents 1,
+1/5, and 5 and saturation cut-off at maximum liquid saturation of 0.95 are
+compared.
+{{< img src="../BishopsEffectiveStress.png" >}}