From 2046c6907e9f84588fdcaf2915e73891a67ff60f Mon Sep 17 00:00:00 2001
From: Norihiro Watanabe <norihiro.watanabe@ufz.de>
Date: Fri, 23 Sep 2016 14:23:28 +0200
Subject: [PATCH] add MohrCoulomb model for fractures
---
.../FractureModels/CreateMohrCoulomb.cpp | 65 ++++++++++
.../FractureModels/CreateMohrCoulomb.h | 28 +++++
MaterialLib/FractureModels/MohrCoulomb.cpp | 115 ++++++++++++++++++
MaterialLib/FractureModels/MohrCoulomb.h | 86 +++++++++++++
4 files changed, 294 insertions(+)
create mode 100644 MaterialLib/FractureModels/CreateMohrCoulomb.cpp
create mode 100644 MaterialLib/FractureModels/CreateMohrCoulomb.h
create mode 100644 MaterialLib/FractureModels/MohrCoulomb.cpp
create mode 100644 MaterialLib/FractureModels/MohrCoulomb.h
diff --git a/MaterialLib/FractureModels/CreateMohrCoulomb.cpp b/MaterialLib/FractureModels/CreateMohrCoulomb.cpp
new file mode 100644
index 00000000000..84342c90089
--- /dev/null
+++ b/MaterialLib/FractureModels/CreateMohrCoulomb.cpp
@@ -0,0 +1,65 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, 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 "CreateMohrCoulomb.h"
+
+#include "MohrCoulomb.h"
+
+namespace MaterialLib
+{
+namespace Fracture
+{
+
+template <int DisplacementDim>
+std::unique_ptr<FractureModelBase<DisplacementDim>>
+createMohrCoulomb(
+ std::vector<std::unique_ptr<ProcessLib::ParameterBase>> const& parameters,
+ BaseLib::ConfigTree const& config)
+{
+ config.checkConfigParameter("type", "MohrCoulomb");
+ DBUG("Create MohrCoulomb material");
+
+ auto& Kn = ProcessLib::findParameter<double>(
+ config, "kn", parameters, 1);
+
+ auto& Ks = ProcessLib::findParameter<double>(
+ config, "ks", parameters, 1);
+
+ auto& friction_angle = ProcessLib::findParameter<double>(
+ config, "friction_angle", parameters, 1);
+
+ auto& dilatancy_angle = ProcessLib::findParameter<double>(
+ config, "dilatancy_angle", parameters, 1);
+
+ auto& cohesion = ProcessLib::findParameter<double>(
+ config, "cohesion", parameters, 1);
+
+
+ typename MohrCoulomb<DisplacementDim>::MaterialProperties mp{
+ Kn, Ks, friction_angle, dilatancy_angle, cohesion};
+
+ return std::unique_ptr<MohrCoulomb<DisplacementDim>>{
+ new MohrCoulomb<DisplacementDim>{mp}};
+}
+
+
+template
+std::unique_ptr<FractureModelBase<2>>
+createMohrCoulomb(
+ std::vector<std::unique_ptr<ProcessLib::ParameterBase>> const& parameters,
+ BaseLib::ConfigTree const& config);
+
+template
+std::unique_ptr<FractureModelBase<3>>
+createMohrCoulomb(
+ std::vector<std::unique_ptr<ProcessLib::ParameterBase>> const& parameters,
+ BaseLib::ConfigTree const& config);
+
+} // namespace Fracture
+} // namespace MaterialLib
diff --git a/MaterialLib/FractureModels/CreateMohrCoulomb.h b/MaterialLib/FractureModels/CreateMohrCoulomb.h
new file mode 100644
index 00000000000..3ee7c7da854
--- /dev/null
+++ b/MaterialLib/FractureModels/CreateMohrCoulomb.h
@@ -0,0 +1,28 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, 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 "ProcessLib/Utils/ProcessUtils.h" // required for findParameter
+#include "FractureModelBase.h"
+
+namespace MaterialLib
+{
+namespace Fracture
+{
+
+template <int DisplacementDim>
+std::unique_ptr<FractureModelBase<DisplacementDim>>
+createMohrCoulomb(
+ std::vector<std::unique_ptr<ProcessLib::ParameterBase>> const& parameters,
+ BaseLib::ConfigTree const& config);
+
+} // namespace Fracture
+} // namespace MaterialLib
+
diff --git a/MaterialLib/FractureModels/MohrCoulomb.cpp b/MaterialLib/FractureModels/MohrCoulomb.cpp
new file mode 100644
index 00000000000..8cfddaaf8ca
--- /dev/null
+++ b/MaterialLib/FractureModels/MohrCoulomb.cpp
@@ -0,0 +1,115 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, 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 "MohrCoulomb.h"
+
+#include "BaseLib/Error.h"
+#include "MathLib/MathTools.h"
+
+namespace MaterialLib
+{
+namespace Fracture
+{
+
+namespace
+{
+
+struct MaterialPropertyValues
+{
+ double Kn = 0.0;
+ double Ks = 0.0;
+ double phi = 0.0; // friction angle
+ double psi = 0.0; // dilation angle
+ double c = 0.0;
+
+ template <typename MaterialProperties>
+ MaterialPropertyValues(
+ MaterialProperties const& mp,
+ double const t,
+ ProcessLib::SpatialPosition const& x)
+ {
+ Kn = mp.normal_stiffness(t,x)[0];
+ Ks = mp.shear_stiffness(t,x)[0];
+ phi = MathLib::to_radians(mp.friction_angle(t,x)[0]);
+ psi = MathLib::to_radians(mp.dilatancy_angle(t,x)[0]);
+ c = mp.cohesion(t,x)[0];
+ }
+};
+
+} // no namespace
+
+template <int DisplacementDim>
+void MohrCoulomb<DisplacementDim>::computeConstitutiveRelation(
+ double const t,
+ ProcessLib::SpatialPosition const& x,
+ Eigen::Ref<Eigen::VectorXd const> w_prev,
+ Eigen::Ref<Eigen::VectorXd const> w,
+ Eigen::Ref<Eigen::VectorXd const> sigma_prev,
+ Eigen::Ref<Eigen::VectorXd> sigma,
+ Eigen::Ref<Eigen::MatrixXd> Kep)
+{
+ if (DisplacementDim == 3)
+ {
+ OGS_FATAL("MohrCoulomb fracture model does not support 3D case.");
+ return;
+ }
+ MaterialPropertyValues const mat(_mp, t, x);
+ Eigen::VectorXd const dw = w - w_prev;
+
+ Eigen::MatrixXd Ke(2,2);
+ Ke.setZero();
+ Ke(0,0) = mat.Ks;
+ Ke(1,1) = mat.Kn;
+
+ sigma.noalias() = sigma_prev + Ke * dw;
+
+ // if opening
+ if (sigma[1] > 0)
+ {
+ Kep.setZero();
+ sigma.setZero();
+ return;
+ }
+
+ // check shear yield function (Fs)
+ double const Fs = std::abs(sigma[0]) + sigma[1] * tan(mat.phi) - mat.c;
+ if (Fs < .0)
+ {
+ Kep = Ke;
+ return;
+ }
+
+ Eigen::VectorXd dFs_dS(2);
+ dFs_dS[0] = MathLib::sgn(sigma[0]);
+ dFs_dS[1] = tan(mat.phi);
+
+ // plastic potential function: Qs = |tau| + Sn * tan da
+ Eigen::VectorXd dQs_dS(2);
+ dQs_dS[0] = MathLib::sgn(sigma[0]);
+ dQs_dS[1] = tan(mat.psi);
+
+ // plastic multiplier
+ Eigen::RowVectorXd const A = dFs_dS.transpose() * Ke / (dFs_dS.transpose() * Ke * dQs_dS);
+ double const d_eta = A * dw;
+
+ // plastic part of the dispalcement
+ Eigen::VectorXd const dwp = dQs_dS * d_eta;
+
+ // correct stress
+ sigma.noalias() = sigma_prev + Ke * (dw - dwp);
+
+ // Kep
+ Kep = Ke - Ke * dQs_dS * A;
+}
+
+template class MohrCoulomb<2>;
+template class MohrCoulomb<3>;
+
+} // namespace Fracture
+} // namespace MaterialLib
diff --git a/MaterialLib/FractureModels/MohrCoulomb.h b/MaterialLib/FractureModels/MohrCoulomb.h
new file mode 100644
index 00000000000..40eb458aa2f
--- /dev/null
+++ b/MaterialLib/FractureModels/MohrCoulomb.h
@@ -0,0 +1,86 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2016, 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 <Eigen/Eigen>
+
+#include "ProcessLib/Parameter/Parameter.h"
+
+#include "FractureModelBase.h"
+
+namespace MaterialLib
+{
+namespace Fracture
+{
+
+/**
+ * Mohr-Coulomb fracture model
+ */
+template <int DisplacementDim>
+class MohrCoulomb final : public FractureModelBase<DisplacementDim>
+{
+public:
+ /// Variables specific to the material model
+ struct MaterialProperties
+ {
+ using P = ProcessLib::Parameter<double>;
+ using X = ProcessLib::SpatialPosition;
+
+ MaterialProperties(
+ P const& normal_stiffness, P const& shear_stiffness,
+ P const& friction_angle, P const& dilatancy_angle,
+ P const& cohesion)
+ : normal_stiffness(normal_stiffness), shear_stiffness(shear_stiffness),
+ friction_angle(friction_angle), dilatancy_angle(dilatancy_angle),
+ cohesion(cohesion)
+ {
+ }
+
+ P const& normal_stiffness;
+ P const& shear_stiffness;
+ P const& friction_angle;
+ P const& dilatancy_angle;
+ P const& cohesion;
+ };
+
+public:
+
+ explicit MohrCoulomb(
+ MaterialProperties const& material_properties)
+ : _mp(material_properties)
+ {
+ }
+
+ void computeConstitutiveRelation(
+ double const t,
+ ProcessLib::SpatialPosition const& x,
+ Eigen::Ref<Eigen::VectorXd const> w_prev,
+ Eigen::Ref<Eigen::VectorXd const> w,
+ Eigen::Ref<Eigen::VectorXd const> sigma_prev,
+ Eigen::Ref<Eigen::VectorXd> sigma,
+ Eigen::Ref<Eigen::MatrixXd> Kep) override;
+
+private:
+
+ MaterialProperties _mp;
+};
+
+} // namespace Fracture
+} // namespace MaterialLib
+
+namespace MaterialLib
+{
+namespace Fracture
+{
+extern template class MohrCoulomb<2>;
+extern template class MohrCoulomb<3>;
+} // namespace Fractrue
+} // namespace MaterialLib
+
--
GitLab