From 666d456e99b4005809063583405f07c6ea68c6ac Mon Sep 17 00:00:00 2001
From: Francesco Parisio <francesco.parisio@protonmail.com>
Date: Fri, 30 Nov 2018 12:31:55 +0100
Subject: [PATCH] [PL] SDN; Damage calculations.
---
ProcessLib/SmallDeformationNonlocal/Damage.h | 78 +++++++++++++++++++
.../SmallDeformationNonlocalFEM.h | 4 +-
2 files changed, 81 insertions(+), 1 deletion(-)
create mode 100644 ProcessLib/SmallDeformationNonlocal/Damage.h
diff --git a/ProcessLib/SmallDeformationNonlocal/Damage.h b/ProcessLib/SmallDeformationNonlocal/Damage.h
new file mode 100644
index 00000000000..487115ccc71
--- /dev/null
+++ b/ProcessLib/SmallDeformationNonlocal/Damage.h
@@ -0,0 +1,78 @@
+/**
+ * \file
+ *
+ * \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/Eigenvalues>
+#include <logog/include/logog.hpp>
+
+#include "MaterialLib/SolidModels/Ehlers.h"
+
+namespace ProcessLib
+{
+namespace SmallDeformationNonlocal
+{
+/// Computes the damage internal material variable explicitly based on the
+/// results obtained from the local stress return algorithm.
+inline double calculateDamage(double const kappa_d, double const alpha_d,
+ double const beta_d)
+{
+ double const damage = (1 - beta_d) * (1 - std::exp(-kappa_d / alpha_d));
+
+ if (damage < 0. || damage >= 1.)
+ {
+ OGS_FATAL("Damage value %g outside of [0,1) interval.", damage);
+ }
+
+ return damage;
+}
+
+template <int DisplacementDim, typename KelvinVectorType>
+double calculateDamageKappaD(
+ double const eps_p_eff_diff,
+ KelvinVectorType const& sigma,
+ double const kappa_d_prev,
+ double const h_d,
+ MaterialLib::Solids::Ehlers::MaterialProperties const& mp)
+{
+ // Default case of the rate problem. Updated below if volumetric plastic
+ // strain rate is positive (dilatancy).
+
+ // non-const for Eigen solver.
+ auto sigma_tensor = MathLib::KelvinVector::kelvinVectorToTensor(sigma);
+
+ Eigen::EigenSolver<decltype(sigma_tensor)> eigen_solver(sigma_tensor);
+ auto const principal_stress = real(eigen_solver.eigenvalues().array());
+ double const prod_stress = std::sqrt(principal_stress.square().sum());
+
+ // Brittleness decrease with confinement for the nonlinear flow rule.
+ // ATTENTION: For linear flow rule -> constant brittleness.
+ double const tensile_strength =
+ std::sqrt(3.0) * mp.kappa / (1 + std::sqrt(3.0) * mp.beta);
+ double const r_s = prod_stress / tensile_strength - 1.;
+
+ // Compute normalizing strain.
+ double const x_s = [](double const h_d, double const r_s) {
+ if (r_s < 0)
+ {
+ return 1.;
+ }
+ if (r_s <= 1)
+ {
+ return 1. + h_d * r_s * r_s;
+ }
+ return 1. - 3 * h_d + 4 * h_d * std::sqrt(r_s);
+ }(h_d, r_s);
+
+ return kappa_d_prev + eps_p_eff_diff / x_s;
+}
+} // namespace SmallDeformationNonlocal
+} // namespace ProcessLib
diff --git a/ProcessLib/SmallDeformationNonlocal/SmallDeformationNonlocalFEM.h b/ProcessLib/SmallDeformationNonlocal/SmallDeformationNonlocalFEM.h
index 5fa904f23a5..0d792388b84 100644
--- a/ProcessLib/SmallDeformationNonlocal/SmallDeformationNonlocalFEM.h
+++ b/ProcessLib/SmallDeformationNonlocal/SmallDeformationNonlocalFEM.h
@@ -362,7 +362,9 @@ public:
// Compute sigma_eff from damage total stress sigma
using KelvinVectorType = typename BMatricesType::KelvinVectorType;
KelvinVectorType const sigma_eff_prev =
- sigma_prev / (1. - damage_prev);
+ sigma_prev /
+ (1. - damage_prev); // damage_prev is in [0,1) range. See
+ // calculateDamage() function.
auto&& solution = _ip_data[ip].solid_material.integrateStress(
t, x_position, _process_data.dt, eps_prev, eps, sigma_eff_prev,
--
GitLab