Inelastic SDOF

6 min read • 1,087 words

Integration of an inelastic single-degree-of-freedom (SDOF) system.

March 2020, By Amir Hossein Namadchi

This is an OpenSeesPy simulation of a simple SDOF system with elastoplastic behavior mentioned in Dynamics of Structures book by Ray W. Clough and J. Penzien. This example has been solved in the book, so the result obtained here can be compared with the reference.

This notebook is adapted from https://github.com/AmirHosseinNamadchi/OpenSeesPy-Examples/blob/master/Elastoplastic%20SDOF%20system.ipynb

import numpy as np
import opensees.openseespy as ops
import matplotlib.pyplot as plt

## Units
inch = 1              # inches
kips = 1              # KiloPounds
sec = 1               # Seconds

lb = kips*(sec**2)/inch    # mass unit (derived)

Model specifications are defined as follows:

m = 0.1*lb              # Mass
k = 5.0*(kips/inch)     # Stiffness
c = 0.2*(kips*sec/inch) # Damping
dy_p = 1.2*inch         # Plastic state displacment
alpha_m = c/m           # Rayleigh damping ratio

# Variation of p(t) in tabular form
load_history = np.array([[0,   0],
                         [0.1, 5],
                         [0.2, 8],
                         [0.3, 7],
                         [0.4, 5],
                         [0.5, 3],
                         [0.6, 2],
                         [0.7, 1],
                         [0.8, 0]])

# Dynamic Analysis Parameters
dt = 0.01
time = 1.0

Analysis

Let’s wrap the whole part in a function so that different material behavior could be passed to the function:

def do_analysis(dt, time, material_params):
    
    model = ops.Model(ndm=1, ndf=1)
    time_domain = np.arange(0, time, dt)

    # Nodes
    model.node(1,0.0,0.0)
    model.node(2,0.0,0.0)

    model.uniaxialMaterial(*material_params)
    model.element('ZeroLength', 1, *[1,2], mat=1, dir=1)
    model.mass(2, m)
    model.rayleigh(alpha_m, 0.0, 0.0, 0.0)

    model.fix(1,1)
    model.timeSeries('Path', 1,
                     values=load_history[:,1],
                     time=load_history[:,0])
    
    model.pattern('Plain', 1, 1)
    model.load(2, 1)
    
    # Analysis
    model.constraints('Plain')
    model.numberer('Plain')
    model.system('ProfileSPD')
    model.test('NormUnbalance', 0.0000001, 100)
    model.algorithm('ModifiedNewton')
    model.integrator('Newmark', 0.5, 0.25)
    model.analysis('Transient')
    
    time_lst =[0]           # list to hold time stations for plotting
    response = [0]          # response params of node 2     
    

    for i in range(len(time_domain)):
        model.analyze(1, dt)
        time_lst.append(model.getTime())
        response.append(model.nodeDisp(2,1))
    
    return {'time_list':np.array(time_lst),
            'd': np.array(response)}

For comparison (similar to the book), elastic analysis is also inculded:

epp = do_analysis(dt, time, ['ElasticPP', 1, k, dy_p])     # Elastic-Perfectly Plastic
els = do_analysis(dt, time, ['Elastic', 1, k])             # Elastic

Visualization

plt.figure(figsize=(7,5))
plt.plot(epp['time_list'], epp['d'], color = '#fe4a49', linewidth=1.75, label = 'Nonlinear (EPP)')
plt.plot(els['time_list'], els['d'], color = '#2ab7ca', linewidth=1.75, label = 'Linear (Elastic)')
plt.ylabel('Displacement (inch)', {'fontstyle':'italic','size':14})
plt.xlabel('Time (sec)', {'fontstyle':'italic','size':14})
plt.xlim([0.0, time])
plt.legend()
plt.grid()
plt.yticks(fontsize = 14)
plt.xticks(fontsize = 14);

Output:

<Figure size 700x500 with 1 Axes>

Closure

Results obtained here with OpenSees perfectly agree with the ones in the book.

References

Clough, R.W. and Penzien, J., 2003. Dynamics of structures. Berkeley. CA: Computers and Structures, Inc.

$head

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175


import sdof
import numpy as np

import opensees.openseespy as op
FREE  = 0
FIXED = 1
X, Y, RZ = 1, 2, 3


def plastic_sdof(material, motion, dt, xi=0.05, r_post=0.0):
    """
    Run seismic analysis of a nonlinear SDOF

    :param mass: mass
    :param k: spring stiffness
    :param f_yield: yield strength
    :param motion: list, acceleration values
    :param dt: float, time step of acceleration values
    :param xi: damping ratio
    :param r_post: post-yield stiffness
    :return:
    """
    mass, k, f_yield = material
    op.wipe()
     # 2 dimensions, 3 dof per node
    op.model('basic', '-ndm', 2, '-ndf', 3) 

    # Establish nodes
    bot_node = 1
    top_node = 2
    op.node(bot_node, 0., 0.)
    op.node(top_node, 0., 0.)

    # Fix bottom node
    op.fix(top_node, FREE,  FIXED, FIXED)
    op.fix(bot_node, FIXED, FIXED, FIXED)

    # Set out-of-plane DOFs to be slaved
    op.equalDOF(1, 2, *[2, 3])

    # nodal mass (weight / g):
    op.mass(top_node, mass, 0., 0.)

    # Define material
    bilinear_mat_tag = 1
    mat_type = "Steel01"
    mat_props = [f_yield, k, r_post]
    op.uniaxialMaterial(mat_type, bilinear_mat_tag, *mat_props)

    # Assign zero length element
    beam_tag = 1
    op.element('zeroLength', beam_tag, bot_node, top_node, "-mat", bilinear_mat_tag, "-dir", 1, '-doRayleigh', 1)

    # Define the dynamic analysis
    load_tag_dynamic = 1
    pattern_tag_dynamic = 1

    values = list(-1 * motion)  # should be negative

