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 »

Code: Select all

# -*- coding: utf-8 -*-
"""
Created on Thu Jul 27 17:22:45 2023

@author: praf_
"""
# height is reguired to calculate drift from displacement. This was for experimental specimen
DriftAmp =     [0.25,-0.25,0.25,-0.25
                 ,0.5,-0.5,0.5,-0.5,
                   0.75,-0.75,0.75,-0.75,
                   1,-1,1,-1,
                   1.5,-1.5,1.5,-1.5,
                    2,-2,2,-2,
                    2.5,-2.5,2.5,-2.5,
                    3,-3, 3,-3,
                    4,-4,5,-0.1] 
       
     
                 
RunCyclic2Converge(controlNode, DriftAmp, fsheight,  controlDOF=1)
# ---------RunCyclic2Converge.tcl----------------------------------------------------------------
# Changed into a procedure from Opensees Web-site. Slightly modified
# to stop analyzing if the loadfactor (pseudo-time) falls below zero.
# Also, it accepts an optional input, "dof" which defaults to 1. 
# this is the dof in the POSITIVE OR NEGATIVE direction of which we push.

# Modified from Vam
import openseespy.opensees as ops
import numpy as np

def RunCyclic2Converge(controlNode, Drift, height,  controlDOF=1):
    # Set some parameters
    # Arrays for plotting
    
    currentDisp = ops.nodeDisp(controlNode, controlDOF)
    
    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()
    # ------------------------------------------------------------------------------
    # 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")
    
    
    for Dmax in Disp:
        
        if Dmax > 0:
            Dincr = 0.001
        elif Dmax < 0:
            Dincr =-0.001
                
        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')
                    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)
          
        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)

    errorFileName = "ncerror.out"
    # output success/failure results to a file
    if errorFileName != "":
        with open(errorFileName, "w") as errorFileID:
            errorFileID.write(str(ok))

    if ok != 0:
        print("DispControl Analysis FAILED")
    else:
        print("DispControl Analysis SUCCESSFUL")


Prafulla Malla, Nepal
Praf_malla@hotmail.com
Post Reply