2D Column - Dynamic EQ Ground MotionΒΆ
Converted to openseespy by: Pavan Chigullapally
University of Auckland
Email: pchi893@aucklanduni.ac.nz
EQ ground motion with gravity- uniform excitation of structure
All units are in kip, inch, second
Note: In this example, all input values for Example 1a are replaced by variables. The objective of this example is to demonstrate the use of variables in defining
The OpenSees input and also to run various tests and algorithms at once to increase the chances of convergence
To run EQ ground-motion analysis (
BM68elc.acc
needs to be downloaded into the same directory)The detailed problem description can be found here (example:2a)
The source code is shown below, which can be downloaded
here
.
1# -*- coding: utf-8 -*-
2"""
3Created on Mon Apr 22 15:12:06 2019
4
5@author: pchi893
6"""
7# Converted to openseespy by: Pavan Chigullapally
8# University of Auckland
9# Email: pchi893@aucklanduni.ac.nz
10
11# EQ ground motion with gravity- uniform excitation of structure
12# all units are in kip, inch, second
13##Note: In this example, all input values for Example 1a are replaced by variables. The objective of this example is to demonstrate the use of variables in defining
14#the OpenSees input and also to run various tests and algorithms at once to increase the chances of convergence
15# Example 2a. 2D cantilever column, dynamic eq ground motion
16#To run EQ ground-motion analysis (BM68elc.acc needs to be downloaded into the same directory)
17#the detailed problem description can be found here: http://opensees.berkeley.edu/wiki/index.php/Examples_Manual (example:2a)
18# --------------------------------------------------------------------------------------------------
19# OpenSees (Tcl) code by: Silvia Mazzoni & Frank McKenna, 2006
20
21#
22# ^Y
23# |
24# 2 __
25# | |
26# | |
27# | |
28# (1) LCol
29# | |
30# | |
31# | |
32# =1= _|_ -------->X
33#
34
35# SET UP ----------------------------------------------------------------------------
36import openseespy.opensees as op
37#import the os module
38import os
39import math
40op.wipe()
41
42#########################################################################################################################################################################
43
44#########################################################################################################################################################################
45op.model('basic', '-ndm', 2, '-ndf', 3)
46
47#to create a directory at specified path with name "Data"
48
49
50#this will create the directory with name 'Data' and will update it when we rerun the analysis, otherwise we have to keep deleting the old 'Data' Folder
51dir = "Data-2a"
52if not os.path.exists(dir):
53 os.makedirs(dir)
54
55#this will create just 'Data' folder
56#os.mkdir("Data")
57
58#detect the current working directory
59#path1 = os.getcwd()
60#print(path1)
61
62LCol = 432.0 # column length
63Weight = 2000.0 # superstructure weight
64
65# define section geometry
66HCol = 60.0 # Column Depth
67BCol = 60.0 # Column Width
68
69PCol =Weight # nodal dead-load weight per column
70g = 386.4
71Mass = PCol/g
72
73ACol = HCol*BCol*1000 # cross-sectional area, make stiff
74IzCol = (BCol*math.pow(HCol,3))/12 # Column moment of inertia
75
76op.node(1, 0.0, 0.0)
77op.node(2, 0.0, LCol)
78
79op.fix(1, 1, 1, 1)
80
81op.mass(2, Mass, 1e-9, 0.0)
82
83ColTransfTag = 1
84op.geomTransf('Linear', ColTransfTag)
85#A = 3600000000.0
86#E = 4227.0
87#Iz = 1080000.0
88
89fc = -4.0 # CONCRETE Compressive Strength (+Tension, -Compression)
90Ec = 57*math.sqrt(-fc*1000) # Concrete Elastic Modulus (the term in sqr root needs to be in psi
91
92op.element('elasticBeamColumn', 1, 1, 2, ACol, Ec, IzCol, ColTransfTag)
93
94op.recorder('Node', '-file', 'Data-2a/DFree.out','-time', '-node', 2, '-dof', 1,2,3, 'disp')
95op.recorder('Node', '-file', 'Data-2a/DBase.out','-time', '-node', 1, '-dof', 1,2,3, 'disp')
96op.recorder('Node', '-file', 'Data-2a/RBase.out','-time', '-node', 1, '-dof', 1,2,3, 'reaction')
97#op.recorder('Drift', '-file', 'Data-2a/Drift.out','-time', '-node', 1, '-dof', 1,2,3, 'disp')
98op.recorder('Element', '-file', 'Data-2a/FCol.out','-time', '-ele', 1, 'globalForce')
99op.recorder('Element', '-file', 'Data-2a/DCol.out','-time', '-ele', 1, 'deformations')
100
101#defining gravity loads
102op.timeSeries('Linear', 1)
103op.pattern('Plain', 1, 1)
104op.load(2, 0.0, -PCol, 0.0)
105
106Tol = 1e-8 # convergence tolerance for test
107NstepGravity = 10
108DGravity = 1/NstepGravity
109op.integrator('LoadControl', DGravity) # determine the next time step for an analysis
110op.numberer('Plain') # renumber dof's to minimize band-width (optimization), if you want to
111op.system('BandGeneral') # how to store and solve the system of equations in the analysis
112op.constraints('Plain') # how it handles boundary conditions
113op.test('NormDispIncr', Tol, 6) # determine if convergence has been achieved at the end of an iteration step
114op.algorithm('Newton') # use Newton's solution algorithm: updates tangent stiffness at every iteration
115op.analysis('Static') # define type of analysis static or transient
116op.analyze(NstepGravity) # apply gravity
117
118op.loadConst('-time', 0.0) #maintain constant gravity loads and reset time to zero
119
120#applying Dynamic Ground motion analysis
121GMdirection = 1
122GMfile = 'BM68elc.acc'
123GMfact = 1.0
124
125
126
127Lambda = op.eigen('-fullGenLapack', 1) # eigenvalue mode 1
128import math
129Omega = math.pow(Lambda[0], 0.5)
130betaKcomm = 2 * (0.02/Omega)
131
132xDamp = 0.02 # 2% damping ratio
133alphaM = 0.0 # M-prop. damping; D = alphaM*M
134betaKcurr = 0.0 # K-proportional damping; +beatKcurr*KCurrent
135betaKinit = 0.0 # initial-stiffness proportional damping +beatKinit*Kini
136
137op.rayleigh(alphaM,betaKcurr, betaKinit, betaKcomm) # RAYLEIGH damping
138
139# Uniform EXCITATION: acceleration input
140IDloadTag = 400 # load tag
141dt = 0.01 # time step for input ground motion
142GMfatt = 1.0 # data in input file is in g Unifts -- ACCELERATION TH
143maxNumIter = 10
144op.timeSeries('Path', 2, '-dt', dt, '-filePath', GMfile, '-factor', GMfact)
145op.pattern('UniformExcitation', IDloadTag, GMdirection, '-accel', 2)
146
147op.wipeAnalysis()
148op.constraints('Transformation')
149op.numberer('Plain')
150op.system('BandGeneral')
151op.test('EnergyIncr', Tol, maxNumIter)
152op.algorithm('ModifiedNewton')
153
154NewmarkGamma = 0.5
155NewmarkBeta = 0.25
156op.integrator('Newmark', NewmarkGamma, NewmarkBeta)
157op.analysis('Transient')
158
159DtAnalysis = 0.01
160TmaxAnalysis = 10.0
161
162Nsteps = int(TmaxAnalysis/ DtAnalysis)
163
164ok = op.analyze(Nsteps, DtAnalysis)
165tCurrent = op.getTime()
166
167# for gravity analysis, load control is fine, 0.1 is the load factor increment (http://opensees.berkeley.edu/wiki/index.php/Load_Control)
168
169test = {1:'NormDispIncr', 2: 'RelativeEnergyIncr', 4: 'RelativeNormUnbalance',5: 'RelativeNormDispIncr', 6: 'NormUnbalance'}
170algorithm = {1:'KrylovNewton', 2: 'SecantNewton' , 4: 'RaphsonNewton',5: 'PeriodicNewton', 6: 'BFGS', 7: 'Broyden', 8: 'NewtonLineSearch'}
171
172for i in test:
173 for j in algorithm:
174
175 if ok != 0:
176 if j < 4:
177 op.algorithm(algorithm[j], '-initial')
178
179 else:
180 op.algorithm(algorithm[j])
181
182 op.test(test[i], Tol, 1000)
183 ok = op.analyze(Nsteps, DtAnalysis)
184 print(test[i], algorithm[j], ok)
185 if ok == 0:
186 break
187 else:
188 continue
189
190u2 = op.nodeDisp(2, 1)
191print("u2 = ", u2)
192
193op.wipe()