Example Eigenvalue analysis of a three-dimensional frame

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

Moderators: silvia, selimgunay, Moderators

Post Reply
Prafullamalla
Posts: 160
Joined: Mon Feb 02, 2015 6:32 pm

Example Eigenvalue analysis of a three-dimensional frame

Post by Prafullamalla » Tue Feb 14, 2023 8:29 pm

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")     
[/code]
Prafulla Malla, Nepal
Praf_malla@hotmail.com

Post Reply