NOTE: The slope-deflection sign convention may seem strange to those used to matirx stiffness analysis, but it makes sense. None of the slope deflection equations explicitly state a member 'direction' and it doesn't matter. For example, whether you consider the column AB as going from A to B or as going from B to A, a +ive shear at end A is still directed toward the left. In matrix analysis, that direction matters.
This solves a close approximation to the first-order example of Kulak and Grondin. The major differences are:
In [1]:
from IPython import display
display.SVG('KG-8.2sd.svg')
Out[1]:
In [2]:
import sympy as sy # use symbolic algebra
sy.init_printing() # print pretty math
from sdutil2 import SD, FEF # slope-deflection library
In [3]:
sy.var('theta_a theta_b theta_c theta_d theta_e theta_f Delta_b Delta_c')
Out[3]:
In [4]:
E = 200000
Ic = 222E6 # W310x97
Ib = 488E6 # W460x106
Hf = 21900 # Horizontal load at F
He = 43400 # Horizontal load at E
Lab = 6500 # Length of column
Lbc = 5500 # Length of column
Lbe = 10500 # Length of beam
In [5]:
Mab,Mba,Vab,Vba = SD(Lab,E*Ic,theta_a,theta_b,Delta_b) # column AB, BC
Mbc,Mcb,Vbc,Vcb = SD(Lbc,E*Ic,theta_b,theta_c,Delta_c-Delta_b)
Mde,Med,Vde,Ved = SD(Lab,E*Ic,theta_d,theta_e,Delta_b) # column DE, EF
Mef,Mfe,Vef,Vfe = SD(Lbc,E*Ic,theta_e,theta_f,Delta_c-Delta_b)
Mbe,Meb,Vbe,Veb = SD(Lbe,E*Ib,theta_b,theta_e) + FEF.udl(Lbe,55) # beams BE, CF
Mcf,Mfc,Vcf,Vfc = SD(Lbe,E*Ib,theta_c,theta_f) + FEF.udl(Lbe,45)
In [6]:
eqns = [ Mba+Mbe+Mbc, # sum of moments at B = 0
Mcb+Mcf, # sum of moments at C = 0
Med+Meb+Mef, # sum of moments at E = 0
Mfe+Mfc, # sum of moments at F = 0
-Vab-Vde+Hf+He, # sum of Fx @ base of storey 1 = 0
-Vbc-Vef+Hf, # sum of Fx @ base of storey 2 = 0
theta_a, # fixed support at A, rotation = 0
theta_d] # fixed support at F, rotation = 0
In [7]:
soln = sy.solve(eqns)
soln
Out[7]:
In [8]:
Mab
Out[8]:
In [9]:
Mab.subs(soln)
Out[9]:
In [10]:
Mab.subs(soln).n(4) * 1E-6
Out[10]:
In [11]:
V = globals() # another way to access global variables
V['Mab']
Out[11]:
In [12]:
Mab is V['Mab']
Out[12]:
In [13]:
V['Mab'].subs(soln).n()
Out[13]:
In [14]:
# collect the end moments in all 6 members
allm = []
V = globals()
for m in 'ab,bc,de,ef,be,cf'.split(','):
mj = V['M'+m].subs(soln).n()*1E-6 # Mxy
mk = V['M'+m[::-1]].subs(soln).n()*1E-6 # Myx
allm.append((m.upper(),mj,mk))
allm
Out[14]:
In [15]:
[(m,round(a,1),round(b,1)) for m,a,b in allm] # display to one decimal place
Out[15]:
In [16]:
# collect the end shears in all 6 members
allv = []
V = globals()
for m in 'ab,bc,de,ef,be,cf'.split(','):
mj = V['V'+m].subs(soln).n()*1E-3 # Mxy
mk = V['V'+m[::-1]].subs(soln).n()*1E-3 # Myx
allv.append((m.upper(),mj,mk))
allv
Out[16]:
In [17]:
[(m,round(a,1),round(b,1)) for m,a,b in allv] # display to one decimal place
Out[17]:
In [18]:
import pandas as pd
dd = '../matrix-methods/frame2d/data/KG82sd.d/all' # location of mm solution
In [19]:
mems = pd.DataFrame(allm,columns=['ID','MZJ','MZK']).set_index('ID')
mems
Out[19]:
In [20]:
mefs = pd.read_csv(dd+'/member_end_forces.csv').set_index('MEMBERID').loc[mems.index]
mems2 = mefs[['MZJ','MZK']] * -1E-6 # convert sign and to kN-m
mems2
Out[20]:
In [21]:
mdiff = (100*(1-mems/mems2)) # calculate % diff
mdiff
Out[21]:
In [22]:
mdiff.abs().max()
Out[22]:
The maximum difference in member end moments is 0.000009% (about 7 or 8 sig figs).
In [23]:
mevs = pd.DataFrame(allv,columns=['ID','FYJ','FYK']).set_index('ID')
mevs
Out[23]:
In [24]:
mevs2 = mefs[['FYJ','FYK']] * 1E-3
mevs2[['FYK']] *= -1 # change sign on end k
mevs2
Out[24]:
In [25]:
vdiff = 100*(1-mevs/mevs2)
vdiff
Out[25]:
In [26]:
vdiff.abs().max()
Out[26]:
The maximum difference is about 0.00004% which is high, but end shears on column AB are very small.
In [27]:
deltw = pd.DataFrame([('B', soln[Delta_b], soln[theta_b]),
('C', soln[Delta_c], soln[theta_c])],columns=['ID','DX','RZ']).set_index('ID')
deltw
Out[27]:
In [28]:
disp = pd.read_csv(dd+'/node_displacements.csv').set_index('NODEID').loc[deltw.index][['DX','RZ']]
disp['RZ'] *= -1
disp
Out[28]:
In [29]:
diffd = (100*(1-deltw/disp))
diffd
Out[29]:
In [30]:
diffd.abs().max()
Out[30]:
Max difference in displacement is 0.000006%.
Note that the matrix method solution was accomplished by setting A very high in order to minimize effects of axial deformation. But there is a limit to how high this can be set, due to numerical instability in equation solving (probably). Setting $A=10^{10}$ seems to be about as high as is possible - larger values lead to worse results.
In [ ]: