I have been trying to model a simple structural system for time history analysis, where the mass is concentrated and the columns are represented by spring elements. I tried using two-node link elements for the columns and defining a hysteresis material to represent the overall structural properties.
I am able to produce results, but they are not the expected results. The maximum inter-story displacement angle of my structure is not the expected typical "C" pattern.
Can someone tell me how to model my structure correctly?
import numpy as np
from openseespy.opensees import *
import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import MCS_model_parameter
import NMFS_model_parameter
import opstool as opst
###### 结果输出位置
dir = "Data-2a"
if not os.path.exists(dir):
os.makedirs(dir)
######## 参数输入
NumOfStories = 5
GMfile='EL.txt'
GMfact = 0.45
FloorArea = [495] * NumOfStories # m2
B = 16.5 ## width， ： m
StructuralType = 'C1' # Hazus table 5.1
StoryHeight = [3.3] * NumOfStories ## 楼层高度 unit ： m
StoryHeight[0] = 4.2
pga = 0.15 ## 设计加速度 单位g
year = 2010 ## 建造时间
changdi = ['二类', '第一分组'] ## 先输入【一零类，一一类，二类，三类，四类】，再输入设计分组【第一分组，第二分组，第三分组】，默认['二类','第二分组']
H = np.sum(StoryHeight) ## 结构总高度 ,单位 ： m)
## 地震动归一化
s1 = pd.read_table(GMfile, sep=' ', header=None)
s1 = s1 / np.max(np.abs(s1.values[:, 0])) * 9.8
s1 = pd.DataFrame(s1)
s1.to_csv(GMfile, sep=' ', index=0, header=0)
########
wipe()
model('basic', '-ndm', 2, '-ndf', 3)
## building model
node(0, 0, 0)
node(1, 17, 0)
fix(0, 1, 1, 1) ##剪切弹簧
fix(1, 1, 1, 1) ##弯曲弹簧
# geomTransf('PDelta',2) ##
####计算自振周期
if StructuralType[0] == 'C': # concrete
T1 = 0.25 + 0.00053 * H ** 2 / B ** (1 / 3) ## 周期：s
elif StructuralType[0] == 'S': # steel
T1 = 0.554 + 0.00061 * H ** 2 / B ** (1 / 3)
else:
print('无法计算自振周期')
print('结构经验自振周期：', T1)
##### 计算骨架曲线参数
m, s1p, e1p, s2p, e2p, s3p, e3p, Ms1p, Me1p, Ms2p, Me2p, Ms3p, Me3p, pinchX, pinchY, damage1, damage2, beta = NMFS_model_parameter.NMFS_model_parameter(
NumOfStories, FloorArea, T1, StructuralType, StoryHeight, pga, year, changdi)
##### 生成节点
m = m / 2
for i in range(NumOfStories):
node(i * 2 + 2 , 0, np.sum(StoryHeight[0:i + 1]).tolist())
node(i * 2 + 3 , 17, np.sum(StoryHeight[0:i + 1]).tolist())
# equalDOF(i * 2 + 2, i * 2 + 3, 1, 2)
ja = m * (StoryHeight) ** 2 / 12
mass(2 * i + 2, m.tolist(), m.tolist(), 0)
mass(2 * i + 3, m.tolist(), m.tolist(), 0)
fix(i*2+2, 0, 1, 1) ## 剪切弹簧
fix(i*2+3, 0, 1, 0) ## 弯曲弹簧
# material
for i in range(NumOfStories):
# definition of uniaxialMaterial Hysteretic
# uniaxialMaterial Hysteretic $matTag $s1p $e1p $s2p $e2p <$s3p $e3p> $s1n $e1n $s2n $e2n <$s3n $e3n> $pinchX $pinchY $damage1 $damage2 <$beta>
uniaxialMaterial('Hysteretic', i + 1,
s1p, e1p, s2p, e2p, s3p[i], e3p[i],
-s1p[i], -e1p[i], -s2p[i], -e2p[i], -s3p[i], -e3p[i], pinchX, pinchY,
damage1, damage2, beta)
uniaxialMaterial('Elastic', i + 1000, 1e20)
uniaxialMaterial('Hysteretic', i + 10000,
Ms1p[i], Me1p[i], Ms2p[i], Me2p[i], Ms3p[i], Me3p[i],
-Ms1p[i], -Me1p[i], -Ms2p[i], -Me2p[i], -Ms3p[i], -Me3p[i], pinchX, pinchY,
damage1, damage2, beta)
uniaxialMaterial('Elastic', i + 5000, 0 )
# element twoNodeLink $eleTag $iNode $jNode -mat $matTags -dir $dirs <-orient <$x1 $x2 $x3> $y1 $y2 $y3> <-pDelta (4 $Mratio)> <-shearDist (2 $sDratios)> <-doRayleigh> <-mass $m>
for i in range(NumOfStories):
element('twoNodeLink', 2*i, 2*i, 2*i + 2, '-mat', i + 1000, i + 1, i + 10000,'-dir', 1, 2, 3, '-pDelta', 0, 0, 0,
0) ### 1是纵向弹簧 2是剪切弹簧 3是转动弹簧 ###偶数编号的单元 模拟剪切弹簧
element('twoNodeLink', 2*i+1, 2*i+1, 2*i + 3, '-mat', i + 1000, i + 5000, i + 10000, '-dir', 1, 2, 3, '-pDelta', 0, 0, 0,
0) ### 1是纵向弹簧 2是剪切弹簧 3是转动弹簧 ###奇数编号的单元 模拟弯曲弹簧
# element('Truss', i + 1000, 2*i+2, 2*i+3, 0.001, i + 1000)
##刚性链杆
rigidLink('beam',2*i+2, 2*i+3)
Tol = 1.0e-6
#####################加载重力
timeSeries('Linear', 1) ## 时间序列 线性 ，时间序列编号为 1
pattern('Plain', 1, 1) ## 加载模式为'Plain'，编号为1， 要在荷载模式中使用的时间序列的标记为 1，用于静力分析
for i in range(NumOfStories):
load(2*i+2, 0.0, (-m[i]*9.
.tolist(), 0.0)constraints('Transformation') # 边界约束方程的处理方式
numberer('RCM') ## 结构自由度编号方式
system('BandGeneral') ## 便表示方程的球结合存储，
test( 'NormDispIncr', Tol, 6) # determine if convergence has been achieved at the end of an iteration step
algorithm( 'Newton') ## 牛顿迭代法
integrator ('LoadControl', 0.1) ## 力加载控制方式与系统时间步长，每次加载0.1倍的外力
analysis ('Static') #静力分析
analyze (10) #分析总步长
wipeAnalysis()
loadConst ('-time' ,0.0 )# “loadConst -time 0.0”是很必需的，这个命令可以使重力负载保持恒定（锁定重力荷载）并使系统时间重置为零，在此基础上然后开始下一阶段的分析。
print("gravity is ok!")
###################### dynamic
######### 定义输出，这部分输出只包含动力分析的
for i in range(NumOfStories):
recorder('Node', '-file', 'Data-2a/DFree'+str(2*i+2)+'.txt', '-time', '-node', 2*i+2, '-dof', 1, 'disp') ## ‘-accel’
# define & apply damping
# RAYLEIGH damping parameters, Where to put M/K-prop damping, switches
# (http://opensees.berkeley.edu/OpenSees/m ... l/1099.htm)
# D=$alphaM*M + $betaKcurr*Kcurrent + $betaKcomm*KlastCommit + $beatKinit*$Kinitial
if NumOfStories > 1:
lambdaN = eigen('-fullGenLapack', 2) ## 振型分析
w1 = lambdaN[0] ** 0.5
w2 = lambdaN[1] ** 0.5
T1 = 2.0 * np.pi / w1
T2 = 2.0 * np.pi / w2
# if ifprint:
# print(f'Eigen Analysis: T1 = {T1:.2f} s; T2 = {T2:.2f} s')
else:
lambdaN = eigen('-fullGenLapack', 1)
w1 = lambdaN[0] ** 0.5
T1 = 2.0 * np.pi / w1
# if ifprint:
# print(f'Eigen Analysis: T1 = {T1:.2f} s')
print('结构自振周期：',T1)
DampingRatio=0.05
if NumOfStories > 1:
xDamp = DampingRatio
MpropSwitch = 1.0 ## 详见崔继东博客
KcurrSwitch = 0.0 ## 当前刚度举证
KinitSwitch = 1.0 ## 初始刚度
KcommSwitch = 0.0 ## 上一步的刚度
nEigenI = 1
nEigenJ = 2
lambdaI = lambdaN[nEigenI - 1]
lambdaJ = lambdaN[nEigenJ - 1]
omegaI = lambdaI ** 0.5
omegaJ = lambdaJ ** 0.5
alphaM = MpropSwitch * xDamp * (2.0 * omegaI * omegaJ) / (omegaI + omegaJ)
betaKcurr = KcurrSwitch * 2. * xDamp / (omegaI + omegaJ) # current-K; +beatKcurr*KCurrent
betaKinit = KinitSwitch * 2. * xDamp / (omegaI + omegaJ) # initial-K; +beatKinit*Kini
betaKcomm = KcommSwitch * 2. * xDamp / (omegaI + omegaJ) # last-committed K; +betaKcomm*KlastCommitt
rayleigh(alphaM, betaKcurr, betaKinit, betaKcomm)
else:
xDamp = DampingRatio
MpropSwitch = 1.0
nEigenI = 1
lambdaI = lambdaN[nEigenI - 1]
omegaI = lambdaI ** 0.5
alphaM = MpropSwitch * xDamp * omegaI * 2 ## 只有一阶振型时，只使用质量比例阻尼
rayleigh(alphaM, 0, 0, 0 )
rayleigh(alphaM, 0, 0, 0 )
#############
IDloadTag=2 ## 荷载编号
GMdirection= 1 ## 荷载方向
dt=0.02 ## 地震动时间间隔
print('地震动放大因子:',GMfact)
timeSeries( 'Path',2, '-dt', dt, '-filePath', GMfile, '-factor', GMfact)
pattern('UniformExcitation', IDloadTag, GMdirection, '-accel', 2) ## 用于动力分析， '-accel', 2代表加速度的路径，'UniformExcitation'代表基底一致激励
# dynamic analysis
constraints('Transformation') # 边界约束方程的处理方式
numberer('RCM') ## 结构自由度编号方式
system('BandGeneral')
test ('EnergyIncr',Tol,10) ########### 用位移增量判断收敛，后面的数据代表精度（1.0e-6），最大迭代数（10）
algorithm('Newton')
integrator ('Newmark', 0.5,0.25) ## 加载方式
analysis ('Transient') ## 瞬态分析法
DtAnalysis = 0.02 ## 每隔0.02s输出一个点
TmaxAnalysis = (np.size(pd.read_table(GMfile))-1)*dt ## 地震总时长
Nsteps = int((TmaxAnalysis) / DtAnalysis)
ok = analyze(Nsteps, DtAnalysis)
########### 用位移增量判断收敛，
test1 = {1: 'NormDispIncr', 2: 'RelativeEnergyIncr', 4: 'RelativeNormUnbalance', 5: 'RelativeNormDispIncr',
6: 'NormUnbalance'}
algorithm1 = {1: 'KrylovNewton', 2: 'SecantNewton', 4: 'RaphsonNewton', 5: 'PeriodicNewton', 6: 'BFGS', 7: 'Broyden',
8: 'NewtonLineSearch'}
for i in test1:
for j in algorithm1:
if ok != 0:
if j < 4:
algorithm(algorithm1[j], '-initial')
else:
algorithm(algorithm1[j])
test(test1[i], Tol, 1000)
ok = analyze(Nsteps, DtAnalysis)
print(test1[i], algorithm1[j], ok)
if ok == 0:
break
else:
continue
wipe()