Code: Select all
[code]# -*- coding: utf-8 -*-
"""
Created on Wed Feb 15 10:26:02 2023
@author: praf_malla@hotmail.com
"""
# Seismostruct Example
# EXAMPLE 12 – Eigenvalue analysis of a three-dimensional frame with rigid floor diaphragm
import openseespy.opensees as ops
import opsvis as opsv
from math import asin, sqrt
ops.wipe()
ops.model('basic','-ndm',3,'-ndf',6)
ops.node( 1 , 0 , 0 , 0 )
ops.node( 2 , 35 , 0 , 0 )
ops.node( 3 , 70 , 0 , 0 )
ops.node( 4 , 0 , 25 , 0 )
ops.node( 5 , 35 , 25 , 0 )
ops.node( 6 , 70 , 25 , 0 )
ops.node( 7 , 0 , 50 , 0 )
ops.node( 8 , 35 , 50 , 0 )
ops.node( 9 , 70 , 50 , 0 )
ops.node( 10 , 0 , 0 , 13 )
ops.node( 11 , 35 , 0 , 13 )
ops.node( 12 , 70 , 0 , 13 )
ops.node( 13 , 0 , 25 , 13 )
ops.node( 14 , 35 , 25 , 13 )
ops.node( 15 , 70 , 25 , 13 )
ops.node( 16 , 0 , 50 , 13 )
ops.node( 17 , 35 , 50 , 13 )
ops.node( 18 , 70 , 50 , 13 )
ops.node( 19 , 0 , 0 , 26 )
ops.node( 20 , 35 , 0 , 26 )
ops.node( 21 , 70 , 0 , 26 )
ops.node( 22 , 0 , 25 , 26 )
ops.node( 23 , 35 , 25 , 26 )
ops.node( 24 , 70 , 25 , 26 )
ops.node( 25 , 0 , 50 , 26 )
ops.node( 26 , 35 , 50 , 26 )
ops.node( 27 , 70 , 50 , 26 )
ops.fix( 1 ,1,1,1,1,1,1)
ops.fix( 2 ,1,1,1,1,1,1)
ops.fix( 3 ,1,1,1,1,1,1)
ops.fix( 4 ,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( 9 ,1,1,1,1,1,1)
# Column
Ac = 4
Ic = 1.25
Jc = 2*Ic
Ec = 350000
v = 0.3
Gc = 0.5*Ec/(1+v)
ops.section('Elastic',1,Ec,Ac,Ic,Ic,Gc,Jc)
ops.beamIntegration('Legendre',1,1,2)
#ops.geomTransf('Linear',1,1,0,0,'-jntOffset',0,0,0,0,0,-d)
ops.geomTransf('Linear',1,1,0,0)
ops.element('forceBeamColumn' , 1 , 1 , 10 ,1,1)
ops.element('forceBeamColumn' , 2 , 2 , 11 ,1,1)
ops.element('forceBeamColumn' , 3 , 3 , 12 ,1,1)
ops.element('forceBeamColumn' , 4 , 4 , 13 ,1,1)
ops.element('forceBeamColumn' , 5 , 5 , 14 ,1,1)
ops.element('forceBeamColumn' , 6 , 6 , 15 ,1,1)
ops.element('forceBeamColumn' , 7 , 7 , 16 ,1,1)
ops.element('forceBeamColumn' , 8 , 8 , 17 ,1,1)
ops.element('forceBeamColumn' , 9 , 9 , 18 ,1,1)
ops.element('forceBeamColumn' , 10 , 10 , 19 ,1,1)
ops.element('forceBeamColumn' , 11 , 11 , 20 ,1,1)
ops.element('forceBeamColumn' , 12 , 12 , 21 ,1,1)
ops.element('forceBeamColumn' , 13 , 13 , 22 ,1,1)
ops.element('forceBeamColumn' , 14 , 14 , 23 ,1,1)
ops.element('forceBeamColumn' , 15 , 15 , 24 ,1,1)
ops.element('forceBeamColumn' , 16 , 16 , 25 ,1,1)
ops.element('forceBeamColumn' , 17 , 17 , 26 ,1,1)
ops.element('forceBeamColumn' , 18, 18 , 27 ,1,1)
# Beam
Ab = 5
Ib = 2.61
Jb = 2*Ib
Eb = 500000
v = 0.3
Gb = 0.5*Ec/(1+v)
ops.section('Elastic',3,Eb,Ab,Ib,Ib,1000*Gb,Jb)
ops.beamIntegration('Legendre',3,3,2)
# Beam X direction
ops.geomTransf('Linear',3,0,1,0)
ops.element('forceBeamColumn' , 19 , 10 , 11 ,3,3)
ops.element('forceBeamColumn' , 20 , 11 , 12 ,3,3)
ops.element('forceBeamColumn' , 21 , 13 , 14 ,3,3)
ops.element('forceBeamColumn' , 22 , 14 , 15 ,3,3)
ops.element('forceBeamColumn' , 23 , 16 , 17 ,3,3)
ops.element('forceBeamColumn' , 24 , 17 , 18 ,3,3)
ops.element('forceBeamColumn' , 31 , 19 , 20 ,3,3)
ops.element('forceBeamColumn' , 32 , 20 , 21 ,3,3)
ops.element('forceBeamColumn' , 33 , 22 , 23 ,3,3)
ops.element('forceBeamColumn' , 34 , 23 , 24 ,3,3)
ops.element('forceBeamColumn' , 35 , 25 , 26 ,3,3)
ops.element('forceBeamColumn' , 36 , 26 , 27 ,3,3)
# Beam Y direction
ops.geomTransf('Linear',4,1,0,0)
ops.element('forceBeamColumn' , 25 , 10 , 13 ,4,3)
ops.element('forceBeamColumn' , 26 , 13 , 16 ,4,3)
ops.element('forceBeamColumn' , 27 , 11 , 14 ,4,3)
ops.element('forceBeamColumn' , 28 , 14 , 17 ,4,3)
ops.element('forceBeamColumn' , 29 , 12 , 15 ,4,3)
ops.element('forceBeamColumn' , 30 , 15 , 18 ,4,3)
ops.element('forceBeamColumn' , 37 , 19 , 22 ,4,3)
ops.element('forceBeamColumn' , 38 , 22 , 25 ,4,3)
ops.element('forceBeamColumn' , 39 , 20 , 23 ,4,3)
ops.element('forceBeamColumn' , 40 , 23 , 26 ,4,3)
ops.element('forceBeamColumn' , 41 , 21 , 24 ,4,3)
ops.element('forceBeamColumn' , 42 , 24 , 27 ,4,3)
ops.node( 28 , 38 , 27 , 13 ) #CM1
ops.node( 29 , 38 , 27 , 26 ) #CM2
ops.fix(28,0,0,1,1,1,0)
ops.fix(29,0,0,1,1,1,0)
ops.rigidDiaphragm( 3,28, 10,11,12,13,14,15,16,17,18 )
ops.rigidDiaphragm( 3,29, 19,20,21,22,23,24,25,26,27 )
# Display Model
m = 6.2112
L1=70
L2=50
mz =0 # m*L1*L2/6 # Rotational mass about vertical
ops.mass(28,6.2112,6.2112,0,0,0,mz)
ops.mass(29,6.2112,6.2112,0,0,0,mz)
# calculate eigenvalues & print results
numEigen = 4
eigenValues = ops.eigen('-fullGenLapack', numEigen)
print(eigenValues)
PI = 2 * asin(1.0)
# Display Model
opsv.plot_model()
# Display specific mode shape
for n in range(numEigen):
opsv.plot_mode_shape(n+1)
# Modal Participation factors
N = numEigen
ndf = 6 # Nodal degrees of freedom
omega2 = eigenValues
for n in range(N):
Mn = 0.0
Ln = 0.0
for nd in ops.getNodeTags():
ndMass = ops.nodeMass(nd)
ndEigen = ops.nodeEigenvector(nd,n+1)
Ln += ndEigen[0]*ndMass[0] # 0 for X, 1 for Y, 2 for Z excitation
for dof in range(ndf):
Mn += (ndEigen[dof]**2)*ndMass[dof]
Gamman = Ln/Mn
Tn = 2*3.14159/omega2[n]**0.5
print(f'Mode {n+1}, Tn = {Tn}, Mn = {Mn}, Gamma = {Gamman}')
# compute the modal properties
ops.modalProperties("-print", "-file", "ModalReport.txt", "-unorm")