Nonlin Canti Col Inelstc Uniaxial Mat in Fiber Sec - Dyn EQΒΆ

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. In this example, the Uniaxial Section of Example 2b is replaced by a fiber section. Inelastic uniaxial materials are used in this example,

  3. Which are assigned to each fiber, or patch of fibers, in the section.

  4. In this example the axial and flexural behavior are coupled, a characteristic of the fiber section.

  5. The nonlinear/inelastic behavior of a fiber section is defined by the stress-strain response of the uniaxial materials used to define it.

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

  7. The problem description can be found here (example:2c)

  8. 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 2c. 2D cantilever column, dynamic eq ground motion
 11# EQ ground motion with gravity- uniform excitation of structure
 12#In this example, the Uniaxial Section of Example 2b is replaced by a fiber section. Inelastic uniaxial materials are used in this example, 
 13#which are assigned to each fiber, or patch of fibers, in the section.
 14#In this example the axial and flexural behavior are coupled, a characteristic of the fiber section.
 15#The nonlinear/inelastic behavior of a fiber section is defined by the stress-strain response of the uniaxial materials used to define it.
 16
 17#To run EQ ground-motion analysis (BM68elc.acc needs to be downloaded into the same directory)
 18#the problem description can be found here: http://opensees.berkeley.edu/wiki/index.php/Examples_Manual(example: 2c)
 19# --------------------------------------------------------------------------------------------------
 20#	OpenSees (Tcl) code by:	Silvia Mazzoni & Frank McKenna, 2006
 21
 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#to create a directory at specified path with name "Data"
 44os.chdir('C:\\Opensees Python\\OpenseesPy examples')
 45
 46#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
 47dir = "C:\\Opensees Python\\OpenseesPy examples\\Data-2c"
 48if not os.path.exists(dir):
 49    os.makedirs(dir)
 50#this will create just 'Data' folder    
 51#os.mkdir("Data")    
 52#detect the current working directory
 53#path1 = os.getcwd()
 54#print(path1)
 55#########################################################################################################################################################################
 56
 57#########################################################################################################################################################################
 58op.model('basic', '-ndm', 2, '-ndf', 3) 
 59LCol = 432.0 # column length
 60Weight = 2000.0 # superstructure weight
 61
 62# define section geometry
 63HCol = 60.0 # Column Depth
 64BCol = 60.0 # Column Width
 65
 66PCol =Weight  # nodal dead-load weight per column
 67g = 386.4
 68Mass =  PCol/g
 69
 70ACol = HCol*BCol*1000  # cross-sectional area, make stiff
 71IzCol = (BCol*math.pow(HCol,3))/12 # Column moment of inertia
 72
 73op.node(1, 0.0, 0.0)
 74op.node(2, 0.0, LCol)
 75
 76op.fix(1, 1, 1, 1)
 77
 78op.mass(2, Mass, 1e-9, 0.0)
 79
 80ColSecTag = 1			# assign a tag number to the column section
 81coverCol = 5.0   # Column cover to reinforcing steel NA.
 82numBarsCol = 16  # number of longitudinal-reinforcement bars in column. (symmetric top & bot)
 83barAreaCol = 2.25  # area of longitudinal-reinforcement bars
 84
 85# MATERIAL parameters
 86IDconcU = 1 			# material ID tag -- unconfined cover concrete (here used for complete section)
 87IDreinf = 2 				# material ID tag -- reinforcement
 88
 89# nominal concrete compressive strength
 90fc = -4.0 				# CONCRETE Compressive Strength (+Tension, -Compression)
 91Ec = 57*math.sqrt(-fc*1000) # Concrete Elastic Modulus (the term in sqr root needs to be in psi
 92
 93# unconfined concrete
 94fc1U = fc			# UNCONFINED concrete (todeschini parabolic model), maximum stress
 95eps1U = -0.003			# strain at maximum strength of unconfined concrete
 96fc2U =  0.2*fc1U		# ultimate stress
 97eps2U = -0.01			# strain at ultimate stress
 98Lambda = 0.1				# ratio between unloading slope at $eps2 and initial slope $Ec
 99
