Page MenuHomec4science
Files F121868522

material_viscoelastic_maxwell.cc
No OneTemporary
Actions

File Metadata

Created
Mon, Jul 14, 11:52

material_viscoelastic_maxwell.cc
View Options

/**
* Copyright (©) 2018-2023 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* This file is part of Akantu
*
* 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 "material_viscoelastic_maxwell.hh"
#include "solid_mechanics_model.hh"
namespace akantu {
/* -------------------------------------------------------------------------- */
template <Int dim>
MaterialViscoelasticMaxwell<dim>::MaterialViscoelasticMaxwell(
SolidMechanicsModel & model, const ID & id)
: MaterialElastic<dim>(model, id),
dissipated_energy(this->registerInternal("dissipated_energy", 1)),
mechanical_work(this->registerInternal("mechanical_work", 1)) {
this->registerParam("Einf", Einf, Real(1.), _pat_parsable | _pat_modifiable,
"Stiffness of the elastic element");
this->registerParam("previous_dt", previous_dt, Real(0.), _pat_readable,
"Time step of previous solveStep");
this->registerParam("Eta", Eta, _pat_parsable | _pat_modifiable,
"Viscosity of a Maxwell element");
this->registerParam("Ev", Ev, _pat_parsable | _pat_modifiable,
"Stiffness of a Maxwell element");
}
/* -------------------------------------------------------------------------- */
template <Int dim> void MaterialViscoelasticMaxwell<dim>::initMaterial() {
AKANTU_DEBUG_IN();
this->E = Einf + Ev.lpNorm<1>();
MaterialElastic<dim>::initMaterial();
AKANTU_DEBUG_ASSERT(this->Eta.size() == this->Ev.size(),
"Eta and Ev have different dimensions! Please correct.");
AKANTU_DEBUG_ASSERT(
!this->finite_deformation,
"Current material works only in infinitesimal deformations.");
auto stress_size = dim * dim;
this->registerInternal("sigma_v", stress_size * this->Ev.size());
this->registerInternal("epsilon_v", stress_size * this->Ev.size());
this->sigma_v = this->getSharedPtrInternal("sigma_v");
this->epsilon_v = this->getSharedPtrInternal("epsilon_v");
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <Int dim>
void MaterialViscoelasticMaxwell<dim>::updateInternalParameters() {
MaterialElastic<dim>::updateInternalParameters();
[[maybe_unused]] auto pre_mult = 1 / (1 + this->nu) / (1 - 2 * this->nu);
[[maybe_unused]] auto Miiii = pre_mult * (1 - this->nu);
[[maybe_unused]] auto Miijj = pre_mult * this->nu;
[[maybe_unused]] auto Mijij = pre_mult * 0.5 * (1 - 2 * this->nu);
[[maybe_unused]] auto Diiii = 1;
[[maybe_unused]] auto Diijj = -this->nu;
[[maybe_unused]] auto Dijij = (2 + 2 * this->nu);
C.zero();
D.zero();
if constexpr (dim == 1) {
C(0, 0) = 1;
D(0, 0) = 1;
} else {
C(0, 0) = Miiii;
D(0, 0) = Diiii;
}
if constexpr (dim >= 2) {
auto n = voigt_h::size;
C(1, 1) = Miiii;
C(0, 1) = Miijj;
C(1, 0) = Miijj;
C(n - 1, n - 1) = Mijij;
D(1, 1) = Diiii;
D(0, 1) = Diijj;
D(1, 0) = Diijj;
D(n - 1, n - 1) = Dijij;
}
if constexpr (dim == 3) {
C(2, 2) = Miiii;
C(0, 2) = Miijj;
C(1, 2) = Miijj;
C(2, 0) = Miijj;
C(2, 1) = Miijj;
C(3, 3) = Mijij;
C(4, 4) = Mijij;
D(2, 2) = Diiii;
D(0, 2) = Diijj;
D(1, 2) = Diijj;
D(2, 0) = Diijj;
D(2, 1) = Diijj;
D(3, 3) = Dijij;
D(4, 4) = Dijij;
}
}
/* -------------------------------------------------------------------------- */
template <Int dim>
void MaterialViscoelasticMaxwell<dim>::computeStress(ElementType el_type,
GhostType ghost_type) {
// NOLINTNEXTLINE(bugprone-parent-virtual-call)
Parent::computeStress(el_type, ghost_type);
for (auto && args : getArguments(el_type, ghost_type)) {
computeStressOnQuad(args);
}
}
/* -------------------------------------------------------------------------- */
template <Int dim>
template <class Args>
void MaterialViscoelasticMaxwell<dim>::computeStressOnQuad(Args && args) {
Real dt = this->getModel().getTimeStep();
// Wikipedia convention:
// 2*eps_ij (i!=j) = voigt_eps_I
// http://en.wikipedia.org/wiki/Voigt_notation
auto voigt_current_strain = Material::strainToVoigt<dim>(
Material::gradUToEpsilon<dim>(args["grad_u"_n]));
auto voigt_previous_strain = Material::strainToVoigt<dim>(
Material::gradUToEpsilon<dim>(args["previous_grad_u"_n]));
Vector<Real, voigt_h::size> voigt_stress =
this->Einf * this->C * voigt_current_strain;
Vector<Real, voigt_h::size> stress =
this->C * (voigt_current_strain - voigt_previous_strain);
for (auto && [Eta, Ev, sigma] : zip(this->Eta, this->Ev, args["sigma_v"_n])) {
auto lambda = Eta / Ev;
auto exp_dt_lambda = exp(-dt / lambda);
Real E_additional{0.};
if (exp_dt_lambda == 1) {
E_additional = Ev;
} else {
E_additional = (1 - exp_dt_lambda) * Ev * lambda / dt;
}
voigt_stress += E_additional * stress +
exp_dt_lambda * Material::stressToVoigt<dim>(sigma);
}
auto && sigma = args["sigma"_n];
for (Int I = 0; I < voigt_h::size; ++I) {
auto && [i, j] = voigt_h::vec[I];
sigma(i, j) = sigma(j, i) =
voigt_stress(I) + Math::kronecker(i, j) * args["sigma_th"_n];
}
}
/* -------------------------------------------------------------------------- */
template <Int dim>
void MaterialViscoelasticMaxwell<dim>::computePotentialEnergy(
ElementType el_type) {
for (auto && [args, epot, epsilon_v] :
zip(getArguments(el_type), this->potential_energy(el_type),
make_view((*this->epsilon_v)(el_type), dim, dim, Eta.size()))) {
this->computePotentialEnergyOnQuad(args["grad_u"_n], epot,
args["sigma_v"_n], epsilon_v);
}
}
/* -------------------------------------------------------------------------- */
template <Int dim>
template <typename D1>
void MaterialViscoelasticMaxwell<dim>::computePotentialEnergyOnQuad(
const Eigen::MatrixBase<D1> & grad_u, Real & epot,
Tensor3Proxy<Real> & sigma_v, Tensor3Proxy<Real> & epsilon_v) {
auto voigt_strain =
Material::strainToVoigt<dim>(Material::gradUToEpsilon<dim>(grad_u));
Vector<Real, voigt_h::size> voigt_stress =
this->Einf * this->C * voigt_strain;
epot = 0.5 * voigt_stress.dot(voigt_strain);
for (Int k = 0; k < this->Eta.size(); ++k) {
epot += sigma_v(k).doubleDot(epsilon_v(k)) / 2.;
}
}
/* -------------------------------------------------------------------------- */
template <Int dim>
void MaterialViscoelasticMaxwell<dim>::afterSolveStep(bool converged) {
Material::afterSolveStep(converged);
if (not converged) {
return;
}
if (this->update_variable_flag) {
updateIntVariables();
}
for (const auto & el_type : this->getElementFilter().elementTypes(
_all_dimensions, _not_ghost, _ek_not_defined)) {
this->updateDissipatedEnergy(el_type);
}
}
/* -------------------------------------------------------------------------- */
template <Int dim>
template <typename D1, typename D2>
void MaterialViscoelasticMaxwell<dim>::updateIntVarOnQuad(
const Eigen::MatrixBase<D1> & grad_u,
const Eigen::MatrixBase<D2> & previous_grad_u, Tensor3Proxy<Real> & sigma_v,
Tensor3Proxy<Real> & epsilon_v) {
Matrix<Real, dim, dim> grad_delta_u = grad_u - previous_grad_u;
Real dt = this->getModel().getTimeStep();
auto voigt_delta_strain =
Material::strainToVoigt<dim>(Material::gradUToEpsilon<dim>(grad_delta_u));
for (Idx k = 0; k < this->Eta.size(); ++k) {
auto lambda = this->Eta(k) / this->Ev(k);
auto exp_dt_lambda = exp(-dt / lambda);
Real E_ef_v = this->Ev(k);
if (exp_dt_lambda != 1) {
E_ef_v *= (1 - exp_dt_lambda) * lambda / dt;
}
auto voigt_sigma_v = Material::stressToVoigt<dim>(sigma_v(k));
Vector<Real, voigt_h::size> voigt_epsilon_v =
exp_dt_lambda * voigt_sigma_v + E_ef_v * this->C * voigt_delta_strain;
voigt_epsilon_v = 1 / Ev(k) * this->D * voigt_sigma_v;
for (Int I = 0; I < voigt_h::size; ++I) {
auto && [i, j] = voigt_h::vec[I];
sigma_v(i, j, k) = sigma_v(j, i, k) = voigt_sigma_v(I);
epsilon_v(i, j, k) = epsilon_v(j, i, k) = voigt_epsilon_v(I);
}
}
}
/* -------------------------------------------------------------------------- */
template <Int dim>
void MaterialViscoelasticMaxwell<dim>::computeTangentModuli(
ElementType /*el_type*/, Array<Real> & tangent_matrix,
GhostType /*ghost_type*/) {
AKANTU_DEBUG_IN();
Real dt = this->getModel().getTimeStep();
Real E_ef = this->Einf;
for (Int k = 0; k < Eta.size(); ++k) {
Real lambda = this->Eta(k) / this->Ev(k);
Real exp_dt_lambda = exp(-dt / lambda);
if (exp_dt_lambda == 1) {
E_ef += this->Ev(k);
} else {
E_ef += (1 - exp_dt_lambda) * this->Ev(k) * lambda / dt;
}
}
this->previous_dt = dt;
const auto tangent_size = Material::getTangentStiffnessVoigtSize(dim);
for (auto && tangent :
make_view<tangent_size, tangent_size>(tangent_matrix)) {
this->computeTangentModuliOnQuad(tangent);
}
tangent_matrix *= E_ef;
this->was_stiffness_assembled = true;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <Int dim>
template <typename D1>
void MaterialViscoelasticMaxwell<dim>::computeTangentModuliOnQuad(
Eigen::MatrixBase<D1> & tangent) {
tangent = C;
}
/* -------------------------------------------------------------------------- */
template <Int dim> void MaterialViscoelasticMaxwell<dim>::updateIntVariables() {
for (const auto & el_type : this->getElementFilter().elementTypes()) {
for (auto && [args, epsilon_v] :
zip(getArguments(el_type),
make_view((*this->epsilon_v)(el_type), dim, dim, Eta.size()))) {
updateIntVarOnQuad(args["grad_u"_n], args["previous_grad_u"_n],
args["sigma_v"_n], epsilon_v);
}
}
}
/* -------------------------------------------------------------------------- */
template <Int dim>
void MaterialViscoelasticMaxwell<dim>::updateDissipatedEnergy(
ElementType el_type) {
this->computePotentialEnergy(el_type);
for (auto && [args, epot, edis, work] :
zip(getArguments(el_type), this->potential_energy(el_type),
this->dissipated_energy(el_type), this->mechanical_work(el_type))) {
updateDissipatedEnergyOnQuad(args["grad_u"_n], args["previous_grad_u"_n],
args["sigma"_n], args["previous_sigma"_n],
edis, work, epot);
}
}
/* -------------------------------------------------------------------------- */
template <Int dim>
template <typename D1, typename D2, typename D3, typename D4>
void MaterialViscoelasticMaxwell<dim>::updateDissipatedEnergyOnQuad(
const Eigen::MatrixBase<D1> & grad_u,
const Eigen::MatrixBase<D2> & previous_grad_u,
const Eigen::MatrixBase<D3> & sigma,
const Eigen::MatrixBase<D4> & previous_sigma, Real & dis_energy,
Real & mech_work, const Real & pot_energy) {
Real dt = this->getModel().getTimeStep();
auto && strain_rate = (grad_u - previous_grad_u) / dt;
auto && av_stress = (sigma + previous_sigma) / 2.;
mech_work += av_stress.doubleDot(strain_rate) * dt;
dis_energy = mech_work - pot_energy;
}
/* -------------------------------------------------------------------------- */
template <Int dim>
Real MaterialViscoelasticMaxwell<dim>::getDissipatedEnergy() const {
AKANTU_DEBUG_IN();
auto & fem = this->getFEEngine();
Real de = 0.;
/// integrate the dissipated energy for each type of elements
for (auto && type : this->getElementFilter().elementTypes(dim, _not_ghost)) {
de += fem.integrate(this->dissipated_energy(type, _not_ghost), type,
_not_ghost, this->getElementFilter(type, _not_ghost));
}
AKANTU_DEBUG_OUT();
return de;
}
/* -------------------------------------------------------------------------- */
template <Int dim>
Real MaterialViscoelasticMaxwell<dim>::getDissipatedEnergy(
const Element & element) const {
auto & fem = this->getFEEngine();
auto nb_quadrature_points = fem.getNbIntegrationPoints(element.type);
auto it =
make_view(this->dissipated_energy(element.type), nb_quadrature_points)
.begin();
auto mat_element = element;
mat_element.element = this->getElementFilter()(element);
return fem.integrate(it[element.element], mat_element);
}
/* -------------------------------------------------------------------------- */
template <Int dim>
Real MaterialViscoelasticMaxwell<dim>::getMechanicalWork() const {
auto & fem = this->getFEEngine();
Real mw = 0.;
/// integrate the dissipated energy for each type of elements
for (auto && type : this->getElementFilter().elementTypes(dim)) {
mw += fem.integrate(this->mechanical_work(type), type, _not_ghost,
this->getElementFilter(type, _not_ghost));
}
return mw;
}
/* -------------------------------------------------------------------------- */
template <Int dim>
Real MaterialViscoelasticMaxwell<dim>::getMechanicalWork(
const Element & element) const {
auto & fem = this->getFEEngine();
auto nb_quadrature_points = fem.getNbIntegrationPoints(element.type);
auto it = make_view(this->mechanical_work(element.type), nb_quadrature_points)
.begin();
auto mat_element = element;
mat_element.element = this->getElementFilter()(element);
return fem.integrate(it[element.element], mat_element);
}
/* -------------------------------------------------------------------------- */
template <Int dim>
Real MaterialViscoelasticMaxwell<dim>::getPotentialEnergy() const {
AKANTU_DEBUG_IN();
auto & fem = this->getFEEngine();
Real epot = 0.;
/// integrate the dissipated energy for each type of elements
for (auto && type : this->getElementFilter().elementTypes(dim, _not_ghost)) {
epot += fem.integrate(this->potential_energy(type, _not_ghost), type,
_not_ghost, this->getElementFilter(type, _not_ghost));
}
AKANTU_DEBUG_OUT();
return epot;
}
/* -------------------------------------------------------------------------- */
template <Int dim>
Real MaterialViscoelasticMaxwell<dim>::getPotentialEnergy(
const Element & element) const {
AKANTU_DEBUG_IN();
auto & fem = this->getFEEngine();
auto nb_quadrature_points = fem.getNbIntegrationPoints(element.type);
auto it =
make_view(this->potential_energy(element.type), nb_quadrature_points)
.begin();
auto mat_element = element;
mat_element.element = this->getElementFilter()(element);
AKANTU_DEBUG_OUT();
return fem.integrate(it[element.element], mat_element);
}
/* -------------------------------------------------------------------------- */
template <Int dim>
Real MaterialViscoelasticMaxwell<dim>::getEnergy(const std::string & type) {
if (type == "dissipated") {
return getDissipatedEnergy();
}
if (type == "potential") {
return getPotentialEnergy();
}
if (type == "work") {
return getMechanicalWork();
}
return MaterialElastic<dim>::getEnergy(type);
}
/* -------------------------------------------------------------------------- */
template <Int dim>
Real MaterialViscoelasticMaxwell<dim>::getEnergy(const std::string & energy_id,
const Element & element) {
if (energy_id == "dissipated") {
return getDissipatedEnergy(element);
}
if (energy_id == "potential") {
return getPotentialEnergy(element);
}
if (energy_id == "work") {
return getMechanicalWork(element);
}
return MaterialElastic<dim>::getEnergy(energy_id, element);
}
/* -------------------------------------------------------------------------- */
template <Int dim>
void MaterialViscoelasticMaxwell<dim>::forceUpdateVariable() {
update_variable_flag = true;
}
/* -------------------------------------------------------------------------- */
template <Int dim>
void MaterialViscoelasticMaxwell<dim>::forceNotUpdateVariable() {
update_variable_flag = false;
}
/* -------------------------------------------------------------------------- */
const bool material_is_allocated_viscoelastic_maxwell [[maybe_unused]] =
instantiateMaterial<MaterialViscoelasticMaxwell>("viscoelastic_maxwell");
} // namespace akantu

Event Timeline

c4science · Help