From 95026d544d5fa3fdea47a7143caca9f41d3578ae Mon Sep 17 00:00:00 2001 From: Sudipta Biswas Date: Tue, 31 Mar 2020 16:18:17 -0600 Subject: [PATCH] Update bondslip constraint jacobian #99 --- include/constraints/RebarBondSlipConstraint.h | 2 + src/constraints/RebarBondSlipConstraint.C | 38 ++++++++++++++----- 2 files changed, 30 insertions(+), 10 deletions(-) diff --git a/include/constraints/RebarBondSlipConstraint.h b/include/constraints/RebarBondSlipConstraint.h index 6bdab89b2..dc22bd1d1 100644 --- a/include/constraints/RebarBondSlipConstraint.h +++ b/include/constraints/RebarBondSlipConstraint.h @@ -81,10 +81,12 @@ class RebarBondSlipConstraint : public EqualValueEmbeddedConstraint std::vector _transitional_slip; /// constraint force needed to enforce the constraint RealVectorValue _constraint_residual; + RealVectorValue _constraint_jacobian_axial; /// penalty force for the current constraint RealVectorValue _pen_force; RealVectorValue _slave_tangent; Real _current_elem_volume; bool _bond; Real _bond_stress; + Real _bond_stress_deriv; }; diff --git a/src/constraints/RebarBondSlipConstraint.C b/src/constraints/RebarBondSlipConstraint.C index 7d94e2d28..8815527a4 100644 --- a/src/constraints/RebarBondSlipConstraint.C +++ b/src/constraints/RebarBondSlipConstraint.C @@ -126,7 +126,7 @@ RebarBondSlipConstraint::shouldApply() void RebarBondSlipConstraint::computeTangent() { - _slave_tangent *= 0.0; + _slave_tangent = 0.0; // get normals // get connected elements of the current node @@ -138,7 +138,7 @@ RebarBondSlipConstraint::computeTangent() for (auto & elem : elems) { Elem * elem_ptr = _mesh.elemPtr(elem); - // _assembly.reinit(elem_ptr, 0); + _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"), @@ -204,6 +204,8 @@ RebarBondSlipConstraint::reinitConstraint() 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; + _bond_stress_deriv = 0.0; + if (slip >= bond_slip.slip_max || slip <= bond_slip.slip_min) { if (std::abs(slip) < _transitional_slip[0]) @@ -215,6 +217,10 @@ RebarBondSlipConstraint::reinitConstraint() _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); + _bond_stress_deriv = + _max_bondstress * MathUtils::sign(slip) * + (5.0 / _transitional_slip[0] - 4.5 * 2.0 * slip_ratio / _transitional_slip[0] + + 1.4 * 3.0 * slip_ratio * slip_ratio / _transitional_slip[0]); } else if (slip >= _transitional_slip[0] && slip < _ultimate_slip) { @@ -250,6 +256,8 @@ RebarBondSlipConstraint::reinitConstraint() _bond_stress = (slip - plastic_slip_max) * bond_slip.bondstress_max / (bond_slip.slip_max - plastic_slip_max); + + _bond_stress_deriv = bond_slip.bondstress_max / (bond_slip.slip_max - plastic_slip_max); } else if (slip < plastic_slip_min && slip > bond_slip.slip_min) { @@ -260,6 +268,7 @@ RebarBondSlipConstraint::reinitConstraint() _bond_stress = (slip - plastic_slip_min) * bond_slip.bondstress_min / (bond_slip.slip_min - plastic_slip_min); + _bond_stress_deriv = bond_slip.bondstress_min / (bond_slip.slip_min - plastic_slip_min); } else _bond_stress = _frictional_bondstress; @@ -269,11 +278,14 @@ RebarBondSlipConstraint::reinitConstraint() std::cout << "Bondstress = " << _bond_stress << "\n"; Real bond_force = 2.0 * libMesh::pi * _bar_radius * _current_elem_volume * _bond_stress; + Real bond_force_deriv = + 2.0 * libMesh::pi * _bar_radius * _current_elem_volume * _bond_stress_deriv; RealVectorValue constraint_force_axial = bond_force * _slave_tangent; RealVectorValue constraint_force_normal = _penalty * slip_normal; _constraint_residual = constraint_force_axial + constraint_force_normal; + _constraint_jacobian_axial = bond_force_deriv * _slave_tangent; if (_debug) // if (_current_node->id() == 144) @@ -315,23 +327,29 @@ RebarBondSlipConstraint::computeQpResidual(Moose::ConstraintType type) Real RebarBondSlipConstraint::computeQpJacobian(Moose::ConstraintJacobianType type) { + Real jac_axial = _constraint_jacobian_axial(_component); + switch (type) { case Moose::SlaveSlave: - return _phi_slave[_j][_qp] * _penalty * _test_slave[_i][_qp] * - (1.0 - _slave_tangent(_component) * _slave_tangent(_component)); + return _phi_slave[_j][_qp] * jac_axial * _slave_tangent(_component) * _test_slave[_i][_qp] + + _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)); + return -_phi_master[_j][_qp] * jac_axial * _slave_tangent(_component) * _test_slave[_i][_qp] - + _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)); + return -_test_master[_i][_qp] * jac_axial * _slave_tangent(_component) * _phi_slave[_j][_qp] - + _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)); + return _test_master[_i][_qp] * jac_axial * _slave_tangent(_component) * _phi_master[_j][_qp] + + _test_master[_i][_qp] * _penalty * _phi_master[_j][_qp] * + (1.0 - _slave_tangent(_component) * _slave_tangent(_component)); default: mooseError("Unsupported type");