2D Portal Frame with Units- Uniform Dynamic EQ -bidirectional-acctimeseries

Converted to openseespy by: Pavan Chigullapally
                      University of Auckland
                      Email: pchi893@aucklanduni.ac.nz
  1. To run Uniaxial Inelastic Material, Fiber Section, Nonlinear Mode, Bidirectional-uniform earthquake ground motion

  2. First import the InelasticFiberSectionPortal2Dframe.py

  3. To run EQ ground-motion analysis (ReadRecord.py, H-E12140.AT2, H-E01140.AT2 needs to be downloaded into the same directory)

  4. Bidirectional-uniform support excitation using acceleration timeseries (different accelerations are input at all support nodes in two directions) – two support nodes here

  5. The problem description can be found here (example:4)

  6. 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# Example4. 2D Portal Frame--  Dynamic EQ input analysis-- Bidirectional-uniform support excitation using acceleration timeseries
 12
 13#To run Uniaxial Inelastic Material, Fiber Section, Nonlinear Mode, Bidirectional-uniform earthquake ground motion:First import the InelasticFiberSectionPortal2Dframe.py
 14#(upto gravity loading is already in this script) and run the current script
 15#To run EQ ground-motion analysis (ReadRecord.py, H-E12140.AT2 and H-E01140.AT2 needs to be downloaded into the same directory)
 16# Bidirectional-uniform support excitation using acceleration timeseries (different accelerations are input at all support nodes in two directions) -- two support nodes here
 17#the problem description can be found here: 
 18#http://opensees.berkeley.edu/wiki/index.php/Examples_Manual and http://opensees.berkeley.edu/wiki/index.php/OpenSees_Example_4._Portal_Frame(example: 4)
 19# --------------------------------------------------------------------------------------------------
 20#	OpenSees (Tcl) code by:	Silvia Mazzoni & Frank McKenna, 2006
 21# --------------------------------------------------------------------------------------------------
 22##########################################################################################################################################################################
 23import openseespy.opensees as op
 24#import the os module
 25#import os
 26import math
 27op.wipe()
 28##########################################################################################################################################################################
 29
 30from InelasticFiberSectionPortal2Dframe import *
 31# execute this file after you have built the model, and after you apply gravity
 32#
 33
 34# MultipleSupport Earthquake ground motion (different displacement input at spec'd support nodes) -- two nodes here
 35
 36#applying Dynamic Ground motion analysis
 37#iSupportNode = [1, 2]
 38iGMfact = [1.5, 0.25]
 39iGMdirection = [1, 2]
 40iGMfile = ['H-E01140', 'H-E12140']
 41DtAnalysis = 0.01*sec # time-step Dt for lateral analysis
 42TmaxAnalysis = 10.0*sec # maximum duration of ground-motion analysis
 43Tol = 1e-8
 44
 45# define DAMPING--------------------------------------------------------------------------------------
 46# apply Rayleigh DAMPING from $xDamp
 47# D=$alphaM*M + $betaKcurr*Kcurrent + $betaKcomm*KlastCommit + $beatKinit*$Kinitial
 48Lambda = op.eigen('-fullGenLapack', 1) # eigenvalue mode 1
 49Omega = math.pow(Lambda, 0.5)
 50betaKcomm = 2 * (0.02/Omega)
 51
 52xDamp = 0.02				# 2% damping ratio
 53alphaM = 0.0				# M-prop. damping; D = alphaM*M	
 54betaKcurr = 0.0		# K-proportional damping;      +beatKcurr*KCurrent
 55betaKinit = 0.0 # initial-stiffness proportional damping      +beatKinit*Kini
 56
 57op.rayleigh(alphaM,betaKcurr, betaKinit, betaKcomm) # RAYLEIGH damping
 58#--------------------------------------------------------------------------------------
 59#  ---------------------------------    perform Dynamic Ground-Motion Analysis
 60# the following commands are unique to the Multiple-Support Earthquake excitation
 61# Set some parameters
 62IDloadTag = 400			# load tag
 63
 64
 65# read a PEER strong motion database file, extracts dt from the header and converts the file 
 66# to the format OpenSees expects for Uniform/multiple-support ground motions
 67record = ['H-E01140', 'H-E12140']
 68#dt =[]
 69#nPts = []
 70
 71import ReadRecord
 72
 73#this is similar to uniform excitation in single direction
 74count = 2
 75for i in range(len(iGMdirection)):
 76    IDloadTag = IDloadTag+count
 77    record_single = record[i]
 78    GMfatt = (iGMfact[i])*g
 79    dt, nPts = ReadRecord.ReadRecord(record_single+'.AT2', record_single+'.dat')
 80    op.timeSeries('Path', count, '-dt', dt, '-filePath', record_single+'.dat', '-factor', GMfatt)
 81    op.pattern('UniformExcitation', IDloadTag, iGMdirection[i], '-accel', 2)
 82    count = count + 1
 83
 84maxNumIter = 10
 85op.wipeAnalysis()
 86op.constraints('Transformation')
 87op.numberer('RCM')
 88op.system('BandGeneral')
 89#op.test('EnergyIncr', Tol, maxNumIter)
 90#op.algorithm('ModifiedNewton')
 91#NewmarkGamma = 0.5
 92#NewmarkBeta = 0.25
 93#op.integrator('Newmark', NewmarkGamma, NewmarkBeta)
 94#op.analysis('Transient')
 95
 96
 97#Nsteps =  int(TmaxAnalysis/ DtAnalysis)
 98#
 99#ok = op.analyze(Nsteps, DtAnalysis)
