Homec4science
Diffusion akantu (features/joss-2023) Browse

rAKA/src/model/solid_mechanics/materials/material_damage/material_anisotropic_damage_tmpl.hh
34b888c7b651features/joss-2023

Last Change
Raw File

material_anisotropic_damage_tmpl.hh
View Options

/**
* @file material_anisotropic_damage_tmpl.hh
*
* @author Nicolas Richart
*
* @date creation mar jun 25 2019
*
* @brief A Documented file.
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu is distributed in the hope that it will be useful, but WITHOUT ANY
* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
* A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "aka_iterators.hh"
#include "material_anisotropic_damage.hh"
/* -------------------------------------------------------------------------- */
#ifndef __AKANTU_MATERIAL_ANISOTROPIC_DAMAGE_TMPL_HH__
#define __AKANTU_MATERIAL_ANISOTROPIC_DAMAGE_TMPL_HH__
namespace akantu {
struct EmptyIteratorContainer {
struct iterator {
auto & operator++() { return *this; }
Real operator*() { return 0; }
bool operator!=(const iterator &) const { return true; }
bool operator==(const iterator &) const { return false; }
};
auto begin() { return iterator(); }
auto end() { return iterator(); }
};
} // namespace akantu
namespace std {
template <> struct iterator_traits<::akantu::EmptyIteratorContainer::iterator> {
using iterator_category = forward_iterator_tag;
using value_type = akantu::Real;
using difference_type = std::ptrdiff_t;
using pointer = akantu::Real *;
using reference = akantu::Real &;
};
} // namespace std
namespace akantu {
namespace {
template <UInt dim, class Op>
void tensorPlus_(const Matrix<Real> & A, Op && oper) {
Vector<Real> A_eigs(dim);
A.eig(A_eigs);
for (auto & ap : A_eigs) {
oper(ap);
}
}
template <UInt dim> auto tensorPlus2(const Matrix<Real> & A) {
Real square = 0;
tensorPlus_<dim>(A, [&](Real eig) {
eig = std::max(eig, 0.);
square += eig * eig;
});
return square;
}
template <UInt dim> auto tensorPlusTrace(const Matrix<Real> & A) {
Real trace_plus = 0;
Real trace_minus = 0;
tensorPlus_<dim>(A, [&](Real eig) {
trace_plus += std::max(eig, 0.);
trace_minus += std::min(eig, 0.);
});
return std::make_pair(trace_plus, trace_minus);
}
template <UInt dim, class Op>
auto tensorPlusOp(const Matrix<Real> & A, Matrix<Real> & A_directions,
Op && oper, bool sorted = false) {
Vector<Real> A_eigs(dim);
Matrix<Real> A_diag(dim, dim);
A.eig(A_eigs, A_directions, sorted);
for (auto && data : enumerate(A_eigs)) {
auto i = std::get<0>(data);
A_diag(i, i) = oper(std::max(std::get<1>(data), 0.), i);
}
return A_directions * A_diag * A_directions.transpose();
}
template <UInt dim, class Op>
auto tensorPlus(const Matrix<Real> & A, Matrix<Real> & A_directions,
bool sorted = false) {
return tensorPlusOp<dim>(A, A_directions, [](Real x, Real) { return x; }, sorted);
}
template <UInt dim, class Op>
auto tensorPlusOp(const Matrix<Real> & A, Op && oper) {
Matrix<Real> A_directions(dim, dim);
return tensorPlusOp<dim>(A, A_directions, std::forward<Op>(oper));
}
template <UInt dim>
auto tensorPlus(const Matrix<Real> & A) {
return tensorPlusOp<dim>(A, [](Real x, Real) { return x; });
}
template <UInt dim> auto tensorSqrt(const Matrix<Real> & A) {
return tensorPlusOp<dim>(A, [](Real x, UInt) { return std::sqrt(x); });
}
} // namespace
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
template <UInt dim, template <UInt> class EquivalentStrain,
template <UInt> class DamageThreshold, template <UInt> class Parent>
MaterialAnisotropicDamage<dim, EquivalentStrain, DamageThreshold, Parent>::
MaterialAnisotropicDamage(SolidMechanicsModel & model, const ID & id)
: Parent<dim>(model, id), damage("damage_tensor", *this),
elastic_stress("elastic_stress", *this),
equivalent_strain("equivalent_strain", *this),
trace_damage("trace_damage", *this), equivalent_strain_function(*this),
damage_threshold_function(*this) {
this->registerParam("Dc", Dc, _pat_parsable, "Critical damage");
this->damage.initialize(dim * dim);
this->elastic_stress.initialize(dim * dim);
this->equivalent_strain.initialize(1);
this->trace_damage.initialize(1);
this->trace_damage.initializeHistory();
}
/* -------------------------------------------------------------------------- */
template <UInt dim, template <UInt> class EquivalentStrain,
template <UInt> class DamageThreshold, template <UInt> class Parent>
void MaterialAnisotropicDamage<dim, EquivalentStrain, DamageThreshold, Parent>::
damageStress(Matrix<Real> & sigma, const Matrix<Real> & sigma_el,
const Matrix<Real> & D, Real TrD) {
// σ_(n + 1) = (1 − D_(n + 1))^(1/2) σ~_(n + 1) (1 − D_(n + 1))^(1 / 2)
// - ((1 − D_(n + 1)) : σ~_(n + 1))/ (3 - Tr(D_(n+1))) (1 − D_(n + 1))
// + 1/3 (1 - Tr(D_(n+1)) <Tr(σ~_(n + 1))>_+ + <Tr(σ~_(n + 1))>_-) I
auto one_D = Matrix<Real>::eye(dim) - D;
auto sqrt_one_D = tensorSqrt<dim>(one_D);
Real Tr_sigma_plus, Tr_sigma_minus;
std::tie(Tr_sigma_plus, Tr_sigma_minus) = tensorPlusTrace<dim>(sigma_el);
auto I = Matrix<Real>::eye(dim);
sigma = sqrt_one_D * sigma_el * sqrt_one_D -
(one_D.doubleDot(sigma_el) / (dim - TrD) * one_D) +
1. / dim * ((1 - TrD) * Tr_sigma_plus - Tr_sigma_minus) * I;
}
/* -------------------------------------------------------------------------- */
template <UInt dim, template <UInt> class EquivalentStrain,
template <UInt> class DamageThreshold, template <UInt> class Parent>
void MaterialAnisotropicDamage<dim, EquivalentStrain, DamageThreshold,
Parent>::computeStress(ElementType type,
GhostType ghost_type) {
for (auto && data :
zip(make_view(this->stress(type, ghost_type), dim, dim),
make_view(this->gradu(type, ghost_type), dim, dim),
make_view(this->sigma_th(type, ghost_type)),
make_view(this->elastic_stress(type, ghost_type), dim, dim),
make_view(this->equivalent_strain(type, ghost_type)),
make_view(this->damage(type, ghost_type), dim, dim),
make_view(this->trace_damage(type, ghost_type)),
make_view(this->trace_damage.previous(type, ghost_type)),
equivalent_strain_function, damage_threshold_function)) {
auto & sigma = std::get<0>(data);
auto & grad_u = std::get<1>(data);
auto & sigma_th = std::get<2>(data);
auto & sigma_el = std::get<3>(data);
auto & epsilon_hat = std::get<4>(data);
auto & D = std::get<5>(data);
auto & TrD_n_1 = std::get<6>(data);
auto & TrD = std::get<7>(data);
auto & equivalent_strain_data = std::get<8>(data);
auto & damage_threshold_data = std::get<9>(data);
Matrix<Real> Dtmp(dim, dim);
Real TrD_n_1_tmp;
Matrix<Real> epsilon(dim, dim);
// yes you read properly this is a label for a goto
auto computeDamage = [&]() {
MaterialElastic<dim>::computeStressOnQuad(grad_u, sigma_el, sigma_th);
this->template gradUToEpsilon<dim>(grad_u, epsilon);
// evaluate the damage criteria
epsilon_hat = equivalent_strain_function(epsilon, equivalent_strain_data);
// evolve the damage if needed
auto K_TrD = damage_threshold_function.K(TrD, damage_threshold_data);
auto f = epsilon_hat - K_TrD;
// if test function > 0 evolve the damage
if (f > 0) {
TrD_n_1_tmp =
damage_threshold_function.K_inv(epsilon_hat, damage_threshold_data);
auto epsilon_plus = tensorPlus<dim>(epsilon);
auto delta_lambda = (TrD_n_1_tmp - TrD) / (epsilon_hat * epsilon_hat);
Dtmp = D + delta_lambda * epsilon_plus;
return true;
}
return false;
};
// compute a temporary version of the new damage
auto is_damage_updated = computeDamage();
if (is_damage_updated) {
/// Check and correct for broken case
if (Dtmp.trace() > Dc) {
if (epsilon.trace() > 0) { // tensile loading
auto kpa = this->kpa;
auto lambda = this->lambda;
// change kappa to Kappa_broken = (1-Dc) Kappa
kpa = (1 - Dc) * kpa;
this->E = 9 * kpa * (kpa - lambda) / (3 * kpa - lambda);
this->nu = lambda / (3 * kpa - lambda);
this->updateInternalParameters();
computeDamage();
} else if (epsilon.trace() < 1e-10) { // deviatoric case
Matrix<Real> n(dim, dim);
std::vector<UInt> ns;
tensorPlusOp<dim>(Dtmp, n,
[&](Real x, UInt i) {
if (x > this->Dc) {
ns.push_back(i);
return this->Dc;
}
return x;
},
true);
}
}
TrD_n_1 = TrD_n_1_tmp;
D = Dtmp;
} else {
TrD_n_1 = TrD;
}
// apply the damage to the stress
damageStress(sigma, sigma_el, D, TrD_n_1);
}
}
/* -------------------------------------------------------------------------- */
/* EquivalentStrain functions */
/* -------------------------------------------------------------------------- */
template <UInt dim>
class EquivalentStrainMazars : public EmptyIteratorContainer {
public:
EquivalentStrainMazars(Material & /*mat*/) {}
template <class... Other>
Real operator()(const Matrix<Real> & epsilon, Other &&... /*other*/) {
Real epsilon_hat = 0.;
std::tie(epsilon_hat, std::ignore) = tensorPlusTrace<dim>(epsilon);
return std::sqrt(epsilon_hat);
}
};
template <UInt dim>
class EquivalentStrainMazarsDruckerPrager : public EquivalentStrainMazars<dim> {
public:
EquivalentStrainMazarsDruckerPrager(Material & mat)
: EquivalentStrainMazars<dim>(mat) {
mat.registerParam("k", k, _pat_parsable, "k");
}
template <class... Other>
Real operator()(const Matrix<Real> & epsilon, Real) {
Real epsilon_hat = EquivalentStrainMazars<dim>::operator()(epsilon);
epsilon_hat += k * epsilon.trace();
return epsilon_hat;
}
protected:
Real k;
};
/* -------------------------------------------------------------------------- */
/* DamageThreshold functions */
/* -------------------------------------------------------------------------- */
template <UInt dim>
class DamageThresholdLinear : public EmptyIteratorContainer {
public:
DamageThresholdLinear(Material & mat) : mat(mat) {
mat.registerParam("A", A, _pat_parsable, "A");
mat.registerParam("K0", K0, _pat_parsable, "K0");
}
template <class... Other> Real K(Real x, Other &&... /*other*/) {
return 1. / A * x + K0;
}
template <class... Other> Real K_inv(Real x, Other &&... /*other*/) {
return A * (x - K0);
}
private:
Material & mat;
Real A;
Real K0;
};
template <UInt dim> class DamageThresholdTan : public EmptyIteratorContainer {
public:
DamageThresholdTan(Material & mat) : mat(mat) {
mat.registerParam("a", a, _pat_parsable, "a");
mat.registerParam("A", A, _pat_parsable, "A");
mat.registerParam("K0", K0, _pat_parsable, "K0");
}
template <class... Other> Real K(Real x, Other &&... /*other*/) {
return a * std::tan(std::atan2(x, a) - std::atan2(K0, a));
}
template <class... Other> Real K_inv(Real x, Other &&... /*other*/) {
return a * A * (std::atan2(x, a) - std::atan2(K0, a));
}
private:
Material & mat;
Real a{2.93e-4};
Real A{5e3};
Real K0{5e-5};
};
} // namespace akantu
#endif /* __AKANTU_MATERIAL_ANISOTROPIC_DAMAGE_TMPL_HH__ */
c4science · Help