Weird Deformed Shape

[Feras]

Hello

I am modeling the famous SPEAR frame in openseespy. I am running gravity analysis for the first story only and I am getting deformed shape totally different from what I am getting in another software for elements 13, 22, 18, 23. However, axial loads and bending moments almost match for both software. I have tried to check everything that is possible to cause this discrepancy, and I am still getting the same deformed shape in OpenSees.
I am posting script here so could someone help me find what is causing this wired shape.

Thanks in advance.



import openseespy.opensees as ops
import numpy as np
import matplotlib.pyplot as plt
import openseespy.postprocessing.Get_Rendering as opsplt
import openseespy.postprocessing.ops_vis as opsv

################################################################################
# ------------------------------Define units------------------------------------
################################################################################
#----------------------Basic Units m, KN, sec, ton----------------------
m= 1.
KN=1.
ton=1.
sec=1.

mm=0.001*m
N=KN*0.001

KPa= KN/m**2
MPa= KPa*1000
GPa=1000*MPa

g = 9.81*m/sec**2
PI= 22/7

GammaConcrete= 24*(KN/m**3)

ops.wipe()
ops.model('basic','-ndm',3,'-ndf',6)


crdZ=np.array([0.0, 3.0])*m
#----------ground level nodes-------------

ops.node(11, 0.0*m,-6.0*m, crdZ[0])
ops.node(12, -5.0*m,-6.0*m, crdZ[0])
ops.node(6, 0.0*m, 0.0*m, crdZ[0])
ops.node(7, -6.0*m, 0.0*m, crdZ[0])
ops.node(10, 3.0*m, -6.0*m, crdZ[0])
ops.node(3, 0.0*m, 4.25*m, crdZ[0])
ops.node(5, -6.0*m, 4.0*m, crdZ[0])
ops.node(1, 3.0*m, 4.5*m, crdZ[0])
ops.node(8, 3.0*m,- 0.5*m, crdZ[0])
#----------First level nodes-------------
ops.node(24, 0.0*m, -6.0*m, crdZ[1])
ops.node(25, -5.0*m,-6.0*m, crdZ[1])
ops.node(18, 0.0*m, 0.0*m, crdZ[1])
ops.node(20, -6.0*m, 0.0*m, crdZ[1])
ops.node(23, 3.0*m, -6.0*m, crdZ[1])
ops.node(15, 0.0*m, 4.25*m, crdZ[1])
ops.node(17, -6.0*m, 4.0*m, crdZ[1])
ops.node(13, 3.0*m, 4.5*m, crdZ[1])
ops.node(21, 3.0*m, -0.5*m, crdZ[1])
ops.node(14, 0.0*m, 4.5*m, crdZ[1])
ops.node(16, 0.0*m, 4.0*m, crdZ[1])
ops.node(19, -5.0*m, 0.0*m, crdZ[1])
ops.node(22, 0.0*m, -0.5*m, crdZ[1])

################################################################################
#------------------------------------Boundary Conditions-----------------------------
################################################################################

ops.fix(1, 1, 1, 1, 1, 1, 1)
ops.fix(3, 1, 1, 1, 1, 1, 1)
ops.fix(5, 1, 1, 1, 1, 1, 1)
ops.fix(6, 1, 1, 1, 1, 1, 1)
ops.fix(7, 1, 1, 1, 1, 1, 1)
ops.fix(8, 1, 1, 1, 1, 1, 1)
ops.fix(10, 1, 1, 1, 1, 1, 1)
ops.fix(11, 1, 1, 1, 1, 1, 1)
ops.fix(12, 1, 1, 1, 1, 1, 1)


###############################################################################
# --------------- Rigid Diaphragm Definition----------------------
###############################################################################
ops.node(100, -1.58*m, -0.85*m, crdZ[1]) # master node at mass center of level 1
perpDirn= 3

ops.fix(100, 0, 0, 1, 1, 1, 0) # free translation in X and Y and rotation about Z
#rigidDiaphragm(perpDirn, rNodeTag, *cNodeTags
ops.rigidDiaphragm(perpDirn,100 ,*[13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25])

