Homec4science
Diffusion tamaas (features/netcdf-mpi) Browse

tamaas/src/model/elasto_plastic/isotropic_hardening.hh
383b32a7e7befeatures/netcdf-mpi

Last Change
Raw File

isotropic_hardening.hh
View Options

/**
* @file
*
* @author Lucas Frérot <lucas.frerot@epfl.ch>
*
* @section LICENSE
*
* Copyright (©) 2017 EPFL (Ecole Polytechnique Fédérale de
* Lausanne) Laboratory (LSMS - Laboratoire de Simulation en Mécanique des
* Solides)
*
* Tamaas 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.
*
* Tamaas 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 Tamaas. If not, see <http://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"
/* -------------------------------------------------------------------------- */
__BEGIN_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;
Real hardening(Real p) const { return sigma_0 + h * p; }
TAMAAS_ACCESSOR(h, Real, HardeningModulus);
TAMAAS_ACCESSOR(sigma_0, Real, YieldStress);
const Grid<Real, dim>& getPlasticStrain() const { return plastic_strain; }
Grid<Real, dim>& getPlasticStrain() { return plastic_strain; }
protected:
Model* model;
Real sigma_0; /// < initial yield stress
Real h; /// < hardening modulus
Grid<Real, dim> 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 = MatrixProxy<Real, dim, dim>;
using cpmatrix = MatrixProxy<const Real, dim, dim>;
Loop::loop(
[&elasticity,
this](MatrixProxy<Real, dim, dim> dep, // < plastic strain increment
MatrixProxy<const Real, dim, dim> epsilon, // < total strain
MatrixProxy<const Real, dim, dim>
delta_epsilon, // < strain increment
MatrixProxy<Real, dim, dim> 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 - hardening(p) > 0) {
// radial return
Real dp = (von_mises - hardening(p)) / (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) {
this->template computePlasticIncrement<update>(stress, strain,
strain_increment);
const influence::ElasticHelper<dim> elasticity(
this->model->getShearModulus(), this->model->getPoissonRatio());
using pmatrix = MatrixProxy<Real, dim, dim>;
using cpmatrix = MatrixProxy<const Real, dim, dim>;
Loop::loop(
[&elasticity](
MatrixProxy<Real, dim, dim> sigma, // < stress
MatrixProxy<const Real, dim, dim> epsilon, // < total strain
MatrixProxy<const Real, dim, dim>
delta_epsilon, // < strain increment
MatrixProxy<Real, dim, dim> ep /* < plastic strain */) {
// !!! sigma contains delta_ep
sigma *= -1;
sigma += epsilon;
sigma -= ep;
sigma += delta_epsilon;
sigma = elasticity(sigma);
},
range<pmatrix>(stress), range<cpmatrix>(strain),
range<cpmatrix>(strain_increment), range<pmatrix>(this->plastic_strain));
}
__END_TAMAAS__
/* -------------------------------------------------------------------------- */
#endif // ISOTROPIC_HARDENING_HH
c4science · Help