diff --git a/packages/solid_mechanics.cmake b/packages/solid_mechanics.cmake index 27a98df70..f56de607e 100644 --- a/packages/solid_mechanics.cmake +++ b/packages/solid_mechanics.cmake @@ -1,129 +1,131 @@ #=============================================================================== # @file solid_mechanics.cmake # # @author Guillaume Anciaux # @author Nicolas Richart # # @date creation: Mon Nov 21 2011 # @date last modification: Mon Jan 18 2016 # # @brief package description for core # # @section LICENSE # # Copyright (©) 2010-2012, 2014, 2015 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 . # #=============================================================================== package_declare(solid_mechanics DEFAULT ON DESCRIPTION "Solid mechanics model" DEPENDS core ) package_declare_sources(solid_mechanics model/solid_mechanics/material.cc model/solid_mechanics/material.hh model/solid_mechanics/material_inline_impl.cc model/solid_mechanics/material_selector.hh model/solid_mechanics/material_selector_tmpl.hh model/solid_mechanics/materials/internal_field.hh model/solid_mechanics/materials/internal_field_tmpl.hh model/solid_mechanics/materials/random_internal_field.hh model/solid_mechanics/materials/random_internal_field_tmpl.hh model/solid_mechanics/solid_mechanics_model.cc model/solid_mechanics/solid_mechanics_model.hh model/solid_mechanics/solid_mechanics_model_inline_impl.cc model/solid_mechanics/solid_mechanics_model_io.cc model/solid_mechanics/solid_mechanics_model_mass.cc model/solid_mechanics/solid_mechanics_model_material.cc model/solid_mechanics/solid_mechanics_model_tmpl.hh model/solid_mechanics/solid_mechanics_model_event_handler.hh model/solid_mechanics/materials/plane_stress_toolbox.hh model/solid_mechanics/materials/plane_stress_toolbox_tmpl.hh model/solid_mechanics/materials/material_core_includes.hh model/solid_mechanics/materials/material_elastic.cc model/solid_mechanics/materials/material_elastic.hh model/solid_mechanics/materials/material_elastic_inline_impl.cc model/solid_mechanics/materials/material_thermal.cc model/solid_mechanics/materials/material_thermal.hh model/solid_mechanics/materials/material_elastic_linear_anisotropic.cc model/solid_mechanics/materials/material_elastic_linear_anisotropic.hh model/solid_mechanics/materials/material_elastic_linear_anisotropic_inline_impl.cc model/solid_mechanics/materials/material_elastic_orthotropic.cc model/solid_mechanics/materials/material_elastic_orthotropic.hh model/solid_mechanics/materials/material_damage/material_damage.hh model/solid_mechanics/materials/material_damage/material_damage_tmpl.hh model/solid_mechanics/materials/material_damage/material_marigo.cc model/solid_mechanics/materials/material_damage/material_marigo.hh model/solid_mechanics/materials/material_damage/material_marigo_inline_impl.cc model/solid_mechanics/materials/material_damage/material_mazars.cc model/solid_mechanics/materials/material_damage/material_mazars.hh model/solid_mechanics/materials/material_damage/material_mazars_inline_impl.cc model/solid_mechanics/materials/material_finite_deformation/material_neohookean.cc model/solid_mechanics/materials/material_finite_deformation/material_neohookean.hh model/solid_mechanics/materials/material_finite_deformation/material_neohookean_inline_impl.cc model/solid_mechanics/materials/material_plastic/material_plastic.cc model/solid_mechanics/materials/material_plastic/material_plastic.hh model/solid_mechanics/materials/material_plastic/material_plastic_inline_impl.cc model/solid_mechanics/materials/material_plastic/material_linear_isotropic_hardening.cc model/solid_mechanics/materials/material_plastic/material_linear_isotropic_hardening.hh model/solid_mechanics/materials/material_plastic/material_linear_isotropic_hardening_inline_impl.cc model/solid_mechanics/materials/material_viscoelastic/material_standard_linear_solid_deviatoric.cc model/solid_mechanics/materials/material_viscoelastic/material_standard_linear_solid_deviatoric.hh + model/solid_mechanics/materials/material_viscoelastic/material_viscoelastic_maxwell.cc + model/solid_mechanics/materials/material_viscoelastic/material_viscoelastic_maxwell.hh model/solid_mechanics/materials/material_non_local.hh model/solid_mechanics/materials/material_non_local_tmpl.hh model/solid_mechanics/materials/material_non_local_includes.hh ) package_declare_material_infos(solid_mechanics LIST AKANTU_CORE_MATERIAL_LIST INCLUDE material_core_includes.hh ) package_declare_documentation_files(solid_mechanics manual-solidmechanicsmodel.tex manual-constitutive-laws.tex manual-lumping.tex manual-appendix-materials.tex figures/dynamic_analysis.png figures/explicit_dynamic.pdf figures/explicit_dynamic.svg figures/static.pdf figures/static.svg figures/hooke_law.pdf figures/implicit_dynamic.pdf figures/implicit_dynamic.svg figures/problemDomain.pdf_tex figures/problemDomain.pdf figures/static_analysis.png figures/stress_strain_el.pdf figures/tangent.pdf figures/tangent.svg figures/stress_strain_neo.pdf figures/visco_elastic_law.pdf figures/isotropic_hardening_plasticity.pdf figures/stress_strain_visco.pdf ) package_declare_extra_files_to_package(solid_mechanics SOURCES model/solid_mechanics/material_list.hh.in ) diff --git a/src/model/solid_mechanics/materials/material_viscoelastic/material_viscoelastic_maxwell.cc b/src/model/solid_mechanics/materials/material_viscoelastic/material_viscoelastic_maxwell.cc new file mode 100644 index 000000000..62e387649 --- /dev/null +++ b/src/model/solid_mechanics/materials/material_viscoelastic/material_viscoelastic_maxwell.cc @@ -0,0 +1,568 @@ +/** + * @file material_viscoelastic_maxwell.hh + * + * @author Emil Gallyamov + * + * @date creation: Tue May 08 2018 + * @date last modification: Tue May 08 2018 + * + * @brief Material Visco-elastic, based on Maxwell chain, + * see + * [] R. de Borst and A.H. van den Boogaard "Finite-element modeling of + * deformation and cracking in early-age concrete", J.Eng.Mech., 1994 + * as well as + * [] Manual of DIANA FEA Theory manual v.10.2 Section 37.6 + * + * @section LICENSE + * + * Copyright (©) 2010-2018 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 . + * + */ + +/* -------------------------------------------------------------------------- */ +#include "material_viscoelastic_maxwell.hh" +#include "solid_mechanics_model.hh" + +namespace akantu { + +/* -------------------------------------------------------------------------- */ +template +MaterialViscoelasticMaxwell::MaterialViscoelasticMaxwell( + SolidMechanicsModel & model, const ID & id) + : MaterialElastic(model, id), + C(voigt_h::size, voigt_h::size), sigma_v("sigma_v", *this), + dissipated_energy("dissipated_energy", *this) { + + AKANTU_DEBUG_IN(); + + this->registerParam("Eta1", Eta1, Real(1.), _pat_parsable | _pat_modifiable, + "Viscosity"); + this->registerParam("Ev1", Ev1, Real(1.), _pat_parsable | _pat_modifiable, + "Stiffness of the viscous element"); + this->registerParam("Einf", Einf, Real(1.), _pat_parsable | _pat_modifiable, + "Stiffness of the elastic element"); + this->registerParam("static_dt", static_dt, Real(0.), + _pat_parsable | _pat_modifiable, + "fictitious time step for static computation"), + this->registerParam("previous_dt", previous_dt, Real(0.), _pat_readable, + "Time step of previous solveStep"); + UInt stress_size = spatial_dimension * spatial_dimension; + + this->use_previous_stress = true; + this->use_previous_gradu = true; + this->use_previous_stress_thermal = true; + + this->sigma_v.initialize(stress_size); + this->dissipated_energy.initialize(1); + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +template +void MaterialViscoelasticMaxwell::initMaterial() { + AKANTU_DEBUG_IN(); + + // this->E = Einf + Ev1; + this->E = std::min(this->Einf, this->Ev1); + MaterialElastic::initMaterial(); + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +template +void MaterialViscoelasticMaxwell< + spatial_dimension>::updateInternalParameters() { + MaterialElastic::updateInternalParameters(); + + Real pre_mult = 1 / (1 + this->nu) / (1 - 2 * this->nu); + UInt n = voigt_h::size; + Real Miiii = pre_mult * (1 - this->nu); + Real Miijj = pre_mult * this->nu; + Real Mijij = pre_mult * 0.5 * (1 - 2 * this->nu); + + if (spatial_dimension == 1) { + C(0, 0) = 1; + } else { + C(0, 0) = Miiii; + } + if (spatial_dimension >= 2) { + C(1, 1) = Miiii; + C(0, 1) = Miijj; + C(1, 0) = Miijj; + C(n - 1, n - 1) = Mijij; + } + + if (spatial_dimension == 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; + } +} + +/* -------------------------------------------------------------------------- */ +template <> void MaterialViscoelasticMaxwell<2>::updateInternalParameters() { + MaterialElastic<2>::updateInternalParameters(); + + Real pre_mult = 1 / (1 + this->nu) / (1 - 2 * this->nu); + UInt n = voigt_h::size; + Real Miiii = pre_mult * (1 - this->nu); + Real Miijj = pre_mult * this->nu; + Real Mijij = pre_mult * 0.5 * (1 - 2 * this->nu); + + C(0, 0) = Miiii; + C(1, 1) = Miiii; + C(0, 1) = Miijj; + C(1, 0) = Miijj; + C(n - 1, n - 1) = Mijij; +} + +/* -------------------------------------------------------------------------- */ +template +void MaterialViscoelasticMaxwell::setToSteadyState( + ElementType el_type, GhostType ghost_type) { + AKANTU_DEBUG_IN(); + + /* Array & stress_dev_vect = stress_dev(el_type, ghost_type); + Array & history_int_vect = history_integral(el_type, ghost_type); + + Array::matrix_iterator stress_d = + stress_dev_vect.begin(spatial_dimension, spatial_dimension); + Array::matrix_iterator history_int = + history_int_vect.begin(spatial_dimension, spatial_dimension); + + /// Loop on all quadrature points + MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); + + Matrix & dev_s = *stress_d; + Matrix & h = *history_int; + + /// Compute the first invariant of strain + Real Theta = grad_u.trace(); + + for (UInt i = 0; i < spatial_dimension; ++i) + for (UInt j = 0; j < spatial_dimension; ++j) { + dev_s(i, j) = 2 * this->mu * (.5 * (grad_u(i, j) + grad_u(j, i)) - + 1. / 3. * Theta * (i == j)); + h(i, j) = 0.; + } + + /// Save the deviator of stress + ++stress_d; + ++history_int; + + MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; + */ + AKANTU_DEBUG_OUT(); +} +/* -------------------------------------------------------------------------- */ +template +void MaterialViscoelasticMaxwell::computeStress( + ElementType el_type, GhostType ghost_type) { + AKANTU_DEBUG_IN(); + + MaterialThermal::computeStress(el_type, ghost_type); + + auto sigma_th_it = this->sigma_th(el_type, ghost_type).begin(); + + auto previous_sigma_th_it = + this->sigma_th.previous(el_type, ghost_type).begin(); + + auto previous_gradu_it = this->gradu.previous(el_type, ghost_type) + .begin(spatial_dimension, spatial_dimension); + + auto previous_stress_it = this->stress.previous(el_type, ghost_type) + .begin(spatial_dimension, spatial_dimension); + + auto sigma_v_it = this->sigma_v(el_type, ghost_type) + .begin(spatial_dimension, spatial_dimension); + + MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); + + auto & previous_grad_u = *previous_gradu_it; + auto & previous_sigma = *previous_stress_it; + + computeStressOnQuad(grad_u, previous_grad_u, sigma, previous_sigma, + *sigma_v_it, *sigma_th_it, *previous_sigma_th_it); + ++sigma_th_it; + ++previous_sigma_th_it; + ++previous_stress_it; + ++previous_gradu_it; + ++sigma_v_it; + + MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +template +void MaterialViscoelasticMaxwell::computeStressOnQuad( + const Matrix & grad_u, const Matrix & previous_grad_u, + Matrix & sigma, const Matrix & previous_sigma, + Matrix & sigma_v, const Real & sigma_th, + const Real & previous_sigma_th) { + + Matrix delta_sigma(spatial_dimension, spatial_dimension); + Matrix grad_delta_u(grad_u); + grad_delta_u -= previous_grad_u; + Real delta_sigma_th = sigma_th - previous_sigma_th; + + Real dt = this->model.getTimeStep(); + Real lambda1 = this->Eta1 / this->Ev1; + Real exp_dt_lambda = exp(-dt / lambda1); + Real E_ef = this->Einf + (1 - exp_dt_lambda) * this->Ev1 * lambda1 / dt; + + // Wikipedia convention: + // 2*eps_ij (i!=j) = voigt_eps_I + // http://en.wikipedia.org/wiki/Voigt_notation + Vector voigt_delta_strain(voigt_h::size); + Vector voigt_delta_stress(voigt_h::size); + Vector voigt_sigma_v(voigt_h::size); + + for (UInt I = 0; I < voigt_h::size; ++I) { + Real voigt_factor = voigt_h::factors[I]; + UInt i = voigt_h::vec[I][0]; + UInt j = voigt_h::vec[I][1]; + + voigt_delta_strain(I) = + voigt_factor * (grad_delta_u(i, j) + grad_delta_u(j, i)) / 2.; + + voigt_sigma_v(I) = sigma_v(i, j); + } + + voigt_delta_stress = + E_ef * this->C * voigt_delta_strain - (1 - exp_dt_lambda) * voigt_sigma_v; + + for (UInt I = 0; I < voigt_h::size; ++I) { + UInt i = voigt_h::vec[I][0]; + UInt j = voigt_h::vec[I][1]; + + delta_sigma(i, j) = delta_sigma(j, i) = + voigt_delta_stress(I) + (i == j) * delta_sigma_th; + } + /* +Real trace = grad_delta_u.trace(); // trace = (\nabla u)_{kk} + +for (UInt i = 0; i < spatial_dimension; ++i) { + for (UInt j = 0; j < spatial_dimension; ++j) { + delta_sigma(i, j) = + (i == j) * E_ef * this->lambda * trace + + E_ef * this->mu * (grad_delta_u(i, j) + grad_delta_u(j, i)) + + (i == j) * sigma_th; + } +} + */ + sigma = previous_sigma + delta_sigma; +} + +/* -------------------------------------------------------------------------- */ +template +void MaterialViscoelasticMaxwell::afterSolveStep() { + + Material::afterSolveStep(); + + for (auto & el_type : this->element_filter.elementTypes( + _all_dimensions, _not_ghost, _ek_not_defined)) { + auto previous_gradu_it = this->gradu.previous(el_type, _not_ghost) + .begin(spatial_dimension, spatial_dimension); + + auto sigma_v_it = this->sigma_v(el_type, _not_ghost) + .begin(spatial_dimension, spatial_dimension); + + MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, _not_ghost); + + auto & previous_grad_u = *previous_gradu_it; + + updateIntVarOnQuad(grad_u, previous_grad_u, *sigma_v_it); + + ++previous_gradu_it; + ++sigma_v_it; + + MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; + this->updateDissipatedEnergy(el_type, _not_ghost); + } +} +/* -------------------------------------------------------------------------- */ +template +void MaterialViscoelasticMaxwell::updateIntVarOnQuad( + const Matrix & grad_u, const Matrix & previous_grad_u, + Matrix & sigma_v) { + + Matrix grad_delta_u(grad_u); + grad_delta_u -= previous_grad_u; + + Real dt = this->model.getTimeStep(); + Real lambda1 = Eta1 / Ev1; + Real exp_dt_lambda = exp(-dt / lambda1); + Real E_ef_v = (1 - exp_dt_lambda) * this->Ev1 * lambda1 / dt; + + // Wikipedia convention: + // 2*eps_ij (i!=j) = voigt_eps_I + // http://en.wikipedia.org/wiki/Voigt_notation + Vector voigt_delta_strain(voigt_h::size); + Vector voigt_sigma_v(voigt_h::size); + + for (UInt I = 0; I < voigt_h::size; ++I) { + Real voigt_factor = voigt_h::factors[I]; + UInt i = voigt_h::vec[I][0]; + UInt j = voigt_h::vec[I][1]; + + voigt_delta_strain(I) = + voigt_factor * (grad_delta_u(i, j) + grad_delta_u(j, i)) / 2.; + + voigt_sigma_v(I) = sigma_v(i, j); + } + + voigt_sigma_v = + exp_dt_lambda * voigt_sigma_v + E_ef_v * this->C * voigt_delta_strain; + + for (UInt I = 0; I < voigt_h::size; ++I) { + UInt i = voigt_h::vec[I][0]; + UInt j = voigt_h::vec[I][1]; + + sigma_v(i, j) = sigma_v(j, i) = voigt_sigma_v(I); + } +} + +/* -------------------------------------------------------------------------- */ +template +void MaterialViscoelasticMaxwell::computeTangentModuli( + const ElementType & el_type, Array & tangent_matrix, + GhostType ghost_type) { + AKANTU_DEBUG_IN(); + + Real dt = this->model.getTimeStep(); + Real lambda1 = this->Eta1 / this->Ev1; + Real exp_dt_lambda = exp(-dt / lambda1); + Real E_ef = this->Einf + (1 - exp_dt_lambda) * this->Ev1 * lambda1 / dt; + + this->previous_dt = dt; + + MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_BEGIN(tangent_matrix); + this->computeTangentModuliOnQuad(tangent); + MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_END; + + tangent_matrix *= E_ef; + + this->was_stiffness_assembled = true; + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +template +void MaterialViscoelasticMaxwell::computeTangentModuliOnQuad( + Matrix & tangent) { + + tangent.copy(C); +} + +/* -------------------------------------------------------------------------- */ +template +void MaterialViscoelasticMaxwell::savePreviousState() { + + for (auto & el_type : this->element_filter.elementTypes( + _all_dimensions, _not_ghost, _ek_not_defined)) { + + auto sigma_th_it = this->sigma_th(el_type, _not_ghost).begin(); + + auto previous_sigma_th_it = + this->sigma_th.previous(el_type, _not_ghost).begin(); + + auto previous_gradu_it = this->gradu.previous(el_type, _not_ghost) + .begin(spatial_dimension, spatial_dimension); + + auto previous_sigma_it = this->stress.previous(el_type, _not_ghost) + .begin(spatial_dimension, spatial_dimension); + + auto sigma_v_it = this->sigma_v(el_type, _not_ghost) + .begin(spatial_dimension, spatial_dimension); + + MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, _not_ghost); + auto & previous_grad_u = *previous_gradu_it; + auto & previous_sigma = *previous_sigma_it; + + previous_grad_u.copy(*gradu_it); + previous_sigma.copy(*stress_it); + *previous_sigma_th_it = *sigma_th_it; + + ++previous_gradu_it, ++previous_sigma_it, ++previous_sigma_th_it, + ++sigma_v_it, ++sigma_th_it; + MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; + } +} + +/* -------------------------------------------------------------------------- */ +template +void MaterialViscoelasticMaxwell::updateIntVariables() { + + for (auto & el_type : this->element_filter.elementTypes( + _all_dimensions, _not_ghost, _ek_not_defined)) { + + auto previous_gradu_it = this->gradu.previous(el_type, _not_ghost) + .begin(spatial_dimension, spatial_dimension); + auto previous_sigma_it = this->stress.previous(el_type, _not_ghost) + .begin(spatial_dimension, spatial_dimension); + + auto sigma_v_it = this->sigma_v(el_type, _not_ghost) + .begin(spatial_dimension, spatial_dimension); + + MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, _not_ghost); + auto & previous_grad_u = *previous_gradu_it; + + this->updateIntVarOnQuad(grad_u, previous_grad_u, *sigma_v_it); + + ++previous_gradu_it, ++sigma_v_it; + MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; + } +} + +/* -------------------------------------------------------------------------- */ +template +void MaterialViscoelasticMaxwell::updateDissipatedEnergy( + ElementType el_type, GhostType ghost_type) { + AKANTU_DEBUG_IN(); + + // if(ghost_type == _ghost) return 0.; + + Real * dis_energy = dissipated_energy(el_type, ghost_type).storage(); + + /* Array & sigma_vec = stress_dev(el_type, ghost_type); + Array & history_int_vect = history_integral(el_type, ghost_type); + + Array::matrix_iterator stress_d = + stress_dev_vect.begin(spatial_dimension, spatial_dimension); + Array::matrix_iterator history_int = + history_int_vect.begin(spatial_dimension, spatial_dimension); + + Matrix q(spatial_dimension, spatial_dimension); + Matrix q_rate(spatial_dimension, spatial_dimension); + Matrix epsilon_d(spatial_dimension, spatial_dimension); + Matrix epsilon_v(spatial_dimension, spatial_dimension); + + Real dt = this->model.getTimeStep(); + + Real gamma_v = Ev / this->E; + Real alpha = 1. / (2. * this->mu * gamma_v); + + /// Loop on all quadrature points + MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); + + Matrix & dev_s = *stress_d; + Matrix & h = *history_int; + + /// Compute the first invariant of strain + this->template gradUToEpsilon(grad_u, epsilon_d); + + Real Theta = epsilon_d.trace(); + epsilon_v.eye(Theta / Real(3.)); + epsilon_d -= epsilon_v; + + q.copy(dev_s); + q -= h; + q *= gamma_v; + + q_rate.copy(dev_s); + q_rate *= gamma_v; + q_rate -= q; + q_rate /= tau; + + for (UInt i = 0; i < spatial_dimension; ++i) + for (UInt j = 0; j < spatial_dimension; ++j) + *dis_energy += (epsilon_d(i, j) - alpha * q(i, j)) * q_rate(i, j) * dt; + + /// Save the deviator of stress + ++stress_d; + ++history_int; + ++dis_energy; + + MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; + */ + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +template +Real MaterialViscoelasticMaxwell::getDissipatedEnergy() + const { + AKANTU_DEBUG_IN(); + + Real de = 0.; + + /// integrate the dissipated energy for each type of elements + for (auto & type : + this->element_filter.elementTypes(spatial_dimension, _not_ghost)) { + de += + this->fem.integrate(dissipated_energy(type, _not_ghost), type, + _not_ghost, this->element_filter(type, _not_ghost)); + } + + AKANTU_DEBUG_OUT(); + return de; +} + +/* -------------------------------------------------------------------------- */ +template +Real MaterialViscoelasticMaxwell::getDissipatedEnergy( + ElementType type, UInt index) const { + AKANTU_DEBUG_IN(); + + UInt nb_quadrature_points = this->fem.getNbIntegrationPoints(type); + auto it = + this->dissipated_energy(type, _not_ghost).begin(nb_quadrature_points); + UInt gindex = (this->element_filter(type, _not_ghost))(index); + + AKANTU_DEBUG_OUT(); + return this->fem.integrate(it[index], type, gindex); +} + +/* -------------------------------------------------------------------------- */ +template +Real MaterialViscoelasticMaxwell::getEnergy( + const std::string & type) { + if (type == "dissipated") + return getDissipatedEnergy(); + else if (type == "viscoelastic_maxwell") + return getDissipatedEnergy(); + else + return MaterialElastic::getEnergy(type); +} + +/* -------------------------------------------------------------------------- */ +template +Real MaterialViscoelasticMaxwell::getEnergy( + const std::string & energy_id, ElementType type, UInt index) { + if (energy_id == "dissipated") + return getDissipatedEnergy(type, index); + else if (energy_id == "viscoelastic_maxwell") + return getDissipatedEnergy(type, index); + else + return MaterialElastic::getEnergy(energy_id, type, + index); +} + +/* -------------------------------------------------------------------------- */ + +INSTANTIATE_MATERIAL(viscoelastic_maxwell, MaterialViscoelasticMaxwell); + +} // namespace akantu diff --git a/src/model/solid_mechanics/materials/material_viscoelastic/material_viscoelastic_maxwell.hh b/src/model/solid_mechanics/materials/material_viscoelastic/material_viscoelastic_maxwell.hh new file mode 100644 index 000000000..a81ff66c3 --- /dev/null +++ b/src/model/solid_mechanics/materials/material_viscoelastic/material_viscoelastic_maxwell.hh @@ -0,0 +1,194 @@ +/** + * @file material_viscoelastic_maxwell.hh + * + * @author Emil Gallyamov + * + * @date creation: Tue May 08 2018 + * @date last modification: Tue May 08 2018 + * + * @brief Material Visco-elastic, based on Maxwell chain, + * see + * [] R. de Borst and A.H. van den Boogaard "Finite-element modeling of + * deformation and cracking in early-age concrete", J.Eng.Mech., 1994 + * as well as + * [] Manual of DIANA FEA Theory manual v.10.2 Section 37.6 + * + * @section LICENSE + * + * Copyright (©) 2010-2018 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 . + * + */ + +/* -------------------------------------------------------------------------- */ +#include "aka_common.hh" +#include "aka_voigthelper.hh" +#include "material_elastic.hh" +/* -------------------------------------------------------------------------- */ + +#ifndef __AKANTU_MATERIAL_VISCOELASTIC_MAXWELL_HH__ +#define __AKANTU_MATERIAL_VISCOELASTIC_MAXWELL_HH__ + +namespace akantu { + +/** + * Material Viscoelastic based on Maxwell chain + * + * + * @verbatim + + E_0 + ------|\/\/\|------- + | | + ---| |--- + | | + ----|\/\/\|--[|----- + | E_v1 \lambda_1| + ---| |--- + | | + ----|\/\/\|--[|----- + | E_v2 \lambda_2 | + ---| |--- + | | + ----|\/\/\|--[|---- + E_vN \lambda_N + + @endverbatim + * + * keyword : viscoelastic_maxwell + * + * parameters in the material files : + * - N : number of Maxwell elements + * - Einf : one spring element stiffness + * - Ev1 : stiffness of the 1st viscous element + * - lambda1: relaxation time of the 1st Maxwell element + * ... + * - Ev : stiffness of the Nst viscous element + * - lambda: relaxation time of the Nst Maxwell element + */ + +template +class MaterialViscoelasticMaxwell : public MaterialElastic { + /* ------------------------------------------------------------------------ */ + /* Constructors/Destructors */ + /* ------------------------------------------------------------------------ */ +public: + MaterialViscoelasticMaxwell(SolidMechanicsModel & model, const ID & id = ""); + ~MaterialViscoelasticMaxwell() override = default; + + /* ------------------------------------------------------------------------ */ + /* Methods */ + /* ------------------------------------------------------------------------ */ +public: + /// initialize the material computed parameter + void initMaterial() override; + + /// recompute the lame coefficient and effective tangent moduli + void updateInternalParameters() override; + + /// update internal variable on a converged Newton + void afterSolveStep() override; + + /// update internal variable based on previous and current strain values + void updateIntVariables(); + + /// update the internal variable sigma_v on quadrature point + void updateIntVarOnQuad(const Matrix & grad_u, + const Matrix & previous_grad_u, + Matrix & sigma_v); + + /// set material to steady state + void setToSteadyState(ElementType el_type, + GhostType ghost_type = _not_ghost) override; + + /// constitutive law for all element of a type + void computeStress(ElementType el_type, + GhostType ghost_type = _not_ghost) override; + + /// compute the tangent stiffness matrix for an element type + void computeTangentModuli(const ElementType & el_type, + Array & tangent_matrix, + GhostType ghost_type = _not_ghost) override; + + /// save previous stress and strain values into "previous" arrays + void savePreviousState() override; + +protected: + /// update the dissipated energy, is called after the stress have been + /// computed + void updateDissipatedEnergy(ElementType el_type, GhostType ghost_type); + + /// compute stresses on a quadrature point + void computeStressOnQuad(const Matrix & grad_u, + const Matrix & previous_grad_u, + Matrix & sigma, + const Matrix & previous_sigma, + Matrix & sigma_v, const Real & sigma_th, + const Real & previous_sigma_th); + + /// compute tangent moduli on a quadrature point + void computeTangentModuliOnQuad(Matrix & tangent); + + bool hasStiffnessMatrixChanged() override { + + Real dt = this->model.getTimeStep(); + + return ((this->previous_dt == dt) + ? (!(this->previous_dt == dt)) * (this->was_stiffness_assembled) + : (!(this->previous_dt == dt))); + // return (!(this->previous_dt == dt)); + } + + /* ------------------------------------------------------------------------ */ + /* Accessors */ + /* ------------------------------------------------------------------------ */ +public: + /// give the dissipated energy for the time step + Real getDissipatedEnergy() const; + Real getDissipatedEnergy(ElementType type, UInt index) const; + + /// get the energy using an energy type string for the time step + Real getEnergy(const std::string & type) override; + Real getEnergy(const std::string & energy_id, ElementType type, + UInt index) override; + /* ------------------------------------------------------------------------ */ + /* Class Members */ + /* ------------------------------------------------------------------------ */ +private: + using voigt_h = VoigtHelper; + + /// viscosity, viscous elastic modulus, one spring element elastic modulus + Real Eta1, Ev1, Einf; + + /// fictitious time step for static computations + Real static_dt; + + /// time step from previous solveStep + Real previous_dt; + + /// Effective viscoelastic stiffness tensor in voigt notation + Matrix C; + + /// Internal variable: viscous_stress + InternalField sigma_v; + + /// Dissipated energy + InternalField dissipated_energy; +}; + +} // namespace akantu + +#endif /* __AKANTU_MATERIAL_VISCOELASTIC_MAXWELL_HH__ */