Created on Thu May 28 16:27:29 2020
@author: praf_malla@hotmail.com, Nepal
@ Corona time, Lockdown
"""
# =============================================================================
# ANALYSIS OF RC COLUMN
# =============================================================================
# Units in N, mm, sec
print("==========================")
from openseespy.postprocessing.Get_Rendering import *
from openseespy.opensees import *
wipe()
import numpy as np
import matplotlib.pyplot as plt
from math import asin, sqrt
print("Starting RC column analysis")
PI = 2 * asin(1.0)
# Create ModelBuilder (with two-dimensions and 3 DOF/node)
model('basic', '-ndm', 2, '-ndf', 3)
# Set parameters for overall model geometry
height = 2000.0 # height of column
# Create nodes
node(1, 0.0, 0.0)
node(2, 0.0, height)
# Fix supports at base of columns
# tag, DX, DY, RZ
fix(1, 1, 1, 1)
# Define cross-section for nonlinear columns
# some parameters
colWidth = 400 ; # direction normal to loading
colDepth = 400 ; # direction parallel to loading
cover = 40 ; # concrete cover
barDia = 20 ; # longitudinal bar dia
As = PI/4*barDia**2 # area of longitudinal bar
# Properties of concrete
# Define parameters of concrete
fc = -30; # Concrete compressive strength (+Tension, -compression)
Ec = 5700*sqrt(-fc); # Modulus of elasticity of concrete
e0 = 0.002; # strain at compressive strength of concrete
# Define parameters of stirrups for confinement
fyh = 415; # yield strength of stirrups
diaStirrup = 8 ; # Diameter of stirrups
(nLY, nLZ) = (2, 3) # no of Legs of stirrup parallel to Y and Z - direction
sh = 150 ; # Spacing of stirrups
# Derived confinement parameters from stirrups
hP = colDepth-2*cover; # Depth of concrete core to oustide of stirrups
bP = colWidth-2*cover; # Width of concrete core to oustide of stirrups
Aso = PI/4*diaStirrup**2; # Area of stirrups
Vstirrup = Aso*(nLY*hP+nLZ*bP)*1000/sh; # Volume of stirrups per 1000 mm
Vcore = hP*bP*1000; # Volume of concrete core in 1000 mm
Rhos = Vstirrup/Vcore; # ratio of volume of hoop steel to the volume of concrete core measured to the outside of hoops.
# Parameters of Confined concrete
Kfc = 1+Rhos*fyh/(-fc); # Ratio of confined to unconfined concrete strength (by Scott)
fc1C = -Kfc*fc; # Compressive strength of confined concrete
eps1C = Kfc*e0; # Strain at confined compressive strength
fc2C = 0.2*fc1C; # stress of confined concrete at ultimate strain
e50uC = (0.0207+eps1C*fc1C)/(fc1C-6.89)
e50hC = 3/4*Rhos*sqrt(hP/sh)
Zc = 0.5/(e50uC+e50hC-eps1C) ; # Descending branch of Confined concrete
eps2C = 0.8*fc1C/Zc+eps1C; # strain at ultimate stress of confined concrete
ftUC = 0.14*fc1C; EtsC = 2*ftUC/eps1C; lambdaC = 0.1
# Parameters of Unconfined concrete
fc1U =-fc ; # Compressive strength of unconfined concrete
eps1U = e0 ; # strain at compressvie strength of unconfined concrete
fc2U = 0.2*fc1U ; # stress of unconfined concrete at ultimate strain
e50uU = (0.0207+eps1U*fc1U)/(fc1U-6.89)
e50hU = 0.0 ; # There is no sirrups i.e., Rhos = 0.0
Zu = 0.5/(e50uU+e50hU-eps1U); # Descending branch of Unconfined concrete
eps2U = 0.8*fc1U/Zu+eps1U ; # Strain at ultimate stress of unconfined concrete
ftUU = 0.14*fc1U; EtsU = 2*ftUU/eps1U; lambdaU = 0.1
# ------------------------------------------
# CONCRETE tag
# Core concrete (confined)
uniaxialMaterial('Concrete02', 1, -fc1C, -eps1C, -fc2C, -eps2C, lambdaC , ftUC, EtsC)
# Cover concrete (unconfined)
uniaxialMaterial('Concrete02', 2, -fc1U, -eps1U, -fc2U, -eps2U, lambdaU , ftUU, EtsU)
print ("Parameters of Confined concrete", -fc1C, -eps1C, -fc2C, -eps2C)
print ("Parameters of UnConfined concrete", -fc1U, -eps1U, -fc2U, -eps2U)
# Plot stress-strain curve of Confined and Unconfined concrete
plt.figure()
sigC = (0, fc1C, fc2C); epsC = (0, eps1C, eps2C)
sigU = (0, fc1U, fc2U); epsU = (0, eps1U, eps2U)
plt.plot(epsC, sigC, linewidth=2.0)
plt.plot(epsU, sigU, linewidth=2.0)
plt.xlabel('strain')
plt.ylabel('compressive stress MPa')
# Longitudinal Reinforcing steel
fy = 400.0; # Yield stress
E = 200000.0; # Young's modulus
Bs = 0.01; # strain hardening ratio
(R0, cR1, cR2) = (18, 0.925, 0.15) # Steel parameters
# Variables derived for formation of section
coverY = colDepth / 2.0
coverZ = colWidth / 2.0
coreY = coverY-cover
coreZ = coverZ-cover
(nfCoreY, nfCoreZ, nfCoverY, nfCoverZ)= (20, 20, 20, 20)
secTag = 1
# tag
uniaxialMaterial('Steel02', 3, fy, E, Bs, R0, cR1, cR2)
section('Fiber', secTag)
# Create the concrete core fibers
patch('quad', 1, nfCoreZ, nfCoreY, -coreY, coreZ, -coreY, -coreZ, coreY, -coreZ, coreY, coreZ,)
# Create the concrete cover patches (right,left,bottom, top)
patch('quad', 2, 2, nfCoverY, -coverY, coverZ, - coreY, coreZ, coreY, coreZ, coverY, coverZ)
patch('quad', 2, 2, nfCoverY, -coreY, -coreZ, -coverY, -coverZ, coverY, -coverZ, coreY, -coreZ)
patch('quad', 2, nfCoverZ, 2, -coverY, coverZ, -coverY, -coverZ, -coreY, -coreZ, -coreY, coreZ)
patch('quad', 2, nfCoverZ, 2, coreY, coreZ, coreY, -coreZ, coverY, -coverZ, coverY, coverZ)
# Create the reinforcing fibers (top, middle, bottom)
layer('straight', 3, 3, As, coreY, coreZ, coreY, -coreZ)
layer('straight', 3, 2, As, 0.0, -coreZ, 0.0, coreZ)
layer('straight', 3, 3, As, -coreY, coreZ, -coreY, -coreZ)
# Define column elements
# ----------------------
geomTag = 1
# Geometry of column elements
# tag
geomTransf('Corotational', geomTag)
np = 4 # Number of integration points along length of element
# Lobatto integratoin
#beamIntegration('Lobatto', 1, 1, np)
# Create the coulumns using Beam-column elements
# e tag ndI ndJ transfTag integrationTag
#eleType = 'forceBeamColumn'
#element(eleType, 1, 1, 2, 1, 1)
# mixedBeamColumn tag iNode jNode numIntgrPts secTag transftag
eleType = 'mixedBeamColumn'
element(eleType, 1, 1, 2, np, secTag, geomTag, '-integration', 'Lobatto')
# Define gravity loads
# a parameter for the axial load
AxialLoad = 0.1*fc1U*colDepth*colWidth; # 10% of axial capacity of columns
print("10 % of axial load ratio", AxialLoad/1000, "kN")
# Create a Plain load pattern with a Linear TimeSeries
timeSeries ('Linear', 1)
pattern ('Plain', 1, 1)
# Create nodal loads at nodes 3 & 4
# nd FX, FY, MZ
load(2, 0.0, -AxialLoad, 0.0)
# ------------------------------
# End of model generation
# ------------------------------
# ------------------------------
# Start of analysis generation
# ------------------------------
# Gravity Analysis
# Create the system of equation, a sparse solver with partial pivoting
Tol = 1.0e-6
system ('BandGeneral')
constraints ('Plain')
numberer ('RCM')
test ('NormDispIncr', Tol, 10, 3)
algorithm ('Newton')
integrator ('LoadControl', 0.1)
analysis ('Static')
analyze (10)
# Display Model
plot_model()
print("==========================")
print("Gravity Analysis Completed")
# ----------------------------------------------------
# End of Model Generation & Gravity Analysis
# ----------------------------------------------------
print("Pushover Analysis Started!")
# Set the gravity loads to be constant & reset the time in the domain
loadConst ('-time', 0.0)
# Set some parameters
LateralLoad = 1000.0 # Reference lateral load of 1kN
# Set lateral load pattern with a Linear TimeSeries
pattern('Plain', 2, 1)
# Create nodal loads at nodes 3 & 4
# nd FX FY MZ
load(2, LateralLoad, 0.0, 0.0)
# ----------------------------------------------------
# Start of push over analysis
# ----------------------------------------------------
# Set some parameters
dU = 1 # Displacement increment
dUMin = 1
dUMax = 1
(Tol, MaxIter) = (1.0e-12, 10) # Tolerance and Maximum iteration
(controlNode, controlDOF) =(2, 1) # Control node and control DOF
system ('BandGeneral')
constraints ('Plain')
numberer ('RCM')
test ('NormDispIncr', Tol, MaxIter)
analysis ('Static');
# Change the integration scheme to be displacement control
# node dof init Jd min max
integrator ('DisplacementControl', controlNode, controlDOF, dU, 1, dUMin,dUMax)
#integrator ('LoadControl', 1)
algorithm ('Newton')
# Set some parameters
maxU = 140.0 # Max displacement
currentDisp = nodeDisp(controlNode, controlDOF)
ok = 0
# perform the analysis
Nsteps=int(maxU/dU)
import numpy as np
Nsteps=int(maxU/dU) # Number of data to be stored
data=np.zeros( (Nsteps+1, 2) ) # Initial variable to store data
j = 0 # Counter to store data
while (ok == 0 and currentDisp < maxU):
j+=1
analyze(1)
# if the analysis fails try initial tangent iteration
if ok != 0:
print("regular newton failed .. lets try an initail stiffness for this step")
test('NormDispIncr', 1.0e-6, 1000)
algorithm('ModifiedNewton', '-initial')
ok = analyze(1)
algorithm('Newton')
currentDisp = nodeDisp(controlNode, controlDOF)
factorLoad =getTime()
data[j,0]=currentDisp
data[j,1]=factorLoad
print( "Horizontal load", factorLoad, 'kN,','displacement at control node',currentDisp,'mm')
plt.plot(data[:,0], data[:,1], linewidth=2.0)
plt.xlabel('Lateral displacement')
plt.ylabel('Lateral Load')
plt.show()
wipe