Analysis of Bolts: Modified Visual Basic Code for Calculating Bolt Coefficient Using a Method Developed by Donald Brandt

Another blog post for Structural Engineers and Designers…

When I wrote and posted the VBA (Visual Basic for Application) Code for the analysis of eccentrically loaded bolt group in my blogpost, “Eccentrically Loaded Bolt Group: Analysis and Design using 'INSTANTANEOUS CENTER of ROTATION METHOD' with VBA for Excel Implementation”, I was aware of the difference between the values returned by the program and that of the table of C coefficients provided by AISC Steel Construction Manual. I never believed, however, that the values returned by the program are wrong. How could it be when the VBA code is based on the work by Crawford & Kulak as presented in AISC Steel Construction Manual 13th Edition, pages 7-6 to 7-8? Despite the differences, though, I still decided to post the VBA Code.

Not long afterwards, engineers who happened to follow my blog shared and expressed the same observations and 'demanded' an explanation. Here are some of their comments:

1. Blog comment from Josh Gionfriddo, EIT, Structural Engineer, Structures Workshop, Inc. on November 19, 2010:

   ---

   Thanks for this great code. One thing though... The values in the AISC manual seem to hover around 98% of the values I get from your code. Any Idea why the differences?
   
   
2. Comment from Gregory Vidgop via AISC SteelTOOLS on October 11, 2010:

   ---

   Hello Redem,
   
   
   Thank you so much for your effort in developing the very nice and useful program as "SteelPro". Just would like to ask you a simple question: why is the C coefficient for the bolt group analysis slightly different from the book value and C' coefficient is exactly right?
   
   
   Please let me know when your program will be finally issued or any revisions are planned in the future.
   
   
   Thank you again,
   
   
   Gregory
   
   
3. Blog comment from Steven Vicha, P.E., Jacobs Engineering Group, Inc. on February 9, 2011:

   ---

   Great code, but I have a couple of questions. In the code line "Rn = 74 * (1 - Exp(-10 * 0.34)) ^ 0.55" what is the 74? Also the answers are a little off from the tables, any idea why?
   
   
   I need to study it some more.
   
   

The three comments above have something in common. All are asking why the values returned by the VBA Code are a little off from that provided by the table in AISC Manual. I was very hesitant to respond to the first two comments above until I received another one from Steven Vicha. Fortunately, he got a reply:

Hi Steven... You've got the same predicament with josh. Actually, I planned not to respond to the question while still looking for the answer -- until i got another question from you.

I did a review on the equation provided for by AISC Steel Construction Manual 13th Edition on pages 7-6 to 7-8 and ran more analyses -- but still the results are a little bit off from the table. Despite the difference, however, the results returned by the code are still correct: the code SATISFIES the in-plane static equilibrium up to 4 decimal places of accuracy. We may throw the same question to AISC why the differences...

Please give me more time to look into this matter.

As to your question, 'what is 74?' -- the number is the ultimate shear strength of 3/4-in. diameter ASTM A325 bolt, 74 kips, which corresponds to its maximum deformation of 0.34 inches. You may refer to pages 7-6 to 7-8 of the AISC Steel Construction Manual, 13th Edition for a thorough discussion.

Shortly after replying to Steven’s question, I found myself indebted to him for his effort in finding the real answer to his question and of the others. He conducted his own research on the topic and pin-pointed out, bull’s-eye, the precise reason. In his February 14, 2011 follow-up comment to the same blogpost, Steven added:

Thanks for the response. I did find the 74, finally. I think that I may have found the reason for the difference in the tables and the results of your VBA code. There is a paper posted on AISC.org that may explain the difference. I have the paper on my other PC at home. I will try and send the excerpt.

Steven is referring to the paper by Larry S. Muir, P.E. and William A. Thornton, P.E., PhD, both of Cives Steel Company, “Exploring the True Geometry of the Inelastic Instantaneous Center Method for Eccentrically Loaded Bolt Groups”. The paper presented the true location of the instantaneous center (I.C.) of a bolt group. In the paper, the authors argued that the I.C. is not necessarily located along the axis perpendicular to the applied load. The assumption that the I.C. of a bolt group lies along the axis perpendicular to the applied load is incorrect for all but a few very specific instances. Finally, the authors concluded that the search for instantaneous center cannot be limited to the one-dimensional line, but must include all points in the two-dimensional plane.

