Nonlinear Canti Col Uniaxial Inelastic Section- Dyn EQ GM

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. The nonlinear beam-column element that replaces the elastic element of Example 2a requires the definition of the element cross section, or its behavior. In this example,

  3. The Uniaxial Section used to define the nonlinear moment-curvature behavior of the element section is “aggregated” to an elastic response for the axial behavior to define

  4. The required characteristics of the column element in the 2D model. In a 3D model, torsional behavior would also have to be aggregated to this section.

  5. Note:In this example, both the axial behavior (typically elastic) and the flexural behavior (moment curvature) are defined indepenently and are then “aggregated” into a section.

  6. This is a characteristic of the uniaxial section: there is no coupling of behaviors.

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

  8. The problem description can be found here (example:2b)

  9. 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# Example 2b. 2D cantilever column, dynamic eq ground motion
 11# EQ ground motion with gravity- uniform excitation of structure
 12#he nonlinear beam-column element that replaces the elastic element of Example 2a requires the definition of the element cross section, or its behavior. In this example, 
 13#the Uniaxial Section used to define the nonlinear moment-curvature behavior of the element section is "aggregated" to an elastic response for the axial behavior to define 
 14#the required characteristics of the column element in the 2D model. In a 3D model, torsional behavior would also have to be aggregated to this section.
 15#Note:In this example, both the axial behavior (typically elastic) and the flexural behavior (moment curvature) are defined indepenently and are then "aggregated" into a section. 
 16#This is a characteristic of the uniaxial section: there is no coupling of behaviors.
 17
 18#To run EQ ground-motion analysis (BM68elc.acc needs to be downloaded into the same directory)
 19#the problem description can be found here: http://opensees.berkeley.edu/wiki/index.php/Examples_Manual(example:2b)
 20# --------------------------------------------------------------------------------------------------
 21#	OpenSees (Tcl) code by:	Silvia Mazzoni & Frank McKenna, 2006
 22#
 23#    ^Y
 24#    |
 25#    2       __ 
 26#    |          | 
 27#    |          |
 28#    |          |
 29#  (1)       LCol
 30#    |          |
 31#    |          |
 32#    |          |
 33#  =1=      _|_  -------->X
 34#
 35
 36# SET UP ----------------------------------------------------------------------------
 37import openseespy.opensees as op
 38#import the os module
 39import os
 40import math
 41op.wipe()
 42#########################################################################################################################################################################
 43
 44#to create a directory at specified path with name "Data"
 45os.chdir('C:\\Opensees Python\\OpenseesPy examples')
 46
 47#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
 48dir = "C:\\Opensees Python\\OpenseesPy examples\\Data-2b"
 49if not os.path.exists(dir):
 50    os.makedirs(dir)
 51#this will create just 'Data' folder    
 52#os.mkdir("Data")    
 53#detect the current working directory
 54#path1 = os.getcwd()
 55#print(path1)
 56#########################################################################################################################################################################
 57
 58#########################################################################################################################################################################
 59op.model('basic', '-ndm', 2, '-ndf', 3) 
 60LCol = 432.0 # column length
 61Weight = 2000.0 # superstructure weight
 62
 63# define section geometry
 64HCol = 60.0 # Column Depth
 65BCol = 60.0 # Column Width
 66
 67PCol =Weight  # nodal dead-load weight per column
 68g = 386.4
 69Mass =  PCol/g
 70
 71ACol = HCol*BCol*1000  # cross-sectional area, make stiff
 72IzCol = (BCol*math.pow(HCol,3))/12 # Column moment of inertia
 73
 74op.node(1, 0.0, 0.0)
 75op.node(2, 0.0, LCol)
 76
 77op.fix(1, 1, 1, 1)
 78
 79op.mass(2, Mass, 1e-9, 0.0)
 80
 81#Define Elements and Sections
 82ColMatTagFlex  = 2
 83ColMatTagAxial = 3
 84ColSecTag = 1
 85BeamSecTag = 2
 86
 87fc = -4.0 # CONCRETE Compressive Strength (+Tension, -Compression)
 88Ec = 57*math.sqrt(-fc*1000) # Concrete Elastic Modulus (the term in sqr root needs to be in psi
 89
 90#Column Section
 91EICol = Ec*IzCol # EI, for moment-curvature relationship
 92EACol = Ec*ACol # EA, for axial-force-strain relationship
 93MyCol = 130000.0 #yield Moment calculated
 94PhiYCol = 0.65e-4	# yield curvature
 95EIColCrack = MyCol/PhiYCol	# cracked section inertia
 96b = 0.01 # strain-hardening ratio (ratio between post-yield tangent and initial elastic tangent)
 97
 98op.uniaxialMaterial('Steel01', ColMatTagFlex, MyCol, EIColCrack, b) #steel moment curvature isused for Mz of the section only, # bilinear behavior for flexure
 99op.uniaxialMaterial('Elastic', ColMatTagAxial, EACol) # this is not used as a material, this is an axial-force-strain response
