diff --git a/cohesive_extrinsic.cc b/cohesive_extrinsic.cc index 9acc5ab..e244f85 100644 --- a/cohesive_extrinsic.cc +++ b/cohesive_extrinsic.cc @@ -1,137 +1,144 @@ /** * @file cohesive_extrinsic.cc * * @author Zineb Fouad * @author Nicolas Richart * @author Seyedeh Mohadeseh Taheri Mousavi * @author Marco Vocialta * * @date creation: Tue May 08 2012 * @date last modification: Wed Feb 06 2019 * * @brief Cohesive element examples in extrinsic * * * @section LICENSE * * Copyright (©) 2015-2021 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "solid_mechanics_model_cohesive.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ using namespace akantu; -int main(int argc, char * argv[]) { +int main(int argc, char *argv[]) { initialize("material.dat", argc, argv); const UInt spatial_dimension = 3; const UInt max_steps = 10000; Mesh mesh(spatial_dimension); mesh.read("cube.msh"); SolidMechanicsModelCohesive model(mesh); /// model initialization model.initFull(_analysis_method = _explicit_lumped_mass, _is_extrinsic = true); Real time_step = model.getStableTimeStep() * 0.05; model.setTimeStep(time_step); std::cout << "Time step: " << time_step << std::endl; - auto & inserter = model.getElementInserter(); + auto &inserter = model.getElementInserter(); inserter.setLimit(_y, 0.30, 0.20); model.updateAutomaticInsertion(); - const auto & position = mesh.getNodes(); - auto & velocity = model.getVelocity(); + const auto &position = mesh.getNodes(); + auto &velocity = model.getVelocity(); /// boundary conditions - model.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "top"); - model.applyBC(BC::Dirichlet::FixedValue(0.0, _z), "top"); - model.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "bottom"); - model.applyBC(BC::Dirichlet::FixedValue(0.0, _z), "bottom"); + // model.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "top"); + // model.applyBC(BC::Dirichlet::FixedValue(0.0, _z), "top"); + // model.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "bottom"); + // model.applyBC(BC::Dirichlet::FixedValue(0.0, _z), "bottom"); + + Matrix stress{{c, 0., 0.}, {0., c, 0.}, {0., 0., c}}; + model.applyBC(BC::Neumann::FromStress(stress, "top")); + model.applyBC(BC::Neumann::FromStress(stress, "bottom")); model.setBaseName("extrinsic"); model.addDumpFieldVector("displacement"); model.addDumpField("velocity"); model.addDumpField("acceleration"); model.addDumpField("internal_force"); model.addDumpField("stress"); model.addDumpField("blocked_dofs"); model.addDumpField("grad_u"); model.addDumpFieldToDumper("cohesive elements", "displacement"); model.addDumpFieldToDumper("cohesive elements", "velocity"); model.addDumpFieldToDumper("cohesive elements", "tractions"); auto lower = mesh.getLowerBounds(); auto upper = mesh.getUpperBounds(); /// initial conditions Real loading_rate = 0.5; Real disp_update = loading_rate * time_step; - for (auto && data : zip(make_view(position, spatial_dimension), make_view(velocity, spatial_dimension))) { - auto && position = std::get<0>(data); - auto && velocity = std::get<1>(data); - auto pos = position(_z);// - (upper(_z) - lower(_z)) / 2; + for (auto &&data : zip(make_view(position, spatial_dimension), + make_view(velocity, spatial_dimension))) { + auto &&position = std::get<0>(data); + auto &&velocity = std::get<1>(data); + auto pos = position(_z); // - (upper(_z) - lower(_z)) / 2; velocity(_x) = loading_rate * pos; - velocity(_z) = loading_rate / 1000 * pos; + // velocity(_z) = loading_rate / 1000 * pos; } model.dump(); model.dump("cohesive elements"); /// Main loop for (auto s : arange(1, max_steps)) { - model.applyBC(BC::Dirichlet::IncrementValue( disp_update, _x), "top"); - //model.applyBC(BC::Dirichlet::IncrementValue( disp_update/1000, _z), "top"); - model.applyBC(BC::Dirichlet::IncrementValue(-disp_update, _x), "bottom"); - //model.applyBC(BC::Dirichlet::IncrementValue(-disp_update/1000, _z), "bottom"); + // model.applyBC(BC::Dirichlet::IncrementValue(disp_update, _x), "top"); + // model.applyBC(BC::Dirichlet::IncrementValue( disp_update/1000, _z), + // "top"); + // model.applyBC(BC::Dirichlet::IncrementValue(-disp_update, _x), "bottom"); + // model.applyBC(BC::Dirichlet::IncrementValue(-disp_update/1000, _z), + // "bottom"); model.checkCohesiveStress(); model.solveStep(); if (s % 100 == 0) { model.dump(); model.dump("cohesive elements"); std::cout << "passing step " << s << "/" << max_steps << std::endl; } } Real Ed = model.getEnergy("dissipated"); Real Edt = 200 * std::sqrt(2); std::cout << Ed << " " << Edt << std::endl; if (Ed < Edt * 0.999 || Ed > Edt * 1.001 || std::isnan(Ed)) { std::cout << "The dissipated energy is incorrect" << std::endl; return EXIT_FAILURE; } finalize(); return EXIT_SUCCESS; } diff --git a/material.dat b/material.dat index 9c1edbd..2b99d04 100644 --- a/material.dat +++ b/material.dat @@ -1,17 +1,21 @@ model solid_mechanics_model_cohesive [ material elastic [ name = body - rho = 1e3 # density - E = 1e9 # young's modulus - nu = 0.2 # poisson's ratio + rho = 2500 # density + E = 70e9 # young's modulus + nu = 0.22 # poisson's ratio finite_deformation = true ] material cohesive_bilinear [ name = insertion - beta = 1 - G_c = 10 - sigma_c = 1e6 - delta_0 = 1e-7 + beta = 3 + G_c = 7.6 + sigma_c = 70e6 + delta_0 = 0.217e-6 + ] + + cohesive_inserter [ + bounding_box = [[-1., 1.],[-1., 1.],[-1., 1.]] ] ]