Trouble with Rigid Diaphragm

Forum for asking and answering questions related to use of the OpenSeesPy module

Moderators: silvia, selimgunay, Moderators

Post Reply
volkanozs
Posts: 2
Joined: Wed Feb 22, 2017 3:47 pm
Location: IUSS

Trouble with Rigid Diaphragm

Post by volkanozs » Thu Oct 15, 2020 3:09 am

Dear All,

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)

volkanozs
Posts: 2
Joined: Wed Feb 22, 2017 3:47 pm
Location: IUSS

Re: Trouble with Rigid Diaphragm

Post by volkanozs » Thu Oct 15, 2020 3:37 am

All right, I resolved the issue.
Other dofs that don't belong to rigid diaphragm control must have be constrained for the retained nodes

Post Reply