diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear.cc b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear.cc index 21ce92fcc..806cb7ef9 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear.cc +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear.cc @@ -1,418 +1,419 @@ /** * @file material_cohesive_linear.cc * * @author Mauro Corrado * @author Marco Vocialta * * @date creation: Wed Feb 22 2012 * @date last modification: Thu Jan 14 2016 * * @brief Linear irreversible cohesive law of mixed mode loading with * random stress definition for extrinsic type * * @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 . * */ /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ #include "dof_synchronizer.hh" #include "material_cohesive_linear.hh" #include "solid_mechanics_model_cohesive.hh" #include "sparse_matrix.hh" namespace akantu { /* -------------------------------------------------------------------------- */ template MaterialCohesiveLinear::MaterialCohesiveLinear( SolidMechanicsModel & model, const ID & id) : MaterialCohesive(model, id), sigma_c_eff("sigma_c_eff", *this), delta_c_eff("delta_c_eff", *this), insertion_stress("insertion_stress", *this) { AKANTU_DEBUG_IN(); this->registerParam("beta", beta, Real(0.), _pat_parsable | _pat_readable, "Beta parameter"); this->registerParam("G_c", G_c, Real(0.), _pat_parsable | _pat_readable, "Mode I fracture energy"); this->registerParam("penalty", penalty, Real(0.), _pat_parsable | _pat_readable, "Penalty coefficient"); this->registerParam("volume_s", volume_s, Real(0.), _pat_parsable | _pat_readable, "Reference volume for sigma_c scaling"); this->registerParam("m_s", m_s, Real(1.), _pat_parsable | _pat_readable, "Weibull exponent for sigma_c scaling"); this->registerParam("kappa", kappa, Real(1.), _pat_parsable | _pat_readable, "Kappa parameter"); this->registerParam( "contact_after_breaking", contact_after_breaking, false, _pat_parsable | _pat_readable, "Activation of contact when the elements are fully damaged"); this->registerParam("max_quad_stress_insertion", max_quad_stress_insertion, false, _pat_parsable | _pat_readable, "Insertion of cohesive element when stress is high " "enough just on one quadrature point"); this->registerParam("recompute", recompute, false, _pat_parsable | _pat_modifiable, "recompute solution"); this->use_previous_delta_max = true; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinear::initMaterial() { AKANTU_DEBUG_IN(); MaterialCohesive::initMaterial(); sigma_c_eff.initialize(1); delta_c_eff.initialize(1); insertion_stress.initialize(spatial_dimension); if (not Math::are_float_equal(delta_c, 0.)) delta_c_eff.setDefaultValue(delta_c); else delta_c_eff.setDefaultValue(2 * G_c / sigma_c); if (model->getIsExtrinsic()) scaleInsertionTraction(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinear::updateInternalParameters() { /// compute scalars beta2_kappa2 = beta * beta / kappa / kappa; beta2_kappa = beta * beta / kappa; if (Math::are_float_equal(beta, 0)) beta2_inv = 0; else beta2_inv = 1. / beta / beta; } + /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinear::scaleInsertionTraction() { AKANTU_DEBUG_IN(); // do nothing if volume_s hasn't been specified by the user if (Math::are_float_equal(volume_s, 0.)) return; const Mesh & mesh_facets = model->getMeshFacets(); const auto & fe_engine = model->getFEEngine(); const auto & fe_engine_facet = model->getFEEngine("FacetsFEEngine"); Real base_sigma_c = sigma_c; for (auto && type_facet : mesh_facets.elementTypes(spatial_dimension - 1)) { const Array> & facet_to_element = mesh_facets.getElementToSubelement(type_facet); UInt nb_facet = facet_to_element.size(); UInt nb_quad_per_facet = fe_engine_facet.getNbIntegrationPoints(type_facet); // iterator to modify sigma_c for all the quadrature points of a facet auto sigma_c_iterator = sigma_c(type_facet).begin_reinterpret(nb_quad_per_facet, nb_facet); for (UInt f = 0; f < nb_facet; ++f, ++sigma_c_iterator) { const std::vector & element_list = facet_to_element(f); // compute bounding volume Real volume = 0; auto elem = element_list.begin(); auto elem_end = element_list.end(); for (; elem != elem_end; ++elem) { if (*elem == ElementNull) continue; // unit vector for integration in order to obtain the volume UInt nb_quadrature_points = fe_engine.getNbIntegrationPoints(elem->type); Vector unit_vector(nb_quadrature_points, 1); volume += fe_engine.integrate(unit_vector, elem->type, elem->element, elem->ghost_type); } // scale sigma_c *sigma_c_iterator -= base_sigma_c; *sigma_c_iterator *= std::pow(volume_s / volume, 1. / m_s); *sigma_c_iterator += base_sigma_c; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinear::checkInsertion( bool check_only) { AKANTU_DEBUG_IN(); const Mesh & mesh_facets = model->getMeshFacets(); CohesiveElementInserter & inserter = model->getElementInserter(); for (auto && type_facet : mesh_facets.elementTypes(spatial_dimension - 1)) { ElementType type_cohesive = FEEngine::getCohesiveElementType(type_facet); const auto & facets_check = inserter.getCheckFacets(type_facet); auto & f_insertion = inserter.getInsertionFacets(type_facet); auto & f_filter = facet_filter(type_facet); auto & sig_c_eff = sigma_c_eff(type_cohesive); auto & del_c = delta_c_eff(type_cohesive); auto & ins_stress = insertion_stress(type_cohesive); auto & trac_old = tractions.previous(type_cohesive); const auto & f_stress = model->getStressOnFacets(type_facet); const auto & sigma_lim = sigma_c(type_facet); UInt nb_quad_facet = model->getFEEngine("FacetsFEEngine").getNbIntegrationPoints(type_facet); UInt nb_quad_cohesive = model->getFEEngine("CohesiveFEEngine") .getNbIntegrationPoints(type_cohesive); AKANTU_DEBUG_ASSERT(nb_quad_cohesive == nb_quad_facet, "The cohesive element and the corresponding facet do " "not have the same numbers of integration points"); UInt nb_facet = f_filter.size(); // if (nb_facet == 0) continue; auto sigma_lim_it = sigma_lim.begin(); Matrix stress_tmp(spatial_dimension, spatial_dimension); Matrix normal_traction(spatial_dimension, nb_quad_facet); Vector stress_check(nb_quad_facet); UInt sp2 = spatial_dimension * spatial_dimension; const auto & tangents = model->getTangents(type_facet); const auto & normals = model->getFEEngine("FacetsFEEngine") .getNormalsOnIntegrationPoints(type_facet); auto normal_begin = normals.begin(spatial_dimension); auto tangent_begin = tangents.begin(tangents.getNbComponent()); auto facet_stress_begin = f_stress.begin(spatial_dimension, spatial_dimension * 2); std::vector new_sigmas; std::vector> new_normal_traction; std::vector new_delta_c; // loop over each facet belonging to this material for (UInt f = 0; f < nb_facet; ++f, ++sigma_lim_it) { UInt facet = f_filter(f); // skip facets where check shouldn't be realized if (!facets_check(facet)) continue; // compute the effective norm on each quadrature point of the facet for (UInt q = 0; q < nb_quad_facet; ++q) { UInt current_quad = facet * nb_quad_facet + q; const Vector & normal = normal_begin[current_quad]; const Vector & tangent = tangent_begin[current_quad]; const Matrix & facet_stress_it = facet_stress_begin[current_quad]; // compute average stress on the current quadrature point Matrix stress_1(facet_stress_it.storage(), spatial_dimension, spatial_dimension); Matrix stress_2(facet_stress_it.storage() + sp2, spatial_dimension, spatial_dimension); stress_tmp.copy(stress_1); stress_tmp += stress_2; stress_tmp /= 2.; Vector normal_traction_vec(normal_traction(q)); // compute normal and effective stress stress_check(q) = computeEffectiveNorm(stress_tmp, normal, tangent, normal_traction_vec); } // verify if the effective stress overcomes the threshold Real final_stress = stress_check.mean(); if (max_quad_stress_insertion) final_stress = *std::max_element( stress_check.storage(), stress_check.storage() + nb_quad_facet); if (final_stress > *sigma_lim_it) { f_insertion(facet) = true; if (check_only) continue; // store the new cohesive material parameters for each quadrature // point for (UInt q = 0; q < nb_quad_facet; ++q) { Real new_sigma = stress_check(q); Vector normal_traction_vec(normal_traction(q)); if (spatial_dimension != 3) normal_traction_vec *= -1.; new_sigmas.push_back(new_sigma); new_normal_traction.push_back(normal_traction_vec); Real new_delta; // set delta_c in function of G_c or a given delta_c value if (Math::are_float_equal(delta_c, 0.)) new_delta = 2 * G_c / new_sigma; else new_delta = (*sigma_lim_it) / new_sigma * delta_c; new_delta_c.push_back(new_delta); } } } // update material data for the new elements UInt old_nb_quad_points = sig_c_eff.size(); UInt new_nb_quad_points = new_sigmas.size(); sig_c_eff.resize(old_nb_quad_points + new_nb_quad_points); ins_stress.resize(old_nb_quad_points + new_nb_quad_points); trac_old.resize(old_nb_quad_points + new_nb_quad_points); del_c.resize(old_nb_quad_points + new_nb_quad_points); for (UInt q = 0; q < new_nb_quad_points; ++q) { sig_c_eff(old_nb_quad_points + q) = new_sigmas[q]; del_c(old_nb_quad_points + q) = new_delta_c[q]; for (UInt dim = 0; dim < spatial_dimension; ++dim) { ins_stress(old_nb_quad_points + q, dim) = new_normal_traction[q](dim); trac_old(old_nb_quad_points + q, dim) = new_normal_traction[q](dim); } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinear::computeTraction( const Array & normal, ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); /// define iterators auto traction_it = tractions(el_type, ghost_type).begin(spatial_dimension); auto opening_it = opening(el_type, ghost_type).begin(spatial_dimension); auto contact_traction_it = contact_tractions(el_type, ghost_type).begin(spatial_dimension); auto contact_opening_it = contact_opening(el_type, ghost_type).begin(spatial_dimension); auto normal_it = normal.begin(spatial_dimension); auto traction_end = tractions(el_type, ghost_type).end(spatial_dimension); auto sigma_c_it = sigma_c_eff(el_type, ghost_type).begin(); auto delta_max_it = delta_max(el_type, ghost_type).begin(); auto delta_c_it = delta_c_eff(el_type, ghost_type).begin(); auto damage_it = damage(el_type, ghost_type).begin(); auto insertion_stress_it = insertion_stress(el_type, ghost_type).begin(spatial_dimension); Vector normal_opening(spatial_dimension); Vector tangential_opening(spatial_dimension); /// loop on each quadrature point for (; traction_it != traction_end; ++traction_it, ++opening_it, ++normal_it, ++sigma_c_it, ++delta_max_it, ++delta_c_it, ++damage_it, ++contact_traction_it, ++insertion_stress_it, ++contact_opening_it) { Real normal_opening_norm, tangential_opening_norm; bool penetration; this->computeTractionOnQuad( *traction_it, *opening_it, *normal_it, *delta_max_it, *delta_c_it, *insertion_stress_it, *sigma_c_it, normal_opening, tangential_opening, normal_opening_norm, tangential_opening_norm, *damage_it, penetration, *contact_traction_it, *contact_opening_it); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinear::computeTangentTraction( const ElementType & el_type, Array & tangent_matrix, const Array & normal, GhostType ghost_type) { AKANTU_DEBUG_IN(); /// define iterators auto tangent_it = tangent_matrix.begin(spatial_dimension, spatial_dimension); auto tangent_end = tangent_matrix.end(spatial_dimension, spatial_dimension); auto normal_it = normal.begin(spatial_dimension); auto opening_it = opening(el_type, ghost_type).begin(spatial_dimension); /// NB: delta_max_it points on delta_max_previous, i.e. the /// delta_max related to the solution of the previous incremental /// step auto delta_max_it = delta_max.previous(el_type, ghost_type).begin(); auto sigma_c_it = sigma_c_eff(el_type, ghost_type).begin(); auto delta_c_it = delta_c_eff(el_type, ghost_type).begin(); auto damage_it = damage(el_type, ghost_type).begin(); auto contact_opening_it = contact_opening(el_type, ghost_type).begin(spatial_dimension); Vector normal_opening(spatial_dimension); Vector tangential_opening(spatial_dimension); for (; tangent_it != tangent_end; ++tangent_it, ++normal_it, ++opening_it, ++delta_max_it, ++sigma_c_it, ++delta_c_it, ++damage_it, ++contact_opening_it) { Real normal_opening_norm, tangential_opening_norm; bool penetration; this->computeTangentTractionOnQuad( *tangent_it, *delta_max_it, *delta_c_it, *sigma_c_it, *opening_it, *normal_it, normal_opening, tangential_opening, normal_opening_norm, tangential_opening_norm, *damage_it, penetration, *contact_opening_it); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ INSTANTIATE_MATERIAL(cohesive_linear, MaterialCohesiveLinear); } // akantu diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear.hh b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear.hh index ceb677555..56dd8570e 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear.hh +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear.hh @@ -1,188 +1,188 @@ /** * @file material_cohesive_linear.hh * * @author Mauro Corrado * @author Marco Vocialta * * @date creation: Fri Jun 18 2010 * @date last modification: Thu Jan 14 2016 * * @brief Linear irreversible cohesive law of mixed mode loading with * random stress definition for extrinsic type * * @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 . * */ /* -------------------------------------------------------------------------- */ #include "material_cohesive.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_COHESIVE_LINEAR_HH__ #define __AKANTU_MATERIAL_COHESIVE_LINEAR_HH__ namespace akantu { /** * Cohesive material linear damage for extrinsic case * * parameters in the material files : * - sigma_c : critical stress sigma_c (default: 0) * - beta : weighting parameter for sliding and normal opening (default: * 0) * - G_cI : fracture energy for mode I (default: 0) * - G_cII : fracture energy for mode II (default: 0) * - penalty : stiffness in compression to prevent penetration */ template class MaterialCohesiveLinear : public MaterialCohesive { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: MaterialCohesiveLinear(SolidMechanicsModel & model, const ID & id = ""); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// initialize the material parameters void initMaterial() override; void updateInternalParameters() override; /// check stress for cohesive elements' insertion void checkInsertion(bool check_only = false) override; /// compute effective stress norm for insertion check Real computeEffectiveNorm(const Matrix & stress, const Vector & normal, const Vector & tangent, Vector & normal_stress) const; protected: /// constitutive law void computeTraction(const Array & normal, ElementType el_type, GhostType ghost_type = _not_ghost) override; /// compute tangent stiffness matrix void computeTangentTraction(const ElementType & el_type, Array & tangent_matrix, const Array & normal, GhostType ghost_type) override; /** * Scale insertion traction sigma_c according to the volume of the * two elements surrounding a facet * * see the article: F. Zhou and J. F. Molinari "Dynamic crack * propagation with cohesive elements: a methodology to address mesh * dependency" International Journal for Numerical Methods in * Engineering (2004) */ void scaleInsertionTraction(); /// compute the traction for a given quadrature point inline void computeTractionOnQuad( - Vector & traction, const Vector & opening, + Vector & traction, Vector & opening, const Vector & normal, Real & delta_max, const Real & delta_c, const Vector & insertion_stress, const Real & sigma_c, Vector & normal_opening, Vector & tangential_opening, Real & normal_opening_norm, Real & tangential_opening_norm, Real & damage, bool & penetration, Vector & contact_traction, Vector & contact_opening); inline void computeTangentTractionOnQuad( Matrix & tangent, Real & delta_max, const Real & delta_c, const Real & sigma_c, Vector & opening, const Vector & normal, Vector & normal_opening, Vector & tangential_opening, Real & normal_opening_norm, Real & tangential_opening_norm, Real & damage, bool & penetration, Vector & contact_opening); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get sigma_c_eff AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(InsertionTraction, sigma_c_eff, Real); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// beta parameter Real beta; /// beta square inverse to compute effective norm Real beta2_inv; /// mode I fracture energy Real G_c; /// kappa parameter Real kappa; /// constitutive law scalar to compute delta Real beta2_kappa2; /// constitutive law scalar to compute traction Real beta2_kappa; /// penalty coefficient Real penalty; /// reference volume used to scale sigma_c Real volume_s; /// weibull exponent used to scale sigma_c Real m_s; /// variable defining if we are recomputing the last loading step /// after load_reduction bool recompute; /// critical effective stress RandomInternalField sigma_c_eff; /// effective critical displacement (each element can have a /// different value) CohesiveInternalField delta_c_eff; /// stress at insertion CohesiveInternalField insertion_stress; /// variable saying if there should be penalty contact also after /// breaking the cohesive elements bool contact_after_breaking; /// insertion of cohesive element when stress is high enough just on /// one quadrature point bool max_quad_stress_insertion; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ } // akantu #include "material_cohesive_linear_inline_impl.cc" #endif /* __AKANTU_MATERIAL_COHESIVE_LINEAR_HH__ */ diff --git a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_inline_impl.cc b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_inline_impl.cc index a6b111fb4..49e9c5a22 100644 --- a/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_inline_impl.cc +++ b/src/model/solid_mechanics/solid_mechanics_model_cohesive/materials/constitutive_laws/material_cohesive_linear_inline_impl.cc @@ -1,276 +1,265 @@ /** * @file material_cohesive_linear_inline_impl.cc * * @author Mauro Corrado * @author Marco Vocialta * * @date creation: Wed Apr 22 2015 * @date last modification: Thu Jan 14 2016 * * @brief Inline functions of the MaterialCohesiveLinear * * @section LICENSE * * Copyright (©) 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 . * */ /* -------------------------------------------------------------------------- */ #include "material_cohesive_linear.hh" #include "solid_mechanics_model_cohesive.hh" /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_COHESIVE_LINEAR_INLINE_IMPL_CC__ #define __AKANTU_MATERIAL_COHESIVE_LINEAR_INLINE_IMPL_CC__ /* -------------------------------------------------------------------------- */ namespace akantu { /* -------------------------------------------------------------------------- */ template inline Real MaterialCohesiveLinear::computeEffectiveNorm( const Matrix & stress, const Vector & normal, const Vector & tangent, Vector & normal_traction) const { normal_traction.mul(stress, normal); Real normal_contrib = normal_traction.dot(normal); /// in 3D tangential components must be summed Real tangent_contrib = 0; if (dim == 2) { Real tangent_contrib_tmp = normal_traction.dot(tangent); tangent_contrib += tangent_contrib_tmp * tangent_contrib_tmp; } else if (dim == 3) { for (UInt s = 0; s < dim - 1; ++s) { const Vector tangent_v(tangent.storage() + s * dim, dim); Real tangent_contrib_tmp = normal_traction.dot(tangent_v); tangent_contrib += tangent_contrib_tmp * tangent_contrib_tmp; } } tangent_contrib = std::sqrt(tangent_contrib); normal_contrib = std::max(Real(0.), normal_contrib); return std::sqrt(normal_contrib * normal_contrib + tangent_contrib * tangent_contrib * beta2_inv); } /* -------------------------------------------------------------------------- */ template inline void MaterialCohesiveLinear::computeTractionOnQuad( - Vector & traction, const Vector & opening, + Vector & traction, Vector & opening, const Vector & normal, Real & delta_max, const Real & delta_c, const Vector & insertion_stress, const Real & sigma_c, Vector & normal_opening, Vector & tangential_opening, Real & normal_opening_norm, Real & tangential_opening_norm, Real & damage, bool & penetration, Vector & contact_traction, Vector & contact_opening) { /// compute normal and tangential opening vectors normal_opening_norm = opening.dot(normal); normal_opening = normal; normal_opening *= normal_opening_norm; tangential_opening = opening; tangential_opening -= normal_opening; tangential_opening_norm = tangential_opening.norm(); /** * compute effective opening displacement * @f$ \delta = \sqrt{ * \frac{\beta^2}{\kappa^2} \Delta_t^2 + \Delta_n^2 } @f$ */ Real delta = tangential_opening_norm * tangential_opening_norm * this->beta2_kappa2; - penetration = normal_opening_norm < 0.0; + penetration = normal_opening_norm / delta_c < -Math::getTolerance(); + //penetration = normal_opening_norm < 0.; if (this->contact_after_breaking == false && Math::are_float_equal(damage, 1.)) penetration = false; - /** - * if during the convergence loop a cohesive element continues to - * jumps from penetration to opening, and convergence is not - * reached, its penalty parameter will be reduced in the - * recomputation of the same incremental step. recompute is set - * equal to true when convergence is not reached in the - * solveStepCohesive function and the execution of the program - * goes back to the main file where the variable load_reduction - * is set equal to true. - */ - if (penetration) { /// use penalty coefficient in case of penetration contact_traction = normal_opening; contact_traction *= this->penalty; contact_opening = normal_opening; /// don't consider penetration contribution for delta - // opening = tangential_opening; + opening = tangential_opening; normal_opening.clear(); - } else { delta += normal_opening_norm * normal_opening_norm; contact_traction.clear(); contact_opening.clear(); } delta = std::sqrt(delta); /// update maximum displacement and damage delta_max = std::max(delta_max, delta); damage = std::min(delta_max / delta_c, Real(1.)); /** * Compute traction @f$ \mathbf{T} = \left( * \frac{\beta^2}{\kappa} \Delta_t \mathbf{t} + \Delta_n * \mathbf{n} \right) \frac{\sigma_c}{\delta} \left( 1- * \frac{\delta}{\delta_c} \right)@f$ */ if (Math::are_float_equal(damage, 1.)) traction.clear(); else if (Math::are_float_equal(damage, 0.)) { if (penetration) traction.clear(); else traction = insertion_stress; } else { traction = tangential_opening; traction *= this->beta2_kappa; traction += normal_opening; AKANTU_DEBUG_ASSERT(delta_max != 0., "Division by zero, tolerance might be too low"); traction *= sigma_c / delta_max * (1. - damage); } } /* -------------------------------------------------------------------------- */ template inline void MaterialCohesiveLinear::computeTangentTractionOnQuad( Matrix & tangent, Real & delta_max, const Real & delta_c, const Real & sigma_c, Vector & opening, const Vector & normal, Vector & normal_opening, Vector & tangential_opening, Real & normal_opening_norm, Real & tangential_opening_norm, Real & damage, bool & penetration, Vector & contact_opening) { /** * During the update of the residual the interpenetrations are * stored in the array "contact_opening", therefore, in the case * of penetration, in the array "opening" there are only the * tangential components. */ opening += contact_opening; /// compute normal and tangential opening vectors normal_opening_norm = opening.dot(normal); normal_opening = normal; normal_opening *= normal_opening_norm; tangential_opening = opening; tangential_opening -= normal_opening; tangential_opening_norm = tangential_opening.norm(); Real delta = tangential_opening_norm * tangential_opening_norm * this->beta2_kappa2; penetration = normal_opening_norm < 0.0; if (this->contact_after_breaking == false && Math::are_float_equal(damage, 1.)) penetration = false; Real derivative = 0; // derivative = d(t/delta)/ddelta Real t = 0; Matrix n_outer_n(spatial_dimension, spatial_dimension); n_outer_n.outerProduct(normal, normal); if (penetration) { /// stiffness in compression given by the penalty parameter tangent += n_outer_n; tangent *= penalty; opening = tangential_opening; normal_opening_norm = opening.dot(normal); normal_opening = normal; normal_opening *= normal_opening_norm; } else { delta += normal_opening_norm * normal_opening_norm; } delta = std::sqrt(delta); /** * Delta has to be different from 0 to have finite values of * tangential stiffness. At the element insertion, delta = * 0. Therefore, a fictictious value is defined, for the * evaluation of the first value of K. */ if (delta < Math::getTolerance()) delta = delta_c / 1000.; if (delta >= delta_max) { if (delta <= delta_c) { derivative = -sigma_c / (delta * delta); t = sigma_c * (1 - delta / delta_c); } else { derivative = 0.; t = 0.; } } else if (delta < delta_max) { Real tmax = sigma_c * (1 - delta_max / delta_c); t = tmax / delta_max * delta; } /// computation of the derivative of the constitutive law (dT/ddelta) Matrix I(spatial_dimension, spatial_dimension); I.eye(this->beta2_kappa); Matrix nn(n_outer_n); nn *= (1. - this->beta2_kappa); nn += I; nn *= t / delta; Vector t_tilde(normal_opening); t_tilde *= (1. - this->beta2_kappa2); Vector mm(opening); mm *= this->beta2_kappa2; t_tilde += mm; Vector t_hat(normal_opening); t_hat += this->beta2_kappa * tangential_opening; Matrix prov(spatial_dimension, spatial_dimension); prov.outerProduct(t_hat, t_tilde); prov *= derivative / delta; prov += nn; Matrix prov_t = prov.transpose(); tangent += prov_t; } /* -------------------------------------------------------------------------- */ } // akantu /* -------------------------------------------------------------------------- */ #endif //__AKANTU_MATERIAL_COHESIVE_LINEAR_INLINE_IMPL_CC__