Page MenuHomec4science
Files F61083069

isotropic_hardening.hh
No OneTemporary
Actions

File Metadata

Created
Sat, May 4, 10:33

isotropic_hardening.hh
View Options

/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
* Copyright (©) 2016-2023 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
* Copyright (©) 2020-2023 Lucas Frérot
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Affero General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program 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 Affero General Public License for more details.
*
* You should have received a copy of the GNU Affero General Public License
* along with this program. If not, see <https://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#ifndef ISOTROPIC_HARDENING_HH
#define ISOTROPIC_HARDENING_HH
/* -------------------------------------------------------------------------- */
#include "grid.hh"
#include "influence.hh"
#include "model.hh"
#include "model_type.hh"
/* -------------------------------------------------------------------------- */
namespace tamaas {
/* -------------------------------------------------------------------------- */
template <model_type type>
class IsotropicHardening {
using trait = model_type_traits<type>;
static constexpr UInt dim = trait::dimension;
public:
/// Constructor
IsotropicHardening(Model* model, Real sigma_0, Real h);
/// Compute plastic strain increment with radial return algorithm
template <bool update>
void computePlasticIncrement(Grid<Real, dim>& increment,
const Grid<Real, dim>& strain,
const Grid<Real, dim>& strain_increment);
/// Compute stress
template <bool update>
void computeStress(Grid<Real, dim>& stress, const Grid<Real, dim>& strain,
const Grid<Real, dim>& strain_increment);
void applyTangentIncrement(Grid<Real, dim>& output,
const Grid<Real, dim>& input,
const Grid<Real, dim>& strain,
const Grid<Real, dim>& strain_increment) const;
static __device__ __host__ Real hardening(Real p, Real h, Real sigma_0) {
return sigma_0 + h * p;
}
Real getHardeningModulus() const { return h; }
Real getYieldStress() const { return sigma_0; }
void setHardeningModulus(Real h_) {
if (h_ < 0)
TAMAAS_EXCEPTION("Hardening modulus should be positive");
h = h_;
}
void setYieldStress(Real sigma_0_) {
if (sigma_0_ < 0)
TAMAAS_EXCEPTION("Yield stress should be positive");
sigma_0 = sigma_0_;
}
const GridBase<Real>& getPlasticStrain() const { return *plastic_strain; }
GridBase<Real>& getPlasticStrain() { return *plastic_strain; }
protected:
Model* model;
Real sigma_0; /// < initial yield stress
Real h; /// < hardening modulus
std::shared_ptr<Grid<Real, 3>> plastic_strain, cumulated_plastic_strain;
};
/* -------------------------------------------------------------------------- */
template <>
template <bool update>
void IsotropicHardening<model_type::volume_2d>::computePlasticIncrement(
Grid<Real, dim>& increment, const Grid<Real, dim>& strain,
const Grid<Real, dim>& strain_increment) {
influence::ElasticHelper<dim> elasticity(this->model->getShearModulus(),
this->model->getPoissonRatio());
using pmatrix = SymMatrixProxy<Real, dim>;
using cpmatrix = SymMatrixProxy<const Real, dim>;
Real sigma_0 = this->sigma_0, h = this->h; // for captures
Loop::loop(
[elasticity, sigma_0,
h] CUDA_LAMBDA(pmatrix dep, // < plastic strain increment
cpmatrix epsilon, // < total strain
cpmatrix delta_epsilon, // < strain increment
pmatrix ep, // < plastic strain
Real & p /* < cumulated plastic strain */) {
// Trial elastic stress
auto sigma_tr = elasticity(epsilon - ep + delta_epsilon);
// Zero increment by default
dep = 0;
// Deviatoric trial stress
decltype(sigma_tr) dev;
dev.deviatoric(sigma_tr);
// Von-Mises stress
Real von_mises = std::sqrt(1.5) * dev.l2norm();
// Plasticity condition
if (von_mises - IsotropicHardening<model_type::volume_2d>::hardening(
p, h, sigma_0) >
0) {
// radial return
Real dp =
(von_mises - IsotropicHardening<model_type::volume_2d>::hardening(
p, h, sigma_0)) /
(3 * elasticity.mu + h);
dev *= 3 * dp / (2 * von_mises); // saving memory (dev is delta ep)
dep = dev;
if (update) { // constexpr if when C++17
p += dp;
ep += dep;
}
}
},
range<pmatrix>(increment), range<cpmatrix>(strain),
range<cpmatrix>(strain_increment), range<pmatrix>(*this->plastic_strain),
*this->cumulated_plastic_strain);
}
/* -------------------------------------------------------------------------- */
template <>
template <bool update>
void IsotropicHardening<model_type::volume_2d>::computeStress(
Grid<Real, dim>& stress, const Grid<Real, dim>& strain,
const Grid<Real, dim>& strain_increment) {
using pmatrix = SymMatrixProxy<Real, dim>;
using cpmatrix = SymMatrixProxy<const Real, dim>;
const influence::ElasticHelper<dim> elasticity(
this->model->getShearModulus(), this->model->getPoissonRatio());
this->template computePlasticIncrement<update>(stress, strain,
strain_increment);
Loop::loop(
[elasticity] CUDA_LAMBDA(pmatrix sigma, // < stress
cpmatrix epsilon, // < total strain
cpmatrix delta_epsilon, // < strain increment
cpmatrix ep // < plastic strain
) {
// !!! sigma contains delta_ep
if (not update) {
// dep is not contained in ep
sigma *= -1;
sigma += epsilon;
} else {
// dep is contained in ep, can discard value of sigma
sigma = epsilon;
}
sigma -= ep;
sigma += delta_epsilon;
sigma = elasticity(sigma);
},
range<pmatrix>(stress), range<cpmatrix>(strain),
range<cpmatrix>(strain_increment),
range<cpmatrix>(*this->plastic_strain));
}
} // namespace tamaas
/* -------------------------------------------------------------------------- */
#endif // ISOTROPIC_HARDENING_HH

Event Timeline

c4science · Help