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 index 0de48b797..8ab5c4c7a 100644 --- a/src/model/solid_mechanics/materials/material_viscoelastic/material_viscoelastic_maxwell.cc +++ b/src/model/solid_mechanics/materials/material_viscoelastic/material_viscoelastic_maxwell.cc @@ -1,576 +1,735 @@ /** * @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) { + C(voigt_h::size, voigt_h::size), D(voigt_h::size, voigt_h::size), + sigma_v("sigma_v", *this), epsilon_v("epsilon_v", *this), + dissipated_energy("dissipated_energy", *this), + mechanical_work("mechanical_work", *this) { AKANTU_DEBUG_IN(); 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"); + this->update_variable_flag = true; this->use_previous_stress = true; this->use_previous_gradu = true; this->use_previous_stress_thermal = true; this->dissipated_energy.initialize(1); + this->mechanical_work.initialize(1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialViscoelasticMaxwell::initMaterial() { AKANTU_DEBUG_IN(); - // this->E = Einf + Ev1; - this->E = std::min(this->Einf, this->Ev(0)); + this->E = Einf + Ev.norm(); +// this->E = std::min(this->Einf, this->Ev(0)); MaterialElastic::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."); UInt stress_size = spatial_dimension * spatial_dimension; this->sigma_v.initialize(stress_size * this->Ev.size()); + this->epsilon_v.initialize(stress_size * this->Ev.size()); 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); + Real Diiii = 1; + Real Diijj = -this->nu; + Real Dijij = (2 + 2 * this->nu); + if (spatial_dimension == 1) { C(0, 0) = 1; + D(0, 0) = 1; } else { C(0, 0) = Miiii; + D(0, 0) = Diiii; } if (spatial_dimension >= 2) { 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 (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; + + 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 <> 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); + Real Diiii = 1; + Real Diijj = -this->nu; + Real Dijij = (2 + 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(); + D(0, 0) = Diiii; + D(1, 1) = Diiii; + D(0, 1) = Diijj; + D(1, 0) = Diijj; + D(n - 1, n - 1) = Dijij; } + /* -------------------------------------------------------------------------- */ 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, this->Eta.size()); 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); + computeStressOnQuad(grad_u, *previous_gradu_it, sigma, *sigma_v_it, + *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, - Tensor3 & sigma_v, const Real & sigma_th, - const Real & previous_sigma_th) { - - Matrix grad_delta_u(grad_u); - grad_delta_u -= previous_grad_u; - Real delta_sigma_th = sigma_th - previous_sigma_th; + Matrix & sigma, Tensor3 & sigma_v, const Real & sigma_th) { // Wikipedia convention: // 2*eps_ij (i!=j) = voigt_eps_I // http://en.wikipedia.org/wiki/Voigt_notation Vector voigt_current_strain(voigt_h::size); Vector voigt_previous_strain(voigt_h::size); Vector voigt_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_current_strain(I) = voigt_factor * (grad_u(i, j) + grad_u(j, i)) / 2.; voigt_previous_strain(I) = voigt_factor * (previous_grad_u(i, j) + previous_grad_u(j, i)) / 2.; } voigt_stress = this->Einf * this->C * voigt_current_strain; Real dt = this->model.getTimeStep(); for (UInt k = 0; k < Eta.size(); ++k) { Real lambda = this->Eta(k) / this->Ev(k); Real exp_dt_lambda = exp(-dt / lambda); - Real E_additional = (1 - exp_dt_lambda) * this->Ev(k) * lambda / dt; + Real E_additional; + + if (exp_dt_lambda == 1) { + E_additional = this->Ev(k); + } else { + E_additional = (1 - exp_dt_lambda) * this->Ev(k) * lambda / dt; + } for (UInt I = 0; I < voigt_h::size; ++I) { UInt i = voigt_h::vec[I][0]; UInt j = voigt_h::vec[I][1]; voigt_sigma_v(I) = sigma_v(i, j, k); } voigt_stress += E_additional * this->C * (voigt_current_strain - voigt_previous_strain) + 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]; sigma(i, j) = sigma(j, i) = voigt_stress(I) + (i == j) * sigma_th; } } +/* -------------------------------------------------------------------------- */ +template +void MaterialViscoelasticMaxwell::computePotentialEnergy( + ElementType el_type, GhostType ghost_type) { + AKANTU_DEBUG_IN(); + + MaterialThermal::computePotentialEnergy(el_type, + ghost_type); + + if (ghost_type != _not_ghost) + return; + + auto epot = this->potential_energy(el_type, ghost_type).begin(); + auto sigma_v_it = this->sigma_v(el_type, ghost_type) + .begin(spatial_dimension, spatial_dimension, this->Eta.size()); + auto epsilon_v_it = this->epsilon_v(el_type, ghost_type) + .begin(spatial_dimension, spatial_dimension, this->Eta.size()); + + + + MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); + + this->computePotentialEnergyOnQuad(grad_u, *epot, *sigma_v_it, *epsilon_v_it); + ++epot; + ++sigma_v_it; + ++epsilon_v_it; + + MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +template +void MaterialViscoelasticMaxwell::computePotentialEnergyOnQuad( + const Matrix & grad_u, Real & epot, + Tensor3 & sigma_v, Tensor3 & epsilon_v) { + + Vector voigt_strain(voigt_h::size); + Vector voigt_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_strain(I) = + voigt_factor * (grad_u(i, j) + grad_u(j, i)) / 2.; + } + + voigt_stress = this->Einf * this->C * voigt_strain; + epot = 0.5 * voigt_stress.dot(voigt_strain); + + for (UInt k = 0; k < this->Eta.size(); ++k) { + Matrix stress_v = sigma_v(k); + Matrix strain_v = epsilon_v(k); + epot += 0.5 * stress_v.doubleDot(strain_v); + } +} + /* -------------------------------------------------------------------------- */ 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); + if (this->update_variable_flag) { + 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, this->Eta.size()); + auto sigma_v_it = this->sigma_v(el_type, _not_ghost) + .begin(spatial_dimension, spatial_dimension, this->Eta.size()); - MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, _not_ghost); + auto epsilon_v_it = this->epsilon_v(el_type, _not_ghost) + .begin(spatial_dimension, spatial_dimension, this->Eta.size()); - auto & previous_grad_u = *previous_gradu_it; + MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, _not_ghost); - updateIntVarOnQuad(grad_u, previous_grad_u, *sigma_v_it); + updateIntVarOnQuad(grad_u, *previous_gradu_it, *sigma_v_it, *epsilon_v_it); - ++previous_gradu_it; - ++sigma_v_it; + ++previous_gradu_it; + ++sigma_v_it; + ++epsilon_v_it; - MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; + 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, - Tensor3 & sigma_v) { + Tensor3 & sigma_v, Tensor3 & epsilon_v) { Matrix grad_delta_u(grad_u); grad_delta_u -= previous_grad_u; Real dt = this->model.getTimeStep(); Vector voigt_delta_strain(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.; } - for (UInt k = 0; k < Eta.size(); ++k) { + for (UInt k = 0; k < this->Eta.size(); ++k) { Real lambda = this->Eta(k) / this->Ev(k); Real exp_dt_lambda = exp(-dt / lambda); - Real E_ef_v = (1 - exp_dt_lambda) * this->Ev(k) * lambda / dt; + Real E_ef_v; + + if (exp_dt_lambda == 1) { + E_ef_v = this->Ev(k); + } else { + E_ef_v = (1 - exp_dt_lambda) * this->Ev(k) * lambda / dt; + } Vector voigt_sigma_v(voigt_h::size); + Vector voigt_epsilon_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_sigma_v(I) = sigma_v(i, j, k); } voigt_sigma_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 (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, 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 void MaterialViscoelasticMaxwell::computeTangentModuli( const ElementType & el_type, Array & tangent_matrix, GhostType ghost_type) { AKANTU_DEBUG_IN(); Real dt = this->model.getTimeStep(); Real E_ef = this->Einf; for (UInt k = 0; k < Eta.size(); ++k) { Real lambda = this->Eta(k) / this->Ev(k); Real exp_dt_lambda = exp(-dt / lambda); - E_ef += (1 - exp_dt_lambda) * this->Ev(k) * lambda / dt; + 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; 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, this->Eta.size()); 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, this->Eta.size()); + auto epsilon_v_it = this->epsilon_v(el_type, _not_ghost) + .begin(spatial_dimension, spatial_dimension, this->Eta.size()); + 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); + updateIntVarOnQuad(grad_u, *previous_gradu_it, *sigma_v_it, *epsilon_v_it); + + ++previous_gradu_it; + ++sigma_v_it; + ++epsilon_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(); + if (ghost_type != _not_ghost) + return; - /* Array & sigma_vec = stress_dev(el_type, ghost_type); - Array & history_int_vect = history_integral(el_type, ghost_type); + this->computePotentialEnergy(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); + auto epot = this->potential_energy(el_type, ghost_type).begin(); + auto dis_energy = this->dissipated_energy(el_type, ghost_type).begin(); + auto mech_work = this->mechanical_work(el_type, ghost_type).begin(); + auto sigma_v_it = this->sigma_v(el_type, ghost_type) + .begin(spatial_dimension, spatial_dimension, this->Eta.size()); + auto epsilon_v_it = this->epsilon_v(el_type, ghost_type) + .begin(spatial_dimension, spatial_dimension, this->Eta.size()); + auto previous_gradu_it = this->gradu.previous(el_type, ghost_type) + .begin(spatial_dimension, spatial_dimension); + auto previous_sigma_it = this->stress.previous(el_type, ghost_type) + .begin(spatial_dimension, spatial_dimension); - Real dt = this->model.getTimeStep(); + /// Loop on all quadrature points + MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); - Real gamma_v = Ev / this->E; - Real alpha = 1. / (2. * this->mu * gamma_v); + updateDissipatedEnergyOnQuad(grad_u, *previous_gradu_it, sigma, + *previous_sigma_it, *dis_energy, *mech_work, *epot); + ++previous_gradu_it; + ++previous_sigma_it; + ++dis_energy; + ++mech_work; + ++epot; - /// Loop on all quadrature points - MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); + MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; - Matrix & dev_s = *stress_d; - Matrix & h = *history_int; + AKANTU_DEBUG_OUT(); +} - /// 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; +/* -------------------------------------------------------------------------- */ +template +void MaterialViscoelasticMaxwell::updateDissipatedEnergyOnQuad( + const Matrix & grad_u, const Matrix & previous_grad_u, + const Matrix & sigma, const Matrix & previous_sigma, + Real & dis_energy, Real & mech_work, const Real & pot_energy) { - q.copy(dev_s); - q -= h; - q *= gamma_v; + Real dt = this->model.getTimeStep(); - q_rate.copy(dev_s); - q_rate *= gamma_v; - q_rate -= q; - q_rate /= tau; + Matrix strain_rate = grad_u; + strain_rate -= previous_grad_u; + strain_rate /= dt; - 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; + Matrix av_stress = sigma; + av_stress += previous_sigma; + av_stress /= 2; - /// Save the deviator of stress - ++stress_d; - ++history_int; - ++dis_energy; + mech_work += av_stress.doubleDot(strain_rate) * dt; - MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; - */ - AKANTU_DEBUG_OUT(); + dis_energy = mech_work - pot_energy; } /* -------------------------------------------------------------------------- */ 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, + this->fem.integrate(this->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::getMechanicalWork() + const { + AKANTU_DEBUG_IN(); + + Real mw = 0.; + + /// integrate the dissipated energy for each type of elements + for (auto & type : + this->element_filter.elementTypes(spatial_dimension, _not_ghost)) { + mw += + this->fem.integrate(this->mechanical_work(type, _not_ghost), type, + _not_ghost, this->element_filter(type, _not_ghost)); + } + + AKANTU_DEBUG_OUT(); + return mw; +} + +/* -------------------------------------------------------------------------- */ +template +Real MaterialViscoelasticMaxwell::getMechanicalWork( + ElementType type, UInt index) const { + AKANTU_DEBUG_IN(); + + UInt nb_quadrature_points = this->fem.getNbIntegrationPoints(type); + auto it = + this->mechanical_work(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::getPotentialEnergy() + const { + AKANTU_DEBUG_IN(); + + Real epot = 0.; + + /// integrate the dissipated energy for each type of elements + for (auto & type : + this->element_filter.elementTypes(spatial_dimension, _not_ghost)) { + epot += + this->fem.integrate(this->potential_energy(type, _not_ghost), type, + _not_ghost, this->element_filter(type, _not_ghost)); + } + + AKANTU_DEBUG_OUT(); + return epot; +} + +/* -------------------------------------------------------------------------- */ +template +Real MaterialViscoelasticMaxwell::getPotentialEnergy( + ElementType type, UInt index) const { + AKANTU_DEBUG_IN(); + + UInt nb_quadrature_points = this->fem.getNbIntegrationPoints(type); + auto it = + this->potential_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 if (type == "potential") + return getPotentialEnergy(); + else if (type == "work") + return getMechanicalWork(); 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 if (energy_id == "potential") + return getPotentialEnergy(type, index); + else if (energy_id == "work") + return getMechanicalWork(type, index); else return MaterialElastic::getEnergy(energy_id, type, index); } +/* -------------------------------------------------------------------------- */ +template +void MaterialViscoelasticMaxwell::forceUpdateVariable() { + update_variable_flag = true; +} + + +/* -------------------------------------------------------------------------- */ +template +void MaterialViscoelasticMaxwell::forceNotUpdateVariable() { + update_variable_flag = false; +} + + /* -------------------------------------------------------------------------- */ 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 index f351d5098..9d8feaaea 100644 --- a/src/model/solid_mechanics/materials/material_viscoelastic/material_viscoelastic_maxwell.hh +++ b/src/model/solid_mechanics/materials/material_viscoelastic/material_viscoelastic_maxwell.hh @@ -1,193 +1,226 @@ /** * @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 \Eta_1| ---| |--- | | ----|\/\/\|--[|----- | E_v2 \Eta_2 | ---| |--- | | ----|\/\/\|--[|---- E_vN \Eta_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 * - Eta1: viscosity of the 1st Maxwell element * ... * - Ev : stiffness of the Nst viscous element * - Eta: viscosity 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, - Tensor3 & sigma_v); - - /// set material to steady state - void setToSteadyState(ElementType el_type, - GhostType ghost_type = _not_ghost) override; + Tensor3 & sigma_v, Tensor3 & epsilon_v); /// 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; + /// change flag of updateIntVar to true + void forceUpdateVariable(); + + /// change flag of updateIntVar to false + void forceNotUpdateVariable(); + + /// compute the elastic potential energy + void computePotentialEnergy(ElementType el_type, + GhostType ghost_type = _not_ghost) override; protected: + + void computePotentialEnergyOnQuad(const Matrix & grad_u, Real & epot, + Tensor3 & sigma_v, Tensor3 & epsilon_v); + + /// update the dissipated energy, is called after the stress have been /// computed void updateDissipatedEnergy(ElementType el_type, GhostType ghost_type); + + void updateDissipatedEnergyOnQuad(const Matrix & grad_u, + const Matrix & previous_grad_u, + const Matrix & sigma, const Matrix & previous_sigma, + Real & dis_energy, Real & mech_work, const Real & pot_energy); + /// compute stresses on a quadrature point void computeStressOnQuad(const Matrix & grad_u, const Matrix & previous_grad_u, Matrix & sigma, - const Matrix & previous_sigma, - Tensor3 & sigma_v, const Real & sigma_th, - const Real & previous_sigma_th); + Tensor3 & sigma_v, const Real & 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 + /// give the dissipated energy Real getDissipatedEnergy() const; Real getDissipatedEnergy(ElementType type, UInt index) const; + /// get the potential energy + Real getPotentialEnergy() const; + Real getPotentialEnergy(ElementType type, UInt index) const; + + /// get the potential energy + Real getMechanicalWork() const; + Real getMechanicalWork(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: +protected: using voigt_h = VoigtHelper; /// Vectors of viscosity, viscous elastic modulus, one spring element elastic modulus Vector Eta; Vector Ev; Real Einf; /// time step from previous solveStep Real previous_dt; - /// Effective viscoelastic stiffness tensor in voigt notation + /// Stiffness matrix template Matrix C; + /// Compliance matrix template + Matrix D; /// Internal variable: viscous_stress InternalField sigma_v; + /// Internal variable: spring strain in Maxwell element + InternalField epsilon_v; + /// Dissipated energy InternalField dissipated_energy; + + /// Mechanical work + InternalField mechanical_work; + + /// Update internal variable after solve step or not + bool update_variable_flag; }; } // namespace akantu #endif /* __AKANTU_MATERIAL_VISCOELASTIC_MAXWELL_HH__ */