From dfaf9b9fc8fd30df2bd16fb3f73c1a7ba015045e Mon Sep 17 00:00:00 2001
From: Francesco Parisio <francesco.parisio@protonmail.com>
Date: Fri, 10 Feb 2017 12:11:53 +0100
Subject: [PATCH] [MatL] Ehlers: Add damage state and updateDamage.
---
MaterialLib/SolidModels/Ehlers-impl.h | 64 ++++++++++++++++++++++++++-
MaterialLib/SolidModels/Ehlers.h | 23 +++++++++-
2 files changed, 83 insertions(+), 4 deletions(-)
diff --git a/MaterialLib/SolidModels/Ehlers-impl.h b/MaterialLib/SolidModels/Ehlers-impl.h
index d5415da9d24..92c401d855c 100644
--- a/MaterialLib/SolidModels/Ehlers-impl.h
+++ b/MaterialLib/SolidModels/Ehlers-impl.h
@@ -464,6 +464,50 @@ void calculatePlasticJacobian(
// G_52, G_53, G_55 are zero
}
+template <int DisplacementDim>
+void SolidEhlers<DisplacementDim>::updateDamage(
+ double const t, ProcessLib::SpatialPosition const& x,
+ typename MechanicsBase<DisplacementDim>::MaterialStateVariables&
+ material_state_variables)
+{
+ assert(dynamic_cast<MaterialStateVariables*>(&material_state_variables) !=
+ nullptr);
+ MaterialStateVariables& state =
+ static_cast<MaterialStateVariables&>(material_state_variables);
+
+ // Default case of the rate problem. Updated below if volumetric plastic
+ // strain rate is positive (dilatancy).
+ state.kappa_d = state.kappa_d_prev;
+
+ // Compute damage current step
+ double const eps_p_V_dot = state.eps_p_V - state.eps_p_V_prev;
+ if (eps_p_V_dot > 0)
+ {
+ double const h_d = _damage_properties->h_d(t, x)[0];
+ double const eps_p_eff_dot = state.eps_p_eff - state.eps_p_eff_prev;
+ double const r_s = eps_p_eff_dot / eps_p_V_dot;
+
+ // Brittleness decrease with confinement for the nonlinear flow rule.
+ // ATTENTION: For linear flow rule -> constant brittleness.
+ double x_s = 0;
+ if (r_s < 1)
+ {
+ x_s = 1 + h_d * r_s * r_s;
+ }
+ else
+ {
+ x_s = 1 - 3 * h_d + 4 * h_d * std::sqrt(r_s);
+ }
+ state.kappa_d += eps_p_eff_dot / x_s;
+ }
+
+ double const alpha_d = _damage_properties->alpha_d(t, x)[0];
+ double const beta_d = _damage_properties->beta_d(t, x)[0];
+
+ // Update internal damage variable.
+ state.damage = (1 - beta_d) * (1 - std::exp(-state.kappa_d / alpha_d));
+}
+
/// Calculates the derivative of the residuals with respect to total
/// strain. Implementation fully implicit only.
template <int DisplacementDim>
@@ -551,8 +595,17 @@ bool SolidEhlers<DisplacementDim>::computeConstitutiveRelation(
double const G = _mp.G(t, x)[0];
double const K = _mp.K(t, x)[0];
+ KelvinVector sigma_eff_prev = sigma_prev; // In case without damage the
+ // effective value is same as the
+ // previous one.
+ if (_damage_properties)
+ {
+ // Compute sigma_eff from damage total stress sigma, which is given by
+ // sigma_eff=sigma_prev / (1-damage)
+ sigma_eff_prev = sigma_prev / (1 - _state.damage_prev);
+ }
KelvinVector sigma =
- predict_sigma<DisplacementDim>(G, K, sigma_prev, eps, eps_prev, eps_V);
+ predict_sigma<DisplacementDim>(G, K, sigma_eff_prev, eps, eps_prev, eps_V);
// update parameter
_mp.calculateIsotropicHardening(t, x, _state.eps_p_eff);
@@ -644,6 +697,9 @@ bool SolidEhlers<DisplacementDim>::computeConstitutiveRelation(
dresidual_deps.template block<KelvinVectorSize, KelvinVectorSize>(0, 0)
.noalias() = calculateDResidualDEps<DisplacementDim>(K, G);
+ if (_damage_properties)
+ updateDamage(t, x, _state);
+
// Extract consistent tangent.
C.noalias() =
_mp.G(t, x)[0] *
@@ -652,7 +708,11 @@ bool SolidEhlers<DisplacementDim>::computeConstitutiveRelation(
}
// Update sigma.
- sigma_final.noalias() = G * sigma;
+ if (_damage_properties)
+ sigma_final = G * sigma * (1 - _state.damage);
+ else
+ sigma_final.noalias() = G * sigma;
+
return true;
}
diff --git a/MaterialLib/SolidModels/Ehlers.h b/MaterialLib/SolidModels/Ehlers.h
index b7e0ec00479..7b60c649c91 100644
--- a/MaterialLib/SolidModels/Ehlers.h
+++ b/MaterialLib/SolidModels/Ehlers.h
@@ -143,6 +143,8 @@ public:
eps_p_D = eps_p_D_prev;
eps_p_V = eps_p_V_prev;
eps_p_eff = eps_p_eff_prev;
+ kappa_d = kappa_d_prev;
+ damage = damage_prev;
lambda = 0;
}
@@ -151,6 +153,8 @@ public:
eps_p_D_prev = eps_p_D;
eps_p_V_prev = eps_p_V;
eps_p_eff_prev = eps_p_eff; // effective part of trace(eps_p)
+ kappa_d_prev = kappa_d;
+ damage_prev = damage;
lambda = 0;
}
@@ -159,13 +163,16 @@ public:
KelvinVector eps_p_D; ///< deviatoric plastic strain
double eps_p_V = 0; ///< volumetric strain
double eps_p_eff = 0; ///< effective plastic strain
+ double kappa_d = 0; ///< damage driving variable
+ double damage = 0; ///< isotropic damage variable
// Initial values from previous timestep
KelvinVector eps_p_D_prev; ///< \copydoc eps_p_D
double eps_p_V_prev = 0; ///< \copydoc eps_p_V
double eps_p_eff_prev = 0; ///< \copydoc eps_p_eff
-
- double lambda = 0; ///< plastic multiplier
+ double kappa_d_prev = 0; ///< \copydoc kappa_d
+ double damage_prev = 0; ///< \copydoc damage
+ double lambda = 0; ///< plastic multiplier
#ifndef NDEBUG
friend std::ostream& operator<<(std::ostream& os,
@@ -174,8 +181,12 @@ public:
os << "State:\n"
<< "eps_p_D: " << m.eps_p_D << "\n"
<< "eps_p_eff: " << m.eps_p_eff << "\n"
+ << "kappa_d: " << m.kappa_d << "\n"
+ << "damage: " << m.damage << "\n"
<< "eps_p_D_prev: " << m.eps_p_D_prev << "\n"
<< "eps_p_eff_prev: " << m.eps_p_eff_prev << "\n"
+ << "kappa_d_prev: " << m.kappa_d_prev << "\n"
+ << "damage_prev: " << m.damage_prev << "\n"
<< "lambda: " << m.lambda << "\n";
return os;
}
@@ -216,6 +227,14 @@ public:
typename MechanicsBase<DisplacementDim>::MaterialStateVariables&
material_state_variables) override;
+private:
+ /// Computes the damage internal material variable explicitly based on the
+ /// results obtained from the local stress return algorithm.
+ void updateDamage(
+ double const t, ProcessLib::SpatialPosition const& x,
+ typename MechanicsBase<DisplacementDim>::MaterialStateVariables&
+ material_state_variables);
+
private:
MaterialProperties _mp;
std::unique_ptr<EhlersDamageProperties> _damage_properties;
--
GitLab