Cyclic Loading Algorithm

Forum for asking and answering questions related to use of the OpenSeesPy module

Moderators: silvia, selimgunay, Moderators

Post Reply
Prafullamalla
Posts: 160
Joined: Mon Feb 02, 2015 6:32 pm

Cyclic Loading Algorithm

Post by Prafullamalla » Thu Dec 22, 2022 5:22 am

Its converted from the available TCL file.

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()


Prafulla Malla, Nepal
Praf_malla@hotmail.com

Post Reply