Modify tests and add documentation idaholab#99

doc/content/bib/blackbear.bib

@@ -75,3 +75,21 @@ @Article{wald_2017
 		Title = {Development and multiaxial distribution of expansions in reinforced concrete elements affected by alkali--silica reaction},
 		Volume = {18},
 		Year = {2017}}
+
+@article{bondlsip_adina_1985,
+    Author = {Mehlhorn, G. and Kollegger, J. and Keuser, M. and Kolmar, W.},
+    Journal = {Computers and Structures},
+    Number = {},
+    Pages = {69--80},
+    Title = {Nonlinear Contact Problems - A Finite Element Approach Implemented in ADINA},
+    Volume = {21},
+    Year = {1985}}
+
+@article{hameed_2013,
+		Author = {Hameed, R. and Turatsinze, A. and Duprat, F. and Sellier A. },
+		Journal = {KSCE:Journal of Civil Engineering},
+		Month = jan,
+		Pages = {1700--1701},
+		Title = {Bond stress-slip Behavior of Reinforcing Bars Embedded in Hybrid Fiber-reinforced Concrete},
+		Volume = {17},
+		Year = {2013}}

doc/content/source/constraints/RebarBondSlipConstraint.md

@@ -0,0 +1,24 @@
+# RebarBondSlipConstraint
+
+!syntax description /Constraints/RebarBondSlipConstraint
+
+`RebarBondSlipConstraint` implements a node-to-element constraint  designed to apply the bond-slip relations between concrete and reinforcement bars. It uses a simplistic bond-slip model that assigns bond stress based on the slip values calculated as the relative displacement between concrete and rebar:
+\begin{equation}
+  \begin{aligned}
+    \sigma_s = & \sigma_{max} \left[ 0.5 \left(\frac{\Delta}{\Delta_1}\right) \, - \, 4.5 \, \left(\frac{\Delta}{\Delta_1}\right)^2 \, + \, 1.4 \, \left(\frac{\Delta}{\Delta_1}\right)^3 \right] \;  \mathrm{for} \, \Delta < \Delta_1 \, \,  \\
+    = & 1.9 \, \sigma_{max} \;  \mathrm{for}  \, \Delta >= \Delta_1
+  \end{aligned}
+\end{equation}
+Here, $\sigma_s$ is the bondstress, $\sigma_{max}$ is the maximum bondstress related to the compressive strength of the concrete, $\Delta$ is the slip calculated as the relative displacement between the concrete and rebar in the axial direction of the rebar, and $\Delta_1$ is the slip magnitude at which the maximum bondstress is reached. This model is similar to what was implemented in [!cite](bondlsip_adina_1985).
+
+## Example Input File Syntax
+
+!listing rebar_bondslip/RCBeam_constraint.i block=Constraints
+
+!syntax parameters /Constraints/RebarBondSlipConstraint
+
+!syntax inputs /Constraints/RebarBondSlipConstraint
+
+!syntax children  /Constraints/RebarBondSlipConstraint
+
+!bibtex bibliography

include/constraints/RebarBondSlipConstraint.h

@@ -17,8 +17,10 @@ class RebarBondSlipConstraint;
 
 template <>
 InputParameters validParams<RebarBondSlipConstraint>();
-
-/// A RebarBondSlipConstraint enforces concrete-rebar constraint
+/**
+ * A RebarBondSlipConstraint enforces the constraint between concrete and
+ * reinforcing bars establishing a slip vs. bondstress relationship
+ */
 class RebarBondSlipConstraint : public EqualValueEmbeddedConstraint
 {
 public:
@@ -40,14 +42,24 @@ class RebarBondSlipConstraint : public EqualValueEmbeddedConstraint
   bool getCoupledVarComponent(unsigned int var_num, unsigned int & component);
 
 protected:
+  /// method to calculate the tangential and the normal direction for the rebars
   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 designed to hold info about the bond-slip history
+   * slip_min miminum slip value at the current step
+   * slip_max maximum slip value at the current step
+   * slip_min_old minimum slip value from the history
+   * slip_max_old maximum slip value from the history
+   * bondstress_min miminum bondstress value at the current step
+   * bondstress_max maximum bondstress value at the current step
+   * bondstress_min_old minimum bondstress value from the history
+   * bondstress_max_old maximum bondstress value from the history
    */
   struct bondSlipData
   {
@@ -60,6 +72,7 @@ class RebarBondSlipConstraint : public EqualValueEmbeddedConstraint
     Real bondstress_min_old;
     Real bondstress_max_old;
 
+    /// initializing the bond-slip data
     bondSlipData()
       : slip_min(0.0),
         slip_max(0.0),
@@ -73,29 +86,55 @@ class RebarBondSlipConstraint : public EqualValueEmbeddedConstraint
     }
   };
 
-  // Bond-slip data
+  /// storing the bond-slip history values for each of the nodes
   std::map<dof_id_type, bondSlipData> _bondslip;
 
+  /// the direction in which the constraint works
   const unsigned _component;
+
+  /// problem dimesion
   const unsigned int _mesh_dimension;
+
+  /// displacement variables
   std::vector<unsigned int> _var_nums;
   std::vector<MooseVariable *> _vars;
 
+  /// flag to turn on printing values for debugging
   const bool _debug;
 
