https://openseespydoc.readthedocs.io/en ... anti2.html
Code: Select all
import math
import os
from openseespy.opensees import *
wipe()
a = os.path.exists('recorder')
b = os.getcwd()
if a != True:
os.mkdir(b + '\\recorder')
else:
pass
model('basic', '-ndm', 3, '-ndf', 6)
PierHeight = 6.5
ElementLength = 0.5
D = 1.2
def model_make(PierHeight, ElementLength, D):
Area = 3.14 * (D ** 2) / 4
g = 9.8
N_PierNode = int((PierHeight/ElementLength)+1)
print(f'node number:{N_PierNode}')
PierNodalMass = 2.5 * Area * ElementLength
for i in range(1, N_PierNode+1):
PierNode = i
PierCoorX = 0
PierCoorY = 0
PierCoorZ = ElementLength * (i - 1)
node(PierNode, PierCoorX, PierCoorY, PierCoorZ)
mass(PierNode, PierNodalMass, PierNodalMass, PierNodalMass, 0.0, 0.0, 0.0)
DeckNode = N_PierNode + 1
DeckCoorX = 0
DeckCoorY = 0
DeckCoorZ = (N_PierNode-1)*ElementLength
node(DeckNode, DeckCoorX, DeckCoorY, DeckCoorZ)
DeckMass = 65000 / g
mass(DeckNode, DeckMass, DeckMass, DeckMass, 0.0, 0.0, 0.0)
fix(14, 1, 1, 1, 1, 1, 1)
print(f'node build is done!')
import fiber
import material
SecTag = 1
D = 1.2
a = 0.02
Cover = material.HRB400[0]
Core = material.C40_gg[0]
fiber.Circle1(SecTag, D, a, Cover, Core)
print(f'fiber build is done!')
AggSecTag1 = 3
section('Aggregator', AggSecTag1, 1001,'T', '-section', SecTag)
ColTransfTag = 1
geomTransf('Linear', ColTransfTag, 2, 1, 1)
numIntgrPts = 4
for i in range(1,N_PierNode):
eleTag = i
element('nonlinearBeamColumn', eleTag, i, i + 1, numIntgrPts, SecTag, ColTransfTag)
equalDOF(14, 15, 1, 2, 3, 4, 5, 6)
print(f'element build is done!')
recorder('Node', '-file', 'recorder/disp1.out', '-time', '-node', 1, '-dof', 1, 2, 3, 'disp')
recorder('Node', '-file', 'recorder/disp5.out', '-time', '-node', 5, '-dof', 1, 2, 3, 'disp')
recorder('Node', '-file', 'recorder/disp12.out', '-time', '-node', 12, '-dof', 1, 2, 3, 'disp')
recorder('Node', '-file', 'recorder/reaction1.out', '-time', '-node', 1, '-dof', 1, 2, 3, 'reaction')
recorder('Node', '-file', 'recorder/reaction5.out', '-time', '-node', 5, '-dof', 1, 2, 3, 'reaction')
recorder('Node', '-file', 'recorder/reaction12.out', '-time', '-node', 12, '-dof', 1, 2, 3, 'reaction')
recorder('Element', '-file', 'recorder/globalForce1.out', '-time', '-ele', 1, 'globalForce')
recorder('Element', '-file', 'recorder/globalForce5.out', '-time', '-ele', 5, 'globalForce')
recorder('Element', '-file', 'recorder/globalForce12.out', '-time', '-ele', 12, 'globalForce')
recorder('Element', '-file', 'recorder/Force1.out', '-time', '-ele', 1, 'Force')
recorder('Element', '-file', 'recorder/Force5.out', '-time', '-ele', 5, 'Force')
recorder('Element', '-file', 'recorder/Force12.out', '-time', '-ele', 12, 'Force')
recorder('Element', '-file', 'recorder/deformations1.out', '-time', '-ele', 1, 'deformations')
recorder('Element', '-file', 'recorder/deformations5.out', '-time', '-ele', 5, 'deformations')
recorder('Element', '-file', 'recorder/deformations12.out', '-time', '-ele', 12, 'deformations')
timeSeries('Linear', 1)
pattern('Plain', 1, 1)#'-fact', fact
PierNodalMass = 2.5*Area*ElementLength
for i in range(1,N_PierNode+1):
load(i, 0.0, 0.0, -PierNodalMass*g, 0.0, 0.0, 0.0)
load(N_PierNode, 0.0, 0.0, DeckMass*g, 0.0, 0.0, 0.0)
print('Gravity analysis is done!')
print('Model Built!')
model_make(PierHeight, ElementLength, D)
# -*- coding: utf-8 -*-
from openseespy.opensees import *
import math
wipe()
import model
# 应用动态地震动分析
Tol = 1e-8
GMdirection = 1
GMfile = 'acc2.txt'
GMfact = 1.0
Lambda = eigen('-fullGenLapack', 1) # eigenvalue mode 1
Omega = math.pow(Lambda[0], 0.5)
betaKcomm = 2 * (0.02/Omega)
xDamp = 0.02 # 阻尼比 2%
alphaM = 0.3785 # M-prop.阻尼;D = alphaM * M
betaKcurr = 0.0 # k -比例阻尼;+ beatKcurr * KCurrent
betaKinit = 9.2e-4 # 初始刚度比例阻尼+beatKinit*Kini
rayleigh(alphaM,betaKcurr, betaKinit, betaKcomm) # 瑞利阻尼
# 均匀激励:加速度输入
IDloadTag = 400 # 加载号
dt = 0.02 # 输入地震动时间步长
GMfatt = 1.0 # 输入文件中的数据于 g Unifts ——加速度 TH
maxNumIter = 100
timeSeries('Path', 2, '-dt', dt, '-filePath', GMfile, '-factor', GMfact) # 离散荷载的输入
pattern('UniformExcitation', IDloadTag, GMdirection, '-accel', 2) # 统一激励模式
wipeAnalysis()
constraints('Transformation') # 此命令用于构造转换约束处理程序,该处理程序使用转换方法强制约束。下面是构造转换约束处理程序的命令
numberer('Plain') # 该命令用于构造一个Plain自由度编号对象,以提供节点的自由度与方程数之间的映射。一个Plain号只接受域给它的节点和它们编号的任何顺序,这个顺序既依赖于节点编号和模型的大小。
system('BandGeneral') # 该命令用于构造一个方程对象的BandGeneralSOE# 的线性系统。顾名思义,该类用于具有带状轮廓的矩阵系统。这个矩阵被存储在如下所示的一维数组中,其大小等于带宽乘以未知数的个数。当需要解决方案时,使用Lapack例程DGBSV和SGBTRS
test('EnergyIncr', Tol, maxNumIter) # 创建一个Energylncr测试,它使用解向量的点积和矩阵方程右边的范数来确定是否已经达到收敛
algorithm('ModifiedNewton') # 创建一个修正牛顿算法。与牛顿的不同之处在于,迭代中使用的是初始猜测处的切线,而不是当前的切线。
NewmarkGamma = 0.5
NewmarkBeta = 0.25
integrator('Newmark', NewmarkGamma, NewmarkBeta)
analysis('Transient')
DtAnalysis = 0.02 # 横向分析时步 dt
TmaxAnalysis = 30.0 # 地动分析的最大持续时间
Nsteps = int(TmaxAnalysis / DtAnalysis) # 步数
ok = analyze(Nsteps, DtAnalysis) # 执行分析,如果成功返回o,如果不成功返回<o
tCurrent = getTime() # 返回域中的当前时间
test = {1:'NormDispIncr', 2: 'RelativeEnergyIncr', 4: 'RelativeNormUnbalance',5: 'RelativeNormDispIncr', 6: 'NormUnbalance'}
algorithm = {1:'KrylovNewton', 2: 'SecantNewton' , 4: 'RaphsonNewton',5: 'PeriodicNewton', 6: 'BFGS', 7: 'Broyden', 8: 'NewtonLineSearch'}
print(test[1])