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
EQ ground motion with gravity- uniform excitation of structure
In this example, the Uniaxial Section of Example 2b is replaced by a fiber section. Inelastic uniaxial materials are used in this example,
Which are assigned to each fiber, or patch of fibers, in the section.
In this example the axial and flexural behavior are coupled, a characteristic of the fiber section.
The nonlinear/inelastic behavior of a fiber section is defined by the stress-strain response of the uniaxial materials used to define it.
To run EQ ground-motion analysis (
BM68elc.acc
needs to be downloaded into the same directory)The problem description can be found here (example:2c)
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()