The authors are actually embracing and employing Donald Brandt’s method in determining the instantaneous center. The Brandt’s method is the very same method being used to calculate the values presented in AISC Manual. The authors have the following remarks as to the calculation of the values presented in the table. Excerpts:

THE AISC COEFFICIENTS C FOR ECCENTRICALLY LOADED BOLT GROUPS

Finding errors in the presentation of the theory underlying the calculation of the instantaneous center calls into question the values in Tables 7-17 through 7-24 (Tables 7-7 through 7-14 in AISC Manual 13^th Edition) presented in the AISC Manual. Fortunately these values were produced by a program based on Brandt’s work. Brandt’s procedure never restricts its search for the instantaneous center to the line perpendicular to the applied load and passing through the centroid of the connection, though from a programming standpoint this would seem to be the most efficient approach. Instead both the location of the instantaneous center and the angle of the resultant are allowed to drift off of the values implied by Crawford and Kulak’s theory until equilibrium is satisfied. Being based solely on the load-deformation relationship and equilibrium, the C-values presented by AISC are correct. The authors have checked numerous cases and are satisfied that both the load-deformation relationship and equilibrium are satisfied using Brandt’s approach.

So, finally the question has been answered through Steven's effort and resourcefulness. Grateful for Steven’s effort, I responded:

Thanks, Steve, for the excerpts...

To complete the reading, I downloaded the same article written by Larry Muir and Thornton from the AISC website. After reading, I now understand why the differences between the values returned by the VBA code and that of the table provided by AISC Steel Construction Manual. It is all about Crawford & Kulak's method against that of Brandt’s approach. The table is derived using Brandt's while the VBA Code is based on Crawford and Kulak's.

Having learned that the values in Tables presented in AISC Manual were produced by a program based on Brandt’s work, I will try to write another VBA Code based on the same approach -- on my subsequent blogposts...

Graphically, the figures below show how the 'Crawford & Kulak's Approach' and 'Donald Brandt's Approach' differ from each other:

Instantaneous Center in Crawford and Kulak's Approach	Instantaneous Center in Donald Brandt's Approach
In Crawford and Kulak's Approach, the instantaneous center is searched and limited only to the one-dimensional line that passes thru the centroid of the bolt group and perpendicular to the applied load. In contrast, the search for instantaneous center in Brandt's Approach includes all points in a two-dimensional plane. The two methods, however, would agree on results when the geometry of the bolt group is symmetrical about its centroidal axis parallel to the applied load.

And finally, here's the Modified VBA Code for calculating the Coefficient C for Eccetrically Loaded Bolt Group. With the VBA Code now modified based on Brandt's Approach, the resulting values for bolt coefficient C now hover between 98% and 99.9% relative to the values provided by the tables in AISC Manual with values returned by the VBA Code being the higher ones.