#   op.timeSeries('Path', load_tag_dynamic, dt=dt, values=values)
    op.timeSeries('Path', load_tag_dynamic, "-dt", dt, "-values", *values)

#   op.pattern('UniformExcitation', pattern_tag_dynamic, X, accel=load_tag_dynamic)
    op.pattern('UniformExcitation', pattern_tag_dynamic, X, "-accel", load_tag_dynamic)

    # set damping based on first eigen mode
    eig = op.eigen('-fullGenLapack', 1)
    try:
        angular_freq = eig**0.5
    except:
        angular_freq = eig[0]**0.5
    alpha_m = 0.0
    beta_k = 2 * xi / angular_freq
    beta_k_comm = 0.0
    beta_k_init = 0.0

    op.rayleigh(alpha_m, beta_k, beta_k_init, beta_k_comm)

    # Run the dynamic analysis

#   op.wipeAnalysis()

    op.algorithm('Newton')
#   op.system('SparseGeneral')
    op.numberer('RCM')
    op.constraints('Transformation')
    op.integrator('Newmark', 0.5, 0.25)
    op.analysis('Transient')

    tol = 1.0e-10
    iterations = 10
    op.test('EnergyIncr', tol, iterations, 0, 2)
    analysis_time = (len(values) - 1) * dt
    analysis_dt = 0.001
    outputs = {
        "time": [],
        "rel_disp": [],
        "rel_accel": [],
        "rel_vel": [],
        "force": []
    }

    while op.getTime() < analysis_time:
        curr_time = op.getTime()
        if op.analyze(1, analysis_dt) != 0:
            print(f"Failed at time {op.getTime()}")
            break

        outputs["time"].append(curr_time)
        outputs["rel_disp"].append(op.nodeDisp(top_node, 1))
        outputs["rel_vel"].append(op.nodeVel(top_node, 1))
        outputs["rel_accel"].append(op.nodeAccel(top_node, 1))
        op.reactions()
        outputs["force"].append(-op.nodeReaction(bot_node, 1))  # Negative since diff node

    op.wipe()
    for item in outputs:
        outputs[item] = np.array(outputs[item])

    return outputs


def main():
    """
    Create a plot of an elastic analysis, nonlinear analysis and closed form elastic

    :return:
    """
    import eqsig
    import matplotlib.pyplot as plt

    record_filename = 'test_motion_dt0p01.txt'
    dt = 0.01
    rec = np.loadtxt(record_filename)
    acc_signal = eqsig.AccSignal(rec, dt)
    period = 1.0
    xi = 0.05
    mass = 1.0
    f_yield = 1.5  # Reduce this to make it nonlinear
    r_post = 0.0

    periods = np.array([period])

    k = 4 * np.pi ** 2 * mass / period ** 2

    outputs = plastic_sdof((mass, k, f_yield), rec, dt, xi=xi, r_post=r_post)
    outputs_elastic = plastic_sdof((mass, k, f_yield * 100), rec, dt, xi=xi, r_post=r_post)
    ux_opensees = outputs["rel_disp"]
    ux_opensees_elastic = outputs_elastic["rel_disp"]
    print(outputs)

    bf, sps = plt.subplots(nrows=2)
    sps[0].plot(outputs["time"], ux_opensees, label="OpenSees fy=%.3gN" % f_yield, ls="-")
    sps[0].plot(outputs["time"], ux_opensees_elastic, label="OpenSees fy=%.3gN" % (f_yield * 100), ls="--")
    time = acc_signal.time
    acc_opensees_elastic = np.interp(time, outputs_elastic["time"], outputs_elastic["rel_accel"]) - rec


    resp_u, resp_v, resp_a = sdof.integrate(rec, dt, k, 2*xi*mass*np.sqrt(k/mass), mass, fy=f_yield)
    sps[0].plot(acc_signal.time, resp_u, label="sdof")
    sps[1].plot(acc_signal.time, resp_a, label="sdof")

#   resp_u, resp_v, resp_a = duhamels.response_series(motion=rec, dt=dt, periods=periods, xi=xi)
#   sps[0].plot(acc_signal.time, resp_u[0], label="Eqsig")
#   sps[1].plot(acc_signal.time, resp_a[0], label="Eqsig")  # Elastic solution
#   print("diff", sum(acc_opensees_elastic - resp_a[0]))
    sps[1].plot(time, acc_opensees_elastic, label="Opensees fy=%.2gN" % (f_yield * 100), ls="--")
    sps[0].legend()
    sps[1].legend()
    plt.show()


if __name__ == '__main__':
    main()

Inelastic Cantilever

Linearized Buckling

On this page: