openseesyPy analysis question

[gojaewook]

Select Topic Area

Question

Body

I'm doing structural analysis with open source.

While developing the algorithm, I completed the analysis model and ran the analysis.
On the way, I applied a load but got the result that the node has no reaction force.
What is the error?

import numpy as np
import numpy.linalg
#import vfo.vfo as vfo
import openseespy.opensees as ops
import opsvis as opsv
import matplotlib.pyplot as plt
import math

ops.wipe()

ops.model("Basic", '-ndm', 3, '-ndf', 6)

Define Material properties

m = 1
mm = 0.001*m
mm2 = mm ** 2
mm3 = mm ** 3
mm4 = mm ** 4
IDconc = 1

Concrete

fc = -37.08
Ec = 8500*(-fc)**(1/3)
fck = -40
epsc = -0.002
epscu = -0.003
fct = 2.62
et = fct/Ec

Material

ops.uniaxialMaterial("Elastic", IDconc, Ec)

node

ops.node(1, 0, 0, 0)
ops.node(2, 500, 0, 0)
ops.node(3, 1000, 0, 0)

fix

ops.fix(1, 1, 1, 1, 1, 0, 0)
ops.fix(3, 1, 0, 1, 1, 0, 0)

Section

v = 0.2
Iz = 121732535459 * mm4
Iy = 13754180616 * mm4
J = Iz+Iy
G = Ec/2*(1+v)
GJ = G*J

Geometry pro

section5 = np.array([-300, -300, 300, -300, -300, 300, -300, 300])
section1 = np.array([0.495, -0.5, 0.555, -0.5, 0.555, 0.5, 0.495, 0.5])*1000
section2 = np.array([0.395, -0.1, 0.495, -0.2, 0.495, 0.2, 0.495, 0.1])*1000
section3 = np.array([-0.205, -0.1, 0.395, -0.1, 0.395, 0.1, -0.205, 0.1])*1000
section4 = np.array([-0.405, -0.225, -0.205, -0.1, -0.405, 0.1, -0.205, -0.225])*1000
section5 = np.array([-0.555, -0.225, -0.405, -0.225, -0.405, 0.225, -0.555, 0.225])*1000
nfY1, nfZ1 = 4,20
nfY2, nfZ2 = 2,8
nfY3, nfZ3 = 8,8
nfY4, nfZ4 = 4,8
nfY5, nfZ5 = 2,8

section

ops.section('Fiber', 1, '-GJ', GJ)
#ops.patch('rect', IDconc, nfY1, nfZ1, 1, 2, *section5)
ops.patch('quad', IDconc, nfY1, nfZ1, 375, -500, 475, -500, 475, 500, 375, 500)
ops.patch('quad', IDconc, nfY2, nfZ2, 275, -100, 375, -200, 375, 200, 275, 100)
ops.patch('quad', IDconc, nfY3, nfZ3, -325, -100, 275, -100, 275, 100, -325, 100)
ops.patch('quad', IDconc, nfY4, nfZ4, -525, -225, -325, -100, -325, 100, -525, 225)
ops.patch('quad', IDconc, nfY5, nfZ5, -675, -225, -525, -225, -525, 225, -675, 225)

ops.section('Fiber', 2, '-GJ', GJ)
#ops.patch('rect', IDconc, nfY1, nfZ1, 1, 2, *section5)
ops.patch('quad', IDconc, nfY1, nfZ1, 375, -500, 475, -500, 475, 500, 375, 500)
ops.patch('quad', IDconc, nfY2, nfZ2, 275, -100, 375, -200, 375, 200, 275, 100)
ops.patch('quad', IDconc, nfY3, nfZ3, -325, -100, 275, -100, 275, 100, -325, 100)
ops.patch('quad', IDconc, nfY4, nfZ4, -525, -225, -325, -100, -325, 100, -525, 225)
ops.patch('quad', IDconc, nfY5, nfZ5, -675, -225, -525, -225, -525, 225, -675, 225)