100
101tCurrent = op.getTime()
102
103# for gravity analysis, load control is fine, 0.1 is the load factor increment (http://opensees.berkeley.edu/wiki/index.php/Load_Control)
104
105test = {1:'NormDispIncr', 2: 'RelativeEnergyIncr', 3:'EnergyIncr', 4: 'RelativeNormUnbalance',5: 'RelativeNormDispIncr', 6: 'NormUnbalance'}
106algorithm = {1:'KrylovNewton', 2: 'SecantNewton' , 3:'ModifiedNewton' , 4: 'RaphsonNewton',5: 'PeriodicNewton', 6: 'BFGS', 7: 'Broyden', 8: 'NewtonLineSearch'}
107
108#tFinal = TmaxAnalysis
109tFinal = nPts*dt
110time = [tCurrent]
111u3 = [0.0]
112u4 = [0.0]
113ok = 0
114
115while tCurrent < tFinal:
116#    ok = op.analyze(1, .01)     
117    for i in test:
118        for j in algorithm: 
119            if j < 4:
120                op.algorithm(algorithm[j], '-initial')
121                
122            else:
123                op.algorithm(algorithm[j])
124            while ok == 0 and tCurrent < tFinal:
125                    
126                op.test(test[i], Tol, maxNumIter)        
127                NewmarkGamma = 0.5
128                NewmarkBeta = 0.25
129                op.integrator('Newmark', NewmarkGamma, NewmarkBeta)
130                op.analysis('Transient')
131                ok = op.analyze(1, .01)
132                
133                if ok == 0 :
134                    tCurrent = op.getTime()                
135                    time.append(tCurrent)
136                    u3.append(op.nodeDisp(3,1))
137                    u4.append(op.nodeDisp(4,1))
138                    print(test[i], algorithm[j], 'tCurrent=', tCurrent)
139       
140import matplotlib.pyplot as plt
141plt.figure(figsize=(8,8))
142plt.plot(time, u3)
143plt.ylabel('Horizontal Displacement of node 3 (in)')
144plt.xlabel('Time (s)')
145plt.savefig('Horizontal Disp at Node 3 vs time-uniform excitation-acctime.jpeg', dpi = 500)
146plt.show()
147
148plt.figure(figsize=(8,8))
149plt.plot(time, u4)
150plt.ylabel('Horizontal Displacement of node 4 (in)')
151plt.xlabel('Time (s)')
152plt.savefig('Horizontal Disp at Node 4 vs time-uniform excitation-acctime.jpeg', dpi = 500)
153plt.show() 
154#
155
156op.wipe()