diff --git a/include/constraints/RebarBondSlipConstraint.h b/include/constraints/RebarBondSlipConstraint.h new file mode 100644 index 000000000..6bdab89b2 --- /dev/null +++ b/include/constraints/RebarBondSlipConstraint.h @@ -0,0 +1,90 @@ +//* This file is part of the MOOSE framework +//* https://www.mooseframework.org +//* +//* All rights reserved, see COPYRIGHT for full restrictions +//* https://github.com/idaholab/moose/blob/master/COPYRIGHT +//* +//* Licensed under LGPL 2.1, please see LICENSE for details +//* https://www.gnu.org/licenses/lgpl-2.1.html + +#pragma once + +// MOOSE includes +#include "EqualValueEmbeddedConstraint.h" + +// Forward Declarations +class RebarBondSlipConstraint; + +template <> +InputParameters validParams(); + +/// A RebarBondSlipConstraint enforces concrete-rebar constraint +class RebarBondSlipConstraint : public EqualValueEmbeddedConstraint +{ +public: + static InputParameters validParams(); + + RebarBondSlipConstraint(const InputParameters & parameters); + virtual void initialSetup() override; + virtual void timestepSetup() override; + bool shouldApply() override; + void reinitConstraint(); + +protected: + virtual void computeTangent(); + virtual Real computeQpResidual(Moose::ConstraintType type) override; + virtual Real computeQpJacobian(Moose::ConstraintJacobianType type) override; + virtual Real computeQpOffDiagJacobian(Moose::ConstraintJacobianType type, + unsigned int jvar) override; + + /** + * Struct designed to hold info about the bonds slip history + */ + struct bondSlipData + { + Real slip_min; + Real slip_max; + Real slip_min_old; + Real slip_max_old; + Real bondstress_min; + Real bondstress_max; + Real bondstress_min_old; + Real bondstress_max_old; + + bondSlipData() + : slip_min(0.0), + slip_max(0.0), + slip_min_old(0.0), + slip_max_old(0.0), + bondstress_min(0.0), + bondstress_max(0.0), + bondstress_min_old(0.0), + bondstress_max_old(0.0) + { + } + }; + + // Bond-slip data + std::map _bondslip; + + const unsigned _component; + const unsigned int _mesh_dimension; + std::vector _var_nums; + std::vector _vars; + + const bool _debug; + + const Real _max_bondstress; + const Real _frictional_bondstress; + const Real _ultimate_slip; + const Real _bar_radius; + std::vector _transitional_slip; + /// constraint force needed to enforce the constraint + RealVectorValue _constraint_residual; + /// penalty force for the current constraint + RealVectorValue _pen_force; + RealVectorValue _slave_tangent; + Real _current_elem_volume; + bool _bond; + Real _bond_stress; +}; diff --git a/src/constraints/RebarBondSlipConstraint.C b/src/constraints/RebarBondSlipConstraint.C new file mode 100644 index 000000000..7d94e2d28 --- /dev/null +++ b/src/constraints/RebarBondSlipConstraint.C @@ -0,0 +1,348 @@ +//* This file is part of the MOOSE framework +//* https://www.mooseframework.org +//* +//* All rights reserved, see COPYRIGHT for full restrictions +//* https://github.com/idaholab/moose/blob/master/COPYRIGHT +//* +//* Licensed under LGPL 2.1, please see LICENSE for details +//* https://www.gnu.org/licenses/lgpl-2.1.html + +// MOOSE includes +#include "RebarBondSlipConstraint.h" +#include "Assembly.h" +#include "SystemBase.h" +#include "FEProblem.h" +#include "MathUtils.h" + +#include "libmesh/string_to_enum.h" + +registerMooseObject("BlackBearApp", RebarBondSlipConstraint); + +defineLegacyParams(RebarBondSlipConstraint); + +InputParameters +RebarBondSlipConstraint::validParams() +{ + InputParameters params = EqualValueEmbeddedConstraint::validParams(); + params.addClassDescription( + "This is a constraint enforcing the bodslip behavior between concrete and rebar"); + params.addRequiredParam("component", + "An integer corresponding to the direction " + "the variable this kernel acts in. (0 for x, " + "1 for y, 2 for z)"); + params.addCoupledVar( + "displacements", + "The displacements appropriate for the simulation geometry and coordinate system"); + params.addParam("debug", false, "whether to print out debug messages"); + params.addParam("max_bondstress", 0.0, "Maximum bond stress"); + params.addParam("frictional_bondstress", 0.0, "Bond stress due to friction"); + + params.addParam>( + "transitional_slip_values", + "Singificant slip values at which the bondstress curve changes pattern/slope or " + "trnsitions to a different function"); + params.addParam( + "ultimate_slip", 0.05, "Ultimate value of slip at which the concrete and rebar debonds"); + params.addParam("rebar_radius", 1.0, "Radius of the rebar"); + return params; +} + +RebarBondSlipConstraint::RebarBondSlipConstraint(const InputParameters & parameters) + : EqualValueEmbeddedConstraint(parameters), + _component(getParam("component")), + _mesh_dimension(_mesh.dimension()), + _var_nums(_mesh_dimension, libMesh::invalid_uint), + _vars(_mesh_dimension, nullptr), + _debug(getParam("debug")), + _max_bondstress(getParam("max_bondstress")), + _frictional_bondstress(getParam("frictional_bondstress")), + _ultimate_slip(getParam("ultimate_slip")), + _bar_radius(getParam("rebar_radius")), + _transitional_slip(getParam>("transitional_slip_values")) +{ + if (_mesh_dimension != coupledComponents("displacements")) + mooseError("In RebarBondSlipConstraint, number of displacements must equal the mesh dimension"); + + for (unsigned int i = 0; i < _mesh_dimension; ++i) + { + _var_nums[i] = coupled("displacements", i); + _vars[i] = getVar("displacements", i); + } +} + +void +RebarBondSlipConstraint::initialSetup() +{ + for (auto it = _slave_to_master_map.begin(); it != _slave_to_master_map.end(); ++it) + if (_bondslip.find(it->first) == _bondslip.end()) + _bondslip.insert(std::pair(it->first, + bondSlipData())); // initialize +} + +void +RebarBondSlipConstraint::timestepSetup() +{ + for (auto iter = _slave_to_master_map.begin(); iter != _slave_to_master_map.end(); ++iter) + { + dof_id_type node_id = iter->first; + auto it = _bondslip.find(node_id); + mooseAssert(it != _bondslip.end(), "Node not found in bond-slip map"); + + _bondslip[node_id].slip_min_old = _bondslip[node_id].slip_min; + _bondslip[node_id].slip_max_old = _bondslip[node_id].slip_max; + _bondslip[node_id].bondstress_min_old = _bondslip[node_id].bondstress_min; + _bondslip[node_id].bondstress_max_old = _bondslip[node_id].bondstress_max; + } +} + +bool +RebarBondSlipConstraint::shouldApply() +{ + if (_debug) + // if (_current_node->id() == 144) + { + std::cout << "===========================================\n"; + std::cout << "node id: " << _current_node->id() << std::endl; + std::cout << "at coord: " << (Point)*_current_node << std::endl; + } + auto it = _slave_to_master_map.find(_current_node->id()); + + if (it != _slave_to_master_map.end()) + { + const Elem * master_elem = _mesh.elemPtr(it->second); + std::vector points = {*_current_node}; + + // reinit variables on the master element at the slave point + _fe_problem.setNeighborSubdomainID(master_elem, 0); + _fe_problem.reinitNeighborPhys(master_elem, points, 0); + + reinitConstraint(); + + return true; + } + return false; +} + +void +RebarBondSlipConstraint::computeTangent() +{ + _slave_tangent *= 0.0; + + // get normals + // get connected elements of the current node + const std::map> & node_to_elem_map = _mesh.nodeToElemMap(); + auto node_to_elem_pair = node_to_elem_map.find(_current_node->id()); + mooseAssert(node_to_elem_pair != node_to_elem_map.end(), "Missing entry in node to elem map"); + const std::vector & elems = node_to_elem_pair->second; + + for (auto & elem : elems) + { + Elem * elem_ptr = _mesh.elemPtr(elem); + // _assembly.reinit(elem_ptr, 0); + _current_elem_volume += _assembly.elemVolume(); + // calculate phi and dphi for this element + FEType fe_type(Utility::string_to_enum("first"), + Utility::string_to_enum("lagrange")); + std::unique_ptr fe(FEBase::build(1, fe_type)); + fe->attach_quadrature_rule(const_cast(_assembly.qRule())); + const std::vector * tangents = &fe->get_dxyzdxi(); + unsigned side = 0; + fe->reinit(elem_ptr, side); + for (unsigned i = 0; i < tangents->size(); i++) + _slave_tangent += (*tangents)[i]; + } + + _slave_tangent /= _slave_tangent.norm(); + _current_elem_volume /= elems.size(); + + if (_debug) + std::cout << "tangent: " << _slave_tangent << std::endl; +} + +void +RebarBondSlipConstraint::reinitConstraint() +{ + computeTangent(); + + // Build up residual vector + + RealVectorValue relative_disp; + for (unsigned int i = 0; i < _mesh_dimension; ++i) + relative_disp(i) = ((_vars[i]->dofValues())[0] - (_vars[i]->slnNeighbor())[0]); + + Real slip = relative_disp * _slave_tangent; + RealVectorValue slip_axial = slip * _slave_tangent; + RealVectorValue slip_normal = relative_disp - slip_axial; + Real slip_ratio = std::abs(slip) / _transitional_slip[0]; + // Real bond_stress; + if (_debug) + // if (_current_node->id() == 144) + std::cout << "Slip = " << slip << ".\n"; + + const Node * node = _current_node; + auto it = _bondslip.find(node->id()); + mooseAssert(it != _bondslip.end(), "Node not found in bond-slip map"); + bondSlipData bond_slip = it->second; + + bond_slip.slip_min = std::min(bond_slip.slip_min_old, slip); + bond_slip.slip_max = std::max(bond_slip.slip_max_old, slip); + + if (_debug) + // if (_current_node->id() == 144) + { + std::cout << "Slip_min = " << bond_slip.slip_min << ".\n"; + std::cout << "Slip_min_old = " << bond_slip.slip_min_old << ".\n"; + std::cout << "Slip_max = " << bond_slip.slip_max << ".\n"; + std::cout << "Slip_max_old = " << bond_slip.slip_max_old << ".\n"; + std::cout << "Bondstress_min = " << bond_slip.bondstress_min << ".\n"; + std::cout << "Bondstress_min_old = " << bond_slip.bondstress_min_old << ".\n"; + std::cout << "Bondstress_max = " << bond_slip.bondstress_max << ".\n"; + std::cout << "Bondstress_max_old = " << bond_slip.bondstress_max_old << ".\n"; + } + + Real slope = 5.0 * _max_bondstress / _transitional_slip[0]; + Real plastic_slip_max = bond_slip.slip_max - bond_slip.bondstress_max / slope; + Real plastic_slip_min = bond_slip.slip_min - bond_slip.bondstress_min / slope; + + if (slip >= bond_slip.slip_max || slip <= bond_slip.slip_min) + { + if (std::abs(slip) < _transitional_slip[0]) + { + if (_debug) + // if (_current_node->id() == 144) + std::cout << "Calculating bondstress for Case Ia" + << ".\n"; + _bond_stress = _max_bondstress * MathUtils::sign(slip) * + (5.0 * slip_ratio - 4.5 * slip_ratio * slip_ratio + + 1.4 * slip_ratio * slip_ratio * slip_ratio); + } + else if (slip >= _transitional_slip[0] && slip < _ultimate_slip) + { + if (_debug) + // if (_current_node->id() == 144) + std::cout << "Calculating bondstress for Case Ib" + << ".\n"; + _bond_stress = 1.9 * _max_bondstress; + } + else if (slip <= -_transitional_slip[0] && slip > -_ultimate_slip) + { + if (_debug) + // if (_current_node->id() == 144) + std::cout << "Calculating bondstress for Case Ic" + << ".\n"; + _bond_stress = -1.9 * _max_bondstress; + } + else + { + if (_debug) + // if (_current_node->id() == 144) + std::cout << "Calculating bondstress for Case Id" + << ".\n"; + _bond_stress = _frictional_bondstress * MathUtils::sign(slip); + } + } + else if (slip > plastic_slip_max && slip < bond_slip.slip_max) + { + if (_debug) + // if (_current_node->id() == 144) + std::cout << "Calculating bondstress for Case II" + << ".\n"; + + _bond_stress = (slip - plastic_slip_max) * bond_slip.bondstress_max / + (bond_slip.slip_max - plastic_slip_max); + } + else if (slip < plastic_slip_min && slip > bond_slip.slip_min) + { + if (_debug) + // if (_current_node->id() == 144) + std::cout << "Calculating bondstress for Case III" + << ".\n"; + + _bond_stress = (slip - plastic_slip_min) * bond_slip.bondstress_min / + (bond_slip.slip_min - plastic_slip_min); + } + else + _bond_stress = _frictional_bondstress; + + if (_debug) + // if (_current_node->id() == 144) + std::cout << "Bondstress = " << _bond_stress << "\n"; + + Real bond_force = 2.0 * libMesh::pi * _bar_radius * _current_elem_volume * _bond_stress; + + RealVectorValue constraint_force_axial = bond_force * _slave_tangent; + RealVectorValue constraint_force_normal = _penalty * slip_normal; + + _constraint_residual = constraint_force_axial + constraint_force_normal; + + if (_debug) + // if (_current_node->id() == 144) + { + std::cout << "Constraint Residual Axial = " << constraint_force_axial << "\n"; + std::cout << "Constraint Residual Normal = " << constraint_force_normal << "\n"; + std::cout << "Constraint Residual = " << _constraint_residual << "\n"; + } + + bond_slip.bondstress_min = std::min(bond_slip.bondstress_min_old, _bond_stress); + bond_slip.bondstress_max = std::max(bond_slip.bondstress_max_old, _bond_stress); + + if (_fe_problem.converged()) + { + _bondslip[node->id()].slip_min = bond_slip.slip_min; + _bondslip[node->id()].slip_max = bond_slip.slip_max; + _bondslip[node->id()].bondstress_min = bond_slip.bondstress_min; + _bondslip[node->id()].bondstress_max = bond_slip.bondstress_max; + } +} + +Real +RebarBondSlipConstraint::computeQpResidual(Moose::ConstraintType type) +{ + Real resid = _constraint_residual(_component); + + switch (type) + { + case Moose::Slave: + return _test_slave[_i][_qp] * resid; + + case Moose::Master: + return _test_master[_i][_qp] * -resid; + } + + return 0.0; +} + +Real +RebarBondSlipConstraint::computeQpJacobian(Moose::ConstraintJacobianType type) +{ + switch (type) + { + case Moose::SlaveSlave: + return _phi_slave[_j][_qp] * _penalty * _test_slave[_i][_qp] * + (1.0 - _slave_tangent(_component) * _slave_tangent(_component)); + + case Moose::SlaveMaster: + return -_phi_master[_j][_qp] * _penalty * _test_slave[_i][_qp] * + (1.0 - _slave_tangent(_component) * _slave_tangent(_component)); + + case Moose::MasterSlave: + return -_test_master[_i][_qp] * _penalty * _phi_slave[_j][_qp] * + (1.0 - _slave_tangent(_component) * _slave_tangent(_component)); + + case Moose::MasterMaster: + return _test_master[_i][_qp] * _penalty * _phi_master[_j][_qp] * + (1.0 - _slave_tangent(_component) * _slave_tangent(_component)); + + default: + mooseError("Unsupported type"); + break; + } + return 0.0; +} + +Real +RebarBondSlipConstraint::computeQpOffDiagJacobian(Moose::ConstraintJacobianType /*type*/, + unsigned int /*jvar*/) +{ + return 0.0; +}