diff --git a/packages/cohesive_element.cmake b/packages/cohesive_element.cmake index 334840d1d..e0323e6ee 100644 --- a/packages/cohesive_element.cmake +++ b/packages/cohesive_element.cmake @@ -1,125 +1,127 @@ #=============================================================================== # @file cohesive_element.cmake # # @author Mauro Corrado # @author Nicolas Richart # @author Marco Vocialta # # @date creation: Tue Oct 16 2012 # @date last modification: Tue Jan 12 2016 # # @brief package description for cohesive elements # # @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(cohesive_element DESCRIPTION "Use cohesive_element package of Akantu" DEPENDS lapack) package_declare_sources(cohesive_element fe_engine/cohesive_element.cc fe_engine/cohesive_element.hh fe_engine/fe_engine_template_cohesive.cc fe_engine/shape_cohesive.hh fe_engine/shape_cohesive_inline_impl.cc mesh_utils/cohesive_element_inserter.cc mesh_utils/cohesive_element_inserter.hh model/solid_mechanics/materials/material_cohesive/cohesive_internal_field.hh model/solid_mechanics/materials/material_cohesive/cohesive_internal_field_tmpl.hh model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_bilinear.cc model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_bilinear.hh model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_exponential.cc model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_exponential.hh model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear.cc model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear.hh model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_fatigue.cc model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_fatigue.hh model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_friction.cc model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_friction.hh + model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_uncoupled.cc + model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_uncoupled.hh model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_inline_impl.cc model/solid_mechanics/materials/material_cohesive/material_cohesive.cc model/solid_mechanics/materials/material_cohesive/material_cohesive.hh model/solid_mechanics/materials/material_cohesive/material_cohesive_inline_impl.cc model/solid_mechanics/materials/material_cohesive_includes.hh model/solid_mechanics/solid_mechanics_model_cohesive/fragment_manager.cc model/solid_mechanics/solid_mechanics_model_cohesive/fragment_manager.hh model/solid_mechanics/solid_mechanics_model_cohesive/material_selector_cohesive.cc model/solid_mechanics/solid_mechanics_model_cohesive/material_selector_cohesive.hh model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive.cc model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive.hh model/solid_mechanics/solid_mechanics_model_cohesive/solid_mechanics_model_cohesive_inline_impl.cc ) package_declare_elements(cohesive_element ELEMENT_TYPES _cohesive_2d_4 _cohesive_2d_6 _cohesive_1d_2 _cohesive_3d_6 _cohesive_3d_12 _cohesive_3d_8 _cohesive_3d_16 KIND cohesive GEOMETRICAL_TYPES _gt_cohesive_2d_4 _gt_cohesive_2d_6 _gt_cohesive_1d_2 _gt_cohesive_3d_6 _gt_cohesive_3d_12 _gt_cohesive_3d_8 _gt_cohesive_3d_16 FE_ENGINE_LISTS gradient_on_integration_points interpolate_on_integration_points compute_normals_on_integration_points inverse_map contains get_shapes_derivatives ) package_declare_material_infos(cohesive_element LIST AKANTU_COHESIVE_MATERIAL_LIST INCLUDE material_cohesive_includes.hh ) package_declare_documentation_files(cohesive_element manual-cohesive_elements.tex manual-cohesive_elements_insertion.tex manual-cohesive_laws.tex manual-appendix-materials-cohesive.tex figures/cohesive2d.pdf figures/cohesive_exponential.pdf figures/linear_cohesive_law.pdf figures/bilinear_cohesive_law.pdf ) package_declare_documentation(cohesive_element "This package activates the cohesive elements engine within Akantu." "It depends on:" "\\begin{itemize}" " \\item A fortran compiler." " \\item An implementation of BLAS/LAPACK." "\\end{itemize}" ) diff --git a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_friction.cc b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_friction.cc index aae0fb37e..94e88ed3a 100644 --- a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_friction.cc +++ b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_friction.cc @@ -1,402 +1,401 @@ /** * @file material_cohesive_linear_friction.cc * * @author Mauro Corrado * * @date creation: Tue Jan 12 2016 * @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 (©) 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_friction.hh" #include "solid_mechanics_model_cohesive.hh" __BEGIN_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(); /// define iterators Array::vector_iterator traction_it = this->tractions(el_type, ghost_type).begin(spatial_dimension); Array::vector_iterator traction_end = this->tractions(el_type, ghost_type).end(spatial_dimension); Array::vector_iterator opening_it = this->opening(el_type, ghost_type).begin(spatial_dimension); /// opening_prec is the opening at the end of the previous step in /// the Newton-Raphson loop Array::vector_iterator opening_prec_it = this->opening_prec(el_type, ghost_type).begin(spatial_dimension); /// opening_prec_prev is the opening (opening + penetration) /// referred to the convergence of the previous incremental step Array::vector_iterator opening_prec_prev_it = this->opening_prec.previous(el_type, ghost_type).begin(spatial_dimension); Array::vector_iterator contact_traction_it = this->contact_tractions(el_type, ghost_type).begin(spatial_dimension); Array::vector_iterator contact_opening_it = this->contact_opening(el_type, ghost_type).begin(spatial_dimension); Array::const_iterator< Vector > normal_it = this->normal.begin(spatial_dimension); Array::iteratorsigma_c_it = this->sigma_c_eff(el_type, ghost_type).begin(); Array::iteratordelta_max_it = this->delta_max(el_type, ghost_type).begin(); Array::iteratordelta_max_prev_it = this->delta_max.previous(el_type, ghost_type).begin(); Array::iteratordelta_c_it = this->delta_c_eff(el_type, ghost_type).begin(); Array::iteratordamage_it = this->damage(el_type, ghost_type).begin(); Array::vector_iterator insertion_stress_it = this->insertion_stress(el_type, ghost_type).begin(spatial_dimension); Array::iteratorres_sliding_it = this->residual_sliding(el_type, ghost_type).begin(); Array::iteratorres_sliding_prev_it = this->residual_sliding.previous(el_type, ghost_type).begin(); Array::vector_iterator friction_force_it = this->friction_force(el_type, ghost_type).begin(spatial_dimension); Array::iterator reduction_penalty_it = this->reduction_penalty(el_type, ghost_type).begin(); Vector normal_opening(spatial_dimension); Vector tangential_opening(spatial_dimension); if (! this->model->isExplicit()) this->delta_max(el_type, ghost_type).copy(this->delta_max.previous(el_type, ghost_type)); - // std::cout<< "NEW LOOOOOOOOOOOOOOOOOOOOOP " << std::endl; - /// loop on each quadrature point for (; traction_it != traction_end; ++traction_it, ++opening_it, ++opening_prec_prev_it, ++opening_prec_it, ++normal_it, ++sigma_c_it, ++delta_max_it, ++delta_c_it, ++damage_it, ++contact_traction_it, ++insertion_stress_it, ++contact_opening_it, ++delta_max_prev_it, ++res_sliding_it, ++res_sliding_prev_it, ++friction_force_it, ++reduction_penalty_it) { Real normal_opening_norm, tangential_opening_norm; bool penetration; Real current_penalty = 0.; this->computeTractionOnQuad(*traction_it, *delta_max_it, *delta_max_prev_it, *delta_c_it, *insertion_stress_it, *sigma_c_it, *opening_it, *opening_prec_it, *normal_it, normal_opening, tangential_opening, normal_opening_norm, tangential_opening_norm, *damage_it, penetration, *reduction_penalty_it, current_penalty, *contact_traction_it, *contact_opening_it); 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; /// the definition of tau_max refers to the opening /// (penetration) of the previous incremental step Real normal_opening_prev_norm = std::min(opening_prec_prev_it->dot(*normal_it), Real(0.)); // Vector normal_opening_prev = (*normal_it); // normal_opening_prev *= normal_opening_prev_norm; Real tau_max = mu * current_penalty * (std::abs(normal_opening_prev_norm)); Real delta_sliding_norm = std::abs(tangential_opening_norm - *res_sliding_prev_it); /// tau is the norm of the friction force, acting tangentially to the surface Real tau = std::min(friction_penalty * delta_sliding_norm, tau_max); if ((tangential_opening_norm - *res_sliding_prev_it) < 0.0) tau = -tau; /// from tau get the x and y components of friction, to be added in the force vector Vector tangent_unit_vector(spatial_dimension); tangent_unit_vector = tangential_opening / tangential_opening_norm; *friction_force_it = tau * tangent_unit_vector; /// update residual_sliding *res_sliding_it = tangential_opening_norm - (std::abs(tau) / friction_penalty); } else { friction_force_it->clear(); } *traction_it += *friction_force_it; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinearFriction::checkDeltaMax(GhostType ghost_type) { AKANTU_DEBUG_IN(); MaterialParent::checkDeltaMax(); Mesh & mesh = this->fem_cohesive->getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type, _ek_cohesive); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type, _ek_cohesive); for(; it != last_type; ++it) { Array & elem_filter = this->element_filter(*it, ghost_type); UInt nb_element = elem_filter.getSize(); if (nb_element == 0) continue; ElementType el_type = *it; /// define iterators Array::iteratordelta_max_it = this->delta_max(el_type, ghost_type).begin(); Array::iteratordelta_max_end = this->delta_max(el_type, ghost_type).end(); Array::iteratorres_sliding_it = this->residual_sliding(el_type, ghost_type).begin(); Array::iteratorres_sliding_prev_it = this->residual_sliding.previous(el_type, ghost_type).begin(); /// loop on each quadrature point for (; delta_max_it != delta_max_end; ++delta_max_it, ++res_sliding_it, ++res_sliding_prev_it) { /** * in case convergence is not reached, set "residual_sliding" * for the friction behaviour to the value referred to the * previous incremental step */ *res_sliding_it = *res_sliding_prev_it; } } } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinearFriction::computeTangentTraction(const ElementType & el_type, Array & tangent_matrix, __attribute__((unused)) const Array & normal, GhostType ghost_type) { AKANTU_DEBUG_IN(); /// define iterators Array::matrix_iterator tangent_it = tangent_matrix.begin(spatial_dimension, spatial_dimension); Array::matrix_iterator tangent_end = tangent_matrix.end(spatial_dimension, spatial_dimension); Array::const_vector_iterator normal_it = this->normal.begin(spatial_dimension); Array::vector_iterator opening_it = this->opening(el_type, ghost_type).begin(spatial_dimension); Array::vector_iterator opening_prec_prev_it = this->opening_prec.previous(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 */ Array::iteratordelta_max_it = this->delta_max.previous(el_type, ghost_type).begin(); Array::iteratorsigma_c_it = this->sigma_c_eff(el_type, ghost_type).begin(); Array::iteratordelta_c_it = this->delta_c_eff(el_type, ghost_type).begin(); Array::iteratordamage_it = this->damage(el_type, ghost_type).begin(); Array::iterator< Vector > contact_opening_it = this->contact_opening(el_type, ghost_type).begin(spatial_dimension); Array::iteratorres_sliding_prev_it = this->residual_sliding.previous(el_type, ghost_type).begin(); Array::iterator reduction_penalty_it = this->reduction_penalty(el_type, ghost_type).begin(); Vector normal_opening(spatial_dimension); Vector tangential_opening(spatial_dimension); for (; tangent_it != tangent_end; ++tangent_it, ++normal_it, ++opening_it, ++opening_prec_prev_it, ++delta_max_it, ++sigma_c_it, ++delta_c_it, ++damage_it, ++contact_opening_it, ++res_sliding_prev_it, ++reduction_penalty_it) { Real normal_opening_norm, tangential_opening_norm; bool penetration; Real current_penalty = 0.; 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, *reduction_penalty_it, current_penalty, *contact_opening_it); if (penetration) { // Real damage = std::min(*delta_max_it / *delta_c_it, Real(1.)); Real mu = mu_max; // * damage; Real normal_opening_prev_norm = std::min(opening_prec_prev_it->dot(*normal_it), Real(0.)); // Vector normal_opening_prev = (*normal_it); // normal_opening_prev *= normal_opening_prev_norm; Real tau_max = mu * current_penalty * (std::abs(normal_opening_prev_norm)); Real delta_sliding_norm = std::abs(tangential_opening_norm - *res_sliding_prev_it); // tau is the norm of the friction force, acting tangentially to the surface Real tau = std::min(friction_penalty * delta_sliding_norm, tau_max); if (tau < tau_max && tau_max > Math::getTolerance()){ Matrix I(spatial_dimension, spatial_dimension); I.eye(1.); Matrix n_outer_n(spatial_dimension, spatial_dimension); n_outer_n.outerProduct(*normal_it, *normal_it); Matrix nn(n_outer_n); I -= nn; *tangent_it += 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(MaterialCohesiveLinearFriction); __END_AKANTU__ diff --git a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_uncoupled.cc b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_uncoupled.cc new file mode 100644 index 000000000..16ff1414b --- /dev/null +++ b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_uncoupled.cc @@ -0,0 +1,588 @@ +/** + * @file material_cohesive_linear_uncoupled.cc + * + * @author Mauro Corrado + * + * @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 "material_cohesive_linear_uncoupled.hh" +#include "solid_mechanics_model_cohesive.hh" + +__BEGIN_AKANTU__ + +/* -------------------------------------------------------------------------- */ +template +MaterialCohesiveLinearUncoupled +::MaterialCohesiveLinearUncoupled(SolidMechanicsModel & model, + const ID & id) : + MaterialCohesiveLinear(model,id), + delta_nt_max("delta_nt_max", *this), + delta_st_max("delta_st_max", *this), + damage_nt("damage_nt", *this), + damage_st("damage_st", *this){ + AKANTU_DEBUG_IN(); + + this->registerParam("roughness", R, Real(1.), + _pat_parsable | _pat_readable, + "Roughness to define coupling between mode II and mode I"); + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +template +void MaterialCohesiveLinearUncoupled::initMaterial() { + AKANTU_DEBUG_IN(); + + MaterialCohesiveLinear::initMaterial(); + + delta_nt_max.initialize(1); + delta_st_max.initialize(1); + damage_nt.initialize(1); + damage_st.initialize(1); + + delta_nt_max.initializeHistory(); + delta_st_max.initializeHistory(); + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +template +void MaterialCohesiveLinearUncoupled::computeTraction(const Array & normal, + ElementType el_type, + GhostType ghost_type) { + AKANTU_DEBUG_IN(); + + delta_nt_max.resize(); + delta_st_max.resize(); + damage_nt.resize(); + damage_st.resize(); + + /// define iterators + Array::vector_iterator traction_it = + this->tractions(el_type, ghost_type).begin(spatial_dimension); + + Array::vector_iterator opening_it = + this->opening(el_type, ghost_type).begin(spatial_dimension); + + /// opening_prec is the opening of the previous step in the + /// Newton-Raphson loop + Array::vector_iterator opening_prec_it = + this->opening_prec(el_type, ghost_type).begin(spatial_dimension); + + Array::vector_iterator contact_traction_it = + this->contact_tractions(el_type, ghost_type).begin(spatial_dimension); + + Array::vector_iterator contact_opening_it = + this->contact_opening(el_type, ghost_type).begin(spatial_dimension); + + Array::const_vector_iterator normal_it = + this->normal.begin(spatial_dimension); + + Array::vector_iterator traction_end = + this->tractions(el_type, ghost_type).end(spatial_dimension); + + Array::scalar_iterator sigma_c_it = + this->sigma_c_eff(el_type, ghost_type).begin(); + + Array::scalar_iterator delta_nt_max_it = + this->delta_nt_max(el_type, ghost_type).begin(); + + Array::scalar_iterator delta_st_max_it = + this->delta_st_max(el_type, ghost_type).begin(); + + Array::scalar_iterator delta_nt_max_prev_it = + this->delta_nt_max.previous(el_type, ghost_type).begin(); + + Array::scalar_iterator delta_st_max_prev_it = + this->delta_st_max.previous(el_type, ghost_type).begin(); + + Array::scalar_iterator delta_c_it = + this->delta_c_eff(el_type, ghost_type).begin(); + + Array::scalar_iterator damage_nt_it = + this->damage_nt(el_type, ghost_type).begin(); + + Array::scalar_iterator damage_st_it = + this->damage_st(el_type, ghost_type).begin(); + + Array::vector_iterator insertion_stress_it = + this->insertion_stress(el_type, ghost_type).begin(spatial_dimension); + + Array::scalar_iterator reduction_penalty_it = + this->reduction_penalty(el_type, ghost_type).begin(); + + Vector normal_opening(spatial_dimension); + Vector tangential_opening(spatial_dimension); + + if (! this->model->isExplicit()){ + this->delta_nt_max(el_type, ghost_type).copy(this->delta_nt_max.previous(el_type, ghost_type)); + this->delta_st_max(el_type, ghost_type).copy(this->delta_st_max.previous(el_type, ghost_type)); + } + + /// loop on each quadrature point + for (; traction_it != traction_end; + ++traction_it, ++opening_it, ++opening_prec_it, + ++normal_it, ++sigma_c_it, ++delta_c_it, ++delta_nt_max_it, + ++delta_st_max_it, ++damage_nt_it, ++damage_st_it, + ++contact_traction_it, ++contact_opening_it, ++insertion_stress_it, + ++reduction_penalty_it) { + + Real normal_opening_norm, tangential_opening_norm; + bool penetration; + Real current_penalty = 0.; + Real delta_c_s = this->kappa / this->beta * (*delta_c_it); + Real delta_c_s2_R2 = delta_c_s * delta_c_s / (R * R); + + /// compute normal and tangential opening vectors + normal_opening_norm = opening_it->dot(*normal_it); + Vector normal_opening = *normal_it; + normal_opening *= normal_opening_norm; + + Vector tangential_opening = *opening_it; + tangential_opening -= normal_opening; + tangential_opening_norm = tangential_opening.norm(); + + /// compute effective opening displacement + Real delta_nt = tangential_opening_norm * tangential_opening_norm * this->beta2_kappa2; + Real delta_st = tangential_opening_norm * tangential_opening_norm; + + penetration = normal_opening_norm < 0.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. + */ + Real normal_opening_prec_norm = opening_prec_it->dot(*normal_it); + *opening_prec_it = *opening_it; + + if (!this->model->isExplicit() && !this->recompute) + if ((normal_opening_prec_norm * normal_opening_norm) < 0.0) { + *reduction_penalty_it = true; + } + + if (penetration) { + if (this->recompute && *reduction_penalty_it){ + /// the penalty parameter is locally reduced + current_penalty = this->penalty / 100.; + } + else + current_penalty = this->penalty; + + /// use penalty coefficient in case of penetration + *contact_traction_it = normal_opening; + *contact_traction_it *= current_penalty; + *contact_opening_it = normal_opening; + + /// don't consider penetration contribution for delta + *opening_it = tangential_opening; + normal_opening.clear(); + + } + else { + delta_nt += normal_opening_norm * normal_opening_norm; + delta_st += normal_opening_norm * normal_opening_norm * delta_c_s2_R2; + contact_traction_it->clear(); + contact_opening_it->clear(); + } + + delta_nt = std::sqrt(delta_nt); + delta_st = std::sqrt(delta_st); + + /// update maximum displacement and damage + *delta_nt_max_it = std::max(*delta_nt_max_it, delta_nt); + *damage_nt_it = std::min(*delta_nt_max_it / *delta_c_it, Real(1.)); + + *delta_st_max_it = std::max(*delta_st_max_it, delta_st); + *damage_st_it = std::min(*delta_st_max_it / delta_c_s, 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$ + */ + + Vector traction_n(spatial_dimension); + Vector traction_s(spatial_dimension); + + /// NORMAL TRACTIONS + if (Math::are_float_equal(*damage_nt_it, 1.)) + traction_n.clear(); + else if (Math::are_float_equal(*damage_nt_it, 0.)) { + if (penetration) + traction_n.clear(); + else + traction_n = *insertion_stress_it; + } + else { + traction_n = normal_opening; + + AKANTU_DEBUG_ASSERT(*delta_nt_max_it != 0., + "Division by zero, tolerance might be too low"); + + traction_n *= *sigma_c_it / (*delta_nt_max_it) * (1. - *damage_nt_it); + } + + /// SHEAR TRACTIONS + if (Math::are_float_equal(*damage_st_it, 1.)) + traction_s.clear(); + else if (Math::are_float_equal(*damage_st_it, 0.)) { + traction_s.clear(); + } + else { + traction_s = tangential_opening; + + AKANTU_DEBUG_ASSERT(*delta_st_max_it != 0., + "Division by zero, tolerance might be too low"); + + traction_s *= *sigma_c_it * this->beta / (*delta_st_max_it) * (1. - *damage_st_it); + } + + *traction_it = traction_n + traction_s; + + } + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +template +void MaterialCohesiveLinearUncoupled::checkDeltaMax(GhostType ghost_type) { + AKANTU_DEBUG_IN(); + + /** + * This function set a predefined value to the parameter + * delta_max_prev of the elements that have been inserted in the + * last loading step for which convergence has not been + * reached. This is done before reducing the loading and re-doing + * the step. Otherwise, the updating of delta_max_prev would be + * done with reference to the non-convergent solution. In this + * function also other variables related to the contact and + * friction behavior are correctly set. + */ + + Mesh & mesh = this->fem_cohesive->getMesh(); + Mesh::type_iterator it = mesh.firstType(spatial_dimension, + ghost_type, _ek_cohesive); + Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, + ghost_type, _ek_cohesive); + + /** + * the variable "recompute" is set to true to activate the + * procedure that reduce the penalty parameter for + * compression. This procedure is set true only during the phase of + * load_reduction, that has to be set in the maiin file. The + * penalty parameter will be reduced only for the elements having + * reduction_penalty = true. + */ + this->recompute = true; + + for(; it != last_type; ++it) { + Array & elem_filter = this->element_filter(*it, ghost_type); + + UInt nb_element = elem_filter.getSize(); + if (nb_element == 0) continue; + + ElementType el_type = *it; + + /// define iterators + Array::scalar_iterator delta_nt_max_it = + delta_nt_max(el_type, ghost_type).begin(); + + Array::scalar_iterator delta_nt_max_end = + delta_nt_max(el_type, ghost_type).end(); + + Array::scalar_iterator delta_nt_max_prev_it = + delta_nt_max.previous(el_type, ghost_type).begin(); + + Array::scalar_iterator delta_st_max_it = + delta_st_max(el_type, ghost_type).begin(); + + Array::scalar_iterator delta_st_max_prev_it = + delta_st_max.previous(el_type, ghost_type).begin(); + + Array::scalar_iterator delta_c_it = + this-> delta_c_eff(el_type, ghost_type).begin(); + + Array::vector_iterator opening_prec_it = + this->opening_prec(el_type, ghost_type).begin(spatial_dimension); + + Array::vector_iterator opening_prec_prev_it = + this->opening_prec.previous(el_type, ghost_type).begin(spatial_dimension); + + /// loop on each quadrature point + for (; delta_nt_max_it != delta_nt_max_end; + ++delta_nt_max_it, ++delta_st_max_it, ++delta_c_it, + ++delta_nt_max_prev_it, ++delta_st_max_prev_it, + ++opening_prec_it, ++opening_prec_prev_it) { + + if (*delta_nt_max_prev_it == 0) + /// elements inserted in the last incremental step, that did + /// not converge + *delta_nt_max_it = *delta_c_it / 1000; + else + /// elements introduced in previous incremental steps, for + /// which a correct value of delta_max_prev already exists + *delta_nt_max_it = *delta_nt_max_prev_it; + + if (*delta_st_max_prev_it == 0) + *delta_st_max_it = *delta_c_it * this->kappa / this->beta / 1000; + else + *delta_st_max_it = *delta_st_max_prev_it; + + /// in case convergence is not reached, set opening_prec to the + /// value referred to the previous incremental step + *opening_prec_it = *opening_prec_prev_it; + } + } +} + +/* -------------------------------------------------------------------------- */ +template +void MaterialCohesiveLinearUncoupled::computeTangentTraction(const ElementType & el_type, + Array & tangent_matrix, + const Array & normal, + GhostType ghost_type) { + AKANTU_DEBUG_IN(); + + /// define iterators + Array::matrix_iterator tangent_it = + tangent_matrix.begin(spatial_dimension, spatial_dimension); + + Array::matrix_iterator tangent_end = + tangent_matrix.end(spatial_dimension, spatial_dimension); + + Array::const_vector_iterator normal_it = + normal.begin(spatial_dimension); + + Array::vector_iterator opening_it = + this->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 + Array::scalar_iterator delta_nt_max_it = + this->delta_nt_max.previous(el_type, ghost_type).begin(); + + Array::scalar_iterator delta_st_max_it = + this->delta_st_max.previous(el_type, ghost_type).begin(); + + Array::scalar_iterator sigma_c_it = + this->sigma_c_eff(el_type, ghost_type).begin(); + + Array::scalar_iterator delta_c_it = + this->delta_c_eff(el_type, ghost_type).begin(); + + // Array::scalar_iterator damage_it = + // damage(el_type, ghost_type).begin(); + + Array::vector_iterator contact_opening_it = + this->contact_opening(el_type, ghost_type).begin(spatial_dimension); + + Array::scalar_iterator reduction_penalty_it = + this->reduction_penalty(el_type, ghost_type).begin(); + + Vector normal_opening(spatial_dimension); + Vector tangential_opening(spatial_dimension); + + for (; tangent_it != tangent_end; + ++tangent_it, ++normal_it, ++opening_it, ++sigma_c_it, + ++delta_c_it, ++delta_nt_max_it, ++delta_st_max_it, + ++contact_opening_it, ++reduction_penalty_it) { + + Real normal_opening_norm, tangential_opening_norm; + bool penetration; + Real current_penalty = 0.; + Real delta_c_s = this->kappa / this->beta * (*delta_c_it); + Real delta_c_s2_R2 = delta_c_s * delta_c_s / (R * R); + + /** + * 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_it += *contact_opening_it; + + /// compute normal and tangential opening vectors + normal_opening_norm = opening_it->dot(*normal_it); + Vector normal_opening = *normal_it; + normal_opening *= normal_opening_norm; + + Vector tangential_opening = *opening_it; + tangential_opening -= normal_opening; + tangential_opening_norm = tangential_opening.norm(); + + Real delta_nt = tangential_opening_norm * tangential_opening_norm * this->beta2_kappa2; + Real delta_st = tangential_opening_norm * tangential_opening_norm; + + 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 tangent_n(spatial_dimension, spatial_dimension); + // Matrix tangent_s(spatial_dimension, spatial_dimension); + + /// TANGENT STIFFNESS FOR NORMAL TRACTIONS + Matrix n_outer_n(spatial_dimension, spatial_dimension); + n_outer_n.outerProduct(*normal_it, *normal_it); + + if (penetration){ + if (this->recompute && *reduction_penalty_it) + current_penalty = this->penalty / 100.; + else + current_penalty = this->penalty; + + /// stiffness in compression given by the penalty parameter + *tangent_it = n_outer_n; + *tangent_it *= current_penalty; + + *opening_it = tangential_opening; + normal_opening.clear(); + } + else{ + delta_nt += normal_opening_norm * normal_opening_norm; + delta_nt = std::sqrt(delta_nt); + + delta_st += normal_opening_norm * normal_opening_norm * delta_c_s2_R2; + + /** + * 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_nt < Math::getTolerance()) + delta_nt = *delta_c_it / 1000.; + + if (delta_nt >= *delta_nt_max_it){ + if (delta_nt <= *delta_c_it){ + derivative = - (*sigma_c_it) / (delta_nt * delta_nt); + t = *sigma_c_it * (1 - delta_nt / *delta_c_it); + } else { + derivative = 0.; + t = 0.; + } + } else if (delta_nt < *delta_nt_max_it){ + Real tmax = *sigma_c_it * (1 - *delta_nt_max_it / *delta_c_it); + t = tmax / *delta_nt_max_it * delta_nt; + } + + /// computation of the derivative of the constitutive law (dT/ddelta) + Matrix nn(n_outer_n); + nn *= t / delta_nt; + + Vector t_tilde(normal_opening); + t_tilde *= (1. - this->beta2_kappa2); + Vector mm(*opening_it); + mm *= this->beta2_kappa2; + t_tilde += mm; + + Vector t_hat(normal_opening); + Matrix prov(spatial_dimension, spatial_dimension); + prov.outerProduct(t_hat, t_tilde); + prov *= derivative / delta_nt; + prov += nn; + + *tangent_it = prov.transpose(); + } + + derivative = 0.; + t = 0.; + /// TANGENT STIFFNESS FOR SHEAR TRACTIONS + delta_st = std::sqrt(delta_st); + + /** + * 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_st < Math::getTolerance()) + delta_st = delta_c_s / 1000.; + + if (delta_st >= *delta_st_max_it){ + if (delta_st <= delta_c_s){ + derivative = - this->beta * (*sigma_c_it) / (delta_st * delta_st); + t = this->beta * (*sigma_c_it) * (1 - delta_st / delta_c_s); + } else { + derivative = 0.; + t = 0.; + } + } else if (delta_st < *delta_st_max_it){ + Real tmax = this->beta * (*sigma_c_it) * (1 - *delta_st_max_it / delta_c_s); + t = tmax / *delta_st_max_it * delta_st; + } + + /// computation of the derivative of the constitutive law (dT/ddelta) + Matrix I(spatial_dimension, spatial_dimension); + I.eye(); + Matrix nn(n_outer_n); + I -= nn; + I *= t / delta_st; + + Vector t_tilde(normal_opening); + t_tilde *= (delta_c_s2_R2 - 1.); + Vector mm(*opening_it); + t_tilde += mm; + + Vector t_hat(tangential_opening); + Matrix prov(spatial_dimension, spatial_dimension); + prov.outerProduct(t_hat, t_tilde); + prov *= derivative / delta_st; + prov += I; + + Matrix prov_t = prov.transpose(); + *tangent_it += prov_t; + } + + AKANTU_DEBUG_OUT(); +} +/* -------------------------------------------------------------------------- */ + + +INSTANTIATE_MATERIAL(MaterialCohesiveLinearUncoupled); + +__END_AKANTU__ diff --git a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_uncoupled.hh b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_uncoupled.hh new file mode 100644 index 000000000..b81d0f7e6 --- /dev/null +++ b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_uncoupled.hh @@ -0,0 +1,115 @@ +/** + * @file material_cohesive_linear.hh + * + * @author Mauro Corrado + * + * @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_linear.hh" + +/* -------------------------------------------------------------------------- */ +#ifndef __AKANTU_MATERIAL_COHESIVE_LINEAR_UNCOUPLED_HH__ +#define __AKANTU_MATERIAL_COHESIVE_LINEAR_UNCOUPLED_HH__ + +/* -------------------------------------------------------------------------- */ + +__BEGIN_AKANTU__ + +/** + * Cohesive material linear with two different laws for mode I and + * mode II, for extrinsic case + * + * parameters in the material files : + * - roughness : define the interaction between mode I and mode II (default: 0) + */ +template +class MaterialCohesiveLinearUncoupled : public MaterialCohesiveLinear { + /* ------------------------------------------------------------------------ */ + /* Constructors/Destructors */ + /* ------------------------------------------------------------------------ */ + // typedef MaterialCohesiveLinear MaterialParent; +public: + + MaterialCohesiveLinearUncoupled(SolidMechanicsModel & model, const ID & id = ""); + + /* ------------------------------------------------------------------------ */ + /* Methods */ + /* ------------------------------------------------------------------------ */ +public: + + /// initialize the material parameters + virtual void initMaterial(); + +protected: + + /// constitutive law + virtual void computeTraction(const Array & normal, + ElementType el_type, + GhostType ghost_type = _not_ghost); + + /// check delta_max for cohesive elements in case of no convergence + /// in the solveStep (only for extrinsic-implicit) + virtual void checkDeltaMax(GhostType ghost_type = _not_ghost); + + /// compute tangent stiffness matrix + virtual void computeTangentTraction(const ElementType & el_type, + Array & tangent_matrix, + const Array & normal, + GhostType ghost_type); + + /* ------------------------------------------------------------------------ */ + /* Accessors */ + /* ------------------------------------------------------------------------ */ +public: + + /* ------------------------------------------------------------------------ */ + /* Class Members */ + /* ------------------------------------------------------------------------ */ +protected: + /// parameter to tune the interaction between mode II and mode I + Real R; + + /// maximum normal displacement + CohesiveInternalField delta_nt_max; + + /// maximum tangential displacement + CohesiveInternalField delta_st_max; + + /// damage associated to normal tractions + CohesiveInternalField damage_nt; + + /// damage associated to shear tractions + CohesiveInternalField damage_st; + +}; + +__END_AKANTU__ + +#endif /* __AKANTU_MATERIAL_COHESIVE_LINEAR_UNCOUPLED_HH__ */ diff --git a/src/model/solid_mechanics/materials/material_cohesive_includes.hh b/src/model/solid_mechanics/materials/material_cohesive_includes.hh index 00bfb1433..4ad8a634d 100644 --- a/src/model/solid_mechanics/materials/material_cohesive_includes.hh +++ b/src/model/solid_mechanics/materials/material_cohesive_includes.hh @@ -1,47 +1,49 @@ /** * @file material_cohesive_includes.hh * * @author Nicolas Richart * * @date creation: Sun Sep 26 2010 * @date last modification: Tue Jan 12 2016 * * @brief List of includes for cohesive elements * * @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 . * */ /* -------------------------------------------------------------------------- */ #ifndef AKANTU_CMAKE_LIST_MATERIALS #include "material_cohesive.hh" #include "material_cohesive_linear.hh" #include "material_cohesive_linear_fatigue.hh" #include "material_cohesive_linear_friction.hh" +#include "material_cohesive_linear_uncoupled.hh" #include "material_cohesive_bilinear.hh" #include "material_cohesive_exponential.hh" #endif #define AKANTU_COHESIVE_MATERIAL_LIST \ ((2, (cohesive_linear, MaterialCohesiveLinear ))) \ ((2, (cohesive_linear_fatigue, MaterialCohesiveLinearFatigue ))) \ ((2, (cohesive_linear_friction, MaterialCohesiveLinearFriction))) \ + ((2, (cohesive_linear_uncoupled, MaterialCohesiveLinearUncoupled))) \ ((2, (cohesive_bilinear , MaterialCohesiveBilinear ))) \ ((2, (cohesive_exponential , MaterialCohesiveExponential )))