Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F80545447
isotropic_hardening.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Sun, Sep 1, 08:01
Size
8 KB
Mime Type
text/x-c
Expires
Tue, Sep 3, 08:01 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
20395098
Attached To
rTAMAAS tamaas
isotropic_hardening.cpp
View Options
/*
* SPDX-License-Indentifier: AGPL-3.0-or-later
*
* Copyright (©) 2016-2024 EPFL (École Polytechnique Fédérale de Lausanne),
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
* Copyright (©) 2020-2024 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/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "isotropic_hardening.hh"
#include "influence.hh"
#include "tamaas.hh"
/* -------------------------------------------------------------------------- */
namespace tamaas {
/* -------------------------------------------------------------------------- */
IsotropicHardening::IsotropicHardening(Model* model, Real sigma_0, Real h)
: Material(model), sigma_0(sigma_0), h(h) {
if (model->getType() != type)
throw model_type_error{TAMAAS_MSG(
"IsotropicHardening law only valid for model_type::volume_2d")};
plastic_strain = model->request<type, false, Real>(
"plastic_strain", model->getDiscretization(), trait::voigt);
cumulated_plastic_strain = model->request<type, false, Real>(
"cumulated_plastic_strain", model->getDiscretization(), 1);
plastic_strain.initialize();
cumulated_plastic_strain.initialize();
}
/* -------------------------------------------------------------------------- */
void IsotropicHardening::computeInelasticDeformationIncrement(
View& increment, const View& strain, const View& strain_increment) {
influence::ElasticHelper<dim> elasticity(this->model->getShearModulus(),
this->model->getPoissonRatio());
plastic_strain.reset();
cumulated_plastic_strain.reset();
Real sigma_0 = this->sigma_0, h = this->h; // for captures
Loop::loop(
[elasticity, sigma_0,
h] CUDA_LAMBDA(Mat dep, // < plastic strain increment
CMat epsilon, // < total strain
CMat delta_epsilon, // < strain increment
Mat 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
const Real von_mises = std::sqrt(1.5) * dev.l2norm();
const Real yield =
von_mises - IsotropicHardening::hardening(p, h, sigma_0);
// Plasticity condition
if (yield > 0) {
// radial return
const Real dp = yield / (3 * elasticity.mu + h);
dev *= 3 * dp / (2 * von_mises); // saving memory (dev is delta ep)
dep = dev;
p += dp;
ep += dep;
}
},
range<Mat>(increment), range<CMat>(strain), range<CMat>(strain_increment),
range<Mat>(*plastic_strain), *cumulated_plastic_strain);
}
/* -------------------------------------------------------------------------- */
void IsotropicHardening::computeStress(Field& stress, const Field& strain,
const Field& strain_increment) {
const auto n = this->model->getDiscretization();
constexpr auto comp = trait::voigt;
Grid<Real, dim> stressv{n, comp, stress.view()};
const Grid<Real, dim> strainv{n, comp, strain.view()},
inc{n, comp, strain_increment.view()};
this->computeInelasticDeformationIncrement(stressv, strainv, inc);
const influence::ElasticHelper<dim> elasticity(
this->model->getShearModulus(), this->model->getPoissonRatio());
Loop::loop(
[elasticity] CUDA_LAMBDA(Mat sigma, // < stress
CMat epsilon, // < total strain
CMat delta_epsilon, // < strain increment
CMat ep // < plastic strain
) {
sigma = epsilon;
sigma += delta_epsilon;
sigma -= ep;
sigma = elasticity(sigma);
},
range<Mat>(stressv), range<CMat>(strainv), range<CMat>(inc),
range<CMat>(*plastic_strain));
}
/* -------------------------------------------------------------------------- */
void IsotropicHardening::computeEigenStress(Field& stressv,
const Field& strainv,
const Field& strain_incrementv) {
const auto n = this->model->getDiscretization();
constexpr auto comp = trait::voigt;
// Ensure grids have the right shape
Grid<Real, dim> stress{n, comp, stressv.view()};
const Grid<Real, dim> strain{n, comp, strainv.view()},
strain_increment{n, comp, strain_incrementv.view()};
this->computeInelasticDeformationIncrement(stress, strain, strain_increment);
const influence::ElasticHelper<dim> elasticity(
this->model->getShearModulus(), this->model->getPoissonRatio());
Loop::loop([elasticity] CUDA_LAMBDA(Mat sigma, // < stress
CMat ep, // < plastic strain
CMat ep_prev // < previous plastic strain
) { sigma = elasticity(ep - ep_prev); },
range<Mat>(stress), range<CMat>(*plastic_strain),
range<CMat>(*plastic_strain.previous));
}
/* -------------------------------------------------------------------------- */
void IsotropicHardening::update() {
plastic_strain.update();
cumulated_plastic_strain.update();
}
/* -------------------------------------------------------------------------- */
void IsotropicHardening::applyTangent(Field& outputv, const Field& inputv,
const Field& strainv,
const Field& strain_incrementv) {
const auto n = this->model->getDiscretization();
constexpr auto comp = trait::voigt;
// Ensure grids have the right shape
Grid<Real, dim> output{n, comp, outputv.view()};
const Grid<Real, dim> input{n, comp, inputv.view()},
strain{n, comp, strainv.view()},
strain_increment{n, comp, strain_incrementv.view()};
const influence::ElasticHelper<dim> elasticity(
this->model->getShearModulus(), this->model->getPoissonRatio());
auto sigma_0 = this->sigma_0, h = this->h; // for captures
plastic_strain.reset();
cumulated_plastic_strain.reset();
Loop::loop(
[elasticity, sigma_0, h] CUDA_LAMBDA(Mat out, CMat in, CMat epsilon,
CMat delta_epsilon, CMat ep,
const Real& p) {
auto sigma_tr = elasticity(epsilon - ep + delta_epsilon);
decltype(sigma_tr) dev;
dev.deviatoric(sigma_tr, 3);
const Real von_mises = std::sqrt(1.5) * dev.l2norm();
const Real yield =
von_mises - IsotropicHardening::hardening(p, h, sigma_0);
if (yield > 0) {
const Real dp = yield / (3 * elasticity.mu + h);
// Applying tangent from Bonnet & Frangi, p.175
const Real beta = 3 * elasticity.mu * dp / von_mises;
const Real gamma = 3 * elasticity.mu / (3 * elasticity.mu + h);
const Real dot = dev.dot(in);
dev *= 3 * elasticity.mu * (gamma - beta) * dot /
(von_mises * von_mises);
out.deviatoric(in);
out *= 2 * elasticity.mu * beta;
out += dev;
} else
out = 0;
},
range<Mat>(output), range<CMat>(input), range<CMat>(strain),
range<CMat>(strain_increment), range<CMat>(*this->plastic_strain),
*this->cumulated_plastic_strain);
}
} // namespace tamaas
Event Timeline
Log In to Comment