截面纤维划分python脚本:SecMesh
该文件已托管至我的 gitee 账号,欢迎下载使用,https://gitee.com/yexiang-yan/ops_utility
使用过OpenSees(注意大小写拼写,其余的都是错的哦,Scott教授在一个视频里强调了好几遍😃)的朋友们肯定都体验过创建纤维截面时的痛苦,OS只提供了几种规则形状的命令(点,线,矩形,四边形,圆(扇)形),当你的截面不规则乃至构件是变截面时,写纤维命令简直太痛苦。而且,如果你想做参数化分析,生成纤维截面更是难上加难。鉴于此,我最近重构了以前写的截面纤维划分的python脚本接口,使用几行命令即可创造出任意截面的纤维划分。
话不多说,先看看效果图
-
环境依赖
环境依赖有 shapely, pygmsh, OpenSeesPy, numpy, matpltlib 这几个包
-
其中 numpy 和 matplotlib 如果你安装了 Anaconda 软件会自带; -
OpenSeesPy 是 OpenSees 的 python 封装包,强烈建议在 python 中使用 OpenSees; -
shapely 是一个几何库,它是对 GEOS 这个库的 python 封装; -
pygmsh 是对网格划分包 gmsh 的一个python封装。其实gmsh本身也提供了python接口,但是非常繁琐。这个pygmsh设计的挺不错的,于是就选择用它,毕竟少写代码,少掉头发。
安装这些包非常简单,使用 pip install shapely pygmsh 即可自动安装包及其依赖
我已经给模块命名为 SecMesh, 所以下面使用时皆从该模块中导入。
该文件已托管至我的 gitee 账号,欢迎下载使用,https://gitee.com/yexiang-yan/ops_utility
目前仅支持如下几个截面类型:
-
单一材料截面:SecPlain 类 -
钢筋混凝土截面:SecRC 类 -
钢骨混凝土截面:SecSRC 类 -
钢管混凝土截面:SecSTC 类
-
创建单一材料截面
所谓单一材料截面是指截面只有一种材料组成,可以是钢截面,也可以是混凝土,可以包含钢筋。单一材料截面被定义在类 SecPlain 中,该类提供了如下几个类方法接口:
-
add_outline(points):为截面增加外轮廓线,参数points是一个列表形状的数据,其中包含了外轮廓线的坐标点, 这些点之间为直线段,因此你只要输入截面轮廓的控制点就可以了,点闭合不闭合都可以, 因为内部设置了自动闭合 -
add_hole(points): 运行一次就会为截面增加一个开口,参数 points 和 add_outline 相同。 -
add_rebar_outline(dia, gap, d): 增加一个与外轮廓线保持相同形状的钢筋环,dia是直径,gap是钢筋之间中心间距,d是钢筋外边缘与截面轮廓线间的距离,也就是说这些钢筋是根据外轮廓线偏移得到的 -
add_rebar_hole(dia, gap, d):根据内部开口的形状偏移得到钢筋 -
add_rebar_line(self, points, dia, gap):增加任意的钢筋形状,points如上描述,可以输入控制点(只要点之间是直线即可),程序会根据gap自动剖分 -
mesh(mesh_size): 根据尺寸 mesh_size 划分截面三角形纤维 -
centring(): 中心化截面,也就是说上述坐标点可以是任意的,程序会自动把中心转移到原点。 -
rotate(theta=0):顺时针转动截面,譬如 theta=90 ,会顺时针转动截面90度 -
sec_props(): 返回截面特性 -
mesh_view(fill=True):可视化截面,fill为True的化会用颜色填充多边形;否则绘制网格
以下是案例,使用非常简单
-
多边形截面
from SecMesh import SecPlain # 导入类
points = [[0, 0], [1, 0], [2, 1], [1, 2], [0, 1]] # 因为是五边形,所以只需要输入五个顶点坐标即可
hole = [(0.5, 0.5), (1.0, 0.5), (1, 1), (0.5, 1)]
sec = SecPlain() # 实例化类
sec.add_outline(points) # 增加外轮廓线
sec.add_hole(hole) # 增加开口,注意开口只能一个个增加,可以试下注释掉该行
sec.add_rebar_outline(dia=0.032, gap=0.15, d=0.1) # 根据外轮廓线偏移钢筋
sec.add_rebar_hole(dia=0.032, gap=0.15, d=0.1) # 根据开口偏移钢筋
sec.mesh(mesh_size=0.1) # 划分纤维
sec.rotate(theta=0) # 转动截面,可以试下其他角度
sec.centring() # 中心化
sec.mesh_view() # 可视化
-
圆形截面
我们可以使用以直代曲来绘制圆形截面,注意曲线的分割点数量不可过多,过密不然会覆盖掉你的纤维尺寸,建议分成20-30段即可;
其他包括直线曲线混合的截面也类似
import numpy as np
sec = SecPlain()
# 如下创建圆形轮廓
points = []
n = 30
r = 1.2 # 圆半径
for angle in np.linspace(0, 2 * np.pi, n, endpoint=False): # endpoint 需要设置为 False
points.append((r * np.cos(angle), r * np.sin(angle)))
# 创建开口点
hole = []
r = 0.6
for angle in np.linspace(0, 2 * np.pi, n, endpoint=False):
hole.append((r * np.cos(angle), r * np.sin(angle)))
sec.add_outline(points)
sec.add_rebar_outline(dia=0.032, gap=0.15, d=0.1)
sec.add_hole(hole)
sec.mesh(mesh_size=0.1) # 可以试下 如果 n 设置太多,mesh_size 不会再控制你的纤维尺寸,而是由n控制, 这个bug以后再想法解决,但是30似乎就够圆了
# h.centring()
sec.mesh_view()
-
钢筋混凝土截面
与单一材料截面不同的是,钢筋混凝土截面多了保护层。对应的类是 SecRC。
因此该类除了具有 SecPlain 所有接口外,多了一个 add_cover(d) 命令,增加一个厚度为d的保护层
from SecMesh import SecRC
points = [(-4, 4), (-2, 4), (-2, 2), (2, 2), (2, 4), (4, 4),
(4, -4), (2, -4), (2, -2), (-2, -2), (-2, -4), (-4, -4)]
holes = [[(-3, 1), (-2, 1), (-2, -1), (-3, -1)],
[(3, 1), (2, 1), (2, -1), (3, -1)]]
sec = SecRC() # 实例化
sec.add_outline(points)
sec.add_cover(d=0.1) # 增加一个厚度为0.1 的保护层
for hole in holes:
sec.add_hole(hole) # 一次增加一个开口
sec.add_rebar_outline(dia=0.06, gap=0.15, d=0.1)
sec.add_rebar_hole(dia=0.06, gap=0.15, d=0.1)
sec.mesh(mesh_size=0.4)
sec.rotate(theta=0)
sec.centring()
sec.mesh_view(fill=True)
import numpy as np
x = np.linspace(-2, 2, 30)
y1 = np.sqrt(1 - (np.abs(x) - 1) * (np.abs(x) - 1))
y2 = np.arccos(1 - np.abs(x)) - np.pi
lines = [(xi, yi) for xi, yi in zip(x, y1)]
lines_ = [(xi, yi) for xi, yi in zip(x[-2:0:-1], y2[-2:0:-1])]
lines.extend(lines_)
sec = SecRC() # 实例化
sec.add_outline(lines)
sec.add_cover(d=0.1) # 增加一个厚度为0.1 的保护层
sec.add_rebar_outline(dia=0.032, gap=0.15, d=0.1)
sec.add_rebar_hole(dia=0.02, gap=0.15, d=0.1)
sec.mesh(mesh_size=0.3)
sec.rotate(theta=0)
sec.centring()
sec.mesh_view(fill=True)
-
钢骨混凝土截面
对应的类为 SecSRC
比 RC 截面类又多了一个 增加钢骨的命令
from SecMesh import SecSRC
points = []
r = 4
for angle in np.linspace(0, 2 * np.pi, 6, endpoint=False):
points.append((r * np.cos(angle), r * np.sin(angle)))
bone = [(-2, 1), (-1, 1), (-1,0.5), (1, 0.5), (1, 1), (2, 1),
(2, -1), (1, -1), (1,-0.5), (-1, -0.5), (-1, -1), (-2, -1)]
h = SecSRC()
h.add_outline(points) # 外轮廓线
h.add_cover(d=0.1) # 保护层
h.add_rebar_outline(dia=0.06, gap=0.15, d=0.1) # 偏移得到钢筋线
h.add_steel_bone(bone) # 增加钢骨,类似开口,一个个增加
h.mesh(mesh_size=0.3)
h.rotate(theta=0)
h.centring()
h.mesh_view()
lines = [(-2, 2), (2, 2), (2, -2), (-2, -2)]
bone = []
r = 1
n = 30
for angle in np.linspace(0, 2 * np.pi, n, endpoint=False):
bone.append((r * np.cos(angle), r * np.sin(angle)))
h = SecSRC()
h.add_outline(lines)
h.add_cover(d=0.1)
h.add_rebar_outline(dia=0.06, gap=0.15, d=0.1)
h.add_steel_bone(bone)
h.mesh(mesh_size=0.25)
h.rotate(theta=0)
h.centring()
h.mesh_view()
-
钢管混凝土截面
from SecMesh import SecSTC
lines = []
r = 1
for angle in np.linspace(0, 2 * np.pi, 24, endpoint=False):
lines.append((r * np.cos(angle), r * np.sin(angle)))
h = SecSTC()
h.add_outline(lines) # 外轮廓线
h.add_steel_tube(d=0.05) # 增加钢管层,厚度为d
h.add_rebar_tube(dia=0.06, gap=0.15, d=0.05) # 根据外轮廓线偏移钢筋
h.mesh(mesh_size=0.2)
# h.rotate(theta=0)
h.centring()
h.mesh_view()
-
生成 OpenSees 命令
以上只是划分了纤维,另外一个重要的东西是如何生成 OpenSees 命令呢?
这里提供了两个接口,都是直接生成点纤维,坐标就是每个三角形的中心,面积就是三角形的面积
-
一个是 直接生成 OpenSeesPy 命令,直接在python中调用:ops_cmds_py(self, secTag: int, matTags: ArrayLike, GJ: float) -
一个是 生成 tcl 命令的 文件,供 tcl版本读取调用:ops_cmds_tcl(self, output_dir: str, secTag: int, matTags: ArrayLike, GJ: float)
secTag: 是 要生成的截面号
matTags: 是 已经定义好的 材料号列表,注意材料号顺序是按照截面从外到内的顺序,譬如钢筋混凝土截面,就是 保护层-核心,钢骨混凝土就是 保护层-核心-钢骨,不管是任何截面,都应该把钢筋的材料号放在最后
GJ: 截面扭转刚度
-
ops_cmds_py
进行命令生成之前,假设创建了一个如下模型,定义材料
import openseespy.opensees as ops
# 定义opensees模型中的材料,依次为保护层,核心区混凝土,以及钢筋
ops.wipe()
ops.model('basic', '-ndm', 3, '-ndf', 6)
# 混凝土材料,采用Cncrete01模型,如不确定下面不用改
vc = 0.2 # 混凝土泊松比
fc = -20.1E3 # 混凝土峰值强度
ec = -2.E-3 # 混凝土峰值应变
fcu = -16.5E3 # 混凝土极限强度
ecu = -4.E-3 # 混凝土极限应变
Ec = 2 * fc / ec # 弹性模量
Gc = Ec / (2 * (1 + vc)) # 剪切模量
# ft40 = 2.39E+3 # 不考虑受拉强度
# et40 = 100.E-6
fccore = -26.8E3 # 约束混凝土峰值强度
eccore = -5.3E-3 # 约束混凝土峰值应变
fcucore = -23. # 约束混凝土极限强度
ecucore = -0.0157 # 约束混凝土极限应变
# 钢筋材料 (Steel01 Material. )
Fys = 300.E3 # 屈服应力
Es = 2.0E8 # 钢筋弹性模量
bs = 0.01 # 刚度硬化比
matTagC = 1
matTagCCore = 2
matTagSteel = 3
ops.uniaxialMaterial('Concrete01', matTagC, fc, ec, fcu, ecu)
ops.uniaxialMaterial('Concrete01', matTagCCore,
fccore, eccore, fcucore, ecucore)
ops.uniaxialMaterial('Steel01', matTagSteel, Fys, Es, bs)
# 建立节点
ops.node(1, 0.0, 0.0, 0.0)
ops.node(2, 0.0, 0.0, 0.0)
ops.fix(1, *([1]*6))
创建一个 RC 截面
from SecMesh import SecRC
points = [(-4, 4), (-2, 4), (-2, 2), (2, 2), (2, 4), (4, 4),
(4, -4), (2, -4), (2, -2), (-2, -2), (-2, -4), (-4, -4)]
holes = [[(-3, 1), (-2, 1), (-2, -1), (-3, -1)],
[(3, 1), (2, 1), (2, -1), (3, -1)]]
sec = SecRC() # 实例化
sec.add_outline(points)
sec.add_cover(d=0.1) # 增加一个厚度为0.1 的保护层
for hole in holes:
sec.add_hole(hole) # 一次增加一个开口
sec.add_rebar_outline(dia=0.06, gap=0.15, d=0.1)
sec.add_rebar_hole(dia=0.06, gap=0.15, d=0.1)
sec.mesh(mesh_size=0.4)
sec.rotate(theta=0)
sec.centring()
sec.mesh_view(fill=True)
# 生成 命令
matTags = [matTagC, matTagCCore, matTagSteel]
sec_Tag = 1
sec.ops_cmds_py(sec_Tag, matTags, GJ=10000) # 自动生成命令
为了验证 生成的命令对不对,我们可以使用以前介绍的函数 ops_sec_vis 来从 OpenSees 内部提取截面纤维可视化
from SecMesh import ops_sec_vis
# ops_sec_vis 必须在单元定义后才可使用
# 定义单元,单元号1
ops.element('zeroLengthSection', 1, 1, 2, sec_Tag)
ops_sec_vis(ele_tag=1, rebar_d=0.06) # 可视化
可以看到,根据OpenSees内部提取的纤维和我们创造的是一致的,说明 代码ok
-
ops_cmds_tcl
from SecMesh import SecRC
points = [(-4, 4), (-2, 4), (-2, 2), (2, 2), (2, 4), (4, 4),
(4, -4), (2, -4), (2, -2), (-2, -2), (-2, -4), (-4, -4)]
holes = [[(-3, 1), (-2, 1), (-2, -1), (-3, -1)],
[(3, 1), (2, 1), (2, -1), (3, -1)]]
sec = SecRC() # 实例化
sec.add_outline(points)
sec.add_cover(d=0.1) # 增加一个厚度为0.1 的保护层
for hole in holes:
sec.add_hole(hole) # 一次增加一个开口
sec.add_rebar_outline(dia=0.06, gap=0.15, d=0.1)
sec.add_rebar_hole(dia=0.06, gap=0.15, d=0.1)
sec.mesh(mesh_size=0.4)
sec.rotate(theta=0)
sec.centring()
sec.mesh_view(fill=True)
output_dir = "sec_mesh.tcl" # 多了这么个参数,输出文件
matTags = [1, 2, 3]
sec_Tag = 1
sec.ops_cmds_tcl(output_dir, sec_Tag, matTags, GJ=10000) # 自动生成命令
看一下生成的 tcl 文件截图,已生成完毕,ok。其中temp_matTag 为临时生成的材料号,方便后续修改
-
截面特性
props = sec.sec_props()
print(props)
{'area': 44.00000000000003, 'center': (5.129735020597595e-15, -9.235036977563796e-16, 0), 'Iy': 190.41472565048613, 'Iz': 294.42365487178506, 'ry': 2.0802902204133464, 'rz': 2.5867830797404476, 'Ip': 484.8383805222712, 'Iyz': 0.00100144291708415}
-
结语
暂时支持这么些功能,日常的使用应该足够了,后续再慢慢增加。
利用上述类代码,你还可以定制任意类型截面的参数化函数或类。譬如你可以定制一个开口矩形截面,利用长、宽、壁厚等参数来生成截面,这里就不再演示了。
如有代码问题,可在gitee中提交,或与本人联系,
创作不易,欢迎关注、点赞、转发、支持。谢谢!