Unexpected: Steel4 and non-zero reactions in unrestrained DOFs

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

Moderators: silvia, selimgunay, Moderators

Post Reply
ioannisVM
Posts: 3
Joined: Tue Jul 07, 2020 11:46 am
Location: Berkeley, CA
Contact:

Unexpected: Steel4 and non-zero reactions in unrestrained DOFs

Post by ioannisVM »

Hi everyone,

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

Post Reply