VBA Code: BOLT COEFFICIENT CALCULATOR Employing Brandt's Approach

 
''' 
'''Redem Legaspi Jr 
'''March 14, 2011 
'''Current Location: Manila, Philippines 
'''Hometown: Pio V. Corpus, Masbate, Philippines  

Option Explicit  

Type BoltInfo    
  Dv As Double    
  Dh As Double 
End Type   

'''Effective Bolt Coefficient 
Function BoltCoefficient(Bolt_Row As Integer, Bolt_Column As Integer, Row_Spacing As Double, Column_Spacing As Double, Eccentricity As Double, Optional Rotation As Double = 0) As Double     
   Dim i, k, n As Integer     
   Dim Fm, Fyy As Double     
   Dim xRo, yRo As Double     
   Dim Mo, Fy As Double     
   Dim Fx, Fxx As Double     
   Dim xi, yi As Double     
   Dim ri, x1 As Double     
   Dim y1, Rot As Double     
   Dim Rn, iRn As Double     
   Dim Rv, Rh As Integer     
   Dim Sh, Sv As Double     
   Dim Ec As Double     
   Dim Delta, rmax As Double     
   Dim BoltLoc() As BoltInfo     
   Dim Stp As Boolean     
   Dim J As Double     
   Dim MapFunc As Double     
   Dim MirrFlag As Integer          
   Rv = Bolt_Row     
   Rh = Bolt_Column     
   Sv = Row_Spacing     
   Sh = Column_Spacing          
   ReDim BoltLoc(Rv * Rh - 1)          
   
   On Error Resume Next          

   Rot = Application.WorksheetFunction.Radians(Rotation)     
   Ec = Abs(Eccentricity * Cos(Rot))     
   MirrFlag = IIf(Sgn(Cos(Rot)) = -1, -1, 1)          
   
  If Ec = 0 Then GoTo ForcedExit          
     n = 0     
     For i = 0 To Rv - 1         
        For k = 0 To Rh - 1             
           y1 = (i * Sv) - (Rv - 1) * Sv / 2             
           x1 = (k * Sh) - (Rh - 1) * Sh / 2             
           With BoltLoc(n)                 
              .Dv =  x1 * Sin(Rot) + y1 * Cos(Rot)              '''Rotate Vertical Coordinate                 
              .Dh = (x1 * Cos(Rot) - y1 * Sin(Rot)) * MirrFlag  '''Rotate Horizontal Coordinate             
           End With             
           n = n + 1         
       Next     
    Next          
    
    xRo = 0     
    yRo = 0     
    Stp = False          
    ''Rn = 74 * (1 - Exp(-10 * 0.34)) ^ 0.55     '''Original Code     (inaccurate)
    Rn = 74                                      '''Per David Fanning (very accurate) 
    Do While Stp = False         
       rmax = 0         
       For i = 0 To Rv * Rh - 1             
          xi = BoltLoc(i).Dh + xRo
          yi = BoltLoc(i).Dv + yRo
          rmax = Application.WorksheetFunction.Max(rmax, Sqr(xi ^ 2 + yi ^ 2))             
          DoEvents         
       Next                  
       Fx = 0:  J = 0
       Fy = 0: Mo = 0
       For i = 0 To Rv * Rh - 1
          xi = BoltLoc(i).Dh + xRo
          yi = BoltLoc(i).Dv + yRo
          ri = Sqr(xi ^ 2 + yi ^ 2)
          Delta = 0.34 * ri / rmax
          iRn = 74 * (1 - Exp(-10 * Delta)) ^ 0.55
          Mo = Mo + (iRn / Rn) * ri              '''Moment
          Fy = Fy + (iRn / Rn) * (xi / ri)       '''Vertical Force            
          Fx = Fx + (iRn / Rn) * (yi / ri)       '''Horizontal Force             
          J = J + 0.015531555 + 0.441786467 * ri ^ 2
          DoEvents
       Next
          
       Fm = Mo / (Ec + xRo)
       Fyy = Fm - Fy                             '''Vertical Unbalanced Force
       Fxx = -Fx                                 '''Horizontal Unbalanced Force
       
       Stp = Abs(Fyy) <= 0.00001 And Abs(Fxx) <= 0.00001
        
       MapFunc = J / (Rv * Rh * Mo)              '''Mapping function
       xRo = xRo + Fyy * MapFunc                 '''I.C. location on x-axis from bolt group C.G.
       yRo = yRo + Fxx * MapFunc                 '''I.C. location on y-axis from bolt group C.G.
             
       DoEvents
    Loop
    BoltCoefficient = (Fm + Fy) / 2
    
    Exit Function
ForcedExit:
    BoltCoefficient = Rv * Rh
End Function

Note:

---

To implement and run the VBA Code in Microsoft Excel Environment, please read my previous blogpost, "Eccentrically Loaded Bolt Group: Analysis and Design using 'INSTANTANEOUS CENTER of ROTATION METHOD' with VBA for Excel Implementation", for instructions.

---

ALTERNATE VBA Code for BOLT COEFFICIENT CALCULATOR (Updated)