+  /// maximum bond stress
   const Real _max_bondstress;
+
+  /// residual bond stress due to friction after joint failure
   const Real _frictional_bondstress;
+
+  /// ultimate slip value attainable before failure
   const Real _ultimate_slip;
+
+  /// radius of the reinforcing bars
   const Real _bar_radius;
+
+  /// slip values at the transition points of the bond-slip curve
   std::vector<Real> _transitional_slip;
+
   /// constraint force needed to enforce the constraint
   RealVectorValue _constraint_residual;
+
+  /// constraint force needed to enforce the constraint
   RealVectorValue _constraint_jacobian_axial;
+
   /// penalty force for the current constraint
   RealVectorValue _pen_force;
+
+  /// tangent direction for the rebars
   RealVectorValue _secondary_tangent;
+
+  /// current element volume/length for the rabar
   Real _current_elem_volume;
-  bool _bond;
+
+  /// bond stress value
   Real _bond_stress;
+
+  /// redivative of the bond stress function w.r.t slip
   Real _bond_stress_deriv;
 };

src/constraints/RebarBondSlipConstraint.C

@@ -25,7 +25,7 @@ RebarBondSlipConstraint::validParams()
 {
   InputParameters params = EqualValueEmbeddedConstraint::validParams();
   params.addClassDescription(
-      "This is a constraint enforcing the bodslip behavior between concrete and rebar");
+      "This is a constraint enforcing the bod-slip behavior between concrete and rebar");
   params.addRequiredParam<unsigned int>("component",
                                         "An integer corresponding to the direction "
                                         "the variable this kernel acts in. (0 for x, "
@@ -100,12 +100,11 @@ 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;
-    }
+  {
+    std::cout << "===========================================\n";
+    std::cout << "node id: " << _current_node->id() << std::endl;
+    std::cout << "at coord: " << (Point)*_current_node << std::endl;
+  }
   auto it = _secondary_to_primary_map.find(_current_node->id());
 
   if (it != _secondary_to_primary_map.end())
@@ -157,8 +156,7 @@ RebarBondSlipConstraint::computeTangent()
   _current_elem_volume /= elems.size();
 
   if (_debug)
-    if (_current_node->id() == 144)
-      std::cout << "tangent: " << _secondary_tangent << std::endl;
+    std::cout << "tangent: " << _secondary_tangent << std::endl;
 }
 
 void
@@ -167,7 +165,6 @@ 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]);
@@ -176,10 +173,9 @@ RebarBondSlipConstraint::reinitConstraint()
   RealVectorValue slip_axial = slip * _secondary_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";
+    std::cout << "Slip = " << slip << ".\n";
 
   const Node * node = _current_node;
   auto it = _bondslip.find(node->id());
@@ -190,17 +186,16 @@ RebarBondSlipConstraint::reinitConstraint()
   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";
-    }
+  {
+    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;
@@ -212,10 +207,6 @@ RebarBondSlipConstraint::reinitConstraint()
   {
     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);
@@ -225,49 +216,21 @@ RebarBondSlipConstraint::reinitConstraint()
            1.4 * 3.0 * slip_ratio * slip_ratio / _transitional_slip[0]);
     }
     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);
 
     _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)
   {
-    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);
     _bond_stress_deriv = bond_slip.bondstress_min / (bond_slip.slip_min - plastic_slip_min);
@@ -276,11 +239,10 @@ RebarBondSlipConstraint::reinitConstraint()
     _bond_stress = _frictional_bondstress;
 
   if (_debug)
-    if (_current_node->id() == 144)
-    {
-      std::cout << "Bondstress = " << _bond_stress << "\n";
-      std::cout << "Bondstress Derivative = " << _bond_stress_deriv << "\n";
-    }
+  {
+    std::cout << "Bondstress = " << _bond_stress << "\n";
+    std::cout << "Bondstress Derivative = " << _bond_stress_deriv << "\n";
+  }
 
   Real bond_force = 2.0 * libMesh::pi * _bar_radius * _current_elem_volume * _bond_stress;
   Real bond_force_deriv =
@@ -292,14 +254,6 @@ RebarBondSlipConstraint::reinitConstraint()
   _constraint_residual = constraint_force_axial + constraint_force_normal;
   _constraint_jacobian_axial = bond_force_deriv * _secondary_tangent;
 
-  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);
 

test/tests/rebar_bondslip/RCBeam_constraint.i

@@ -76,29 +76,29 @@
     type = RebarBondSlipConstraint
     secondary = 2
     primary = 1
-    penalty = 1e6
+    penalty = 1e12
     variable = 'disp_x'
     primary_variable = 'disp_x'
     component = 0
     max_bondstress = 100
     transitional_slip_values = 0.0005
     ultimate_slip = 0.1
     rebar_radius = 7.98e-3
-    debug = true
+    # debug = true
   []
   [rebar_y]
     type = RebarBondSlipConstraint
     secondary = 2
     primary = 1
-    penalty = 1e6
+    penalty = 1e12
     variable = 'disp_y'
     primary_variable = 'disp_y'
     component = 1
     max_bondstress = 100
     transitional_slip_values = 0.0005
     ultimate_slip = 0.1
     rebar_radius = 7.98e-3
-    debug = true
+    # debug = true
   []
 []