2D Column - Dynamic EQ Ground MotionΒΆ

Converted to openseespy by: Pavan Chigullapally
                      University of Auckland
                      Email: pchi893@aucklanduni.ac.nz
  1. EQ ground motion with gravity- uniform excitation of structure

  2. All units are in kip, inch, second

  3. 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

  4. The OpenSees input and also to run various tests and algorithms at once to increase the chances of convergence

  5. To run EQ ground-motion analysis (BM68elc.acc needs to be downloaded into the same directory)

  6. The detailed problem description can be found here (example:2a)

  7. 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()