Merge pull request #1434 from norihiro-w/fracture-model

Material model for fracture deformation

MaterialLib/CMakeLists.txt

@@ -2,6 +2,7 @@
 get_source_files(SOURCES)
 append_source_files(SOURCES Adsorption)
 append_source_files(SOURCES SolidModels)
+append_source_files(SOURCES FractureModels)
 
 append_source_files(SOURCES Fluid)
 append_source_files(SOURCES Fluid/Density)

MaterialLib/FractureModels/CreateLinearElasticIsotropic.cpp

+/**
+ * \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 "CreateLinearElasticIsotropic.h"
+
+#include "ProcessLib/Utils/ProcessUtils.h"  // required for findParameter
+#include "LinearElasticIsotropic.h"
+
+namespace MaterialLib
+{
+namespace Fracture
+{
+
+template <int DisplacementDim>
+std::unique_ptr<FractureModelBase<DisplacementDim>>
+createLinearElasticIsotropic(
+    std::vector<std::unique_ptr<ProcessLib::ParameterBase>> const& parameters,
+    BaseLib::ConfigTree const& config)
+{
+    config.checkConfigParameter("type", "LinearElasticIsotropic");
+    DBUG("Create LinearElasticIsotropic material");
+
+    auto& Kn = ProcessLib::findParameter<double>(
+        config, "normal_stiffness", parameters, 1);
+
+    auto& Ks = ProcessLib::findParameter<double>(
+        config, "shear_stiffness", parameters, 1);
+
+    typename LinearElasticIsotropic<DisplacementDim>::MaterialProperties mp{
+        Kn, Ks};
+
+    return std::unique_ptr<LinearElasticIsotropic<DisplacementDim>>{
+        new LinearElasticIsotropic<DisplacementDim>{mp}};
+}
+
+
+template
+std::unique_ptr<FractureModelBase<2>>
+createLinearElasticIsotropic(
+    std::vector<std::unique_ptr<ProcessLib::ParameterBase>> const& parameters,
+    BaseLib::ConfigTree const& config);
+
+template
+std::unique_ptr<FractureModelBase<3>>
+createLinearElasticIsotropic(
+    std::vector<std::unique_ptr<ProcessLib::ParameterBase>> const& parameters,
+    BaseLib::ConfigTree const& config);
+
+}  // namespace Fracture
+}  // namespace MaterialLib

MaterialLib/FractureModels/CreateLinearElasticIsotropic.h

+/**
+ * \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
+ *
+ */
+
+#ifndef MATERIALLIB_FRACTURE_CREATELINEARELASTICISOTROPIC_H_
+#define MATERIALLIB_FRACTURE_CREATELINEARELASTICISOTROPIC_H_
+
+#include "FractureModelBase.h"
+
+namespace  MaterialLib
+{
+namespace Fracture
+{
+
+template <int DisplacementDim>
+std::unique_ptr<FractureModelBase<DisplacementDim>>
+createLinearElasticIsotropic(
+    std::vector<std::unique_ptr<ProcessLib::ParameterBase>> const& parameters,
+    BaseLib::ConfigTree const& config);
+
+}  // namespace Fracture
+}  // namespace MaterialLib
+
+#endif

MaterialLib/FractureModels/CreateMohrCoulomb.cpp

+/**
+ * \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 "ProcessLib/Utils/ProcessUtils.h"  // required for findParameter
+#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, "normal_stiffness", parameters, 1);
+
+    auto& Ks = ProcessLib::findParameter<double>(
+        config, "shear_stiffness", 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

MaterialLib/FractureModels/CreateMohrCoulomb.h

+/**
+ * \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
+ *
+ */
+
+#ifndef MATERIALLIB_FRACTURE_CREATEMOHRCOULOMB_H_
+#define MATERIALLIB_FRACTURE_CREATEMOHRCOULOMB_H_
+
+#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
+
+#endif

MaterialLib/FractureModels/FractureModelBase.h

+/**
+ * \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
+ *
+ */
+
+#ifndef MATERIALLIB_FRACTURE_FRACTUREMODELBASE_H_
+#define MATERIALLIB_FRACTURE_FRACTUREMODELBASE_H_
+
+#include <Eigen/Eigen>
+
+#include "ProcessLib/Parameter/Parameter.h"
+
+namespace MaterialLib
+{
+namespace Fracture
+{
+
+/**
+ * Interface for mechanical fracture material models. Provides updates of the
+ * stress for a given current state and also a tangent at that position.
+ */
+template <int DisplacementDim>
+class FractureModelBase
+{
+public:
+    virtual ~FractureModelBase() {}
+
+    /**
+     * Computation of the constitutive relation for specific material model.
+     * This should be implemented in the derived model.
+     *
+     * @param t           current time
+     * @param x           current position in space
+     * @param w_prev      fracture displacement at previous time step
+     * @param w           fracture displacement at current time step
+     * @param sigma_prev  stress at previous time step
+     * @param sigma       stress at current time step
+     * @param C           tangent matrix for stress and fracture displacements
+     */
+    virtual 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> C)  = 0;
+
+};
+
+}  // namespace Fracture
+}  // namespace MaterialLib
+
+#endif

MaterialLib/FractureModels/LinearElasticIsotropic.cpp

+/**
+ * \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 "LinearElasticIsotropic.h"
+
+namespace MaterialLib
+{
+namespace Fracture
+{
+
+template <int DisplacementDim>
+void LinearElasticIsotropic<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> C)
+{
+    const int index_ns = DisplacementDim - 1;
+    C.setZero();
+    for (int i=0; i<index_ns; i++)
+        C(i,i) = _mp.shear_stiffness(t, x)[0];
+    C(index_ns, index_ns) = _mp.normal_stiffness(t, x)[0];
+
+    sigma.noalias() = sigma_prev + C * (w - w_prev);
+
+    // correct if fracture is opening
+    if (sigma[index_ns] > 0)
+    {
+        C.setZero();
+        sigma.setZero();
+    }
+}
+
+template class LinearElasticIsotropic<2>;
+template class LinearElasticIsotropic<3>;
+
+
+}   // namespace Fracture
+}   // namespace MaterialLib

MaterialLib/FractureModels/LinearElasticIsotropic.h

+/**
+ * \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
+ *
+ */
+
+#ifndef MATERIALLIB_FRACTURE_LINEARELASTICISOTROPIC_H_
+#define MATERIALLIB_FRACTURE_LINEARELASTICISOTROPIC_H_
+
+#include <Eigen/Eigen>
+
+#include "ProcessLib/Parameter/Parameter.h"
+
+#include "FractureModelBase.h"
+
+namespace MaterialLib
+{
+namespace Fracture
+{
+
+template <int DisplacementDim>
+class LinearElasticIsotropic 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_)
+            : normal_stiffness(normal_stiffness_), shear_stiffness(shear_stiffness_)
+        {
+        }
+
+        P const& normal_stiffness;
+        P const& shear_stiffness;
+    };
+
+public:
+    explicit LinearElasticIsotropic(
+        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> C) override;
+
+private:
+    MaterialProperties _mp;
+};
+
+}  // namespace Fracture
+}  // namespace MaterialLib
+
+namespace MaterialLib
+{
+namespace Fracture
+{
+extern template class LinearElasticIsotropic<2>;
+extern template class LinearElasticIsotropic<3>;
+}  // namespace Fractrue
+}  // namespace MaterialLib
+
+#endif

MaterialLib/FractureModels/MohrCoulomb.cpp

