From 83a9820fe6f591573f86fc30f9821ca95980b276 Mon Sep 17 00:00:00 2001
From: Dmitri Naumov <dmitri.naumov@ufz.de>
Date: Thu, 18 Jan 2018 10:51:14 +0100
Subject: [PATCH] [MatL] LIE/Fracture. Cohesive zone mode I model.
---
.../FractureModels/CohesiveZoneModeI.cpp | 132 ++++++++++++
.../FractureModels/CohesiveZoneModeI.h | 193 ++++++++++++++++++
.../CreateCohesiveZoneModeI.cpp | 70 +++++++
.../FractureModels/CreateCohesiveZoneModeI.h | 27 +++
4 files changed, 422 insertions(+)
create mode 100644 MaterialLib/FractureModels/CohesiveZoneModeI.cpp
create mode 100644 MaterialLib/FractureModels/CohesiveZoneModeI.h
create mode 100644 MaterialLib/FractureModels/CreateCohesiveZoneModeI.cpp
create mode 100644 MaterialLib/FractureModels/CreateCohesiveZoneModeI.h
diff --git a/MaterialLib/FractureModels/CohesiveZoneModeI.cpp b/MaterialLib/FractureModels/CohesiveZoneModeI.cpp
new file mode 100644
index 00000000000..91189306b08
--- /dev/null
+++ b/MaterialLib/FractureModels/CohesiveZoneModeI.cpp
@@ -0,0 +1,132 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2018, 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 <iostream>
+
+#include "CohesiveZoneModeI.h"
+#include "LogPenalty.h"
+
+#include "BaseLib/Error.h"
+#include "MathLib/MathTools.h"
+
+namespace MaterialLib
+{
+namespace Fracture
+{
+namespace CohesiveZoneModeI
+{
+namespace
+{
+double computeDamage(double const damage_prev,
+ double const w_n,
+ double const w_np,
+ double const w_nf)
+{
+ return std::min(
+ 1.0,
+ std::max(damage_prev, std::max(0.0, (w_n - w_np)) / (w_nf - w_np)));
+}
+
+} // namespace
+
+template <int DisplacementDim>
+void CohesiveZoneModeI<DisplacementDim>::computeConstitutiveRelation(
+ double const t,
+ ProcessLib::SpatialPosition const& x,
+ double const aperture0,
+ Eigen::Ref<Eigen::VectorXd const>
+ sigma0,
+ 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,
+ typename FractureModelBase<DisplacementDim>::MaterialStateVariables&
+ material_state_variables)
+{
+ assert(dynamic_cast<StateVariables<DisplacementDim> const*>(
+ &material_state_variables) != nullptr);
+
+ StateVariables<DisplacementDim>& state =
+ static_cast<StateVariables<DisplacementDim> &>(
+ material_state_variables);
+ //reset damage in each iteration
+ state.setInitialConditions();
+
+ auto const mp = evaluatedMaterialProperties(t, x);
+
+ C.setZero();
+
+ // Separately compute shear and normal stresses because of the penalty for
+ // the normal component.
+ const int index_ns = DisplacementDim - 1;
+ double const w_n = w[index_ns];
+ for (int i = 0; i < index_ns; i++)
+ C(i, i) = mp.Ks;
+
+ sigma.noalias() = C * w;
+
+ double const aperture = w_n + aperture0;
+
+ sigma.coeffRef(index_ns) =
+ mp.Kn * w_n * logPenalty(aperture0, aperture, _penalty_aperture_cutoff);
+
+ C(index_ns, index_ns) =
+ mp.Kn *
+ logPenaltyDerivative(aperture0, aperture, _penalty_aperture_cutoff);
+
+ // Exit if fracture is closing
+ // TODO (nagel) to be based on the stress state when initial stress effects
+ // are included.
+ if (w_n < 0)
+ {
+ return; /// Undamaged stiffness used in compression.
+ }
+
+ //
+ // Continue with fracture opening.
+ //
+
+ state.damage = computeDamage(state.damage_prev, w_n, mp.w_np, mp.w_nf);
+ const double degradation = ((1 - state.damage) * mp.w_np) /
+ (mp.w_np + state.damage * (mp.w_nf - mp.w_np));
+
+ // Degrade stiffness tensor in tension.
+ C = C * degradation;
+ sigma = C * w;
+
+ if (state.damage > state.damage_prev)
+ {
+ // If damage is increasing, provide extension to consistent tangent.
+
+ Eigen::Matrix<double, DisplacementDim, 1> dd_dw =
+ Eigen::Matrix<double, DisplacementDim, 1>::Zero();
+
+ double const tmp = mp.w_np + state.damage * (mp.w_nf - mp.w_np);
+ dd_dw[index_ns] =
+ (mp.w_np * mp.w_nf) / ((mp.w_nf - mp.w_np) * (tmp * tmp));
+
+ C -= C * w * (dd_dw).transpose();
+ }
+
+ // TODO (nagel) Initial stress not considered, yet.
+ // sigma.noalias() += sigma0;
+}
+
+template class CohesiveZoneModeI<2>;
+template class CohesiveZoneModeI<3>;
+
+} // namespace CohesiveZoneModeI
+} // namespace Fracture
+} // namespace MaterialLib
diff --git a/MaterialLib/FractureModels/CohesiveZoneModeI.h b/MaterialLib/FractureModels/CohesiveZoneModeI.h
new file mode 100644
index 00000000000..8bbfc83acae
--- /dev/null
+++ b/MaterialLib/FractureModels/CohesiveZoneModeI.h
@@ -0,0 +1,193 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2018, 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 <utility>
+
+#include "ProcessLib/Parameter/Parameter.h"
+
+#include "FractureModelBase.h"
+
+namespace MaterialLib
+{
+namespace Fracture
+{
+namespace CohesiveZoneModeI
+{
+/// Variables specific to the material model
+struct MaterialPropertiesParameters
+{
+ using P = ProcessLib::Parameter<double>;
+ using X = ProcessLib::SpatialPosition;
+
+ MaterialPropertiesParameters(P const& normal_stiffness_,
+ P const& shear_stiffness_,
+ P const& fracture_toughness_,
+ P const& peak_normal_traction_)
+ : normal_stiffness(normal_stiffness_),
+ shear_stiffness(shear_stiffness_),
+ fracture_toughness(fracture_toughness_),
+ peak_normal_traction(peak_normal_traction_)
+ {
+ }
+
+ /// Assuming initially stress-free state.
+ double fracture_opening_at_peak_traction(double const t, X const& x) const
+ {
+ return peak_normal_traction(t, x)[0] / normal_stiffness(t, x)[0];
+ }
+
+ /// Assuming initially stress-free state.
+ double fracture_opening_at_residual_traction(double const t,
+ X const& x) const
+ {
+ return 2 * fracture_toughness(t, x)[0] / peak_normal_traction(t, x)[0];
+ }
+
+ /// Normal stiffness given in units of stress per length.
+ P const& normal_stiffness;
+ /// Shear stiffness given in units of stress per length.
+ P const& shear_stiffness;
+
+ /// Fracture toughness/critical energy release rate given in of stress
+ /// times lengths.
+ P const& fracture_toughness;
+ /// Peak normal traction given in units of stress.
+ P const& peak_normal_traction;
+};
+
+/// Evaluated MaterialPropertiesParameters container, see its documentation for
+/// details.
+struct MaterialProperties final
+{
+ MaterialProperties(double const t, ProcessLib::SpatialPosition const& x,
+ MaterialPropertiesParameters const& mp)
+ : Kn(mp.normal_stiffness(t, x)[0]),
+ Ks(mp.shear_stiffness(t, x)[0]),
+ Gc(mp.fracture_toughness(t, x)[0]),
+ t_np(mp.peak_normal_traction(t, x)[0]),
+ w_np(mp.fracture_opening_at_peak_traction(t, x)),
+ w_nf(mp.fracture_opening_at_residual_traction(t, x))
+ {
+ }
+
+ double const Kn;
+ double const Ks;
+ double const Gc;
+ double const t_np;
+ double const w_np;
+ double const w_nf;
+};
+
+template <int DisplacementDim>
+struct StateVariables
+ : public FractureModelBase<DisplacementDim>::MaterialStateVariables
+{
+ void setInitialConditions() { damage = damage_prev; }
+
+ void pushBackState() override { damage_prev = damage; }
+
+ double damage; ///< damage part of the state.
+
+ // Initial values from previous timestep
+ double damage_prev; ///< \copydoc damage
+};
+
+template <int DisplacementDim>
+class CohesiveZoneModeI final : public FractureModelBase<DisplacementDim>
+{
+public:
+ std::unique_ptr<
+ typename FractureModelBase<DisplacementDim>::MaterialStateVariables>
+ createMaterialStateVariables() override
+ {
+ return std::make_unique<StateVariables<DisplacementDim>>();
+ }
+
+public:
+public:
+ explicit CohesiveZoneModeI(double const penalty_aperture_cutoff,
+ bool const tension_cutoff,
+ MaterialPropertiesParameters material_properties)
+ : _penalty_aperture_cutoff(penalty_aperture_cutoff),
+ _tension_cutoff(tension_cutoff),
+ _mp(std::move(material_properties))
+ {
+ }
+
+ /**
+ * Computation of the constitutive relation for the Mohr-Coulomb model.
+ *
+ * @param t current time
+ * @param x current position in space
+ * @param aperture0 initial fracture's aperture
+ * @param sigma0 initial stress
+ * @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 Kep tangent matrix for stress and fracture displacements
+ * @param material_state_variables material state variables
+ */
+ void computeConstitutiveRelation(
+ double const t,
+ ProcessLib::SpatialPosition const& x,
+ double const aperture0,
+ Eigen::Ref<Eigen::VectorXd const>
+ sigma0,
+ 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,
+ typename FractureModelBase<DisplacementDim>::MaterialStateVariables&
+ material_state_variables) override;
+
+ MaterialProperties evaluatedMaterialProperties(
+ double const t, ProcessLib::SpatialPosition const& x) const
+ {
+ return MaterialProperties(t, x, _mp);
+ }
+
+private:
+ /// Compressive normal displacements above this value will not enter the
+ /// computation of the normal stiffness modulus of the fracture.
+ /// \note Setting this to the initial aperture value allows negative
+ /// apertures.
+ double const _penalty_aperture_cutoff;
+
+ /// If set no resistance to open the fracture over the initial aperture is
+ /// opposed.
+ bool const _tension_cutoff;
+
+ MaterialPropertiesParameters _mp;
+};
+
+} // namespace CohesiveZoneModeI
+} // namespace Fracture
+} // namespace MaterialLib
+
+namespace MaterialLib
+{
+namespace Fracture
+{
+namespace CohesiveZoneModeI
+{
+extern template class CohesiveZoneModeI<2>;
+extern template class CohesiveZoneModeI<3>;
+} // namespace CohesiveZoneModeI
+} // namespace Fracture
+} // namespace MaterialLib
diff --git a/MaterialLib/FractureModels/CreateCohesiveZoneModeI.cpp b/MaterialLib/FractureModels/CreateCohesiveZoneModeI.cpp
new file mode 100644
index 00000000000..8d4fcced10f
--- /dev/null
+++ b/MaterialLib/FractureModels/CreateCohesiveZoneModeI.cpp
@@ -0,0 +1,70 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2018, 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 "CreateCohesiveZoneModeI.h"
+
+#include "CohesiveZoneModeI.h"
+#include "ProcessLib/Utils/ProcessUtils.h" // required for findParameter
+
+namespace MaterialLib
+{
+namespace Fracture
+{
+namespace CohesiveZoneModeI
+{
+template <int DisplacementDim>
+std::unique_ptr<FractureModelBase<DisplacementDim>> createCohesiveZoneModeI(
+ std::vector<std::unique_ptr<ProcessLib::ParameterBase>> const& parameters,
+ BaseLib::ConfigTree const& config)
+{
+ //! \ogs_file_param{material__fracture_model__type}
+ config.checkConfigParameter("type", "CohesiveZoneModeI");
+ DBUG("Create CohesiveZoneModeI material");
+
+ auto& Kn = ProcessLib::findParameter<double>(
+ //! \ogs_file_param_special{material__fracture_model__CohesiveZoneModeI__normal_stiffness}
+ config, "normal_stiffness", parameters, 1);
+
+ auto& Ks = ProcessLib::findParameter<double>(
+ //! \ogs_file_param_special{material__fracture_model__CohesiveZoneModeI__shear_stiffness}
+ config, "shear_stiffness", parameters, 1);
+
+ auto& Gc = ProcessLib::findParameter<double>(
+ //! \ogs_file_param_special{material__fracture_model__CohesiveZoneModeI__fracture_toughness}
+ config, "fracture_toughness", parameters, 1);
+
+ auto& t_np = ProcessLib::findParameter<double>(
+ //! \ogs_file_param_special{material__fracture_model__CohesiveZoneModeI__peak_normal_traction}
+ config, "peak_normal_traction", parameters, 1);
+
+ auto const penalty_aperture_cutoff =
+ //! \ogs_file_param{material__fracture_model__CohesiveZoneModeI__penalty_aperture_cutoff}
+ config.getConfigParameter<double>("penalty_aperture_cutoff");
+
+ auto const tension_cutoff =
+ //! \ogs_file_param{material__fracture_model__CohesiveZoneModeI__tension_cutoff}
+ config.getConfigParameter<bool>("tension_cutoff");
+
+ MaterialPropertiesParameters mp{Kn, Ks, Gc, t_np};
+
+ return std::make_unique<CohesiveZoneModeI<DisplacementDim>>(
+ penalty_aperture_cutoff, tension_cutoff, mp);
+}
+
+template std::unique_ptr<FractureModelBase<2>> createCohesiveZoneModeI(
+ std::vector<std::unique_ptr<ProcessLib::ParameterBase>> const& parameters,
+ BaseLib::ConfigTree const& config);
+
+template std::unique_ptr<FractureModelBase<3>> createCohesiveZoneModeI(
+ std::vector<std::unique_ptr<ProcessLib::ParameterBase>> const& parameters,
+ BaseLib::ConfigTree const& config);
+
+} // namespace CohesiveZoneModeI
+} // namespace Fracture
+} // namespace MaterialLib
diff --git a/MaterialLib/FractureModels/CreateCohesiveZoneModeI.h b/MaterialLib/FractureModels/CreateCohesiveZoneModeI.h
new file mode 100644
index 00000000000..bcb6620aa08
--- /dev/null
+++ b/MaterialLib/FractureModels/CreateCohesiveZoneModeI.h
@@ -0,0 +1,27 @@
+/**
+ * \copyright
+ * Copyright (c) 2012-2018, 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 "FractureModelBase.h"
+
+namespace MaterialLib
+{
+namespace Fracture
+{
+namespace CohesiveZoneModeI
+{
+template <int DisplacementDim>
+std::unique_ptr<FractureModelBase<DisplacementDim>> createCohesiveZoneModeI(
+ std::vector<std::unique_ptr<ProcessLib::ParameterBase>> const& parameters,
+ BaseLib::ConfigTree const& config);
+
+} // namespace CohesiveZoneModeI
+} // namespace Fracture
+} // namespace MaterialLib
--
GitLab