opensees吧 关注:986贴子:3,239
  • 0回复贴,共1

请问有大佬知道建模过程哪里有问题吗,第一次使用openseespy建模

只看楼主收藏回复

import openseespy.opensees as ops
from typing import Literal
import numpy as np
print("开始钢筋混凝土简支T梁桥地震分析示例")
ops.wipe()
# 单位规定:kN, m, s
# 创建ModelBuilder(三维模型,每个节点6个自由度)
ops.model('basic', '-ndm', 3, '-ndf', 6)
#材料参数
fc = 40.0 # 混凝土抗压强度(MPa)
E_concrete = 30000.0 # 混凝土弹性模量
density=25 # 混凝土密度(kN/m^3)
fy = 400.0 # 钢筋屈服强度(MPa)
Es = 200e6 # 钢筋弹性模量
b = 0.01 # 钢筋硬化率
g = 9.81 # 重力加速度(m/s^2)
#T梁信息
L = 30 # 跨度长度(m)
num=5 # T梁数目
wing_width = 3 # 翼板宽度和T梁间隔
wing_thickness = 0.2 # 翼板厚度
web_height = 2.3 # 腹板总高度
web_width_top = 0.2 # 上部腹板宽度
web_width_bottom = 0.6 # 底部腹板宽度(加宽部分)
web_transition_height = 0.2 # 过渡区域的高度
web_thickness_bottom = 0.2 # 底部加宽的高度
#墩柱信息
H = 15 # 墩高(m)
r=1 # 墩柱的直径
space=10 # 两墩的间距
#分段信息
nL = 10 #长度方向的分段数
nH = 5 #高度方向的分段数
# 创建墩的节点并添加质量
y = 0.0 # 两墩连线中点y坐标
M_pier = np.pi * r * r * H * density # 单根墩柱总质量
# 左墩1、2
for i in range(2): # i=0 左墩1,i=1 左墩2
y = -5 if i == 0 else 5
for j in range(nH + 1):
z = j * H / nH
node_tag = (i + 1) * 100 + j + 1
ops.node(node_tag, 0, y, z)
# 节点质量分配
if j == 0 or j == nH:
m = M_pier / (2 * nH) # 顶部和底部节点质量为总质量的一半
else:
m = M_pier / nH # 中间节点质量
ops.mass(node_tag, m, m, m, 0.0, 0.0, 0.0)
# 右墩1、2
for i in range(2): # i=0 右墩1,i=1 右墩2
y = -5 if i == 0 else 5
for j in range(nH + 1):
z = j * H / nH
node_tag = (i + 3) * 100 + j + 1
ops.node(node_tag, L, y, z)
# 节点质量分配
if j == 0 or j == nH:
m = M_pier / (2 * nH) # 顶部和底部节点质量为总质量的一半
else:
m = M_pier / nH # 中间节点质量
ops.mass(node_tag, m, m, m, 0.0, 0.0, 0.0)
# 创建主梁的节点
s1=wing_width*wing_thickness
s2=web_width_top*(web_height-web_transition_height-web_thickness_bottom)
s3=0.5*(web_width_top+web_width_bottom)*web_transition_height
s4=web_width_bottom*web_thickness_bottom
s=s1+s2+s3+s4
M_girder=s*L*density
for n in range(num):
y=-6+3*n
for i in range(nL+1):
x = i * L / nL
z = H
node_tag = 1000 + n * (nL + 1) + i
ops.node(node_tag, x, y, z)
m = M_girder/nL/2 if (i == 0 or i == nL) else M_girder/nL
ops.mass(node_tag, m, m, m, 0.0, 0.0, 0.0)
# 固定墩柱底部
ops.fix(101, 1, 1, 1, 1, 1, 1)
ops.fix(201, 1, 1, 1, 1, 1, 1)
ops.fix(301, 1, 1, 1, 1, 1, 1)
ops.fix(401, 1, 1, 1, 1, 1, 1)
# 定义截面
# 定义非线性柱的材料
# 混凝土
# 约束核心混凝土
ops.uniaxialMaterial("Concrete01", 1, -41000, -0.004, -34470, -0.014)
# 非约束保护层混凝土
ops.uniaxialMaterial("Concrete01", 2, -34470, -0.002, -25000, -0.006)
# 钢筋材料
ops.uniaxialMaterial("Steel01", 3, fy*1e3, Es, b)
# 墩柱截面
# 墩柱尺寸
D = 1.0 # 墩柱总直径
core_diameter = 0.8 # 核心混凝土直径
cover_thickness = (D - core_diameter) / 2 # 保护层厚度
steel_area = 0.02 # 单根钢筋面积
num_bars = 8 # 周围均匀分布的钢筋数量
bar_radius = core_diameter / 2 # 钢筋半径位置
# 创建墩柱纤维截面
ops.section("Fiber", 991, "-GJ", 1e10)
# 定义混凝土纤维
core_area = np.pi * (core_diameter / 2) ** 2
ops.patch("circ", 1, 16, 16, 0.0, 0.0, 0.0, core_diameter / 2, 0.0, 2.0 * np.pi) # 材料Tag 1 (核心混凝土)
ops.patch("circ", 2, 16, 16, 0.0, 0.0, core_diameter / 2, D / 2, 0.0, 2.0 * np.pi) # 材料Tag 2 (保护层混凝土)
# 定义钢筋纤维
for i in range(num_bars):
angle = 2 * np.pi * i / num_bars
x = bar_radius * np.cos(angle)
y = bar_radius * np.sin(angle)
ops.fiber(x, y, steel_area, 3) # 钢筋材料Tag 3
#创建墩柱截面
pier_sec = 1
ops.uniaxialMaterial("Elastic", 103, 1e10)
ops.section("Aggregator", pier_sec, *[103, "T"], "-section", 991)
wing_width = 3 # 翼板宽度和T梁间隔
wing_thickness = 0.2 # 翼板厚度
web_height = 2.3 # 腹板总高度
web_width_top = 0.2 # 上部腹板宽度
web_width_bottom = 0.6 # 底部腹板宽度(加宽部分)
web_transition_height = 0.2 # 过渡区域的高度
web_thickness_bottom = 0.2 # 底部加宽的高度
# 梁截面
# 混凝土材料
ops.uniaxialMaterial("Concrete01", 4, -40000, -0.004, -34470, -0.014)
# 创建T梁纤维截面
ops.section("Fiber", 992, "-GJ", 1e10)
# T梁翼板混凝土部分
ops.patch("rect", 4, 10, 4, -wing_width / 2, - wing_thickness, wing_width / 2, 0)
# 上腹板
ops.patch("rect", 4, 4, 10, -web_width_top / 2, -(wing_thickness + web_height - web_transition_height - web_thickness_bottom), web_width_top / 2, wing_thickness )
# 过渡区
ops.patch("rect", 4, 4, 4, -0.5 * (web_width_top + web_width_bottom), -(wing_thickness + web_height - web_thickness_bottom), 0.5 * (web_width_top + web_width_bottom), -(wing_thickness + web_height - web_transition_height - web_thickness_bottom))
# 底腹板
ops.patch("rect", 4, 4, 4, -web_width_bottom / 2, -(wing_thickness + web_height), web_width_bottom / 2, -(wing_thickness + web_height - web_thickness_bottom))
# 翼板钢筋布置(沿翼板宽度方向均匀分布)
num_steel_bars_wing = 10 # 翼板钢筋数量
for i in range(num_steel_bars_wing):
# 钢筋在翼板宽度方向的位置
x = -wing_width / 2 + i * (wing_width-0.05)/ (num_steel_bars_wing - 1)+0.025
# 钢筋在翼板厚度方向的位置
y = -wing_thickness / 2
ops.fiber(x, y, steel_area, 3)
# 腹板钢筋布置(在宽度方向两侧均匀分布)
num_steel_bars_web = 6 # 腹板钢筋数量
for i in range(num_steel_bars_web):
# 钢筋在腹板高度方向的位置
y = -(wing_thickness + i * (web_height - web_transition_height-web_thickness_bottom) / (num_steel_bars_web - 1))
# 钢筋在腹板宽度方向的位置,两侧对称分布
ops.fiber(-web_width_top / 2+0.025, y, steel_area, 3)
ops.fiber(web_width_top / 2-0.025, y, steel_area, 3)
# 底部抗拉钢筋布置
num_steel_bars_bottom = 6 # 假设底部钢筋数量为6根
y_position = -(wing_thickness + web_height) + 0.025 # 底腹板的底部位置加上保护层
steel_area = 0.02 # 每根钢筋的面积(假设值,可以根据实际需求调整)
# 沿着底腹板宽度布置钢筋
for n in range(2):
for i in range(num_steel_bars_bottom):
# 钢筋在底腹板宽度方向均匀分布
x = -web_width_bottom / 2 + i * ((web_width_bottom-0.05) / (num_steel_bars_bottom - 1))+0.025
if n==0:
y = y_position # 底腹板的底部位置加上保护层
else:
y = y_position+0.15 # 底腹板的底部位置加上保护层
# 将钢筋布置在底腹板区域(钢筋Tag 6)
ops.fiber(x, y, steel_area, 3)
#定义单元
# 桥墩单元
# 定义单元之前先要定义几何变换
ops.geomTransf('Linear', 1, 0, 1, 0)
# 设置元素长度方向的积分点数量(dispBeamColumn需要)
NP = 5
# 使用Lobatto积分,id为2
ops.beamIntegration('Lobatto', 2, pier_sec, NP)
# 使用塑性梁柱单元创建桥墩
# 左墩1、2
for i in range(2): # i=0 左墩1,i=1 左墩2
for j in range(nH): # 每个墩分成nH段
# 确定两个节点的标签
node_tag_start = (i + 1) * 100 + j + 1
node_tag_end = (i + 1) * 100 + j + 2
# 定义单元,使用先前定义的截面和几何变换
ops.element("dispBeamColumn", node_tag_start * 10 + j, node_tag_start, node_tag_end, 1, 2)
# 右墩1、2
for i in range(2): # i=0 右墩1,i=1 右墩2
for j in range(nH): # 每个墩分成nH段
# 确定两个节点的标签
node_tag_start = (i + 3) * 100 + j + 1
node_tag_end = (i + 3) * 100 + j + 2
# 定义单元,使用先前定义的截面和几何变换
ops.element("dispBeamColumn", node_tag_start * 10 + j, node_tag_start, node_tag_end, 1, 2)
#梁单元
ops.geomTransf('Linear', 2, 0, 1, 0)
NP = 5
# 使用Lobatto积分,id为3
ops.beamIntegration('Lobatto', 3, 992, NP)
for n in range(num): # num 为墩柱的数量
y = -6 + 3 * n # 计算y坐标的位置
for i in range(nL): # nL 为每一层的节点数量
# 获取相邻两个节点的标签
node_tag_start = 1000 + n * (nL + 1) + i # 起始节点标签
node_tag_end = 1000 + n * (nL + 1) + i + 1 # 结束节点标签
ops.element("dispBeamColumn", node_tag_start * 10 + i, node_tag_start, node_tag_end,2, 3)
#定义横向连接
# 定义横向连接的弹性材料(调整横向刚度值kx和短梁单元的刚度特性)
# 定义五片梁的横向连接
for i in range(num - 1): # num 是梁的数量
for j in range(nL + 1): # nL 是每片梁的节点数
node_tag1 = 1000 + i * (nL + 1) + j
node_tag2 = 1000 + (i + 1) * (nL + 1) + j
eleTag = 5000 + i * (nL + 1) + j
ops.element('elasticBeamColumn', eleTag, node_tag1, node_tag2, 0.86, 210e6, 23.2, 2.32, 81e6, 3.13, 3)
#定义支座
def define_bearings(type: Literal['Rubber', 'Frame'], num_beams: int = 5, nH: int = 0, nL: int = 0, M_girder: float = 0.0, rigid_tag: int = 1, free_tag: int = 2):
if type == 'Rubber':
for i in range(1, num_beams + 1):
# 使用 equalDOF 将梁和桥墩节点在所有自由度上进行约束
ops.equalDOF(100 + nH +1 , 1000 + i, *[1, 2, 3, 4, 5, 6])
ops.equalDOF(200 + nH +1 , 1000 + i, *[1, 2, 3, 4, 5, 6])
ops.equalDOF(300 + nH +1 , 1000 + nL + i, *[1, 2, 3, 4, 5, 6])
ops.equalDOF(400 + nH +1 , 1000 + nL + i, *[1, 2, 3, 4, 5, 6])
elif type == 'Frame':
for i in range(1, num_beams + 1):
# 使用 equalDOF 将梁和桥墩节点在所有自由度上进行约束
ops.equalDOF(100 + nH +1, 1000 + i, *[1, 2, 3, 4, 5, 6])
ops.equalDOF(200 + nH +1, 1000 + i, *[1, 2, 3, 4, 5, 6])
ops.equalDOF(300 + nH +1, 1000 + nL + i, *[1, 2, 3, 4, 5, 6])
ops.equalDOF(400 + nH +1, 1000 + nL + i, *[1, 2, 3, 4, 5, 6])
else:
raise NotImplementedError(f"支座类型{type}未实现")
#定义重力荷载
ops.timeSeries('Linear', 1)
ops.pattern('Plain', 1, 1)
# 为每个点创建荷载
for node in ops.getNodeTags():
# 读取节点质量,这里读取的是z向质量(dir=3)
P = ops.nodeMass(node,3)*g # 节点z向质量*重力加速度
if P > 0:
ops.load(node, 0.0, 0.0, -P, 0.0, 0.0, 0.0)


IP属地:上海1楼2024-11-14 15:09回复