Update bondslip constraint jacobian #99

idaholab · Mar 31, 2020 · 95026d5 · 95026d5
1 parent 44709c8
commit 95026d5
Show file tree

Hide file tree

Showing 2 changed files with 30 additions and 10 deletions.
diff --git a/include/constraints/RebarBondSlipConstraint.h b/include/constraints/RebarBondSlipConstraint.h
@@ -81,10 +81,12 @@ class RebarBondSlipConstraint : public EqualValueEmbeddedConstraint
   std::vector<Real> _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
@@ -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<Order>("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");