################################################################################
# ---------------------------- Elastic test sections---------------------------
################################################################################
elsSecTag= np.array([400, 500, 600])
fc = -26.5*MPa
Ec= 4700*(np.sqrt(-fc/MPa))*MPa
v=0.3 # Poisson ratio
#Input parameters and elastic section for sq columns
colWidth = 250*mm # Width along local z direction
colDepth = 250*mm # Depth along local y direction
Iz_sq= (colWidth*colWidth**3/12)
Iy_sq= Iz_sq
G_sq= ((0.5*Ec/MPa)/(1+v))*MPa
J_sq= 2*Iz_sq
E_Sq= 4700*(np.sqrt(-fc/MPa))*MPa
A_Sq= colWidth*colWidth
#section('Elastic', secTag, E_mod, A, Iz, Iy, G_mod, Jxx)
ops.section('Elastic',int(elsSecTag[0]), E_Sq, A_Sq, Iz_sq, Iy_sq, G_sq, J_sq) #elastic section for square columns
#Input parameters and elastic section for C6 columns
colC6Width = 250*mm
colC6Depth = 750*mm
Iz_C6= (colC6Width*colC6Depth**3/12)
Iy_C6= (colC6Depth*colC6Width**3/12)
G_C6 = ((0.5*Ec/MPa)/(1+v))*MPa
J_C6= 1/12*colC6Width*colC6Depth*(colC6Width**2+colC6Depth**2)
E_C6= 4700*(np.sqrt(-fc/MPa))*MPa
A_C6= colC6Width*colC6Depth
ops.section('Elastic', int(elsSecTag[1]), E_C6, A_C6, Iz_C6, Iy_C6, G_C6, J_C6) # elastic section for C6 columns
#Input parameters and elastic section for Beams
beamWidth = 250*mm
beamDepth = 500*mm
Iz_b= (beamWidth*beamDepth**3/12)
Iy_b= (beamDepth*beamWidth**3/12)
G_b= ((0.5*Ec/MPa)/(1+v))*MPa
J_b= 1/12 *beamWidth* beamDepth*(beamDepth**2+beamWidth**2)
E_b= 4700*(np.sqrt(-fc/MPa))*MPa
A_b= beamWidth*beamDepth
ops.section('Elastic', int(elsSecTag[2]), E_b, A_b, Iz_b, Iy_b, G_b, J_b) # elastic section for beams

################################################################################
# -------------------------------Define Elements--------------------------------
################################################################################


transfTag= np.array([1, 2]) # geomTransf tags 1=columns, 2= beams
# geomTransf for columns
vecxzCol = np.array ([1.0, 0.0, 0.0])
ops.geomTransf('Linear', int(transfTag[0]), *vecxzCol ) #local z axis is aligned with global X axis

intTag = np.array([1, 2, 3]) # beamIntegration tags

# #beamIntegration('Lobatto', tag, secTag, N)
ops.beamIntegration('Lobatto',int(intTag[0]), int(elsSecTag[0]), 3) # beamIntegration for square columns
ops.beamIntegration('Lobatto',int(intTag[1]), int(elsSecTag[1]), 3) #beamIntegration for C6 column

################################################################################
#----------------------------Generating columns for ground stories--------------------------
################################################################################
#element('forceBeamColumn', eleTag, *eleNodes, transfTag, integrationTag, '-mass', mass=0.0)
ops.element('forceBeamColumn', 8, *[11, 24], int(transfTag[0]), int(intTag[0]) )
ops.element('forceBeamColumn', 6, *[12, 25], int(transfTag[0]), int(intTag[0]))
ops.element('forceBeamColumn', 7, *[6, 18], int(transfTag[0]), int(intTag[0]) )
ops.element('forceBeamColumn', 9, *[7, 20], int(transfTag[0]), int(intTag[0]))
ops.element('forceBeamColumn', 3, *[10, 23], int(transfTag[0]), int(intTag[0]))
ops.element('forceBeamColumn', 4, *[3, 15], int(transfTag[0]), int(intTag[1]) )
ops.element('forceBeamColumn', 5, *[5, 17], int(transfTag[0]),int(intTag[0]) )
ops.element('forceBeamColumn', 1, *[1, 13], int(transfTag[0]), int(intTag[0]) )
ops.element('forceBeamColumn', 2, *[8, 21], int(transfTag[0]), int(intTag[0]))

################################################################################
#----------------------------Generating Beam Elements--------------------------
################################################################################

ops.beamIntegration('Lobatto',int(intTag[2]), int(elsSecTag[2]), 3)

#geomTransf for beams
vecxzBeamx = np.array ([0.0, 0.0, 1.0])
ops.geomTransf('Linear', int(transfTag[1]), *vecxzBeamx ) # Local z is aligned with global Z

