sj464

Posted on Jan 22

I am doing pushover analysis using Opensees in Spyder.

#python #opensees #beginners #pushoveranalysis

I only get a linear performance curve. Could you help me to find out where the problem is?

print("==========================")

from openseespy.opensees import *

import os, sys
import numpy as np
import math as math
import matplotlib.pyplot as plt

print("Starting RCFrameGravity example")

wipe()

Create ModelBuilder (with two-dimensions and 3 DOF/node)

model('basic', '-ndm', 2, '-ndf', 3)

Create nodes

------------

Set parameters for overall model geometry

height = 1300.0

Create nodes

tag, X, Y

node(1, 0.0, 0.0)
node(2, 0.0, 0.0)
node(3, 0.0, height)

Fix supports at base of columns

tag, DX, DY, RZ

fix(1, 1, 1, 1)
fix(2, 1, 1, 0)

%%

Define materials

uniaxialMaterial('Hysteretic', 1, 0.006, 75514606.25, 0.03, 75514606.25*1.1,
-0.006, -75514606.25, -0.03, -75514606.25*1.1, 1.0, 1.0, 0.0, 0.0, 0.7)

%%

Define elements

----------------------

Geometry of column elements

tag

geomTransf('PDelta', 1)

Number of integration points along length of element

np = 10

Lobatto integratoin

beamIntegration('Lobatto', 1, 1, np)

Create the coulumns

element('zeroLength', 1, 1, 2, '-mat',1, '-dir',6)
bc = 250.0
dc = 250.0
A = bc*dc
Iy = bc*dc**3/12
Ec = 26692.696
element('elasticBeamColumn', 2, 2, 3, A, Ec, Iy, 1)

import vfo rendering module

import vfo.vfo as vfo

render the model after defining all the nodes and elements

vfo.plot_model(show_eletags="yes")

plot mode shape

vfo.plot_modeshape(modenumber=3)

%%

Define gravity loads

--------------------

a parameter for the axial load

P = 200000.0

Create a Plain load pattern with a Linear TimeSeries

timeSeries('Linear', 1)
pattern('Plain', 1, 1)

Create nodal loads at nodes 3

nd FX, FY, MZ

load(3, -P, 0.0, 0.0)

load(3, 0.0, -P, 0.0)

------------------------------

End of model generation

------------------------------

Start of analysis generation

------------------------------

Create the system of equation, a sparse solver with partial pivoting

system('ProfileSPD')

Create the constraint handler, the transformation method

constraints('Transformation')

Create the DOF numberer, the reverse Cuthill-McKee algorithm

numberer('RCM')

Create the convergence test, the norm of the residual with a tolerance of

1e-12 and a max number of iterations of 50

test('NormDispIncr', 1.0e-12, 50, 3)

Create the solution algorithm, a Newton-Raphson algorithm

algorithm('Newton')

Create the integration scheme, the LoadControl scheme using steps of 0.1

integrator('LoadControl', 0.1)

Create the analysis object

analysis('Static')

------------------------------

End of analysis generation

------------------------------

Finally perform the analysis

------------------------------

perform the gravity load analysis, requires 50 steps to reach the load level

analyze(50)

Print out the state of nodes 3

print node 3

Print out the state of element 1 and 2

print ele 1 2

u3 = nodeDisp(3, 1)

results = open('results.out', 'a+')

results.close()

print("==========================")

%%

----------------------------------------------------

Start of Model Generation & Initial Gravity Analysis

----------------------------------------------------

Set the gravity loads to be constant & reset the time in the domain

loadConst('-time', 0.0)

----------------------------------------------------

End of Model Generation & Initial Gravity Analysis

----------------------------------------------------

Start of additional modelling for lateral loads

----------------------------------------------------

Define lateral loads

--------------------

Set some parameters

H = 200000.0 # Reference lateral load

Set lateral load pattern with a Linear TimeSeries

pattern('Plain', 2, 1)

Create nodal loads at nodes 3

nd FX FY MZ

load(3, H, 0.0, 0.0)

----------------------------------------------------

End of additional modelling for lateral loads

----------------------------------------------------

Start of modifications to analysis for push over

----------------------------------------------------

Set some parameters

dU = 1.0 # Displacement increment

Change the integration scheme to be displacement control

node dof init Jd min max

integrator('DisplacementControl', 3, 1, dU, 1, dU, dU)

----------------------------------------------------

End of modifications to analysis for push over

----------------------------------------------------

------------------------------

Start of recorder generation

------------------------------

Stop the old recorders by destroying them

remove recorders

Create a recorder to monitor nodal displacements

recorder Node -file node3.out -time -node 3 -dof 1 2 3 disp

Create a recorder to monitor element forces in columns

recorder EnvelopeElement -file ele12.out -time -ele 1 2 forces

--------------------------------

End of recorder generation

---------------------------------

------------------------------

Finally perform the analysis

------------------------------

Set some parameters

maxU = 250.0 # Max displacement
currentDisp = 0.0
ok = 0

test('NormDispIncr', 1.0e-12, 1000)
algorithm('ModifiedNewton', '-initial')

topdispl = [] #empty list to record top node displacement
ndispl = {1: [], 2: [], 3: []}
eforce = {(1,1): [], (1,2):[], (1,3):[], (1,4):[], (1,5):[], (1,6):[],
(2,1): [], (2,2):[], (2,3):[], (2,4):[], (2,5):[], (2,6):[]}

while ok == 0 and currentDisp < maxU:

ok = analyze(1)
print('currentDisp =',currentDisp)

# if the analysis fails try initial tangent iteration
if ok != 0:
    print("modified newton failed")
    break
# print "regular newton failed .. lets try an initail stiffness for this step"
# test('NormDispIncr', 1.0e-12,  1000)
# # algorithm('ModifiedNewton', '-initial')
# ok = analyze(1)
# if ok == 0:
#     print "that worked .. back to regular newton"

# test('NormDispIncr', 1.0e-12,  10)
# algorithm('Newton')

currentDisp = nodeDisp(3, 1)

#record each step
topdispl.append(nodeDisp(3,1))
for ntag in [1,2,3]:
    ndispl[ntag].append(nodeDisp(ntag,1))

for etag in [1,2]:
    eforce[etag,1].append(eleForce(etag,1))
    eforce[etag,2].append(eleForce(etag,2))
    eforce[etag,3].append(eleForce(etag,3))
    eforce[etag,4].append(eleForce(etag,4))
    eforce[etag,5].append(eleForce(etag,5))
    eforce[etag,6].append(eleForce(etag,6))