I am a Ph.D. candidate working on structural analysis using OpenSees, and I have encountered an issue that I am struggling to understand. While simulating some BRBF archetypes, I noticed that using the Steel4 material sometimes resulted in non-zero reactions for DOFs that were not restrained.
I am hoping that someone with more experience using OpenSees might be able to shed some light on why this is happening and how to fix it. Any guidance would be greatly appreciated. Thank you in advance for your help.
Here is a python script that demonstrates this behavior:
Code: Select all
"""
Demonstrate issue with ops.reactions() when using Steel4.
Description
-----------
I always expected reactions to have a nonzero value only when a DOF is
restrained, but it looks like in this example we get a reaction for a
DOF that is unrestrained (node 1, DOF 1). Although the reactions of
the restrained are reported correctly, this can be a problem when all
reactions are accumulated to calculate quantities such as building
base shear and the like.
In its current form, the example is for a 1D problem, but I have
confirmed that the same unexpected behavior is observed with
ops.model('basic', '-ndm', 3,'-ndf', 6).
Also, the same issue appears when using a zeroLength element instead
of a truss, indicating that the issue is unrelated to the truss
element. When using an elastic uniaxial material, I don't get a
reaction at the unrestrained DOF, as expected. Another interesting
observation is that when I reduce the displacement increment, this
strange behavior goes away.
Please run this code and modify the investigative parameters to
reproduce my observations.
Is this a bug? If not, can you help me understand why we get a
reaction there?
---
Date Created: Sun Jan 8 20:33:08 PST 2023
Author: John Vouvakis Manousakis
email: ioannis_vm@berkeley.edu
"""
import openseespy.opensees as ops
import numpy as np
import matplotlib.pyplot as plt
# ~~~ investigative parameters ~~~
# material to use: change to True to switch to Elastic
use_elastic = False
# smaller displacement increment: change to True and the reaction goes
# away
use_small_step = False
# use zeroLength instead of Truss: makes no difference in the observed
# behavior
use_zerolength = False
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# define peak displacements for cycling loading
peaks = np.array(
[
0.286, -0.253, 0.337, -0.244, 0.320, -0.236, 0.328, -0.228,
0.328, -0.228, 0.328, -0.228, 0.758, -0.699, 0.749, -0.682,
0.749, -0.682, 0.749, -0.682, 1.583, -1.491, 1.575, -1.482,
1.566, -1.474, 1.575, -1.474, 2.391, -2.282, 2.391, -2.282,
1.575, -1.465, 3.309, -3.133, 3.276, -3.124, 3.276, -3.116,
3.267, -3.107, 3.267, -3.107, -0.152
]) # inch
# define parameters
# (I am simulating a BRB)
young_eff = 34000000. # lb/in2
fyield = 50000.0 # lb
elm_len = 200.00 # in
area = 3.00 # in2
# initialize
ops.wipe()
ops.model('basic', '-ndm', 1, '-ndf', 1)
# define nodes
ops.node(0, 0.00)
if use_zerolength:
ops.node(1, 0.00)
else:
ops.node(1, elm_len)
# define restraints
ops.fix(0, True)
ops.fix(1, False)
# note that DOF 1 of node 1 is *not* restrained
if use_elastic:
ops.uniaxialMaterial(
'Elastic', 2, young_eff)
else:
ops.uniaxialMaterial(
'Steel4', 2, fyield, young_eff, '-asym',
'-kin', 0.003, 25.0,
0.9, 0.15, 0.023, 25.0, 0.9, 0.15,
'-iso', 0.0025, 1.34,
0.004, 1.0, 1.0, 0.0045, 0.77, 0.004, 1.0
)
if use_zerolength:
# connect the two nodes with a zerolength element
ops.element('zeroLength', 1, 0, 1,
'-mat', 2, '-dir', 1)
else:
# connect the two nodes with a truss element
ops.element(
'Truss', 1, 0, 1, area, 2,
'-rho', 0.0, '-cMass', 0, '-doRayleigh', 0)
# run cyclic pushover analysis
ops.timeSeries('Linear', 1)
ops.pattern('Plain', 1, 1)
ops.load(1, 1.00)
ops.numberer('Plain')
ops.constraints('Plain')
# maximum displacement increment
incr = 0.1
if use_small_step:
incr *= 0.1
if use_zerolength:
# convert to strain
incr /= elm_len
peaks /= elm_len
# instantiate lists to append results
forces_i = []
forces_j = []
deformations = []
# run
curr_defo = ops.nodeDisp(1, 1)
# for each of the required displacement peaks
for target_defo in peaks:
# determine push direction
if curr_defo < target_defo:
incr = abs(incr)
sign = +1.00
else:
incr = -abs(incr)
sign = -1.00
# move towards the peak
while curr_defo * sign < target_defo * sign:
# determine increment
if abs(curr_defo - target_defo) < \
abs(incr):
# close to the peak, use what's remaining
incr_anl = sign * abs(curr_defo - target_defo)
else:
# away from peak, use max
incr_anl = incr
# assign analysis settings
ops.test('NormDispIncr', 1.0e-8, 25, 0)
ops.algorithm('Newton')
ops.integrator("DisplacementControl", 1, 1, incr_anl)
ops.system('FullGeneral')
ops.analysis("Static")
# execute and check
flag = ops.analyze(1)
if flag != 0:
# failed. report, stop and debug.
print('Analysis failed')
break
else:
# successful.
curr_defo = ops.nodeDisp(1, 1)
if use_zerolength:
# record reaction forces
ops.reactions()
forces_i.append(-ops.nodeReaction(0)[0]*area)
forces_j.append(-ops.nodeReaction(1)[0]*area)
# record displacement of node (1)
deformations.append(curr_defo*elm_len)
else:
# record reaction forces
ops.reactions()
forces_i.append(-ops.nodeReaction(0)[0])
forces_j.append(-ops.nodeReaction(1)[0])
# record displacement of node (1)
deformations.append(curr_defo)
force_vec_i = np.array(forces_i)/1e3
force_vec_j = np.array(forces_j)/1e3
displacement = np.array(deformations)
plt.plot(displacement, force_vec_i, label='Free node reaction')
plt.plot(displacement, force_vec_j, label='Restrained node reaction')
plt.grid()
plt.legend()
plt.xlabel('Displacement [in]')
plt.ylabel('Force [kips]')
plt.tight_layout()
plt.show()