ops.element('forceBeamColumn', 10, *[23, 24], int(transfTag[1]), int(intTag[2]) )
ops.element('forceBeamColumn', 11, *[24, 25], int(transfTag[1]), int(intTag[2]))
ops.element('forceBeamColumn', 12, *[21, 22], int(transfTag[1]), int(intTag[2]))
ops.element('forceBeamColumn', 13, *[18, 19], int(transfTag[1]), int(intTag[2]))
ops.element('forceBeamColumn', 14, *[13, 14], int(transfTag[1]), int(intTag[2]) )
ops.element('forceBeamColumn', 15, *[16, 17],int(transfTag[1]), int(intTag[2]) )
ops.element('forceBeamColumn', 16, *[19, 25], int(transfTag[1]),int(intTag[2]) )
ops.element('forceBeamColumn', 17, *[17, 20], int(transfTag[1]), int(intTag[2]) )
ops.element('forceBeamColumn', 18, *[22, 24], int(transfTag[1]), int(intTag[2]) )
ops.element('forceBeamColumn', 19, *[15, 18], int(transfTag[1]), int(intTag[2]))
ops.element('forceBeamColumn', 20, *[21, 23], int(transfTag[1]), int(intTag[2]) )
ops.element('forceBeamColumn', 21, *[13, 21], int(transfTag[1]), int(intTag[2]) )
ops.element('forceBeamColumn', 22, *[19, 20], int(transfTag[1]), int(intTag[2]) )
ops.element('forceBeamColumn', 23, *[18, 22], int(transfTag[1]), int(intTag[2]) )

##################################################################################
# ----------------------- Add Rigid Link at C6 Column-----------------------------
##################################################################################
E_mod= E_b
G_mod= 0.4*E_mod
Iy= (0.5*0.25**3)/12
Iz= (0.25*0.5**3)/12
Jxx= Iy+Iz
A= 0.25*0.5
#element('elasticBeamColumn', eleTag, *eleNodes, Area, E_mod, G_mod, Jxx, Iy, Iz, transfTag)
ops.element('elasticBeamColumn',25 ,*[14,15], A, E_mod*1000, G_mod, Jxx, Iy, Iz,int(transfTag[1]))
ops.element('elasticBeamColumn',26 ,*[15,16], A, E_mod*1000, G_mod, Jxx, Iy, Iz,int(transfTag[1]))

################################################################################
#---------------Gravity Loads & Analysis----------------------------
################################################################################
opsplt.plot_model('nodes')
opsplt.createODB("3DFrame", "Gravity")
wY= np.array([0.00099834, 0.00106014, 0.00161634, 0.00155454, 0.00093654, 0.00081295, 0.00158639, 0.00132772, 0.00130734, 0.00143093, 0.00068935, 0.00068935, 0.00081295, 0.00130734])*g*1000
ops.timeSeries('Constant',1)
ops.pattern('Plain',1,1)
#('-ele', *eleTags,'-type', '-beamUniform', Wy, <Wz>, Wx=0.0
ops.eleLoad('-ele', 10,'-type','-beamUniform',0.0, -wY[0], 0.0)
ops.eleLoad('-ele', 11,'-type','-beamUniform',0.0, -wY[1], 0.0)
ops.eleLoad('-ele', 12,'-type','-beamUniform',0.0, -wY[2], 0.0)
ops.eleLoad('-ele', 13,'-type','-beamUniform',0.0, -wY[3], 0.0)
ops.eleLoad('-ele', 14,'-type','-beamUniform',0.0, -wY[4], 0.0)
ops.eleLoad('-ele', 15,'-type','-beamUniform',0.0, -wY[5], 0.0)
ops.eleLoad('-ele', 16,'-type','-beamUniform',0.0, -wY[6], 0.0)
ops.eleLoad('-ele', 17,'-type','-beamUniform',0.0, -wY[7], 0.0)
ops.eleLoad('-ele', 18,'-type','-beamUniform',0.0, -wY[8], 0.0)
ops.eleLoad('-ele', 19,'-type','-beamUniform',0.0, -wY[9], 0.0)
ops.eleLoad('-ele', 20,'-type','-beamUniform',0.0, -wY[10], 0.0)
ops.eleLoad('-ele', 21,'-type','-beamUniform',0.0, -wY[11], 0.0)
ops.eleLoad('-ele', 22,'-type','-beamUniform',0.0, -wY[12], 0.0)
ops.eleLoad('-ele', 23,'-type','-beamUniform',0.0, -wY[13], 0.0)

NstepGravity= 10
DGravity=1/NstepGravity
#integrator('LoadControl', incr)
ops.integrator('LoadControl',DGravity)
ops.system('BandGen')
ops.test('NormDispIncr', 1e-8, 10)
ops.numberer('RCM')
ops.constraints('Transformation')
ops.algorithm('Newton','-initial')
ops.analysis('Static')
ops.initialize()
ops.analyze(NstepGravity)

ops.wipe()

opsplt.plot_deformedshape(Model="3DFrame", LoadCase="Gravity", scale = 100.0)

[mhscott]

I think your beams are oriented for weak axis bending.

https://portwooddigital.com/2020/08/08/ ... x-z-plane/