diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_friction.cc b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_friction.cc index e633a403b..f2edfc994 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_friction.cc +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_friction.cc @@ -1,423 +1,428 @@ /** * @file material_cohesive_linear_friction.cc * * @author Mauro Corrado * * @date creation: Tue Jan 12 2016 * @date last modification: Fri Dec 11 2020 * * @brief Linear irreversible cohesive law of mixed mode loading with * random stress definition for extrinsic type * * * @section LICENSE * * Copyright (©) 2015-2021 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_cohesive_linear_friction.hh" #include "solid_mechanics_model_cohesive.hh" namespace akantu { /* -------------------------------------------------------------------------- */ template MaterialCohesiveLinearFriction:: MaterialCohesiveLinearFriction(SolidMechanicsModel & model, const ID & id) : MaterialParent(model, id), residual_sliding("residual_sliding", *this), friction_force("friction_force", *this) { AKANTU_DEBUG_IN(); this->registerParam("mu", mu_max, Real(0.), _pat_parsable | _pat_readable, "Maximum value of the friction coefficient"); this->registerParam("penalty_for_friction", friction_penalty, Real(0.), _pat_parsable | _pat_readable, "Penalty parameter for the friction behavior"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinearFriction::initMaterial() { AKANTU_DEBUG_IN(); MaterialParent::initMaterial(); friction_force.initialize(spatial_dimension); residual_sliding.initialize(1); residual_sliding.initializeHistory(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ // template // void MaterialCohesiveLinearFriction::computeTraction( // __attribute__((unused)) const Array & normal, ElementType el_type, // GhostType ghost_type) { // AKANTU_DEBUG_IN(); // residual_sliding.resize(); // friction_force.resize(); // auto & traction = this->tractions(el_type, ghost_type); // auto & opening = this->opening(el_type, ghost_type); // auto & opening_prev = this->opening.previous(el_type, ghost_type); // auto & contact_traction = this->contact_tractions(el_type, ghost_type); // auto & contact_opening = this->contact_opening(el_type, ghost_type); // auto & sigma_c = this->sigma_c_eff(el_type, ghost_type); // auto & delta_max = this->delta_max(el_type, ghost_type); // auto & delta_max_prev = this->delta_max.previous(el_type, ghost_type); // auto & delta_c = this->delta_c_eff(el_type, ghost_type); // auto & damage = this->damage(el_type, ghost_type); // auto & insertion_stress = this->insertion_stress(el_type, ghost_type); // auto & res_sliding = this->residual_sliding(el_type, ghost_type); // auto & res_sliding_prev = // this->residual_sliding.previous(el_type, ghost_type); // auto & friction_force = this->friction_force(el_type, ghost_type); // auto & penetration = this->penetration(el_type, ghost_type); // Vector normal_opening(spatial_dimension); // Vector tangential_opening(spatial_dimension); // if (not this->model->isDefaultSolverExplicit()) { // this->delta_max(el_type, ghost_type) // .copy(this->delta_max.previous(el_type, ghost_type)); // } // /// loop on each quadrature point // for (auto && data : // zip(make_view(traction, spatial_dimension), // make_view(opening, spatial_dimension), // make_view(opening_prev, spatial_dimension), // make_view(normal, spatial_dimension), sigma_c, delta_max, // delta_max_prev, delta_c, damage, // make_view(insertion_stress, spatial_dimension), res_sliding, // res_sliding_prev, make_view(friction_force, spatial_dimension), // penetration, make_view(contact_opening, spatial_dimension), // make_view(contact_traction, spatial_dimension))) { // auto & _traction = std::get<0>(data); // auto & _opening = std::get<1>(data); // auto & _opening_prev = std::get<2>(data); // auto & _normal = std::get<3>(data); // auto & _sigma_c = std::get<4>(data); // auto & _delta_max = std::get<5>(data); // auto & _delta_max_prev = std::get<6>(data); // auto & _delta_c = std::get<7>(data); // auto & _damage = std::get<8>(data); // auto & _insertion_stress = std::get<9>(data); // auto & _res_sliding = std::get<10>(data); // auto & _res_sliding_prev = std::get<11>(data); // auto & _friction_force = std::get<12>(data); // auto & _penetration = std::get<13>(data); // auto & _contact_opening = std::get<14>(data); // auto & _contact_traction = std::get<15>(data); // Real normal_opening_norm; // Real tangential_opening_norm; // this->computeTractionOnQuad( // _traction, _opening, _normal, _delta_max, _delta_c, // _insertion_stress, _sigma_c, normal_opening, tangential_opening, // normal_opening_norm, tangential_opening_norm, _damage, _penetration, // _contact_traction, _contact_opening); // if (_penetration) { // /// the friction coefficient mu increases with the damage. It // /// equals the maximum value when damage = 1. // // Real damage = std::min(*delta_max_prev_it / *delta_c_it, // // Real(1.)); // Real mu = mu_max; // * damage; // // Real contact_opening_norm = // // std::min(_contact_opening.dot(_normal), Real(0.)); // Real contact_opening_norm = _contact_opening.dot(_normal); // /// the definition of tau_max refers to the opening // /// (penetration) of the previous incremental step // Real tau_max = mu * this->penalty * (std::abs(contact_opening_norm)); // Real trial_elastic_slip = tangential_opening_norm - _res_sliding_prev; // /// tau is the norm of the friction force, acting tangentially to the // /// surface // Real tau = // std::min(std::abs(friction_penalty * trial_elastic_slip), tau_max); // if (trial_elastic_slip < 0.0) { // tau = -tau; // } // /// from tau get the x and y components of friction, to be added in the // /// force vector // if (tangential_opening_norm > Math::getTolerance()) { // Vector tangent_unit_vector(spatial_dimension); // tangent_unit_vector = tangential_opening / tangential_opening_norm; // _friction_force = tau * tangent_unit_vector; // } else { // _friction_force.zero(); // } // /// update residual_sliding // if (friction_penalty == 0.) { // _res_sliding = tangential_opening_norm; // } else { // _res_sliding = // tangential_opening_norm - (std::abs(tau) / friction_penalty); // } // } else { // _friction_force.zero(); // } // _traction += _friction_force; // } // AKANTU_DEBUG_OUT(); // } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinearFriction::computeTraction( __attribute__((unused)) const Array & normal, ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); residual_sliding.resize(); friction_force.resize(); auto & traction = this->tractions(el_type, ghost_type); auto & opening = this->opening(el_type, ghost_type); auto & opening_prev = this->opening.previous(el_type, ghost_type); auto & contact_traction = this->contact_tractions(el_type, ghost_type); auto & contact_opening = this->contact_opening(el_type, ghost_type); auto & sigma_c = this->sigma_c_eff(el_type, ghost_type); auto & delta_max = this->delta_max(el_type, ghost_type); auto & delta_max_prev = this->delta_max.previous(el_type, ghost_type); auto & delta_c = this->delta_c_eff(el_type, ghost_type); auto & damage = this->damage(el_type, ghost_type); auto & insertion_stress = this->insertion_stress(el_type, ghost_type); auto & res_sliding = this->residual_sliding(el_type, ghost_type); auto & res_sliding_prev = this->residual_sliding.previous(el_type, ghost_type); auto & friction_force = this->friction_force(el_type, ghost_type); auto & penetration = this->penetration(el_type, ghost_type); Vector normal_opening(spatial_dimension); Vector tangential_opening(spatial_dimension); if (not this->model->isDefaultSolverExplicit()) { this->delta_max(el_type, ghost_type) .copy(this->delta_max.previous(el_type, ghost_type)); } /// loop on each quadrature point for (auto && data : zip(make_view(traction, spatial_dimension), make_view(opening, spatial_dimension), make_view(opening_prev, spatial_dimension), make_view(normal, spatial_dimension), sigma_c, delta_max, delta_max_prev, delta_c, damage, make_view(insertion_stress, spatial_dimension), res_sliding, res_sliding_prev, make_view(friction_force, spatial_dimension), penetration, make_view(contact_opening, spatial_dimension), make_view(contact_traction, spatial_dimension))) { auto & _traction = std::get<0>(data); auto & _opening = std::get<1>(data); auto & _opening_prev = std::get<2>(data); auto & _normal = std::get<3>(data); auto & _sigma_c = std::get<4>(data); auto & _delta_max = std::get<5>(data); auto & _delta_max_prev = std::get<6>(data); auto & _delta_c = std::get<7>(data); auto & _damage = std::get<8>(data); auto & _insertion_stress = std::get<9>(data); auto & _res_sliding = std::get<10>(data); auto & _res_sliding_prev = std::get<11>(data); auto & _friction_force = std::get<12>(data); auto & _penetration = std::get<13>(data); auto & _contact_opening = std::get<14>(data); auto & _contact_traction = std::get<15>(data); Real normal_opening_norm; Real tangential_opening_norm; this->computeTractionOnQuad( _traction, _opening, _normal, _delta_max, _delta_c, _insertion_stress, _sigma_c, normal_opening, tangential_opening, normal_opening_norm, tangential_opening_norm, _damage, _penetration, _contact_traction, _contact_opening); if (_penetration) { /// the friction coefficient mu increases with the damage. It /// equals the maximum value when damage = 1. // Real damage = std::min(*delta_max_prev_it / *delta_c_it, // Real(1.)); Real mu = mu_max; // * damage; // Real contact_opening_norm = // std::min(_contact_opening.dot(_normal), Real(0.)); Real contact_opening_norm = _contact_opening.dot(_normal); /// the definition of tau_max refers to the opening /// (penetration) of the previous incremental step auto tau_max = mu * this->penalty * (std::abs(contact_opening_norm)); auto trial_elastic_slip = tangential_opening_norm - _res_sliding_prev; auto tau_trial = std::abs(friction_penalty * trial_elastic_slip); auto tau = std::min(tau_trial, tau_max); if (trial_elastic_slip < 0.0) { tau = -tau; } /// from tau get the x and y components of friction, to be added in the /// force vector if (tangential_opening_norm > Math::getTolerance()) { Vector tangent_unit_vector(spatial_dimension); tangent_unit_vector = tangential_opening / tangential_opening_norm; _friction_force = tau * tangent_unit_vector; } else { _friction_force.zero(); } /// update residual_sliding - if (friction_penalty == 0.) { - _res_sliding = tangential_opening_norm; + if (tau_trial > tau_max) { + if (friction_penalty == 0.) { + _res_sliding = tangential_opening_norm; + } else { + auto lambda = (tau_trial - tau_max) / friction_penalty; + _res_sliding = _res_sliding_prev + lambda; + } } else { - auto lambda = (tau_trial - tau_max) / friction_penalty; - _res_sliding = _res_sliding_prev + lambda; + _res_sliding = _res_sliding_prev; } + } else { _friction_force.zero(); } _traction += _friction_force; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinearFriction::computeTangentTraction( ElementType el_type, Array & tangent_matrix, __attribute__((unused)) const Array & normal, GhostType ghost_type) { AKANTU_DEBUG_IN(); auto & opening = this->opening(el_type, ghost_type); auto & opening_prev = this->opening.previous(el_type, ghost_type); auto & delta_max_prev = this->delta_max.previous(el_type, ghost_type); auto & sigma_c = this->sigma_c_eff(el_type, ghost_type); auto & delta_c = this->delta_c_eff(el_type, ghost_type); auto & damage = this->damage(el_type, ghost_type); auto & contact_opening = this->contact_opening(el_type, ghost_type); auto & contact_opening_prev = this->contact_opening.previous(el_type, ghost_type); auto & res_sliding_prev = this->residual_sliding.previous(el_type, ghost_type); auto & penetration = this->penetration(el_type, ghost_type); Vector normal_opening(spatial_dimension); Vector tangential_opening(spatial_dimension); for (auto && data : zip(make_view(tangent_matrix, spatial_dimension, spatial_dimension), make_view(opening, spatial_dimension), make_view(opening_prev, spatial_dimension), make_view(contact_opening, spatial_dimension), make_view(contact_opening_prev, spatial_dimension), make_view(normal, spatial_dimension), sigma_c, delta_max_prev, delta_c, damage, res_sliding_prev, penetration)) { auto & _tangent = std::get<0>(data); auto & _opening = std::get<1>(data); auto & _opening_prev = std::get<2>(data); auto & _contact_opening = std::get<3>(data); auto & _contact_opening_prev = std::get<4>(data); auto & _normal = std::get<5>(data); auto & _sigma_c = std::get<6>(data); auto & _delta_max_prev = std::get<7>(data); auto & _delta_c = std::get<8>(data); auto & _damage = std::get<9>(data); auto & _res_sliding_prev = std::get<10>(data); auto & _penetration = std::get<11>(data); Real normal_opening_norm; Real tangential_opening_norm; this->computeTangentTractionOnQuad( _tangent, _delta_max_prev, _delta_c, _sigma_c, _opening, _normal, normal_opening, tangential_opening, normal_opening_norm, tangential_opening_norm, _damage, _penetration, _contact_opening); if (_penetration) { // Real damage = std::min(*delta_max_it / *delta_c_it, Real(1.)); Real mu = mu_max; // * damage; Real contact_opening_norm = std::min(_contact_opening.dot(_normal), Real(0.)); // Vector normal_opening_prev = (*normal_it); // normal_opening_prev *= normal_opening_prev_norm; Real tau_max = mu * this->penalty * (std::abs(contact_opening_norm)); Real trial_elastic_slip = tangential_opening_norm - _res_sliding_prev; // tau is the norm of the friction force, acting tangentially to the // surface Real tau = std::min(std::abs(friction_penalty * trial_elastic_slip), tau_max); if ((tau < tau_max && tau_max > Math::getTolerance()) or Math::are_float_equal(tau_max, 0.)) { Matrix I(spatial_dimension, spatial_dimension); I.eye(1.); Matrix n_outer_n(spatial_dimension, spatial_dimension); n_outer_n.outerProduct(_normal, _normal); Matrix nn(n_outer_n); I -= nn; _tangent += I * friction_penalty; } } // check if the tangential stiffness matrix is symmetric // for (UInt h = 0; h < spatial_dimension; ++h){ // for (UInt l = h; l < spatial_dimension; ++l){ // if (l > h){ // Real k_ls = (*tangent_it)[spatial_dimension*h+l]; // Real k_us = (*tangent_it)[spatial_dimension*l+h]; // // std::cout << "k_ls = " << k_ls << std::endl; // // std::cout << "k_us = " << k_us << std::endl; // if (std::abs(k_ls) > 1e-13 && std::abs(k_us) > 1e-13){ // Real error = std::abs((k_ls - k_us) / k_us); // if (error > 1e-10){ // std::cout << "non symmetric cohesive matrix" << std::endl; // // std::cout << "error " << error << std::endl; // } // } // } // } // } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ INSTANTIATE_MATERIAL(cohesive_linear_friction, MaterialCohesiveLinearFriction); } // namespace akantu