SDOF with nonlinear stiffness, vertical excitation

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

Moderators: silvia, selimgunay, Moderators

Post Reply
harrieweijs
Posts: 1
Joined: Tue Nov 01, 2022 6:53 am

SDOF with nonlinear stiffness, vertical excitation

Post by harrieweijs » Thu Dec 15, 2022 1:11 am

Hi all,

For my master's thesis, I want to solve an SDOF system with a non-linear stiffness relationship, which is exited vertically.
After googling, I came across an example by Maxim Millen, which is very similar to what I want to do. I cannot share the link as it is a direct link to a download file, but if you google 'example_inelastic_sdof - ResearchGate', you should find it.

The significant difference between that script and mine is that his script is excited horizontally, whereas mine is excited vertically.
So, besides changing the material input, I tried to change the relevant points from X to Y.
However, if I go and change the points of interest from X to Y, the calculation seems to diverge (I get a displacement of over 90 meters; my expected output should be in the magnitude of 10-40 mm).

My approach thus seems not to be the right one to obtain the solution.
I have also tried to keep the system orientation in the X-direction and apply an additional lateral acceleration of g/2 = 9.81/2 m/s^2 (since the nodal mass is 0.5 * mass). This does return more realistic values, albeit on the high side.

Since my experience with OpenSEES is very limited, I am not sure which approach is the better one. I was hoping one of you would be willing to look alongside me!

Below you find the script.

Kind regards,

Harrie Weijs

Code: Select all

import numpy as np
import os
import pandas as pd
import o3seespy as o3
import matplotlib.pyplot as plt


def get_inelastic_response(mass, k_spring, f_yield, motion, dt, xi, r_post):

    osi = o3.OpenSeesInstance(ndm=2, state=0)

    # Establish nodes.
    bot_node = o3.node.Node(osi, x=0., y=0.)
    top_node = o3.node.Node(osi, x=0., y=0.)

    # Fix bottom node.
    o3.Fix3DOF(osi, node=top_node, x=o3.cc.FIXED, y=o3.cc.FREE, z_rot=o3.cc.FIXED)
    o3.Fix3DOF(osi, node=bot_node, x=o3.cc.FIXED, y=o3.cc.FIXED, z_rot=o3.cc.FIXED)

    # Set out-of-plane DOFs to be constrained.
    o3.EqualDOF(osi, r_node=top_node, c_node=bot_node, dofs=[o3.cc.Y, o3.cc.DOF2D_ROTZ])

    # Nodal mass.
    o3.Mass(osi, node=top_node, x_mass=0., y_mass=mass, rot_mass=0.)

    # Define material.
    bilinear_mat = o3.uniaxial_material.BoucWen(osi, alpha=0.09, ko=41500000, n=1.0, gamma=-0.30, beta=0.10, ao=0.0, delta_a=0.0, delta_nu=0.0, delta_eta=0.0)

    # Assign zero length element. Pass actual node and material objects into element.
    o3.element.ZeroLength(osi, ele_nodes=[top_node, bot_node], mats=[bilinear_mat], dirs=[o3.cc.DOF2D_Y], r_flag=0)

    # Define the dynamic analysis
    acc_series = o3.time_series.Path(osi, dt=dt, values=-motion)
    o3.pattern.UniformExcitation(osi, dir=o3.cc.Y, accel_series=acc_series)


    # Run the dynamic analysis
    o3.wipe_analysis(osi)
    o3.algorithm.Newton(osi)
    o3.numberer.RCM(osi)
    o3.constraints.Transformation(osi)
    o3.integrator.Newmark(osi, gamma=0.5, beta=0.25)
    o3.analysis.Transient(osi)
    o3.test_check.EnergyIncr(osi, tol=1.0e-10, max_iter=10)
    analysis_time = (len(motion) - 1) * dt
    analysis_dt = dt
    outputs = {
        "time": [],
        "displ": [],
        "acc": [],
        "vel": [],
        "force": []
    }

    while o3.get_time(osi) < analysis_time:
        o3.analyze(osi, num_inc=1, dt=analysis_dt)
        curr_time = o3.get_time(osi)
        outputs["time"].append(curr_time)
        outputs["displ"].append(o3.get_node_disp(osi, node=top_node, dof=o3.cc.Y))
        outputs["vel"].append(o3.get_node_vel(osi, node=top_node, dof=o3.cc.Y))
        outputs["acc"].append(o3.get_node_accel(osi, node=top_node, dof=o3.cc.Y))
        o3.gen_reactions(osi)
        outputs["force"].append(-o3.get_node_reaction(osi, node=bot_node, dof=o3.cc.Y))

    o3.wipe(osi)

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

    return outputs


def calc_sdof(record_filename, show=0):

    # Load acceleration signal.
    record_filename = record_filename
    rec = pd.read_csv(record_filename, header=None)[:][0].to_numpy()
    # Input calculation parameters
    dt = 0.01
    xi = 0.05
    mass = 14227
    f_yield = 249 * 10**3
    r_post = ((376 - 249) / (40 - 6)) * 10**6
    k_spring = 4.15 * 10**7
    outputs = get_inelastic_response(mass=mass, k_spring=k_spring, f_yield=f_yield, motion=rec, dt=dt, xi=xi, r_post=r_post)

    if show:
        plt.figure()
        plt.grid(which='major', color='#666666', linestyle='-')
        plt.minorticks_on()
        plt.grid(which='minor', color='#999999', linestyle='-', alpha=0.2)
        plt.plot(outputs["time"], outputs["displ"] * 1000, label='Inelastic response', lw=0.7, c='r')
        plt.ylabel('Disp. [mm]')
        plt.legend(loc='upper right')
        print("")
        print("max displacement = {0:.3f} mm".format(max(outputs["displ"]) * 1000))
        print("max force = {0:.3f} kN".format(max(outputs["force"])))
        plt.show()


if __name__ == '__main__':
    calc_sdof(record_filename=r'\filepath', show=1)

Post Reply