+/**
+ * \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"
+
+#ifndef Q_MOC_RUN // to avoid Qt4 bug, https://bugreports.qt.io/browse/QTBUG-22829
+#include <boost/math/special_functions/sign.hpp>
+#endif
+
+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] * std::tan(mat.phi) - mat.c;
+    if (Fs < .0)
+    {
+        Kep = Ke;
+        return;
+    }
+
+    Eigen::VectorXd dFs_dS(2);
+    dFs_dS[0] = boost::math::sign(sigma[0]);
+    dFs_dS[1] = std::tan(mat.phi);
+
+    // plastic potential function: Qs = |tau| + Sn * tan da
+    Eigen::VectorXd dQs_dS(2);
+    dQs_dS[0] = boost::math::sign(sigma[0]);
+    dQs_dS[1] = std::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

MaterialLib/FractureModels/MohrCoulomb.h

+/**
+ * \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
+ *
+ */
+
+#ifndef MATERIALLIB_FRACTURE_MOHRCOULOMB_H_
+#define MATERIALLIB_FRACTURE_MOHRCOULOMB_H_
+
+#include <Eigen/Eigen>
+
+#include "ProcessLib/Parameter/Parameter.h"
+
+#include "FractureModelBase.h"
+
+namespace MaterialLib
+{
+namespace Fracture
+{
+
+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
+
+#endif

MathLib/MathTools.h

@@ -17,6 +17,10 @@
 #include <omp.h>
 #endif
 
+#ifndef Q_MOC_RUN // to avoid Qt4 bug, https://bugreports.qt.io/browse/QTBUG-22829
+#include <boost/math/constants/constants.hpp>
+#endif
+
 namespace MathLib
 {
 /**
@@ -111,6 +115,11 @@ double sqrDist(const double* p0, const double* p1)
  */
 double getAngle (const double p0[3], const double p1[3], const double p2[3]);
 
+/// converts the given degrees to radians
+inline double to_radians(double degrees) {
+    return degrees*boost::math::constants::pi<double>()/180.;
+}
+
 } // namespace
 
 #endif /* MATHTOOLS_H_ */

Tests/MaterialLib/TestFractureModels.cpp

+/**
+ * \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/LICENSE.txt
+ */
+
+
+#include <limits>
+
+#include <gtest/gtest.h>
+
+#include <Eigen/Eigen>
+
+#include "BaseLib/ConfigTree.h"
+
+#include "MaterialLib/FractureModels/LinearElasticIsotropic.h"
+#include "MaterialLib/FractureModels/MohrCoulomb.h"
+
+#include "ProcessLib/Parameter/ConstantParameter.h"
+
+using namespace MaterialLib::Fracture;
+
+static const double eps_sigma = 1e6*1e-5;
+static const double eps_C = 1e10*1e-5;
+
+TEST(MaterialLib_Fracture, LinearElasticIsotropic)
+{
+    ProcessLib::ConstantParameter<double> const kn(1e11);
+    ProcessLib::ConstantParameter<double> const ks(1e9);
+    LinearElasticIsotropic<2>::MaterialProperties const mp{kn, ks};
+
+    LinearElasticIsotropic<2> fractureModel{mp};
+
+    Eigen::Vector2d w_prev, w, sigma_prev, sigma;
+    Eigen::Matrix2d C;
+
+    ProcessLib::SpatialPosition x;
+    w << -1e-5, -1e-5;
+    fractureModel.computeConstitutiveRelation(0, x, w_prev, w, sigma_prev, sigma, C);
+
+    ASSERT_NEAR(-1e4, sigma[0], eps_sigma);
+    ASSERT_NEAR(-1e6, sigma[1], eps_sigma);
+    ASSERT_NEAR(1e9, C(0,0), eps_C);
+    ASSERT_NEAR(0, C(0,1), eps_C);
+    ASSERT_NEAR(0, C(1,0), eps_C);
+    ASSERT_NEAR(1e11, C(1,1), eps_C);
+
+
+    w << -1e-5, 1e-5;
+    fractureModel.computeConstitutiveRelation(0, x, w_prev, w, sigma_prev, sigma, C);
+
+    ASSERT_NEAR(0, sigma[0], eps_sigma);
+    ASSERT_NEAR(0, sigma[1], eps_sigma);
+    ASSERT_NEAR(0, C(0,0), eps_C);
+    ASSERT_NEAR(0, C(0,1), eps_C);
+    ASSERT_NEAR(0, C(1,0), eps_C);
+    ASSERT_NEAR(0, C(1,1), eps_C);
+}
+
+TEST(MaterialLib_Fracture, MohrCoulomb)
+{
+    ProcessLib::ConstantParameter<double> const kn(50e9);
+    ProcessLib::ConstantParameter<double> const ks(20e9);
+    ProcessLib::ConstantParameter<double> const phi(15);
+    ProcessLib::ConstantParameter<double> const psi(5);
+    ProcessLib::ConstantParameter<double> const c(3e6);
+    MohrCoulomb<2>::MaterialProperties const mp{kn, ks, phi, psi, c};
+
+    MohrCoulomb<2> fractureModel{mp};
+
+    Eigen::Vector2d w_prev, w, sigma_prev, sigma;
+    Eigen::Matrix2d C;
+
+    ProcessLib::SpatialPosition x;
+    sigma_prev << -3.46e6, -2e6;
+    w << -1.08e-5, -0.25e-5;
+    fractureModel.computeConstitutiveRelation(0, x, w_prev, w, sigma_prev, sigma, C);
+
+    ASSERT_NEAR(-3.50360e6, sigma[0], eps_sigma);
+    ASSERT_NEAR(-2.16271e6, sigma[1], eps_sigma);
+    ASSERT_NEAR(1.10723e+09, C(0,0), eps_C);
+    ASSERT_NEAR(1.26558e+10, C(0,1), eps_C);
+    ASSERT_NEAR(4.13226e+09, C(1,0), eps_C);
+    ASSERT_NEAR(4.72319e+10, C(1,1), eps_C);
+}
+