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)