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:
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!
Canam-Heavy, a division of Canam Group Inc.
P 450-641-4000; 3555 / M 514-702-4343
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