Code: Select all
# -*- coding: utf-8 -*-
"""
Created on Wed Nov 30 10:53:29 2022
@author: praf_
"""
import openseespy.opensees as ops
import numpy as np
ops.wipe()
ops.model('basic','-ndm',1,'-ndf',1)
m = 1
g = 981
ops.node(1,0); ops.fix(1,1)
ops.node(2,0); ops.mass(2,m)
F=10000
k = 1000;#24*36000*250**4/12/3000**3
#ops.uniaxialMaterial('Elastic',1,k)
ops.uniaxialMaterial('Steel01 ',1,F,k,-2)
#ops.uniaxialMaterial('Steel02', 1, F,k,0.005, 18, 0.925, 0.15)
ops.element('zeroLength',1,1,2,'-mat',1,'-dir',1)
print("Cyclic Analysis Started!")
(controlNode, controlDOF) =(2, 1)
# Set the gravity loads to be constant & reset the time in the domain
ops.loadConst ('-time', 0.0)
# Set some parameters
LateralLoad = 1000.0 # Reference lateral load of 1kN in N
# Set lateral load pattern with a Linear TimeSeries
# Create a Plain load pattern with a Linear TimeSeries
ops.timeSeries ('Linear', 1)
ops.pattern('Plain', 1, 1)
# Create nodal loads at nodes
# nd FX FY MZ
ops.load(controlNode, LateralLoad)
# Set some parameters
currentDisp = ops.nodeDisp(controlNode, controlDOF)
# Arrays for plotting
outputs = {
"Load": [],
"Displacement": [],
"Drift":["Drift"]
}
# ------------------------
# Displacement steps
# ------------------------
#Drift = [15,0]
height=1
#arr=np.array(Drift) # converting to array to multipy to calculate Displacement from Drift
#arr =height/100*arr # Displacement-Drift*height/100 max for each cycle
#Disp = arr.tolist()
Disp = [5,0,-5,40,0,-40,15,0,-15,35,0,-35,0,10]#,15,0,-15,0,20,0,-20,35,0,-35,35,0,-35,40,0,-40,]
# Set some parameters
currentDisp = ops.nodeDisp(controlNode, controlDOF)
ok = 0 # Counter to store data
currentDisp = ops.nodeDisp(controlNode, controlDOF)
factorLoad =ops.getTime()
outputs["Load"].append(factorLoad)
outputs["Displacement"].append(currentDisp)
outputs["Drift"].append(currentDisp/height*100)
# ------------------------------------------------------------------------------
# Main analysis - performed as a number of analyses equal to the number of steps
# ------------------------------------------------------------------------------
ops.system("BandGeneral") # Overkill, but may need the pivoting!
ops.test("NormDispIncr", 1.0e-7, 10, 0)
ops.numberer("Plain")
ops.constraints("Penalty",1.0e14, 1.0e14)
ops.algorithm("Newton")
ops.analysis("Static")
cycle=0
for Dmax in Disp:
cycle=cycle+1
print('Cycle Number',cycle)
if Dmax > 0:
Dincr = 1
elif Dmax < 0:
Dincr =-1
ops.integrator("DisplacementControl", controlNode, controlDOF, Dincr)
ops.analysis("Static")
if Dincr > 0:
ok = 0
while (ok == 0 and currentDisp < Dmax):
ok = ops.analyze(1)
if ok != 0:
print( " ")
print ("Trying Newton with Initial Tangent ..")
ops.test("NormDispIncr", 1.0e-6, 2000, 0)
ops.algorithm('Newton', '-initial')
ops.analysis("Static")
ok = ops.analyze(1)
ops.test("NormDispIncr", 1.0e-6, 6, 2)
ops.algorithm('Newton')
if ok != 0:
print( " ")
print ("Trying Broyden ..")
ops.test("NormDispIncr", 1.0e-6, 2000, 0)
ops.algorithm('Broyden', 8)
ops.analysis("Static")
ok = ops.analyze(1)
ops.algorithm('Newton')
if ok != 0:
print( " ")
print ("Trying NewtonWithLineSearch ..")
ops.algorithm('NewtonLineSearch', 0.8)
ok = ops.analyze(1)
ops.algorithm('Newton')
currentDisp = ops.nodeDisp(controlNode, controlDOF)
factorLoad =ops.getTime()
outputs["Load"].append(factorLoad)
outputs["Displacement"].append(currentDisp)
outputs["Drift"].append(currentDisp/height*100)
elif Dincr < 0 :
ok = 0
while (ok == 0 and currentDisp > Dmax):
ok = ops.analyze(1)
if ok != 0:
print( " ")
print ("Trying Newton with Initial Tangent ..")
ops.test("NormDispIncr", 1.0e-6, 2000, 0)
ops.algorithm('Newton', '-initial')
ok = ops.analyze(1)
ops.test("NormDispIncr", 1.0e-6, 6, 2)
ops.algorithm('Newton')
if ok != 0:
print( " ")
print ("Trying Broyden ..")
ops.test("NormDispIncr", 1.0e-6, 2000, 0)
ops.algorithm('Broyden', 8)
ok = ops.analyze(1)
ops.algorithm('Newton')
if ok != 0:
print( " ")
print ("Trying NewtonWithLineSearch ..")
ops.algorithm('NewtonLineSearch', 0.8)
ok = ops.analyze(1)
ops.algorithm('Newton')
currentDisp = ops.nodeDisp(controlNode, controlDOF)
factorLoad =ops.getTime()
outputs["Load"].append(factorLoad)
outputs["Displacement"].append(currentDisp)
outputs["Drift"].append(currentDisp/height*100)
if ok ==0:
print ('Cyclic analysis completed SUCCESSFULLY ')
else:
print ("Cyclic analysis FAILED");
for item in outputs:
outputs[item] = np.array(outputs[item])
#P PLOTTING CURVE
# x axis values
x = outputs["Displacement"]
# corresponding y axis values
y = outputs["Load"]
import matplotlib.pyplot as plt
import numpy as np
ops.wipe
# Deflection
## Set parameters for the plot
plt.rcParams.update({'font.size': 7})
plt.figure(figsize=(6,5), dpi=200)
plt.rc('font', family='serif')
plt.plot(x,y, color="red", linewidth=0.8, linestyle="-", label='Cyclic')
plt.axhline(0, color='black', linewidth=0.4)
plt.axvline(0, color='black', linewidth=0.4)
#plt.xlim(-15, 15)
#plt.xticks(np.linspace(-9,9,10,endpoint=True))
plt.grid(linestyle='dotted')
plt.xlabel('Displacement (mm)')
plt.ylabel('Base Shear (kN)')
plt.legend()
#plt.savefig("RCC1_PushoverCurve.png",dpi=1200)
plt.show()