I am trying to understand how to implement rigid diaphragm in an opensees model. Yet I am having troubles in the following example. Model seems to be working without rigid diaphragm, I verified the analysis results with another software. Could you please help me with the problem? I know that constraints handler should be Lagrange in case of rigid diaphragm, yet the issue persists.
Thank you
Code: Select all
import openseespy.opensees as op
from openseespy.postprocessing.Get_Rendering import *
# REFERENCES:
# used in verification by SAP2000 and SeismoStruct:
# SAP2000 Integrated Finite Element Analysis and Design of Structures, Verification Manual,
# Computers and Structures, 2009. Example 1-024.
# SeismoStruct, Verification Report, 2020. Example 12.
# Basic Units
m = 1.0
kN = 1.0
sec = 1.0
LunitTXT = 'meter'
FunitTXT = 'kN'
TunitTXT = 'sec'
# Constants
pi = np.pi
g = 9.81*m/sec**2
# Length
mm = m/1000.0
cm = m/100.0
inch = 25.4*mm
ft = 12.0*inch
# Area
m2 = m**2
cm2 = cm**2
mm2 = mm**2
inch2 = inch**2
ft2 = ft**2
# First Moment of Area
m3 = m**3
cm3 = cm**3
mm3 = mm**3
inch3 = inch**3
ft3 = ft**3
# Second Moment of Area
m4 = m**4
cm4 = cm**4
mm4 = mm**4
inch4 = inch**4
ft4 = ft**4
# Force
N = kN/1000.0
kip = kN*4.448221615
# Moment
kNm = kN*m
# Stress (kN/m2 or kPa)
Pa = N/(m2)
kPa = Pa*1.0e3
MPa = Pa*1.0e6
GPa = Pa*1.0e9
ksi = 6.8947573*MPa
psi = 1e-3*ksi
# Angles
degrees = pi/180.0
op.wipe()
op.model('basic', '-ndm', 3, '-ndf', 6)
# Frame grid
Xs = [0,35*ft,70*ft]
Ys = [0,25*ft,50*ft]
Zs = [0,13*ft,26*ft]
# Center of mass at each floor
Xcm = 38*ft
Ycm = 27*ft
# Lumped floor masses
massX = 6.2112*kip*sec**2/ft
massY = 6.2112*kip*sec**2/ft
# Distributed Loading
wload = -10*kip/ft
# Beam section properties
v = 0.2
Eb = 500000*kip/ft2
Gb = Eb/(2*(1+v))
Ab = 5*ft2
Iyb = 2.61*ft4
Izb = 1.67*ft4
Jb = 0
# Column section properties
v = 0.2
Ec = 350000*kip/ft2
Gc = Ec/(2*(1+v))
Ac = 4*ft2
Izc = 1.25*ft4
Iyc = 1.25*ft4
Jc = 0
# Define Transformation Tags
ColTransf = 1
BeamXTransf = 2
BeamYTransf = 3
op.geomTransf('Linear', ColTransf, 0, 1, 0)
op.geomTransf('Linear', BeamXTransf, 0, 0, 1)
op.geomTransf('Linear', BeamYTransf, 0, 0, 1)
# define NODAL COORDINATES
storey = 0
for k in range(len(Zs)):
no = 1
constrained = []
for i in range(len(Xs)):
for j in range(len(Ys)):
nodeID = int(str(no)+'00'+str(storey))
op.node(nodeID,Xs[i],Ys[j],Zs[k])
if k == 0:
op.fix(nodeID,1,1,1,1,1,1)
else:
constrained.append(nodeID)
no += 1
if k != 0: # Center of mass
nodeID = int(str(no)+'00'+str(storey))
op.node(nodeID,Xcm,Ycm,Zs[k])
op.mass(nodeID,massX,massY,0,0,0,0)
op.rigidDiaphragm(3,nodeID,*constrained)
storey += 1
# Define Columns
colTag = '00'
for no in range(1,(len(Xs))*(len(Ys))+1):
for storey in range(1,len(Zs)):
nodeI = int(str(no)+'00'+str(storey-1))
nodeJ = int(str(no)+'00'+str(storey))
eleTag = int(str(no)+colTag+str(storey))
op.element('elasticBeamColumn',eleTag,nodeI,nodeJ,Ac, Ec, Gc, Jc, Iyc, Izc, ColTransf)
beamEles = []
# Define Beams in X axis
beamXtag = '01'
for storey in range(1,len(Zs)):
no = 1
for i in range(len(Xs)-1):
for j in range(1,len(Ys)+1):
nodeI = int(str(i*len(Ys)+j)+'00'+str(storey))
nodeJ = int(str((i+1)*len(Ys)+j)+'00'+str(storey))
eleTag = int(str(no)+beamXtag+str(storey))
beamEles.append(eleTag)
op.element('elasticBeamColumn',eleTag,nodeI,nodeJ,Ab, Eb, Gb, Jb, Iyb, Izb, BeamXTransf)
no+=1
# Define Beams in Y axis
beamYtag = '02'
for storey in range(1,len(Zs)):
no = 1
for i in range(len(Xs)):
for j in range(1,len(Ys)):
nodeI = int(str(i*len(Ys)+j)+'00'+str(storey))
nodeJ = int(str(i*len(Ys)+(j+1))+'00'+str(storey))
eleTag = int(str(no)+beamYtag+str(storey))
beamEles.append(eleTag)
op.element('elasticBeamColumn',eleTag,nodeI,nodeJ,Ab, Eb, Gb, Jb, Iyb, Izb, BeamYTransf)
no+=1
plot_model('nodes','elements')
# ----------------------------------------------------------------------------
# Gravity Analysis
# ----------------------------------------------------------------------------
# APPLY GRAVITY LOADING
# create TimeSeries
op.timeSeries("Constant", 1)
# create a plain load pattern
op.pattern('Plain', 1, 1)
for ele in beamEles:
op.eleLoad('-ele', ele, '-type', '-beamUniform', 0,wload,0)
# SET ANALYSIS PARAMETERS
op.wipeAnalysis()
op.constraints('Lagrange')
op.numberer('RCM')
op.system('BandGeneral')
op.test('EnergyIncr', 1e-8, 6)
op.algorithm('Newton')
nG = 100
op.integrator('LoadControl', 1/nG)
op.analysis('Static')
# DO THE ANALYSIS
op.analyze(nG)
# maintain constant gravity loads and reset time to zero
op.loadConst('-time', 0.0)
#ì Print reactions
op.reactions('-static')
reactions = 0
no = 1
storey = 0
for i in range(len(Xs)):
for j in range(len(Ys)):
nodeID = int(str(no)+'00'+str(storey))
no+=1
reactions += np.array(op.nodeReaction(nodeID))
print(nodeID,op.nodeReaction(nodeID))
print('Tot ',reactions)
# plot_modeshape(1)