Eccentrically Loaded Bolt Group: BOLT COEFFICIENT CALCULATOR ENHANCED VISUAL BASIC CODE

I just want to share an enhancement to the BOLT COEFFICIENT CALCULATOR VBA Code made by an engineer friend from Canada.

Martin Boisclair, a project engineer from Canam-Heavy of Canam Group Inc., sent me an email sometime ago requesting for a copy of an updated version of the BOLT COEFFICIENT CALCULATOR VBA code. I gave him the code and then an exchange of emails ensued. Fortunately, our exchange of emails turned out to be more than just formalities -- but an exchange of bright ideas. I was so lucky to find that Martin is himself a good programmer. It was good to know that we are on the same level of thinking; basically on the 'same brainwave frequency'. In effect, he helped me improve the code.

Martin found out that the issue about magic numbers has not been fully resolved. The problem in convergence to equilibrium still exists for some input combinations regardless of the modifications previously made to the code. Magic numbers are not really magical as it is just a term being used to describe an input combination that causes the force system not to converge to equilibrium -- that brings about the problem that we have been trying to solve for quite a long time. Martin's enhancement to the code finally solved the problem as he explains in one of his lengthy emails below:

---

Redem,

First, I want to thank you again for sharing your work and I hope you find my work useful. I hope your mother is better. I also want to add that in no way I was expecting that you would respond during that time or during weekends. I only count on your kindness.

My first thought when I looked at the results (with the code revised with the Mark Lasby add)... and after running and comparing with all AISC value, was that I found some non-conservative results and this was not acceptable. I didn't like that. Also, at Canam, we're often challenged by Engineer of Record.

Ok... so... it is not acceptable that a C factor is returned if there's no convergence. The first thing I did was to remove the Stp3 code because this was the reason why values did not converge. Stp3 makes that... in some instances... we were returned with non-conservative values. On the other hand, removing Stp3, made that the function did not converge and stopped and CntMax. Roughly, ~95% of the AISC table worked but we were stuck with some rebellious values (~5% of the table). Magic numbers had to be taken care of.

It took about 2 weeks to come up with the revised VBA code. To make it easier for you to understand the revised code, I attached a pdf file with some comments for you. Note that most of this code revision is based on observation/experiment of the numerical behavior of MapFunc and Pp.

I looked at the behavior of some variables with Debug.Print and figured Magic1 (Bouncing Pp and Bouncing MapFunc). To those values, I wanted out of "MapFunc = J / (Rv * Rh * M)" to try playing with MapFunc by other means. I tried linear interpolation between MapFunc, Pp and MapFuncLast, PpLast but it didn't work. Then, I came up with the idea of manually incrementing MapFunc and changing the direction of the increment as required. It worked but it had to be tweaked and I figured that Incr should not go below a certain value. I figured 0.01 as a minimum value of Incr. Also, sometimes Mapfunc went below 0 and I also had to tweak the VBA code to prevent that from happening.

I ran all the AISC tables and found out that Magic1 solves ~4.5% of the rebellious 5%. I was still stucked with ~0.5% of values that didn't converge. So, I looked at the numerical behavior of the last 0.5% rebellious rebels and figured a way to detect them. After some additional tweaking of the Alternate method of iterating MapFunc, all AISC tables ran and converged.

Iterating for the C factor is always a fragile balance act. It's very sensitive and small differences may make that it does not converge. For instance, with this revised code, if I change all "Dim As Double" to "Dim As Single"... some values don't converge!!

I spoke with the guy that programmed our own Canam code (not as a function) and he told me the same. So even today, our own Canam code made in 1996, even though this code solves all AISC table values, it still don't converge on some combination of inputs. Here's an example where our Canam code don’t converge (2,1,3,0,9.5,78.7). On the other hand the revised function works for this example.

Now, I think it's fair to say that the revision of VBA code I came up with, will not work ALL the time... but it's still an improvement.

Thanks for your support. Awaiting your comments!

Best regards,

Martin Boisclair,P.Eng.
Project engineer
Canam-Heavy, a division of Canam Group Inc.
P 450-641-4000; 3555 / M 514-702-4343
martin.boisclair@canamgroupinc.com