100op.section('Aggregator', ColSecTag, ColMatTagAxial, 'P', ColMatTagFlex, 'Mz')  # combine axial and flexural behavior into one section (no P-M interaction here)
101
102ColTransfTag = 1
103op.geomTransf('Linear', ColTransfTag)
104numIntgrPts = 5
105eleTag = 1
106op.element('nonlinearBeamColumn', eleTag, 1, 2, numIntgrPts, ColSecTag, ColTransfTag)
107
108op.recorder('Node', '-file', 'Data-2b/DFree.out','-time', '-node', 2, '-dof', 1,2,3, 'disp')
109op.recorder('Node', '-file', 'Data-2b/DBase.out','-time', '-node', 1, '-dof', 1,2,3, 'disp')
110op.recorder('Node', '-file', 'Data-2b/RBase.out','-time', '-node', 1, '-dof', 1,2,3, 'reaction')
111#op.recorder('Drift', '-file', 'Data-2b/Drift.out','-time', '-node', 1, '-dof', 1,2,3, 'disp')
112op.recorder('Element', '-file', 'Data-2b/FCol.out','-time', '-ele', 1, 'globalForce')
113op.recorder('Element', '-file', 'Data-2b/ForceColSec1.out','-time', '-ele', 1, 'section', 1, 'force')
114#op.recorder('Element', '-file', 'Data-2b/DCol.out','-time', '-ele', 1, 'deformations')
115
116#defining gravity loads
117op.timeSeries('Linear', 1)
118op.pattern('Plain', 1, 1)
119op.load(2, 0.0, -PCol, 0.0)
120
121Tol = 1e-8 # convergence tolerance for test
122NstepGravity = 10
123DGravity = 1/NstepGravity
124op.integrator('LoadControl', DGravity) # determine the next time step for an analysis
125op.numberer('Plain') # renumber dof's to minimize band-width (optimization), if you want to
126op.system('BandGeneral') # how to store and solve the system of equations in the analysis
127op.constraints('Plain') # how it handles boundary conditions
128op.test('NormDispIncr', Tol, 6) # determine if convergence has been achieved at the end of an iteration step
129op.algorithm('Newton') # use Newton's solution algorithm: updates tangent stiffness at every iteration
130op.analysis('Static') # define type of analysis static or transient
131op.analyze(NstepGravity) # apply gravity
132
133op.loadConst('-time', 0.0) #maintain constant gravity loads and reset time to zero
134 
135#applying Dynamic Ground motion analysis
136GMdirection = 1
137GMfile = 'BM68elc.acc'
138GMfact = 1.0
139
140
141
142Lambda = op.eigen('-fullGenLapack', 1) # eigenvalue mode 1
143import math
144Omega = math.pow(Lambda, 0.5)
145betaKcomm = 2 * (0.02/Omega)
146
147xDamp = 0.02				# 2% damping ratio
148alphaM = 0.0				# M-prop. damping; D = alphaM*M	
149betaKcurr = 0.0		# K-proportional damping;      +beatKcurr*KCurrent
150betaKinit = 0.0 # initial-stiffness proportional damping      +beatKinit*Kini
151
152op.rayleigh(alphaM,betaKcurr, betaKinit, betaKcomm) # RAYLEIGH damping
153
154# Uniform EXCITATION: acceleration input
155IDloadTag = 400			# load tag
156dt = 0.01			# time step for input ground motion
157GMfatt = 1.0			# data in input file is in g Unifts -- ACCELERATION TH
158maxNumIter = 10
159op.timeSeries('Path', 2, '-dt', dt, '-filePath', GMfile, '-factor', GMfact)
160op.pattern('UniformExcitation', IDloadTag, GMdirection, '-accel', 2) 
161
162op.wipeAnalysis()
163op.constraints('Transformation')
164op.numberer('Plain')
165op.system('BandGeneral')
166op.test('EnergyIncr', Tol, maxNumIter)
167op.algorithm('ModifiedNewton')
168
169NewmarkGamma = 0.5
170NewmarkBeta = 0.25
171op.integrator('Newmark', NewmarkGamma, NewmarkBeta)
172op.analysis('Transient')
173
174DtAnalysis = 0.01
175TmaxAnalysis = 10.0
176
177Nsteps =  int(TmaxAnalysis/ DtAnalysis)
178
179ok = op.analyze(Nsteps, DtAnalysis)
180
181tCurrent = op.getTime()
182
183# for gravity analysis, load control is fine, 0.1 is the load factor increment (http://opensees.berkeley.edu/wiki/index.php/Load_Control)
184
185test = {1:'NormDispIncr', 2: 'RelativeEnergyIncr', 4: 'RelativeNormUnbalance',5: 'RelativeNormDispIncr', 6: 'NormUnbalance'}
186algorithm = {1:'KrylovNewton', 2: 'SecantNewton' , 4: 'RaphsonNewton',5: 'PeriodicNewton', 6: 'BFGS', 7: 'Broyden', 8: 'NewtonLineSearch'}
187
188for i in test:
189    for j in algorithm:
190
191        if ok != 0:
192            if j < 4:
193                op.algorithm(algorithm[j], '-initial')
194                
195            else:
196                op.algorithm(algorithm[j])
197                
198            op.test(test[i], Tol, 1000)
199            ok = op.analyze(Nsteps, DtAnalysis)                            
200            print(test[i], algorithm[j], ok)             
201            if ok == 0:
202                break
203        else:
204            continue
205
206u2 = op.nodeDisp(2, 1)
207print("u2 = ", u2)
208
209op.wipe()