ops.section('Fiber', 3, '-GJ', GJ)
#ops.patch('rect', IDconc, nfY1, nfZ1, 1, 2, *section5)
ops.patch('quad', IDconc, nfY1, nfZ1, 375, -500, 475, -500, 475, 500, 375, 500)
ops.patch('quad', IDconc, nfY2, nfZ2, 275, -100, 375, -200, 375, 200, 275, 100)
ops.patch('quad', IDconc, nfY3, nfZ3, -325, -100, 275, -100, 275, 100, -325, 100)
ops.patch('quad', IDconc, nfY4, nfZ4, -525, -225, -325, -100, -325, 100, -525, 225)
ops.patch('quad', IDconc, nfY5, nfZ5, -675, -225, -525, -225, -525, 225, -675, 225)

TransfTag = 1
ops.geomTransf("Linear", TransfTag, 0, -1, 0)

Integrator

#op.beamIntegration('Legendre', int(nodes[a][i-2][0]), int(nodes[a][i-2][0]), numIntgrPts)
#op.element('dispBeamColumn', int(nodes[a][i-2][0]), int(nodes[a][i-2][0]), int(nodes[a][i-1][0]), BeamTransfTag, int(nodes[a][i-2][0]))
numIntegrationPts = 5

ops.beamIntegration("Legendre", 1, 1, numIntegrationPts)
ops.beamIntegration("Legendre", 2, 2, numIntegrationPts)
ops.element('dispBeamColumn', 1, 1, 2, TransfTag, 1)
ops.element('dispBeamColumn', 2, 2,3, TransfTag, 1)

Time Series

tsTag = 1
ops.timeSeries("Linear", tsTag)

9.1 Gravity

9.1.1 Pattern

patternTag1 = 1
ops.pattern('Plain', patternTag1, tsTag)

ops.load(2, 0, 0, -260000, 0, 0, 0)

9.1.2 Constraint

ops.constraints('Plain') # How it handles boundary conditions

9.1.3 Numberer

ops.numberer('Plain') # Renumber DoF's to minimize band-width (optimization)

9.1.4 System

ops.system('BandGen') # How to strore and solve the system of equations in the analysis

9.1.5 Test

Tol = 1.3e-8 # Convergence tolerance for test
iteration = 5 # The max number of iterations to check before returning failure condition
NormType = 2 # Type of norm (2-norm)
ops.test('NormDispIncr', Tol, iteration, 0, NormType) # Determin if convergence has been acheved at the end of an iteration steps

9.1.6 Algorithm

ops.algorithm('Newton') # Create a Linear algorithm to solve the system equations.

9.1.7 Integrator

numSteps = 20
ops.integrator("LoadControl", 1/numSteps) # Determine the next time step for an analysis

9.1.8 Analysis

ops.analysis("Static") # Define type of analysis.

9.1.9 Analyze

ops.analyze(numSteps)

zz = ops.nodeReaction(1)
zz1 = ops.nodeReaction(2)
zz2 = ops.nodeReaction(3)
#zz7 = op.nodeDisp(1)
print(zz, zz1, zz2)

Print the model

ops.printModel()

Get and print all node tags

#node_tags = op.getNodeTags()
#print("Node tags:", node_tags)

Get and print all element tags

#ele_tags = op.getEleTags()
#print("Element tags:", ele_tags)

[Sarveshero3]

The reactions are zero because ops.reactions() was never called before ops.nodeReaction().

In OpenSeesPy, nodeReaction() doesn't compute reactions on the fly, it just reads a cached value that only gets populated when you explicitly call ops.reactions() after the analysis converges. Without that call, every node returns zero reaction regardless of the actual state.

Your displacement data confirms the analysis itself worked correctly and Node 2 shows -0.0031275 in z, meaning the load did propagate through the structure. The solve is fine i think only the reaction reporting step is missing.

Fix:

ops.analyze(numSteps)

ops.reactions()          # <-- add this before reading reactions
zz = ops.nodeReaction(1)
zz1 = ops.nodeReaction(2)
zz2 = ops.nodeReaction(3)
print(zz, zz1, zz2)

If you want reactions including inertial/damping contributions in a dynamic analysis, use ops.reactions('-dynamic') instead, but for your static LoadControl case, the plain ops.reactions() call is what you need.