---

BOLT COEFFICIENT CALCULATOR ENHANCED VISUAL BASIC CODE

(Code update by Martin Boisclair, P.Eng.)

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

'''''' MODIFICATION BY Martin Boisclair
'''''' CANAM Group - Montreal Canada
'''''' Warning - Use at your own risk
'''''' I added some minor correction to Dim definition per Geoff comment 10/02/2015
'''''' AND Rx and Ry calculation sometimes returned an overflow result. I added a few lines to avoid division by 0 (see below)
'''''' I removed Stp3 condition. With Stp3, the function stepped out prematurely without actual Pp converging
'''''' I modified the routine so that magic numbers DO converge by bypassing MapFunc = J / (Rv * Rh * M) and imposing another method of iterating MapFunc.

Option Explicit

Public Cnt As Double   'MARTIN - delete - WARNING - Declared "public" so that value Cnt passed thru the filling table macro to print this value in the results. Delete when done and remove Dim Cnt lower below

Type BoltInfo
    Dv As Double
    Dh As Double
End Type

'''Effective Bolt Coefficient
Function BoltCoefficientTest(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 As Integer, k As Integer, n As Integer
    Dim P As Double, uFy As Double
    Dim xRo As Double, yRo As Double
    Dim M As Double, Ry As Double
    Dim Rx As Double, uFx As Double
    Dim xi As Double, yi As Double
    Dim ri As Double, x1 As Double
    Dim y1 As Double, Rot As Double
    Dim Ru As Integer, iRn As Double
    Dim Rv As Integer, Rh As Integer
    Dim Sh As Double, Sv As Double
    Dim Ec As Double
    Dim Delta As Double, LiMax As Double
    Dim BoltLoc() As BoltInfo
    
    Dim Stp As Boolean
    Dim Stp1 As Boolean
    Dim Stp2 As Boolean
    Dim Magic As Boolean, Magic1 As Boolean, Magic2 As Boolean              'MARTIN
    Magic = False                       'MARTIN
    Magic1 = False                      'MARTIN
    Magic2 = False                      'MARTIN
    Dim Incr As Double, Dir As Integer  'MARTIN
    
    Dim J As Double
    Dim MapFunc As Double, MapFuncLast As Double, MapFuncMin As Double      'MARTIN
    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 = 4999             'MARTIN Safety exit - Will return "Not Converging Error"
'    Dim Cnt As Double              'MARTIN - This variable was declared "Public". Put this line back when done
    Dim CntTemp As Double           'MARTIN For type 2 Magic numbers
    Dim Pp As Double, PpLast As Double, PpLast2 As Double, PpMin As Double      'MARTIN
    PpMin = 9999999                 'MARTIN
    
    On Error Resume Next
    
    Rot = Application.WorksheetFunction.Radians(Rotation)
    Ec = Abs(Eccentricity)
    Cnt = 0                     'MARTIN - This is required since Cnt is "Public" at this time. Remove later
    CntTemp = 0                 'MARTIN
   
'**************** Calc Core START ****************
    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
                If ri = 0 Then                  'MARTIN to mitigate the division by zero in Ry and Rx calc when ri=0
                    Ry = Ry                     'MARTIN to mitigate the division by zero in Ry and Rx calc when ri=0
                    Rx = Rx                     'MARTIN to mitigate the division by zero in Ry and Rx calc when ri=0
                Else
                    Ry = Ry + (iRn / Ru) * (xi / ri)        '''Vertical Force
                    Rx = Rx + (iRn / Ru) * (yi / ri)        '''Horizontal Force
                End If
                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
            
            PpLast2 = PpLast                    'MARTIN
            PpLast = Pp                         'MARTIN
            Pp = (uFy ^ 2 + uFx ^ 2) ^ 0.5
            
'**************** Calc Core END ****************
            
Debug.Print Cnt & " P=" & Round(P, 3) & " " & CntTemp & " " & Magic & " ", Round(Incr * Dir, 6), " With MapFunc=" & Round(MapFunc, 4), " We get Pp=" & Round(Pp, 4)
            
            Stp1 = Abs(uFy) <= 0.001 And Abs(uFx) <= 0.001  'MARTIN reduced from 0.00001 to 0.001 - Verified that All (Fonction value)/(AISC14th table) ratios where between 0.983 and 1.013. Less than that and the results become less conservative.
            Stp2 = (Cnt > CntMax)
            Stp = (Stp1 Or Stp2)                            'MARTIN

'********************************************************************************
'MARTIN:
'Magic numbers are not magic. It's just a term used here that means that
'convergence is just not possible with the use of MapFunc = J / (Rv * Rh * M)
'because of some numerical instabilities. With the modified code below, all AISC
'table were ran and all converge but that doesn't mean that all combination will.
'It is expected that some other values will not converge. In this event the
'function will return a result of 0.01
'I found 2 types of numerical instabilities:
'Magic1, MapFunc and Pp bounce between to identical numbers and don't converge
'Magic2, MapFunc bounces erratically and Pp does not converge
'********************************************************************************

            If Magic = False Then                           'MARTIN
                If Pp > PpLast And Round(Pp, 10) = Round(PpLast2, 10) Then Magic1 = True      'MARTIN - Magic Type 1 detection
                If Cnt > 25 And (Pp > PpMin And PpLast > PpMin And PpLast2 > PpMin) And (MapFunc > MapFuncLast) Then CntTemp = CntTemp + 1  'MARTIN - Magic Type 2 detection
                If CntTemp > 25 Then Magic2 = True          'MARTIN - Magic Type 2 detection
                If Magic1 Or Magic2 Then                    'MARTIN
                    Magic = True                            'MARTIN - To step out of the Mapping Fonction MapFunc = J / (Rv * Rh * M)
                    Incr = MapFuncMin / 10                  'MARTIN
                    Dir = -1                                'MARTIN
                    MapFuncLast = MapFunc                   'MARTIN
                    MapFunc = MapFuncMin + Dir * Incr       'MARTIN
                Else                                        'MARTIN
                    MapFuncLast = MapFunc                   'MARTIN
                    MapFunc = J / (Rv * Rh * M)             '******** PRIMARY Mapping Function *********
                    If Cnt > 0 And PpMin >= Pp Then         'MARTIN
                        MapFuncMin = MapFunc                'MARTIN
                        PpMin = Pp                          'MARTIN
                    End If                                  'MARTIN
                End If                                      'MARTIN
            Else                                            'MARTIN
                If Pp > PpLast Then                         'MARTIN
                    Dir = Dir * -1                          'MARTIN
                    Incr = Application.Max(0.01, Incr / 10) 'MARTIN
                End If                                      'MARTIN
                    If Round(MapFunc, 4) = 0 Then           'MARTIN
                        Incr = Application.Max(0.01, Incr / 10)     'MARTIN
                        Dir = Dir * -1                      'MARTIN
                    End If                                  'MARTIN
                MapFuncLast = MapFunc                       'MARTIN
                MapFunc = MapFuncLast + Incr * Dir          'MARTIN
            End If                                          'MARTIN
            
            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
        
        If Stp2 = True Then                 'MARTIN Stp2 is true if VBA does not converge, Cnt = CntMax
            BoltCoefficientTest = 0.01      'MARTIN This is to return 0.01 if VBA does not converge, if Cnt = CntMax
        Else                                'MARTIN
            BoltCoefficientTest = P         'MARTIN
        End If                              'MARTIN
        If Magic = True Then                'MARTIN - To Delete - To report Magic
            If CntTemp > 25 Then            'MARTIN - To Delete - To report Magic2
                Cnt = Round(Cnt, 0) + 0.2   'MARTIN - To Delete - To report Magic2
            Else                            'MARTIN - To Delete - To report Magic
                Cnt = Round(Cnt, 0) + 0.1   'MARTIN - To Delete - To report Magic1
            End If                          'MARTIN - To Delete - To report Magic
        End If                              'MARTIN - To Delete - To report Magic
        
        Exit Function
ForcedExit:
    BoltCoefficientTest = Rv * Rh

End Function

Note:
Please refer to related post, Analysis of Bolts: Modified Visual Basic Code for Calculating Bolt Coefficient Using a Method Developed by Donald Brandt, for the original version of the above VBA Code.