'''   
'''Redem Legaspi Jr
'''June 24, 2015
'''Current Location: Manila, Philippines
'''Hometown: Pio V. Corpus, Masbate, Philippines  

Option Explicit

Type BoltInfo
  Dv As Double
  Dh As Double
End Type

'''Effective Bolt Coefficient
Function BoltCoefficient(Bolt_Row As Integer, Bolt_Column As Integer, Row_Spacing As Double, Column_Spacing As Double, Eccentricity As Double, Optional Rotation As Double = 0) As Double
   Dim i, k, n As Integer
   Dim P, uFy As Double
   Dim xRo, yRo As Double
   Dim M, Ry As Double
   Dim Rx, uFx As Double
   Dim xi, yi As Double
   Dim ri, x1 As Double
   Dim y1, Rot As Double
   Dim Ru, iRn As Double
   Dim Rv, Rh As Integer
   Dim Sh, Sv As Double
   Dim Ec As Double
   Dim Delta, LiMax As Double
   Dim BoltLoc() As BoltInfo
   Dim Stp As Boolean
   Dim J As Double
   Dim MapFunc As Double
   Dim ro As Double
   Dim Py As Double
   Dim Px As Double
   Rv = Bolt_Row
   Rh = Bolt_Column
   Sv = Row_Spacing
   Sh = Column_Spacing
   ReDim BoltLoc(Rv * Rh - 1)
   
   On Error Resume Next

   Rot = Application.WorksheetFunction.Radians(Rotation)
   Ec = Abs(Eccentricity)
   
  If Ec = 0 Then GoTo ForcedExit
     n = 0
     For i = 0 To Rv - 1
        For k = 0 To Rh - 1
           y1 = (i * Sv) - (Rv - 1) * Sv / 2
           x1 = (k * Sh) - (Rh - 1) * Sh / 2
           With BoltLoc(n)
              .Dv = y1
              .Dh = x1
           End With
           n = n + 1
       Next
    Next
    
    xRo = 0
    yRo = 0
    Stp = False
    Ru = 74
    Do While Stp = False
       LiMax = 0
       For i = 0 To Rv * Rh - 1
          xi = BoltLoc(i).Dh + xRo
          yi = BoltLoc(i).Dv + yRo
          LiMax = Application.WorksheetFunction.Max(LiMax, Sqr(xi ^ 2 + yi ^ 2))
          DoEvents
       Next
       Rx = 0:  J = 0
       Ry = 0:  M = 0
       For i = 0 To Rv * Rh - 1
          xi = BoltLoc(i).Dh + xRo
          yi = BoltLoc(i).Dv + yRo
          ri = Sqr(xi ^ 2 + yi ^ 2)
          Delta = 0.34 * ri / LiMax
          iRn = 74 * (1 - Exp(-10 * Delta)) ^ 0.55
          M = M + (iRn / Ru) * ri                       '''Moment
          Ry = Ry + (iRn / Ru) * (xi / ri)              '''Vertical Force
          Rx = Rx + (iRn / Ru) * (yi / ri)              '''Horizontal Force
          J = J + ri ^ 2
          DoEvents
       Next
       ro = (Ec + xRo) * Cos(Rot) + yRo * Sin(Rot)
       P = M / ro
       Py = P * Cos(Rot)
       Px = P * Sin(Rot)
       
       uFy = Py - Ry                                    '''Vertical Unbalanced Force
       uFx = Px - Rx                                    '''Horizontal Unbalanced Force
       
       Stp = Abs(uFy) <= 0.00001 And Abs(uFx) <= 0.00001
        
       MapFunc = J / (Rv * Rh * M)                      '''Mapping function
       xRo = xRo + uFy * MapFunc                        '''I.C. location on x-axis from bolt group C.G.
       yRo = yRo + uFx * MapFunc                        '''I.C. location on y-axis from bolt group C.G.
             
       DoEvents
    Loop
    BoltCoefficient = P
    
    Exit Function
ForcedExit:
    BoltCoefficient = Rv * Rh
End Function

ALTERNATE VBA Code for BOLT COEFFICIENT CALCULATOR