100# tensile-strength properties
101ftU = -0.14* fc1U		# tensile strength +tension
102Ets = ftU/0.002			# tension softening stiffness
103
104Fy = 66.8			# STEEL yield stress
105Es = 29000.0				# modulus of steel
106Bs = 0.01				# strain-hardening ratio 
107R0 = 18.0				# control the transition from elastic to plastic branches
108cR1 = 0.925				# control the transition from elastic to plastic branches
109cR2 = 0.15				# control the transition from elastic to plastic branches
110
111op.uniaxialMaterial('Concrete02', IDconcU, fc1U, eps1U, fc2U, eps2U, Lambda, ftU, Ets) # build cover concrete (unconfined)
112op.uniaxialMaterial('Steel02', IDreinf, Fy, Es, Bs, R0,cR1,cR2) # build reinforcement material
113# FIBER SECTION properties -------------------------------------------------------------
114# symmetric section
115#                        y
116#                        ^
117#                        |     
118#             ---------------------     --   --
119#             |   o     o     o    |     |    -- cover
120#             |                       |     |
121#             |                       |     |
122#    z <--- |          +           |     H
123#             |                       |     |
124#             |                       |     |
125#             |   o     o     o    |     |    -- cover
126#             ---------------------     --   --
127#             |-------- B --------|
128#
129# RC section: 
130
131coverY = HCol/2.0	# The distance from the section z-axis to the edge of the cover concrete -- outer edge of cover concrete
132coverZ = BCol/2.0	# The distance from the section y-axis to the edge of the cover concrete -- outer edge of cover concrete
133coreY = coverY-coverCol
134coreZ = coverZ-coverCol
135nfY = 16  # number of fibers for concrete in y-direction
136nfZ = 4			# number of fibers for concrete in z-direction
137
138op.section('Fiber', ColSecTag)
139op.patch('quad', IDconcU, nfZ, nfY, -coverY,coverZ, -coverY,-coverZ, coverY,-coverZ, coverY,coverZ) # Define the concrete patch
140op.layer('straight', IDreinf, numBarsCol, barAreaCol, -coreY,coreZ,-coreY,-coreZ)
141op.layer('straight', IDreinf, numBarsCol, barAreaCol, coreY,coreZ, coreY,-coreZ)
142ColTransfTag = 1
143op.geomTransf('Linear', ColTransfTag)
144numIntgrPts = 5
145eleTag = 1
146
147#import InelasticFiberSection
148
149op.element('nonlinearBeamColumn', eleTag, 1, 2, numIntgrPts, ColSecTag, ColTransfTag)
150
151op.recorder('Node', '-file', 'Data-2c/DFree.out','-time', '-node', 2, '-dof', 1,2,3, 'disp')
152op.recorder('Node', '-file', 'Data-2c/DBase.out','-time', '-node', 1, '-dof', 1,2,3, 'disp')
153op.recorder('Node', '-file', 'Data-2c/RBase.out','-time', '-node', 1, '-dof', 1,2,3, 'reaction')
154#op.recorder('Drift', '-file', 'Data-2c/Drift.out','-time', '-node', 1, '-dof', 1,2,3, 'disp')
155op.recorder('Element', '-file', 'Data-2c/FCol.out','-time', '-ele', 1, 'globalForce')
156op.recorder('Element', '-file', 'Data-2c/ForceColSec1.out','-time', '-ele', 1, 'section', 1, 'force')
157#op.recorder('Element', '-file', 'Data-2c/DCol.out','-time', '-ele', 1, 'deformations')
158
159#defining gravity loads
160op.timeSeries('Linear', 1)
161op.pattern('Plain', 1, 1)
162op.load(2, 0.0, -PCol, 0.0)
163
164Tol = 1e-8 # convergence tolerance for test
165NstepGravity = 10
166DGravity = 1/NstepGravity
167op.integrator('LoadControl', DGravity) # determine the next time step for an analysis
168op.numberer('Plain') # renumber dof's to minimize band-width (optimization), if you want to
169op.system('BandGeneral') # how to store and solve the system of equations in the analysis
170op.constraints('Plain') # how it handles boundary conditions
171op.test('NormDispIncr', Tol, 6) # determine if convergence has been achieved at the end of an iteration step
172op.algorithm('Newton') # use Newton's solution algorithm: updates tangent stiffness at every iteration
173op.analysis('Static') # define type of analysis static or transient
174op.analyze(NstepGravity) # apply gravity
175
176op.loadConst('-time', 0.0) #maintain constant gravity loads and reset time to zero
177 
178#applying Dynamic Ground motion analysis
179GMdirection = 1
180GMfile = 'BM68elc.acc'
181GMfact = 1.0
182
183
184
185Lambda = op.eigen('-fullGenLapack', 1) # eigenvalue mode 1
186import math
187Omega = math.pow(Lambda, 0.5)
188betaKcomm = 2 * (0.02/Omega)
189
190xDamp = 0.02				# 2% damping ratio
191alphaM = 0.0				# M-prop. damping; D = alphaM*M	
192betaKcurr = 0.0		# K-proportional damping;      +beatKcurr*KCurrent
193betaKinit = 0.0 # initial-stiffness proportional damping      +beatKinit*Kini
194
195op.rayleigh(alphaM,betaKcurr, betaKinit, betaKcomm) # RAYLEIGH damping
196
197# Uniform EXCITATION: acceleration input
198IDloadTag = 400			# load tag
199dt = 0.01			# time step for input ground motion
200GMfatt = 1.0			# data in input file is in g Unifts -- ACCELERATION TH
201maxNumIter = 10
202op.timeSeries('Path', 2, '-dt', dt, '-filePath', GMfile, '-factor', GMfact)
203op.pattern('UniformExcitation', IDloadTag, GMdirection, '-accel', 2) 
204
205op.wipeAnalysis()
206op.constraints('Transformation')
207op.numberer('Plain')
208op.system('BandGeneral')
209op.test('EnergyIncr', Tol, maxNumIter)
210op.algorithm('ModifiedNewton')
211
212NewmarkGamma = 0.5
213NewmarkBeta = 0.25
214op.integrator('Newmark', NewmarkGamma, NewmarkBeta)
215op.analysis('Transient')
216
217DtAnalysis = 0.01
218TmaxAnalysis = 10.0
219
220Nsteps =  int(TmaxAnalysis/ DtAnalysis)
221
222ok = op.analyze(Nsteps, DtAnalysis)
223
224tCurrent = op.getTime()
225
226# for gravity analysis, load control is fine, 0.1 is the load factor increment (http://opensees.berkeley.edu/wiki/index.php/Load_Control)
227
228test = {1:'NormDispIncr', 2: 'RelativeEnergyIncr', 4: 'RelativeNormUnbalance',5: 'RelativeNormDispIncr', 6: 'NormUnbalance'}
229algorithm = {1:'KrylovNewton', 2: 'SecantNewton' , 4: 'RaphsonNewton',5: 'PeriodicNewton', 6: 'BFGS', 7: 'Broyden', 8: 'NewtonLineSearch'}
230
231for i in test:
232    for j in algorithm:
233
234        if ok != 0:
235            if j < 4:
236                op.algorithm(algorithm[j], '-initial')
237                
238            else:
239                op.algorithm(algorithm[j])
240                
241            op.test(test[i], Tol, 1000)
242            ok = op.analyze(Nsteps, DtAnalysis)                            
243            print(test[i], algorithm[j], ok)             
244            if ok == 0:
245                break
246        else:
247            continue
248
249u2 = op.nodeDisp(2, 1)
250print("u2 = ", u2)
251
252op.wipe()