The OpenSees Community

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