(Updated Code based on Mark Lasby's suggestion)

 
'''
'''Redem Legaspi Jr
'''Aug 8, 2015
'''Current Location: Manila, Philippines
'''Hometown: Pio V. Corpus, Masbate, Philippines

Option Explicit

Type BoltInfo
  Dv As Double
  Dh As Double
End Type

'''Effective Bolt Coefficient
Function BoltCoefficient(Bolt_Row As Integer, Bolt_Column As Integer, Row_Spacing As Double, Column_Spacing As Double, Eccentricity As Double, Optional Rotation As Double = 0) As Double
   Dim i, k, n As Integer
   Dim P, uFy As Double
   Dim xRo, yRo As Double
   Dim M, Ry As Double
   Dim Rx, uFx As Double
   Dim xi, yi As Double
   Dim ri, x1 As Double
   Dim y1, Rot As Double
   Dim Ru, iRn As Double
   Dim Rv, Rh As Integer
   Dim Sh, Sv As Double
   Dim Ec As Double
   Dim Delta, LiMax As Double
   Dim BoltLoc() As BoltInfo
   
   Dim Stp As Boolean
   Dim Stp1 As Boolean
   Dim Stp2 As Boolean
   Dim Stp3 As Boolean
   
   Dim J As Double
   Dim MapFunc As Double
   Dim ro As Double
   Dim Py As Double
   Dim Px As Double
   Rv = Bolt_Row
   Rh = Bolt_Column
   Sv = Row_Spacing
   Sh = Column_Spacing
   ReDim BoltLoc(Rv * Rh - 1)
   
   Const CntMax = 5000
   Dim Cnt As Double
   Dim Pprev As Double
   Dim Pp As Double
   
   On Error Resume Next

   Rot = Application.WorksheetFunction.Radians(Rotation)
   Ec = Abs(Eccentricity)
   
  If Ec = 0 Then GoTo ForcedExit
     n = 0
     For i = 0 To Rv - 1
        For k = 0 To Rh - 1
           y1 = (i * Sv) - (Rv - 1) * Sv / 2
           x1 = (k * Sh) - (Rh - 1) * Sh / 2
           With BoltLoc(n)
              .Dv = y1
              .Dh = x1
           End With
           n = n + 1
       Next
    Next
    
    xRo = 0
    yRo = 0
    Stp = False
    Ru = 74
    Do While Stp = False
       Cnt = Cnt + 1
       LiMax = 0
       For i = 0 To Rv * Rh - 1
          xi = BoltLoc(i).Dh + xRo
          yi = BoltLoc(i).Dv + yRo
          LiMax = Application.WorksheetFunction.Max(LiMax, Sqr(xi ^ 2 + yi ^ 2))
          DoEvents
       Next
       Rx = 0:  J = 0
       Ry = 0:  M = 0
       For i = 0 To Rv * Rh - 1
          xi = BoltLoc(i).Dh + xRo
          yi = BoltLoc(i).Dv + yRo
          ri = Sqr(xi ^ 2 + yi ^ 2)
          Delta = 0.34 * ri / LiMax
          iRn = 74 * (1 - Exp(-10 * Delta)) ^ 0.55
          M = M + (iRn / Ru) * ri                       '''Moment
          Ry = Ry + (iRn / Ru) * (xi / ri)              '''Vertical Force
          Rx = Rx + (iRn / Ru) * (yi / ri)              '''Horizontal Force
          J = J + ri ^ 2
          DoEvents
       Next
       ro = (Ec + xRo) * Cos(Rot) + yRo * Sin(Rot)
       P = M / ro
       Py = P * Cos(Rot)
       Px = P * Sin(Rot)
       
       uFy = Py - Ry                                    '''Vertical Unbalanced Force
       uFx = Px - Rx                                    '''Horizontal Unbalanced Force
       
       Pp = (uFy ^ 2 + uFx ^ 2) ^ 0.5                       ''<--Code courtesy of Mark Lasby
       
       Stp1 = Abs(uFy) <= 0.00001 And Abs(uFx) <= 0.00001
       Stp2 = (Cnt > CntMax)                                ''<--Code courtesy of Mark Lasby
       Stp3 = IIf(Cnt > 1, Pp > Pprev, False)               ''<--Code courtesy of Mark Lasby
       Stp = (Stp1 Or Stp2 Or Stp3)                         ''<--Code courtesy of Mark Lasby
        
       Pprev = Pp                                           ''<--Code courtesy of Mark Lasby
         
       MapFunc = J / (Rv * Rh * M)                      '''Mapping function
       xRo = xRo + uFy * MapFunc                        '''I.C. location on x-axis from bolt group C.G.
       yRo = yRo + uFx * MapFunc                        '''I.C. location on y-axis from bolt group C.G.
             
       DoEvents
    Loop
    BoltCoefficient = P
    
    Exit Function
ForcedExit:
    BoltCoefficient = Rv * Rh
End Function