diff --git a/packages/cohesive_element.cmake b/packages/cohesive_element.cmake index 14b397cee..b3b99a8bd 100644 --- a/packages/cohesive_element.cmake +++ b/packages/cohesive_element.cmake @@ -1,121 +1,123 @@ #=============================================================================== # @file 20_cohesive_element.cmake # # @author Marco Vocialta # @author Nicolas Richart # # @date creation: Tue Oct 16 2012 # @date last modification: Tue Sep 02 2014 # # @brief package description for cohesive elements # # @section LICENSE # # Copyright (©) 2014 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_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.cc b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear.cc index 9ddf21d73..325b064e9 100644 --- a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear.cc +++ b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear.cc @@ -1,984 +1,753 @@ /** * @file material_cohesive_linear.cc * * @author Marco Vocialta * @author Mauro Corrado * * @date creation: Tue May 08 2012 * @date last modification: Thu Nov 19 2015 * * @brief Linear irreversible cohesive law of mixed mode loading with * random stress definition for extrinsic type * * @section LICENSE * * Copyright (©) 2014 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.hh" #include "solid_mechanics_model_cohesive.hh" #include "sparse_matrix.hh" #include "dof_synchronizer.hh" __BEGIN_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), opening_prec("opening_prec", *this), residual_sliding("residual_sliding", *this), friction_force("friction_force", *this), reduction_penalty("reduction_penalty", *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("friction", friction, false, _pat_parsable | _pat_readable, "Activation of friction in case of contact"); 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"); this->registerParam("recompute", recompute, false, _pat_parsable | _pat_modifiable, "recompute solution"); - use_previous_delta_max = true; + this->use_previous_delta_max = true; + this->use_previous_opening = true; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinear::initMaterial() { AKANTU_DEBUG_IN(); MaterialCohesive::initMaterial(); /// 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; sigma_c_eff.initialize(1); delta_c_eff.initialize(1); insertion_stress.initialize(spatial_dimension); opening_prec.initialize(spatial_dimension); friction_force.initialize(spatial_dimension); residual_sliding.initialize(1); reduction_penalty.initialize(1); if (!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(); residual_sliding.initializeHistory(); opening_prec.initializeHistory(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ 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 FEEngine & fe_engine = model->getFEEngine(); const FEEngine & fe_engine_facet = model->getFEEngine("FacetsFEEngine"); // loop over facet type Mesh::type_iterator first = mesh_facets.firstType(spatial_dimension - 1); Mesh::type_iterator last = mesh_facets.lastType(spatial_dimension - 1); Real base_sigma_c = sigma_c; for(;first != last; ++first) { ElementType type_facet = *first; const Array< std::vector > & facet_to_element = mesh_facets.getElementToSubelement(type_facet); UInt nb_facet = facet_to_element.getSize(); UInt nb_quad_per_facet = fe_engine_facet.getNbIntegrationPoints(type_facet); // iterator to modify sigma_c for all the quadrature points of a facet Array::vector_iterator 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; std::vector::const_iterator elem = element_list.begin(); std::vector::const_iterator 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(); Real tolerance = Math::getTolerance(); Mesh::type_iterator it = mesh_facets.firstType(spatial_dimension - 1); Mesh::type_iterator last = mesh_facets.lastType(spatial_dimension - 1); for (; it != last; ++it) { ElementType type_facet = *it; ElementType type_cohesive = FEEngine::getCohesiveElementType(type_facet); const Array & facets_check = inserter.getCheckFacets(type_facet); Array & f_insertion = inserter.getInsertionFacets(type_facet); Array & f_filter = facet_filter(type_facet); Array & sig_c_eff = sigma_c_eff(type_cohesive); Array & del_c = delta_c_eff(type_cohesive); Array & ins_stress = insertion_stress(type_cohesive); Array & trac_old = tractions_old(type_cohesive); Array & open_prec = opening_prec(type_cohesive); - Array & res_sliding = residual_sliding(type_cohesive); Array & red_penalty = reduction_penalty(type_cohesive); const Array & f_stress = model->getStressOnFacets(type_facet); const Array & sigma_lim = sigma_c(type_facet); Real max_ratio = 0.; UInt index_f = 0; UInt index_filter = 0; UInt nn = 0; - UInt nb_quad_facet = model->getFEEngine("FacetsFEEngine").getNbIntegrationPoints(type_facet); UInt nb_facet = f_filter.getSize(); // if (nb_facet == 0) continue; Array::const_iterator 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 Array & tangents = model->getTangents(type_facet); const Array & normals = model->getFEEngine("FacetsFEEngine").getNormalsOnIntegrationPoints(type_facet); Array::const_vector_iterator normal_begin = normals.begin(spatial_dimension); Array::const_vector_iterator tangent_begin = tangents.begin(tangents.getNbComponent()); Array::const_matrix_iterator facet_stress_begin = f_stress.begin(spatial_dimension, spatial_dimension * 2); std::vector new_sigmas; std::vector< 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 - tolerance)) { if (model->isExplicit()){ f_insertion(facet) = true; if (!check_only) { // 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); } } }else{ Real ratio = final_stress/(*sigma_lim_it); if (ratio > max_ratio){ ++nn; max_ratio = ratio; index_f = f; index_filter = f_filter(f); } } } } /// Insertion of only 1 cohesive element in case of implicit approach. The one subjected to the highest stress. if (!model->isExplicit()){ StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); Array abs_max(comm.getNbProc()); abs_max(comm.whoAmI()) = max_ratio; comm.allGather(abs_max.storage(), 1); Array::scalar_iterator it = std::max_element(abs_max.begin(), abs_max.end()); Int pos = it - abs_max.begin(); if (pos != comm.whoAmI()) { AKANTU_DEBUG_OUT(); return; } if (nn) { f_insertion(index_filter) = true; if (!check_only) { // Array::iterator > normal_traction_it = // normal_traction.begin_reinterpret(nb_quad_facet, spatial_dimension, nb_facet); Array::const_iterator sigma_lim_it = sigma_lim.begin(); for (UInt q = 0; q < nb_quad_facet; ++q) { // Vector ins_s(normal_traction_it[index_f].storage() + q * spatial_dimension, // spatial_dimension); Real new_sigma = (sigma_lim_it[index_f]); Vector normal_traction_vec(spatial_dimension, 0.0); 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 = delta_c; else new_delta = 2 * G_c / (new_sigma); new_delta_c.push_back(new_delta); } } } } // update material data for the new elements UInt old_nb_quad_points = sig_c_eff.getSize(); 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); open_prec.resize(old_nb_quad_points + new_nb_quad_points); - res_sliding.resize(old_nb_quad_points + new_nb_quad_points); red_penalty.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]; - res_sliding(old_nb_quad_points + q) = 0.; red_penalty(old_nb_quad_points + q) = false; 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); open_prec(old_nb_quad_points + q, dim) = 0.; } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinear::computeTraction(const Array & normal, ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); /// define iterators Array::vector_iterator traction_it = tractions(el_type, ghost_type).begin(spatial_dimension); Array::vector_iterator opening_it = opening(el_type, ghost_type).begin(spatial_dimension); + Array::vector_iterator opening_prev_it = + opening.previous(el_type, ghost_type).begin(spatial_dimension); + Array::vector_iterator opening_prec_it = opening_prec(el_type, ghost_type).begin(spatial_dimension); - Array::vector_iterator opening_prec_prev_it = - opening_prec.previous(el_type, ghost_type).begin(spatial_dimension); - Array::vector_iterator contact_traction_it = contact_tractions(el_type, ghost_type).begin(spatial_dimension); Array::vector_iterator contact_opening_it = contact_opening(el_type, ghost_type).begin(spatial_dimension); - Array::const_iterator< Vector > normal_it = + Array::const_vector_iterator normal_it = normal.begin(spatial_dimension); Array::vector_iterator traction_end = tractions(el_type, ghost_type).end(spatial_dimension); - Array::iteratorsigma_c_it = + Array::scalar_iterator sigma_c_it = sigma_c_eff(el_type, ghost_type).begin(); - Array::iteratordelta_max_it = + Array::scalar_iterator delta_max_it = delta_max(el_type, ghost_type).begin(); - Array::iteratordelta_max_prev_it = + Array::scalar_iterator delta_max_prev_it = delta_max.previous(el_type, ghost_type).begin(); - Array::iteratordelta_c_it = + Array::scalar_iterator delta_c_it = delta_c_eff(el_type, ghost_type).begin(); - Array::iteratordamage_it = + Array::scalar_iterator damage_it = damage(el_type, ghost_type).begin(); Array::vector_iterator insertion_stress_it = insertion_stress(el_type, ghost_type).begin(spatial_dimension); - Array::iteratorres_sliding_it = + Array::scalar_iterator res_sliding_it = residual_sliding(el_type, ghost_type).begin(); - Array::iteratorres_sliding_prev_it = + Array::scalar_iterator res_sliding_prev_it = residual_sliding.previous(el_type, ghost_type).begin(); Array::vector_iterator friction_force_it = friction_force(el_type, ghost_type).begin(spatial_dimension); - Array::iterator reduction_penalty_it = + Array::scalar_iterator reduction_penalty_it = reduction_penalty(el_type, ghost_type).begin(); - Real * memory_space = new Real[2*spatial_dimension]; - Vector normal_opening(memory_space, spatial_dimension); - Vector tangential_opening(memory_space + spatial_dimension, - spatial_dimension); - - // Vector open_prec(spatial_dimension); - // Vector open(spatial_dimension); + 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)); + /// loop on each quadrature point for (; traction_it != traction_end; - ++traction_it, ++opening_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, ++opening_prec_prev_it, - ++friction_force_it, ++reduction_penalty_it) { - - if (!model->isExplicit()) - *delta_max_it = *delta_max_prev_it; - - // check if the local displacement increment is too large. If so, - // reduce arbitrary it in order to prevent problems of - // instability - - Real normal_opening_prec_norm = opening_prec_it->dot(*normal_it); - -// open_prec = *opening_prec_it; -// Real open_prec_norm = open_prec.norm(); -// open = *opening_it; -// Real open_norm = open.norm(); -// Real delta_open = std::abs(open_norm - open_prec_norm); -// -// // if (delta_open > 1.0e-01){ -// // *opening_it = *opening_prec_it + ((*opening_it - *opening_prec_it) / 50.0); -// // } -// *opening_prec_it = *opening_it; - - /// compute normal and tangential opening vectors - Real normal_opening_norm = opening_it->dot(*normal_it); - normal_opening = (*normal_it); - normal_opening *= normal_opening_norm; - tangential_opening = *opening_it; - tangential_opening -= normal_opening; - - Real 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 * beta2_kappa2; - - bool penetration = normal_opening_norm < -Math::getTolerance(); - if (contact_after_breaking == false && Math::are_float_equal(*damage_it, 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 (!model->isExplicit() && !recompute) - if ((normal_opening_prec_norm * normal_opening_norm) < -Math::getTolerance()){ - *reduction_penalty_it = true; - } - - if (penetration) { - Real current_penalty = 0.; - if (recompute && *reduction_penalty_it){ - /// the penalty parameter is locally reduced - current_penalty = penalty / 1000.; - } - else - current_penalty = penalty; - - /// use penalty coefficient in case of penetration - *contact_traction_it = normal_opening; - *contact_traction_it *= current_penalty; - *contact_opening_it = normal_opening; - - /* ----------------------------------------- */ - /// ADD FRICTION - if (friction){ - Real damage = std::min(*delta_max_prev_it / *delta_c_it, Real(1.)); - Real mu = mu_max * std::sqrt(damage); - Real normal_opening_prec_norm = opening_prec_prev_it->dot(*normal_it); - Vector normal_opening_prec = (*normal_it); - normal_opening_prec *= normal_opening_prec_norm; - Real tau_max = mu * current_penalty * (normal_opening_prec.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_it) < -Math::getTolerance()) - tau = -tau; - // from tau get the x and y components of friction, to be added in the force vector - Vector tangent(spatial_dimension); - tangent = tangential_opening / tangential_opening_norm; - *friction_force_it = tau * tangent; - - //update residual_sliding - *res_sliding_it = tangential_opening_norm - (tau / friction_penalty); - } - /// end add friction - /* ----------------------------------------- */ - - /// don't consider penetration contribution for delta - *opening_it = tangential_opening; - normal_opening.clear(); - - } - else { - delta += normal_opening_norm * normal_opening_norm; - contact_traction_it->clear(); - contact_opening_it->clear(); - friction_force_it->clear(); - } - - delta = std::sqrt(delta); - - /// update maximum displacement and damage - *delta_max_it = std::max(*delta_max_it, delta); - *damage_it = std::min(*delta_max_it / *delta_c_it, 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_it, 1.)) - traction_it->clear(); - else if (Math::are_float_equal(*damage_it, 0.)) { - if (penetration) - traction_it->clear(); - else - *traction_it = *insertion_stress_it; - } - else { - *traction_it = tangential_opening; - *traction_it *= beta2_kappa; - *traction_it += normal_opening; - - AKANTU_DEBUG_ASSERT(*delta_max_it != 0., - "Division by zero, tolerance might be too low"); - - *traction_it *= *sigma_c_it / *delta_max_it * (1. - *damage_it); - } - - if (friction) - *traction_it += *friction_force_it; - + ++traction_it, ++opening_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, ++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); } - delete [] memory_space; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinear::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 = 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. recompute = true; for(; it != last_type; ++it) { Array & elem_filter = 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 = + Array::scalar_iterator delta_max_it = delta_max(el_type, ghost_type).begin(); - Array::iteratordelta_max_end = + Array::scalar_iterator delta_max_end = delta_max(el_type, ghost_type).end(); - Array::iteratordelta_max_prev_it = + Array::scalar_iterator delta_max_prev_it = delta_max.previous(el_type, ghost_type).begin(); - Array::iteratordelta_c_it = + Array::scalar_iterator delta_c_it = delta_c_eff(el_type, ghost_type).begin(); -// Array::vector_iterator opening_prec_it = -// opening_prec(el_type, ghost_type).begin(spatial_dimension); -// -// Array::vector_iterator opening_prec_prev_it = -// opening_prec.previous(el_type, ghost_type).begin(spatial_dimension); + Array::vector_iterator opening_prec_it = + opening_prec(el_type, ghost_type).begin(spatial_dimension); - Array::iteratorres_sliding_it = + Array::vector_iterator opening_prec_prev_it = + opening_prec.previous(el_type, ghost_type).begin(spatial_dimension); + + Array::scalar_iterator res_sliding_it = residual_sliding(el_type, ghost_type).begin(); - Array::iteratorres_sliding_prev_it = + Array::scalar_iterator res_sliding_prev_it = residual_sliding.previous(el_type, ghost_type).begin(); /// loop on each quadrature point for (; delta_max_it != delta_max_end; ++delta_max_it, ++delta_max_prev_it, ++delta_c_it, - // ++opening_prec_it, ++opening_prec_prev_it, + ++opening_prec_it, ++opening_prec_prev_it, ++res_sliding_it, ++res_sliding_prev_it) { if (*delta_max_prev_it == 0) /// elements inserted in the last incremental step, that did /// not converge *delta_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_max_it = *delta_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; + *opening_prec_it = *opening_prec_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 MaterialCohesiveLinear::resetVariables(GhostType ghost_type) { AKANTU_DEBUG_IN(); /// This function set the variables "recompute" and /// "reduction_penalty" to false. It is called by solveStepCohesive /// when convergence is reached. Such variables, in fact, have to be /// false at the beginning of a new incremental step. Mesh & mesh = 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); recompute = false; // std::cout << "RESET VARIABLE" << std::endl; for(; it != last_type; ++it) { Array & elem_filter = element_filter(*it, ghost_type); UInt nb_element = elem_filter.getSize(); if (nb_element == 0) continue; ElementType el_type = *it; - Array::iterator reduction_penalty_it = + Array::scalar_iterator reduction_penalty_it = reduction_penalty(el_type, ghost_type).begin(); - Array::iterator reduction_penalty_end = + Array::scalar_iterator reduction_penalty_end = reduction_penalty(el_type, ghost_type).end(); /// loop on each quadrature point for (; reduction_penalty_it != reduction_penalty_end; ++reduction_penalty_it) { *reduction_penalty_it = false; } } } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinear::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 = opening(el_type, ghost_type).begin(spatial_dimension); - Array::vector_iterator opening_prec_it = - 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 = + Array::scalar_iterator delta_max_it = delta_max.previous(el_type, ghost_type).begin(); - Array::iteratorsigma_c_it = + Array::scalar_iterator sigma_c_it = sigma_c_eff(el_type, ghost_type).begin(); - Array::iteratordelta_c_it = + Array::scalar_iterator delta_c_it = delta_c_eff(el_type, ghost_type).begin(); - Array::iteratordamage_it = + Array::scalar_iterator damage_it = damage(el_type, ghost_type).begin(); - Array::iterator< Vector > contact_opening_it = + Array::vector_iterator contact_opening_it = contact_opening(el_type, ghost_type).begin(spatial_dimension); - Array::iteratorres_sliding_prev_it = + Array::scalar_iterator res_sliding_prev_it = residual_sliding.previous(el_type, ghost_type).begin(); - Array::iterator reduction_penalty_it = + Array::scalar_iterator reduction_penalty_it = 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_it, + ++tangent_it, ++normal_it, ++opening_it, ++delta_max_it, ++sigma_c_it, ++delta_c_it, ++damage_it, - ++contact_opening_it, ++res_sliding_prev_it, ++reduction_penalty_it) { - - /// 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 - Real normal_opening_norm = opening_it->dot(*normal_it); - normal_opening = (*normal_it); - normal_opening *= normal_opening_norm; - - tangential_opening = *opening_it; - tangential_opening -= normal_opening; - - Real tangential_opening_norm = tangential_opening.norm(); - bool penetration = normal_opening_norm < -Math::getTolerance(); - if (contact_after_breaking == false && Math::are_float_equal(*damage_it, 1.)) - penetration = false; - - Real derivative = 0; // derivative = d(t/delta)/ddelta - Real t = 0; - - Real delta = tangential_opening_norm * tangential_opening_norm * beta2_kappa2; - - Matrix n_outer_n(spatial_dimension, spatial_dimension); - n_outer_n.outerProduct(*normal_it, *normal_it); - - if (penetration){ - Real current_penalty = 0.; - if (recompute && *reduction_penalty_it) - current_penalty = penalty / 1000.; - else - current_penalty = penalty; - - /// stiffness in compression given by the penalty parameter - *tangent_it += n_outer_n; - *tangent_it *= current_penalty; - - /* ------------------------------------- */ - /// ADD FRICTION - if (friction){ - Real damage = std::min(*delta_max_it / *delta_c_it, Real(1.)); - // the friction coefficient mu is increasing with the - // damage. It equals the maximum value when damage = 1. - Real mu = mu_max * damage; - Real normal_opening_prec_norm = opening_prec_it->dot(*normal_it); - Vector normal_opening_prec = (*normal_it); - normal_opening_prec *= normal_opening_prec_norm; - Real tau_max = mu * current_penalty * normal_opening_prec.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 nn(n_outer_n); - I -= nn; - *tangent_it += I * friction_penalty; - } - } - /// end add friction - /* ------------------------------------- */ - - *opening_it = tangential_opening; - normal_opening_norm = opening_it->dot(*normal_it); - normal_opening = (*normal_it); - 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_it)/1000.; - - if (delta >= *delta_max_it){ - if (delta <= *delta_c_it){ - derivative = -*sigma_c_it/(delta * delta); - t = *sigma_c_it * (1 - delta / *delta_c_it); - } else { - derivative = 0.; - t = 0.; - } - } else if (delta < *delta_max_it){ - Real tmax = *sigma_c_it * (1 - *delta_max_it / *delta_c_it); - t = tmax / *delta_max_it * delta; - } - - - /// computation of the derivative of the constitutive law (dT/ddelta) - Matrix I(spatial_dimension, spatial_dimension); - I.eye(beta2_kappa); - Matrix nn(n_outer_n); - nn *= (1 - beta2_kappa); - nn += I; - nn *= t/delta; - Vector t_tilde(normal_opening); - t_tilde *= (1 - beta2_kappa2); - Vector mm(*opening_it); - mm *= beta2_kappa2; - t_tilde += mm; - Vector t_hat(normal_opening); - t_hat += beta2_kappa * tangential_opening; - Matrix prov(spatial_dimension, spatial_dimension); - prov.outerProduct(t_hat, t_tilde); - prov *= derivative/delta; - prov += nn; - Matrix tmp(spatial_dimension, spatial_dimension); - for (UInt i = 0; i < spatial_dimension; ++i) { - for (UInt j = 0; j < spatial_dimension; ++j) { - tmp(j,i) = prov(i,j); - } - } - *tangent_it += tmp; + ++contact_opening_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); /// 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(MaterialCohesiveLinear); __END_AKANTU__ diff --git a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear.hh b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear.hh index fbe80397b..34b417771 100644 --- a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear.hh +++ b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear.hh @@ -1,208 +1,227 @@ /** * @file material_cohesive_linear.hh * * @author Marco Vocialta * * @date creation: Tue May 08 2012 * @date last modification: Tue Jul 29 2014 * * @brief Linear irreversible cohesive law of mixed mode loading with * random stress definition for extrinsic type * * @section LICENSE * * Copyright (©) 2014 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__ /* -------------------------------------------------------------------------- */ __BEGIN_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 virtual void initMaterial(); /// check stress for cohesive elements' insertion virtual void checkInsertion(bool check_only = false); /// 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 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) - void checkDeltaMax(GhostType ghost_type = _not_ghost); + virtual void checkDeltaMax(GhostType ghost_type = _not_ghost); /// reset variables when convergence is reached (only for /// extrinsic-implicit) void resetVariables(GhostType ghost_type = _not_ghost); /// compute tangent stiffness matrix - void computeTangentTraction(const ElementType & el_type, - Array & tangent_matrix, - const Array & normal, - GhostType ghost_type); + virtual void computeTangentTraction(const ElementType & el_type, + Array & tangent_matrix, + const Array & normal, + GhostType ghost_type); /** * 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, Real & delta_max, + const Real & delta_max_prev, const Real & delta_c, + const Vector & insertion_stress, const Real & sigma_c, + Vector & opening, Vector & opening_prec, const Vector & normal, + Vector & normal_opening, Vector & tangential_opening, + Real & normal_opening_norm, Real & tangential_opening_norm, Real & damage, + bool & penetration, bool & reduction_penalty, Real & current_penalty, + 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, bool & reduction_penalty, Real & current_penalty, + 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; /// maximum value of the friction coefficient Real mu_max; /// penalty parameter for the friction law Real friction_penalty; /// 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; /// delta of the previous step CohesiveInternalField opening_prec; /// history parameter for the friction law CohesiveInternalField residual_sliding; /// friction force CohesiveInternalField friction_force; /// variable that defines if the penalty parameter for compression /// has to be decreased for problems of convergence in the solution /// loops CohesiveInternalField reduction_penalty; /// 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; /// variable saying if friction has to be added to the cohesive behavior bool friction; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ -#include "material_cohesive_linear_inline_impl.cc" - __END_AKANTU__ +#include "material_cohesive_linear_inline_impl.cc" + #endif /* __AKANTU_MATERIAL_COHESIVE_LINEAR_HH__ */ diff --git a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_fatigue.cc b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_fatigue.cc index 7634f8b39..cfa3635c4 100644 --- a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_fatigue.cc +++ b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_fatigue.cc @@ -1,293 +1,293 @@ /** * @file material_cohesive_linear_fatigue.cc * @author Marco Vocialta * @date Thu Feb 19 14:40:57 2015 * * @brief See material_cohesive_linear_fatigue.hh for information * * @section LICENSE * * Copyright (©) 2010-2011 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_fatigue.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template MaterialCohesiveLinearFatigue ::MaterialCohesiveLinearFatigue(SolidMechanicsModel & model, const ID & id) : MaterialCohesiveLinear(model, id), delta_prec("delta_prec", *this), K_plus("K_plus", *this), K_minus("K_minus", *this), T_1d("T_1d", *this), switches("switches", *this), delta_dot_prec("delta_dot_prec", *this), normal_regime("normal_regime", *this) { this->registerParam("delta_f", delta_f, Real(-1.) , _pat_parsable | _pat_readable, "delta_f"); this->registerParam("progressive_delta_f", progressive_delta_f, false, _pat_parsable | _pat_readable, "Whether or not delta_f is equal to delta_max"); this->registerParam("count_switches", count_switches, false, _pat_parsable | _pat_readable, "Count the opening/closing switches per element"); this->registerParam("fatigue_ratio", fatigue_ratio, Real(1.), _pat_parsable | _pat_readable, "What portion of the cohesive law is subjected to fatigue"); } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinearFatigue::initMaterial() { MaterialCohesiveLinear::initMaterial(); // check that delta_f has a proper value or assign a defaul value if (delta_f < 0) delta_f = this->delta_c_eff; else if (delta_f < this->delta_c_eff) AKANTU_DEBUG_ERROR("Delta_f must be greater or equal to delta_c"); delta_prec.initialize(1); K_plus.initialize(1); K_minus.initialize(1); T_1d.initialize(1); normal_regime.initialize(1); if (count_switches) { switches.initialize(1); delta_dot_prec.initialize(1); } } /* -------------------------------------------------------------------------- */ template void MaterialCohesiveLinearFatigue ::computeTraction(const Array & normal, ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); /// 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); 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 = normal.begin(spatial_dimension); Array::vector_iterator traction_end = this->tractions(el_type, ghost_type).end(spatial_dimension); const Array & sigma_c_array = this->sigma_c_eff(el_type, ghost_type); Array & delta_max_array = this->delta_max(el_type, ghost_type); const Array & delta_c_array = this->delta_c_eff(el_type, ghost_type); Array & damage_array = this->damage(el_type, ghost_type); Array::vector_iterator insertion_stress_it = this->insertion_stress(el_type, ghost_type).begin(spatial_dimension); Array & delta_prec_array = delta_prec(el_type, ghost_type); Array & K_plus_array = K_plus(el_type, ghost_type); Array & K_minus_array = K_minus(el_type, ghost_type); Array & T_1d_array = T_1d(el_type, ghost_type); Array & normal_regime_array = normal_regime(el_type, ghost_type); - Array * switches_array; - Array * delta_dot_prec_array; + Array * switches_array = NULL; + Array * delta_dot_prec_array = NULL; if (count_switches) { switches_array = &switches(el_type, ghost_type); delta_dot_prec_array = &delta_dot_prec(el_type, ghost_type); } Real * memory_space = new Real[2*spatial_dimension]; Vector normal_opening(memory_space, spatial_dimension); Vector tangential_opening(memory_space + spatial_dimension, spatial_dimension); Real tolerance = Math::getTolerance(); /// loop on each quadrature point for (UInt q = 0; traction_it != traction_end; ++traction_it, ++opening_it, ++normal_it, ++contact_traction_it, ++insertion_stress_it, ++contact_opening_it, ++q) { /// compute normal and tangential opening vectors Real normal_opening_norm = opening_it->dot(*normal_it); normal_opening = (*normal_it); normal_opening *= normal_opening_norm; tangential_opening = *opening_it; tangential_opening -= normal_opening; Real 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; bool penetration = normal_opening_norm < -tolerance; if (this->contact_after_breaking == false && Math::are_float_equal(damage_array(q), 1.)) penetration = false; if (penetration) { /// use penalty coefficient in case of penetration *contact_traction_it = normal_opening; *contact_traction_it *= this->penalty; *contact_opening_it = normal_opening; /// don't consider penetration contribution for delta *opening_it = tangential_opening; normal_opening.clear(); } else { delta += normal_opening_norm * normal_opening_norm; contact_traction_it->clear(); contact_opening_it->clear(); } delta = std::sqrt(delta); /** * 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$ */ // update maximum displacement and damage delta_max_array(q) = std::max(delta, delta_max_array(q)); damage_array(q) = std::min(delta_max_array(q) / delta_c_array(q), Real(1.)); Real delta_dot = delta - delta_prec_array(q); // count switches if asked if (count_switches) { if ( (delta_dot > 0. && (*delta_dot_prec_array)(q) <= 0.) || (delta_dot < 0. && (*delta_dot_prec_array)(q) >= 0.) ) ++((*switches_array)(q)); (*delta_dot_prec_array)(q) = delta_dot; } // set delta_f equal to delta_max if desired if (progressive_delta_f) delta_f = delta_max_array(q); // broken element case if (Math::are_float_equal(damage_array(q), 1.)) traction_it->clear(); // just inserted element case else if (Math::are_float_equal(damage_array(q), 0.)) { if (penetration) traction_it->clear(); else *traction_it = *insertion_stress_it; // initialize the 1d traction to sigma_c T_1d_array(q) = sigma_c_array(q); } // normal case else { // if element is closed then there are zero tractions if (delta <= tolerance) traction_it->clear(); // otherwise compute new tractions if the new delta is different // than the previous one else if (std::abs(delta_dot) > tolerance) { // loading case if (delta_dot > 0.) { if (!normal_regime_array(q)) { // equation (4) of the article K_plus_array(q) *= 1. - delta_dot / delta_f; // equivalent to equation (2) of the article T_1d_array(q) += K_plus_array(q) * delta_dot; // in case of reloading, traction must not exceed that of the // envelop of the cohesive law Real max_traction = sigma_c_array(q) * (1 - delta / delta_c_array(q)); bool max_traction_exceeded = T_1d_array(q) > max_traction; if (max_traction_exceeded) T_1d_array(q) = max_traction; // switch to standard linear cohesive law if (delta_max_array(q) > fatigue_ratio * delta_c_array(q)) { // reset delta_max to avoid big jumps in the traction delta_max_array(q) = sigma_c_array(q) / (T_1d_array(q) / delta + sigma_c_array(q) / delta_c_array(q)); damage_array(q) = std::min(delta_max_array(q) / delta_c_array(q), Real(1.)); K_minus_array(q) = sigma_c_array(q) / delta_max_array(q) * (1. - damage_array(q)); normal_regime_array(q) = true; } else { // equation (3) of the article K_minus_array(q) = T_1d_array(q) / delta; // if the traction is following the cohesive envelop, then // K_plus has to be reset if (max_traction_exceeded) K_plus_array(q) = K_minus_array(q); } } else { // compute stiffness according to the standard law K_minus_array(q) = sigma_c_array(q) / delta_max_array(q) * (1. - damage_array(q)); } } // unloading case else if (!normal_regime_array(q)) { // equation (4) of the article K_plus_array(q) += (K_plus_array(q) - K_minus_array(q)) * delta_dot / delta_f; // equivalent to equation (2) of the article T_1d_array(q) = K_minus_array(q) * delta; } // applying the actual stiffness *traction_it = tangential_opening; *traction_it *= this->beta2_kappa; *traction_it += normal_opening; *traction_it *= K_minus_array(q); } } // update precendent delta delta_prec_array(q) = delta; } delete [] memory_space; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ INSTANTIATE_MATERIAL(MaterialCohesiveLinearFatigue); __END_AKANTU__ 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 new file mode 100644 index 000000000..55ca17900 --- /dev/null +++ b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_friction.cc @@ -0,0 +1,382 @@ +/** + * @file material_cohesive_linear_friction.cc + * + * @author Mauro Corrado + * + * @date creation: Tue May 08 2012 + * @date last modification: Thu Nov 19 2015 + * + * @brief Linear irreversible cohesive law of mixed mode loading with + * random stress definition for extrinsic type + * + * @section LICENSE + * + * Copyright (©) 2014 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(const Array & normal, + ElementType el_type, + GhostType ghost_type) { + AKANTU_DEBUG_IN(); + + residual_sliding.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); + + Array::vector_iterator opening_prev_it = + this->opening.previous(el_type, ghost_type).begin(spatial_dimension); + + 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_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)); + + /// loop on each quadrature point + for (; traction_it != traction_end; + ++traction_it, ++opening_it, ++opening_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) { + /* ----------------------------------------- */ + Real damage = std::min(*delta_max_prev_it / *delta_c_it, Real(1.)); + Real mu = mu_max * std::sqrt(damage); + /// the definition of tau_max refers to the opening + /// (penetration) of the previous incremental step + Real normal_opening_prev_norm = opening_prev_it->dot(*normal_it); + Vector normal_opening_prev = (*normal_it); + normal_opening_prev *= normal_opening_prev_norm; + Real tau_max = mu * current_penalty * (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_it) < -Math::getTolerance()) + tau = -tau; + // from tau get the x and y components of friction, to be added in the force vector + Vector tangent(spatial_dimension); + tangent = tangential_opening / tangential_opening_norm; + *friction_force_it = tau * tangent; + + //update residual_sliding + *res_sliding_it = tangential_opening_norm - (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, + 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_prev_it = + this->opening.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_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.)); + // the friction coefficient mu is increasing with the + // damage. It equals the maximum value when damage = 1. + Real mu = mu_max * damage; + Real normal_opening_prev_norm = opening_prev_it->dot(*normal_it); + Vector normal_opening_prev = (*normal_it); + normal_opening_prev *= normal_opening_prev_norm; + Real tau_max = mu * current_penalty * 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_friction.hh b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_friction.hh new file mode 100644 index 000000000..d01b97737 --- /dev/null +++ b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_friction.hh @@ -0,0 +1,110 @@ +/** + * @file material_cohesive_linear_friction.hh + * + * @author Mauro Corrado + * + * @date creation: Tue May 08 2012 + * @date last modification: Tue Jul 29 2014 + * + * @brief Linear irreversible cohesive law of mixed mode loading with + * random stress definition for extrinsic type + * + * @section LICENSE + * + * Copyright (©) 2014 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_FRICTION_HH__ +#define __AKANTU_MATERIAL_COHESIVE_LINEAR_FRICTION_HH__ + +/* -------------------------------------------------------------------------- */ + +__BEGIN_AKANTU__ + +/** + * Cohesive material linear with friction force + * + * parameters in the material files : + * - mu : friction coefficient + * - penalty_for_friction : Penalty parameter for the friction behavior + */ +template +class MaterialCohesiveLinearFriction : public MaterialCohesiveLinear { + /* ------------------------------------------------------------------------ */ + /* Constructors/Destructors */ + /* ------------------------------------------------------------------------ */ + typedef MaterialCohesiveLinear MaterialParent; +public: + + MaterialCohesiveLinearFriction(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: + /// maximum value of the friction coefficient + Real mu_max; + + /// penalty parameter for the friction law + Real friction_penalty; + + /// history parameter for the friction law + CohesiveInternalField residual_sliding; + + /// friction force + CohesiveInternalField friction_force; +}; + +__END_AKANTU__ + +#endif /* __AKANTU_MATERIAL_COHESIVE_LINEAR_FRICTION_HH__ */ diff --git a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_inline_impl.cc b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_inline_impl.cc index 207729906..2aba0b544 100644 --- a/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_inline_impl.cc +++ b/src/model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_inline_impl.cc @@ -1,58 +1,313 @@ /** * @file material_cohesive_linear_inline_impl.cc * @author Marco Vocialta * @date Tue Apr 21 14:49:56 2015 * * @brief Inline functions of the MaterialCohesiveLinear * * @section LICENSE * * Copyright (©) 2010-2011 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_MATERIAL_COHESIVE_LINEAR_INLINE_IMPL_CC__ +#define __AKANTU_MATERIAL_COHESIVE_LINEAR_INLINE_IMPL_CC__ +/* -------------------------------------------------------------------------- */ + +__BEGIN_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, + Real & delta_max, + const Real & delta_max_prev, + const Real & delta_c, + const Vector & insertion_stress, + const Real & sigma_c, + Vector & opening, + Vector & opening_prec, + const Vector & normal, + Vector & normal_opening, + Vector & tangential_opening, + Real & normal_opening_norm, + Real & tangential_opening_norm, + Real & damage, + bool & penetration, + bool & reduction_penalty, + Real & current_penalty, + 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 < -Math::getTolerance(); + 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. + */ + + // opening_prec is the opening of the previous convergence loop + // (NB: it is different from opening_prev, that refers to the solution + // of the previous incremental step) + Real normal_opening_prec_norm = opening_prec.dot(normal); + opening_prec = opening; + + if (!this->model->isExplicit() && !this->recompute) + if ((normal_opening_prec_norm * normal_opening_norm) < 0) { + reduction_penalty = true; + } + + if (penetration) { + if (this->recompute && reduction_penalty){ + /// the penalty parameter is locally reduced + current_penalty = this->penalty / 1000.; + } + else + current_penalty = this->penalty; + + /// use penalty coefficient in case of penetration + contact_traction = normal_opening; + contact_traction *= current_penalty; + contact_opening = normal_opening; + + /// don't consider penetration contribution for delta + 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 *= 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, + bool & reduction_penalty, + Real & current_penalty, + 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(); + penetration = normal_opening_norm < -Math::getTolerance(); + if (contact_after_breaking == false && Math::are_float_equal(damage, 1.)) + penetration = false; + + Real derivative = 0; // derivative = d(t/delta)/ddelta + Real t = 0; + + Real delta = tangential_opening_norm * tangential_opening_norm * this->beta2_kappa2; + + Matrix n_outer_n(spatial_dimension, spatial_dimension); + n_outer_n.outerProduct(normal, normal); + + if (penetration){ + if (recompute && reduction_penalty) + current_penalty = this->penalty / 1000.; + else + current_penalty = this->penalty; + + /// stiffness in compression given by the penalty parameter + tangent += n_outer_n; + tangent *= current_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; +} + +/* -------------------------------------------------------------------------- */ +__END_AKANTU__ + +/* -------------------------------------------------------------------------- */ +#endif //__AKANTU_MATERIAL_COHESIVE_LINEAR_INLINE_IMPL_CC__ diff --git a/src/model/solid_mechanics/materials/material_cohesive/material_cohesive.cc b/src/model/solid_mechanics/materials/material_cohesive/material_cohesive.cc index b8193bb85..140d1f284 100644 --- a/src/model/solid_mechanics/materials/material_cohesive/material_cohesive.cc +++ b/src/model/solid_mechanics/materials/material_cohesive/material_cohesive.cc @@ -1,622 +1,624 @@ /** * @file material_cohesive.cc * * @author Seyedeh Mohadeseh Taheri Mousavi * @author Marco Vocialta * @author Nicolas Richart * * @date creation: Wed Feb 22 2012 * @date last modification: Tue Jul 29 2014 * * @brief Specialization of the material class for cohesive elements * * @section LICENSE * * Copyright (©) 2014 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" #include "solid_mechanics_model_cohesive.hh" #include "sparse_matrix.hh" #include "dof_synchronizer.hh" #include "aka_random_generator.hh" #include "shape_cohesive.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ MaterialCohesive::MaterialCohesive(SolidMechanicsModel & model, const ID & id) : Material(model, id), facet_filter("facet_filter", id, this->getMemoryID()), fem_cohesive(&(model.getFEEngineClass("CohesiveFEEngine"))), reversible_energy("reversible_energy", *this), total_energy("total_energy", *this), opening("opening", *this), opening_old("opening (old)", *this), tractions("tractions", *this), tractions_old("tractions (old)", *this), contact_tractions("contact_tractions", *this), contact_opening("contact_opening", *this), delta_max("delta max", *this), use_previous_delta_max(false), + use_previous_opening(false), damage("damage", *this), sigma_c("sigma_c", *this), normal(0, spatial_dimension, "normal") { AKANTU_DEBUG_IN(); this->model = dynamic_cast(&model); this->registerParam("sigma_c", sigma_c, _pat_parsable | _pat_readable, "Critical stress"); this->registerParam("delta_c", delta_c, Real(0.), _pat_parsable | _pat_readable, "Critical displacement"); this->model->getMesh().initElementTypeMapArray(this->element_filter, 1, spatial_dimension, false, _ek_cohesive); if (this->model->getIsExtrinsic()) this->model->getMeshFacets().initElementTypeMapArray(facet_filter, 1, spatial_dimension - 1); this->reversible_energy.initialize(1 ); this->total_energy .initialize(1 ); this->tractions_old .initialize(spatial_dimension); this->tractions .initialize(spatial_dimension); this->opening_old .initialize(spatial_dimension); this->contact_tractions.initialize(spatial_dimension); this->contact_opening .initialize(spatial_dimension); this->opening .initialize(spatial_dimension); this->delta_max .initialize(1 ); this->damage .initialize(1 ); if (this->model->getIsExtrinsic()) this->sigma_c.initialize(1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ MaterialCohesive::~MaterialCohesive() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::initMaterial() { AKANTU_DEBUG_IN(); Material::initMaterial(); if (this->use_previous_delta_max) this->delta_max.initializeHistory(); + if (this->use_previous_opening) this->opening.initializeHistory(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::assembleResidual(GhostType ghost_type) { AKANTU_DEBUG_IN(); #if defined(AKANTU_DEBUG_TOOLS) debug::element_manager.printData(debug::_dm_material_cohesive, "Cohesive Tractions", tractions); #endif Array & residual = const_cast &>(model->getResidual()); Mesh & mesh = 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 = element_filter(*it, ghost_type); UInt nb_element = elem_filter.getSize(); if (nb_element == 0) continue; const Array & shapes = fem_cohesive->getShapes(*it, ghost_type); Array & traction = tractions(*it, ghost_type); Array & contact_traction = contact_tractions(*it, ghost_type); UInt size_of_shapes = shapes.getNbComponent(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(*it); UInt nb_quadrature_points = fem_cohesive->getNbIntegrationPoints(*it, ghost_type); /// compute @f$t_i N_a@f$ Array * traction_cpy = new Array(nb_element * nb_quadrature_points, spatial_dimension * size_of_shapes); Array::iterator > traction_it = traction.begin(spatial_dimension, 1); Array::iterator > contact_traction_it = contact_traction.begin(spatial_dimension, 1); Array::const_iterator > shapes_filtered_begin = shapes.begin(1, size_of_shapes); Array::iterator > traction_cpy_it = traction_cpy->begin(spatial_dimension, size_of_shapes); Matrix traction_tmp(spatial_dimension, 1); for (UInt el = 0; el < nb_element; ++el) { UInt current_quad = elem_filter(el) * nb_quadrature_points; for (UInt q = 0; q < nb_quadrature_points; ++q, ++traction_it, ++contact_traction_it, ++current_quad, ++traction_cpy_it) { const Matrix & shapes_filtered = shapes_filtered_begin[current_quad]; traction_tmp.copy(*traction_it); traction_tmp += *contact_traction_it; traction_cpy_it->mul(traction_tmp, shapes_filtered); } } /** * compute @f$\int t \cdot N\, dS@f$ by @f$ \sum_q \mathbf{N}^t * \mathbf{t}_q \overline w_q J_q@f$ */ Array * int_t_N = new Array(nb_element, spatial_dimension*size_of_shapes, "int_t_N"); fem_cohesive->integrate(*traction_cpy, *int_t_N, spatial_dimension * size_of_shapes, *it, ghost_type, elem_filter); delete traction_cpy; int_t_N->extendComponentsInterlaced(2, int_t_N->getNbComponent()); Real * int_t_N_val = int_t_N->storage(); for (UInt el = 0; el < nb_element; ++el) { for (UInt n = 0; n < size_of_shapes * spatial_dimension; ++n) int_t_N_val[n] *= -1.; int_t_N_val += nb_nodes_per_element * spatial_dimension; } /// assemble model->getFEEngineBoundary().assembleArray(*int_t_N, residual, model->getDOFSynchronizer().getLocalDOFEquationNumbers(), residual.getNbComponent(), *it, ghost_type, elem_filter, 1); delete int_t_N; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::assembleStiffnessMatrix(GhostType ghost_type) { AKANTU_DEBUG_IN(); SparseMatrix & K = const_cast(model->getStiffnessMatrix()); Mesh & mesh = 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) { UInt nb_quadrature_points = fem_cohesive->getNbIntegrationPoints(*it, ghost_type); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(*it); const Array & shapes = fem_cohesive->getShapes(*it, ghost_type); Array & elem_filter = element_filter(*it, ghost_type); UInt nb_element = elem_filter.getSize(); if(!nb_element) continue; UInt size_of_shapes = shapes.getNbComponent(); Array * shapes_filtered = new Array(nb_element*nb_quadrature_points, size_of_shapes, "filtered shapes"); Real * shapes_val = shapes.storage(); Real * shapes_filtered_val = shapes_filtered->storage(); UInt * elem_filter_val = elem_filter.storage(); for (UInt el = 0; el < nb_element; ++el) { shapes_val = shapes.storage() + elem_filter_val[el] * size_of_shapes * nb_quadrature_points; memcpy(shapes_filtered_val, shapes_val, size_of_shapes * nb_quadrature_points * sizeof(Real)); shapes_filtered_val += size_of_shapes * nb_quadrature_points; } /** * compute A matrix @f$ \mathbf{A} = \left[\begin{array}{c c c c c c c c c c c c} * 1 & 0 & 0 & 0 & 0 & 0 & -1 & 0 & 0 & 0 & 0 & 0 \\ * 0 & 1 & 0 & 0 & 0 & 0 & 0 & -1 & 0 & 0 & 0 & 0 \\ * 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & -1 & 0 & 0 & 0 \\ * 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & -1 & 0 & 0 \\ * 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & -1 & 0 \\ * 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & -1 * \end{array} \right]@f$ **/ // UInt size_of_A = spatial_dimension*size_of_shapes*spatial_dimension*nb_nodes_per_element; // Real * A = new Real[size_of_A]; // memset(A, 0, size_of_A*sizeof(Real)); Matrix A(spatial_dimension*size_of_shapes, spatial_dimension*nb_nodes_per_element); for ( UInt i = 0; i < spatial_dimension*size_of_shapes; ++i) { A(i, i) = 1; A(i, i + spatial_dimension*size_of_shapes) = -1; } /// compute traction. This call is not necessary for the linear /// cohesive law that, currently, is the only one used for the /// extrinsic approach. if (!model->getIsExtrinsic()){ computeTraction(ghost_type); } /// get the tangent matrix @f$\frac{\partial{(t/\delta)}}{\partial{\delta}} @f$ Array * tangent_stiffness_matrix = new Array(nb_element * nb_quadrature_points, spatial_dimension * spatial_dimension, "tangent_stiffness_matrix"); // Array * normal = new Array(nb_element * nb_quadrature_points, spatial_dimension, "normal"); normal.resize(nb_quadrature_points); computeNormal(model->getCurrentPosition(), normal, *it, ghost_type); /// compute openings @f$\mathbf{\delta}@f$ //computeOpening(model->getDisplacement(), opening(*it, ghost_type), *it, ghost_type); tangent_stiffness_matrix->clear(); computeTangentTraction(*it, *tangent_stiffness_matrix, normal, ghost_type); // delete normal; UInt size_at_nt_d_n_a = spatial_dimension*nb_nodes_per_element*spatial_dimension*nb_nodes_per_element; Array * at_nt_d_n_a = new Array (nb_element*nb_quadrature_points, size_at_nt_d_n_a, "A^t*N^t*D*N*A"); Array::iterator > shapes_filt_it = shapes_filtered->begin(size_of_shapes); Array::matrix_iterator D_it = tangent_stiffness_matrix->begin(spatial_dimension, spatial_dimension); Array::matrix_iterator At_Nt_D_N_A_it = at_nt_d_n_a->begin(spatial_dimension * nb_nodes_per_element, spatial_dimension * nb_nodes_per_element); Array::matrix_iterator At_Nt_D_N_A_end = at_nt_d_n_a->end(spatial_dimension * nb_nodes_per_element, spatial_dimension * nb_nodes_per_element); Matrix N (spatial_dimension, spatial_dimension * size_of_shapes); Matrix N_A (spatial_dimension, spatial_dimension * nb_nodes_per_element); Matrix D_N_A(spatial_dimension, spatial_dimension * nb_nodes_per_element); for(; At_Nt_D_N_A_it != At_Nt_D_N_A_end; ++At_Nt_D_N_A_it, ++D_it, ++shapes_filt_it) { N.clear(); /** * store the shapes in voigt notations matrix @f$\mathbf{N} = * \begin{array}{cccccc} N_0(\xi) & 0 & N_1(\xi) &0 & N_2(\xi) & 0 \\ * 0 & * N_0(\xi)& 0 &N_1(\xi)& 0 & N_2(\xi) \end{array} @f$ **/ for (UInt i = 0; i < spatial_dimension ; ++i) for (UInt n = 0; n < size_of_shapes; ++n) N(i, i + spatial_dimension * n) = (*shapes_filt_it)(n); /** * compute stiffness matrix @f$ \mathbf{K} = \delta \mathbf{U}^T * \int_{\Gamma_c} {\mathbf{P}^t \frac{\partial{\mathbf{t}}} {\partial{\delta}} * \mathbf{P} d\Gamma \Delta \mathbf{U}} @f$ **/ N_A.mul(N, A); D_N_A.mul(*D_it, N_A); (*At_Nt_D_N_A_it).mul(D_N_A, N_A); } delete tangent_stiffness_matrix; delete shapes_filtered; Array * K_e = new Array(nb_element, size_at_nt_d_n_a, "K_e"); fem_cohesive->integrate(*at_nt_d_n_a, *K_e, size_at_nt_d_n_a, *it, ghost_type, elem_filter); delete at_nt_d_n_a; model->getFEEngine().assembleMatrix(*K_e, K, spatial_dimension, *it, ghost_type, elem_filter); delete K_e; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- * * Compute traction from displacements * * @param[in] ghost_type compute the residual for _ghost or _not_ghost element */ void MaterialCohesive::computeTraction(GhostType ghost_type) { AKANTU_DEBUG_IN(); #if defined(AKANTU_DEBUG_TOOLS) debug::element_manager.printData(debug::_dm_material_cohesive, "Cohesive Openings", opening); #endif Mesh & mesh = 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 = element_filter(*it, ghost_type); UInt nb_element = elem_filter.getSize(); if (nb_element == 0) continue; UInt nb_quadrature_points = nb_element*fem_cohesive->getNbIntegrationPoints(*it, ghost_type); normal.resize(nb_quadrature_points); /// compute normals @f$\mathbf{n}@f$ computeNormal(model->getCurrentPosition(), normal, *it, ghost_type); /// compute openings @f$\mathbf{\delta}@f$ computeOpening(model->getDisplacement(), opening(*it, ghost_type), *it, ghost_type); /// compute traction @f$\mathbf{t}@f$ computeTraction(normal, *it, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::computeNormal(const Array & position, Array & normal, ElementType type, GhostType ghost_type) { AKANTU_DEBUG_IN(); if (type == _cohesive_1d_2) fem_cohesive->computeNormalsOnIntegrationPoints(position, normal, type, ghost_type); else { #define COMPUTE_NORMAL(type) \ fem_cohesive->getShapeFunctions(). \ computeNormalsOnIntegrationPoints(position, \ normal, \ ghost_type, \ element_filter(type, ghost_type)); AKANTU_BOOST_COHESIVE_ELEMENT_SWITCH(COMPUTE_NORMAL); #undef COMPUTE_NORMAL } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::computeOpening(const Array & displacement, Array & opening, ElementType type, GhostType ghost_type) { AKANTU_DEBUG_IN(); #define COMPUTE_OPENING(type) \ fem_cohesive->getShapeFunctions(). \ interpolateOnIntegrationPoints(displacement, \ opening, \ spatial_dimension, \ ghost_type, \ element_filter(type, ghost_type)); AKANTU_BOOST_COHESIVE_ELEMENT_SWITCH(COMPUTE_OPENING); #undef COMPUTE_OPENING AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::computeEnergies() { AKANTU_DEBUG_IN(); Mesh & mesh = fem_cohesive->getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, _not_ghost, _ek_cohesive); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, _not_ghost, _ek_cohesive); Real * memory_space = new Real[2*spatial_dimension]; Vector b(memory_space, spatial_dimension); Vector h(memory_space + spatial_dimension, spatial_dimension); for(; it != last_type; ++it) { Array::iterator erev = reversible_energy(*it, _not_ghost).begin(); Array::iterator etot = total_energy(*it, _not_ghost).begin(); Array::vector_iterator traction_it = tractions(*it, _not_ghost).begin(spatial_dimension); Array::vector_iterator traction_old_it = tractions_old(*it, _not_ghost).begin(spatial_dimension); Array::vector_iterator opening_it = opening(*it, _not_ghost).begin(spatial_dimension); Array::vector_iterator opening_old_it = opening_old(*it, _not_ghost).begin(spatial_dimension); Array::vector_iterator traction_end = tractions(*it, _not_ghost).end(spatial_dimension); /// loop on each quadrature point for (; traction_it != traction_end; ++traction_it, ++traction_old_it, ++opening_it, ++opening_old_it, ++erev, ++etot) { /// trapezoidal integration b = *opening_it; b -= *opening_old_it; h = *traction_old_it; h += *traction_it; *etot += .5 * b.dot(h); *erev = .5 * traction_it->dot(*opening_it); } } delete [] memory_space; /// update old values it = mesh.firstType(spatial_dimension, _not_ghost, _ek_cohesive); GhostType ghost_type = _not_ghost; for(; it != last_type; ++it) { tractions_old(*it, ghost_type).copy(tractions(*it, ghost_type)); opening_old(*it, ghost_type).copy(opening(*it, ghost_type)); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real MaterialCohesive::getReversibleEnergy() { AKANTU_DEBUG_IN(); Real erev = 0.; /// integrate reversible energy for each type of elements Mesh & mesh = fem_cohesive->getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, _not_ghost, _ek_cohesive); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, _not_ghost, _ek_cohesive); for(; it != last_type; ++it) { erev += fem_cohesive->integrate(reversible_energy(*it, _not_ghost), *it, _not_ghost, element_filter(*it, _not_ghost)); } AKANTU_DEBUG_OUT(); return erev; } /* -------------------------------------------------------------------------- */ Real MaterialCohesive::getDissipatedEnergy() { AKANTU_DEBUG_IN(); Real edis = 0.; /// integrate dissipated energy for each type of elements Mesh & mesh = fem_cohesive->getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, _not_ghost, _ek_cohesive); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, _not_ghost, _ek_cohesive); for(; it != last_type; ++it) { Array dissipated_energy(total_energy(*it, _not_ghost)); dissipated_energy -= reversible_energy(*it, _not_ghost); edis += fem_cohesive->integrate(dissipated_energy, *it, _not_ghost, element_filter(*it, _not_ghost)); } AKANTU_DEBUG_OUT(); return edis; } /* -------------------------------------------------------------------------- */ Real MaterialCohesive::getContactEnergy() { AKANTU_DEBUG_IN(); Real econ = 0.; /// integrate contact energy for each type of elements Mesh & mesh = fem_cohesive->getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, _not_ghost, _ek_cohesive); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, _not_ghost, _ek_cohesive); for(; it != last_type; ++it) { Array & el_filter = element_filter(*it, _not_ghost); UInt nb_quad_per_el = fem_cohesive->getNbIntegrationPoints(*it, _not_ghost); UInt nb_quad_points = el_filter.getSize() * nb_quad_per_el; Array contact_energy(nb_quad_points); Array::vector_iterator contact_traction_it = contact_tractions(*it, _not_ghost).begin(spatial_dimension); Array::vector_iterator contact_opening_it = contact_opening(*it, _not_ghost).begin(spatial_dimension); /// loop on each quadrature point for (UInt el = 0; el < nb_quad_points; ++contact_traction_it, ++contact_opening_it, ++el) { contact_energy(el) = .5 * contact_traction_it->dot(*contact_opening_it); } econ += fem_cohesive->integrate(contact_energy, *it, _not_ghost, el_filter); } AKANTU_DEBUG_OUT(); return econ; } /* -------------------------------------------------------------------------- */ Real MaterialCohesive::getEnergy(std::string type) { AKANTU_DEBUG_IN(); if (type == "reversible") return getReversibleEnergy(); else if (type == "dissipated") return getDissipatedEnergy(); else if (type == "cohesive contact") return getContactEnergy(); AKANTU_DEBUG_OUT(); return 0.; } /* -------------------------------------------------------------------------- */ __END_AKANTU__ diff --git a/src/model/solid_mechanics/materials/material_cohesive/material_cohesive.hh b/src/model/solid_mechanics/materials/material_cohesive/material_cohesive.hh index 61c489626..0eade6996 100644 --- a/src/model/solid_mechanics/materials/material_cohesive/material_cohesive.hh +++ b/src/model/solid_mechanics/materials/material_cohesive/material_cohesive.hh @@ -1,264 +1,267 @@ /** * @file material_cohesive.hh * * @author Seyedeh Mohadeseh Taheri Mousavi * @author Marco Vocialta * @author Nicolas Richart * * @date creation: Wed Feb 22 2012 * @date last modification: Tue Jul 29 2014 * * @brief Specialization of the material class for cohesive elements * * @section LICENSE * * Copyright (©) 2014 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 "material.hh" #include "fe_engine_template.hh" #include "aka_common.hh" #include "cohesive_internal_field.hh" #include "cohesive_element_inserter.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_COHESIVE_HH__ #define __AKANTU_MATERIAL_COHESIVE_HH__ /* -------------------------------------------------------------------------- */ namespace akantu { class SolidMechanicsModelCohesive; } __BEGIN_AKANTU__ class MaterialCohesive : public Material { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: typedef FEEngineTemplate MyFEEngineCohesiveType; public: MaterialCohesive(SolidMechanicsModel& model, const ID & id = ""); virtual ~MaterialCohesive(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// initialize the material computed parameter virtual void initMaterial(); /// compute tractions (including normals and openings) void computeTraction(GhostType ghost_type = _not_ghost); /// assemble residual void assembleResidual(GhostType ghost_type = _not_ghost); /// compute reversible and total energies by element void computeEnergies(); /// check stress for cohesive elements' insertion, by default it /// also updates the cohesive elements' data virtual void checkInsertion(bool check_only = false) { AKANTU_DEBUG_TO_IMPLEMENT(); } /// 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) { AKANTU_DEBUG_TO_IMPLEMENT(); } /// reset variables when convergence is reached (only for /// extrinsic-implicit) virtual void resetVariables(GhostType ghost_type = _not_ghost) { AKANTU_DEBUG_TO_IMPLEMENT(); } /// interpolate stress on given positions for each element (empty /// implemantation to avoid the generic call to be done on cohesive elements) virtual void interpolateStress(__attribute__((unused)) const ElementType type, __attribute__((unused)) Array & result) { }; /// compute the stresses virtual void computeAllStresses(__attribute__((unused)) GhostType ghost_type = _not_ghost) { }; // add the facet to be handled by the material UInt addFacet(const Element & element); protected: virtual void computeTangentTraction(__attribute__((unused)) const ElementType & el_type, __attribute__((unused)) Array & tangent_matrix, __attribute__((unused)) const Array & normal, __attribute__((unused)) GhostType ghost_type = _not_ghost) { AKANTU_DEBUG_TO_IMPLEMENT(); } /// compute the normal void computeNormal(const Array & position, Array & normal, ElementType type, GhostType ghost_type); /// compute the opening void computeOpening(const Array & displacement, Array & normal, ElementType type, GhostType ghost_type); template void computeNormal(const Array & position, Array & normal, GhostType ghost_type); /// assemble stiffness void assembleStiffnessMatrix(GhostType ghost_type); /// constitutive law virtual void computeTraction(const Array & normal, ElementType el_type, GhostType ghost_type = _not_ghost) = 0; /// parallelism functions inline UInt getNbDataForElements(const Array & elements, SynchronizationTag tag) const; inline void packElementData(CommunicationBuffer & buffer, const Array & elements, SynchronizationTag tag) const; inline void unpackElementData(CommunicationBuffer & buffer, const Array & elements, SynchronizationTag tag); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get the opening AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Opening, opening, Real); /// get the traction AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Traction, tractions, Real); /// get damage AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Damage, damage, Real); /// get facet filter AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(FacetFilter, facet_filter, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(FacetFilter, facet_filter, UInt); AKANTU_GET_MACRO(FacetFilter, facet_filter, const ElementTypeMapArray &); // AKANTU_GET_MACRO(ElementFilter, element_filter, const ElementTypeMapArray &); /// compute reversible energy Real getReversibleEnergy(); /// compute dissipated energy Real getDissipatedEnergy(); /// compute contact energy Real getContactEnergy(); /// get energy virtual Real getEnergy(std::string type); /// return the energy (identified by id) for the provided element virtual Real getEnergy(std::string energy_id, ElementType type, UInt index) { return Material::getEnergy(energy_id, type, index); } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// list of facets assigned to this material ElementTypeMapArray facet_filter; /// Link to the cohesive fem object in the model MyFEEngineCohesiveType * fem_cohesive; private: /// reversible energy by quadrature point CohesiveInternalField reversible_energy; /// total energy by quadrature point CohesiveInternalField total_energy; protected: /// opening in all elements and quadrature points CohesiveInternalField opening; /// opening in all elements and quadrature points (previous time step) CohesiveInternalField opening_old; /// traction in all elements and quadrature points CohesiveInternalField tractions; /// traction in all elements and quadrature points (previous time step) CohesiveInternalField tractions_old; /// traction due to contact CohesiveInternalField contact_tractions; /// normal openings for contact tractions CohesiveInternalField contact_opening; /// maximum displacement CohesiveInternalField delta_max; /// tell if the previous delta_max state is needed (in iterative schemes) bool use_previous_delta_max; + /// tell if the previous opening state is needed (in iterative schemes) + bool use_previous_opening; + /// damage CohesiveInternalField damage; /// pointer to the solid mechanics model for cohesive elements SolidMechanicsModelCohesive * model; /// critical stress RandomInternalField sigma_c; /// critical displacement Real delta_c; /// array to temporarily store the normals Array normal; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material_cohesive_inline_impl.cc" __END_AKANTU__ #include "cohesive_internal_field_tmpl.hh" #endif /* __AKANTU_MATERIAL_COHESIVE_HH__ */ diff --git a/src/model/solid_mechanics/materials/material_cohesive_includes.hh b/src/model/solid_mechanics/materials/material_cohesive_includes.hh index 62286a70e..70ede7cf9 100644 --- a/src/model/solid_mechanics/materials/material_cohesive_includes.hh +++ b/src/model/solid_mechanics/materials/material_cohesive_includes.hh @@ -1,44 +1,46 @@ /** * @file material_cohesive_includes.hh * * @author Nicolas Richart * * @date creation: Wed Oct 31 2012 * @date last modification: Fri Mar 21 2014 * * @brief List of includes for cohesive elements * * @section LICENSE * * Copyright (©) 2014 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_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_bilinear , MaterialCohesiveBilinear ))) \ ((2, (cohesive_exponential , MaterialCohesiveExponential )))