openseespy建模

Question

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和短梁单元的刚度特性）
kx = 1e6  # 横向连接刚度值，单位为kN/m
A = 2.5  # 截面积，需根据实际情况设定
E = kx  # 弹性模量设为横向刚度值
I = 1e-6  # 惯性矩，根据需要设置小值
ops.geomTransf('Linear', 3, 1, 0, 0)
# 定义五片梁的横向连接
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

#定义支座
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)

请问有人知道建模问题出在哪吗？

Answer 1

本答案参考 多次询问 ChatGPT-3.5 后选择的最佳回答，若是满意，还请采纳，谢谢

以下是对这段代码可能存在问题的分析：

一、材料定义部分

1. 单位一致性
   • 在定义材料属性时，虽然在代码开头指定了单位（kN, m, s），但是在定义材料属性时，如混凝土抗压强度fc和钢筋屈服强度fy是以MPa为单位给出的，而在后续使用uniaxialMaterial定义材料时，将fy乘以1e3转换为kN/m²，这种单位转换容易出错。如果在其他地方不小心使用了错误的单位量级，可能会导致计算结果错误。
   • 例如，如果在后续计算中，有部分公式假设输入的材料强度是MPa而没有进行正确转换，就会出现问题。
2. 材料模型参数合理性
   • 对于Concrete01材料模型，所给定的应力 - 应变参数（如ops.uniaxialMaterial("Concrete01", 1, -41000, -0.004, -34470, -0.014)）需要根据实际的混凝土材料特性进行仔细校准。这些参数可能因混凝土的等级、配合比等因素而有很大差异，如果使用不合理的参数，可能导致模拟结果与实际结构行为不符。

二、节点和单元定义部分

1. 节点编号和索引管理
   • 在创建节点和单元时，节点编号的分配方式较为复杂，容易出现编号冲突或逻辑错误。例如，在定义梁单元时，node_tag_start = 1000 + n * (nL + 1)+ i和node_tag_end = 1000 + n * (nL + 1)+ i + 1这种编号方式需要确保在整个模型中节点编号是唯一的，并且在后续操作（如施加荷载、约束等）中正确引用。
   • 如果在模型构建过程中，由于节点编号错误导致单元连接到错误的节点上，会影响结构的拓扑结构，进而导致分析结果错误。
2. 单元类型和积分设置
   • 使用dispBeamColumn单元时，积分设置（如Lobatto积分）需要根据具体的结构行为和精度要求进行调整。NP = 5这个积分点数量可能对于某些复杂的结构行为（如大变形、强非线性等）不够精确。不同的积分方案和积分点数量会影响单元的刚度矩阵计算，从而影响整个结构的分析结果。

三、荷载和约束部分

1. 支座定义函数define_bearings
   • 在define_bearings函数中，对于两种支座类型（Rubber和Frame），内部的操作完全相同，都是使用equalDOF进行节点自由度的约束。这可能是代码编写过程中的失误，正常情况下两种不同类型的支座应该有不同的约束方式或者参数设置。
2. 重力荷载施加
   • 在施加重力荷载时，通过P = ops.nodeMass(node,3)*g计算节点荷载，这里假设节点质量在z方向的分量乘以重力加速度就是重力荷载。然而，如果节点质量的计算或者存储方式存在问题（例如，质量矩阵的定义与预期不一致），可能会导致施加的重力荷载不正确。

四、截面定义部分

1. 纤维截面定义
   • 在定义墩柱和梁的纤维截面时，对于纤维的布置（如ops.patch和ops.fiber的使用），需要确保纤维的数量、位置和面积等参数准确反映实际结构的截面特性。例如，在梁的钢筋布置中，如果钢筋数量、间距等设置不符合实际设计要求，会影响梁的受力性能模拟。
   • 特别是在定义墩柱截面时，计算核心混凝土和保护层混凝土的面积以及钢筋位置等操作，如果存在计算错误，会导致截面的力学性能（如惯性矩、抵抗矩等）计算错误，进而影响整个结构的分析结果。