本文经作者同意转载
以下文章来源于潜心画栋亦雕梁 ，作者Yan
许久没有更新公众号了，这次介绍如何在 OpenSeesPy (以下简称OPS) 中快速指定梁单元的局部坐标系。
三维梁单元需要你手动指定局部坐标系，OPS 中关于梁单元的局部坐标系指定方式是这样的：
geomTransf(transfType, transfTag, *vecxz)
其中 transfType 为坐标变换类型，可选的有 ‘Linear’、’PDelta’ 和 ‘Corotational’ 三种，其中 ‘Linear’ 为线性变换，’PDelta’ 考虑了 P-Delta效应， ‘Corotational’意思是使用共旋法的变换，主要用于梁的大转动小应变几何非线性问题。transfTag是坐标变换整数身份标签，vecxz 是局部坐标系定位向量，下面会详解。前面的星号 * 意思是解包符，为什么要加这个，是因为 OpenSeesPy 的封装方式是只使用了一系列位置参数，且这些参数只支持整数、浮点数以及字符串类型，不支持列表、元组等结构，也不支持关键字传参。直观点就是 fun(*[1,2,3]) 等价于 fun(1,2,3)。
为什么需要 vecxz？在 OPS 中，梁单元的局部x轴方向由节点I端指向J端确定，剩下的就是确定局部y轴和z轴，如何确定就需要 vecxz 这个东西。vecxz 用来确定局部坐标系的 x-z 平面，用这个向量叉乘局部x轴就是局部y轴（右手法则，方向从vecxz末端指向局部x轴末端，大拇指方向就是局部y轴）， 局部y轴确定了，再用下右手法则就可以确定局部z轴了。注意，vecxz 向量不是局部z轴，它平行于梁的局部x-z平面，且与局部z轴夹角一定小于90度。
根据这个原理，我们就可以很高效率的进行坐标变换了，当然，不是一个梁一个梁的指定，而是全局指定。我们设想，是否可以让所有的梁单元局部z轴都在局部x轴和整体z轴组成的平面内呢，当然可以，我们只需要让 vecxz = [0, 0, 1] 就行了。但是一个情况是竖直单元，它的局部x轴与整体z轴是共线的，这时 vecxz = [0, 0, 1] 无法与局部x轴叉乘，这时又怎么办呢？很简单，我们换个方向，对于竖向单元，我们让 vecxz = [-1, 0, 0]，这时局部z轴与整体x轴是平行的。综上，我们把梁单元分为两类，竖向单元与非竖向单元，竖向单元 vecxz = [-1, 0, 0]， 非竖向单元 vecxz = [0, 0, 1]。
简单粗暴！下面我们就开始编程看看。
import numpy as np
import openseespy.opensees as ops
import opstool as opst
ops.wipe()
ops.model('basic', '-ndm', 3, '-ndf', 6)
下面我们提前创建好两类梁单元的坐标变换命令，并创造一个函数 beam_geom_transf 用来根据梁两端的节点坐标与变换类型transfType输出要使用的变换标签。
这个代码是通用的，可以用于你的任何模型中！
ops.geomTransf('Linear', 1, *[-1, 0, 0])
ops.geomTransf('PDelta', 2, *[-1, 0, 0])
ops.geomTransf('Corotational', 3, *[-1, 0, 0])
ops.geomTransf('Linear', 4, *[0, 0, 1])
ops.geomTransf('PDelta', 5, *[0, 0, 1])
ops.geomTransf('Corotational', 6, *[0, 0, 1])
geomTransf_ver = {'Linear':1, 'PDelta':2, 'Corotational':3} # 竖向单元标签字典
geomTransf_other = {'Linear':4, 'PDelta':5, 'Corotational':6} # 非竖向单元标签字典
def beam_geom_transf(nodeI, nodeJ, transfType): # 输出相应的变换标签
nodeI = np.array(nodeI)
nodeJ = np.array(nodeJ)
xaxis = nodeJ - nodeI
global_z = [0, 0, 1]
cos_angle = xaxis.dot(global_z)/(np.linalg.norm(xaxis) * np.linalg.norm(global_z))
if np.abs(1 - cos_angle ** 2) < 1e-10: # 判断是否竖向单元
tag = geomTransf_ver[transfType]
else:
tag = geomTransf_other[transfType]
return tag
创建节点，立方体的八个顶点
coords = [(0, 0, 0),
(1, 0, 0),
(1, 1, 0),
(0, 1, 0),
(0, 0, 1),
(1, 0, 1),
(1, 1, 1),
(0, 1, 1)]
beam_link = [(0)]
nodeTags = []
for i, coord in enumerate(coords):
ops.node(i, *coord)
nodeTags.append(i)
ops.fix(0, 1, 1, 1, 1, 1, 1)
ops.fix(1, 1, 1, 1, 1, 1, 1)
ops.fix(2, 1, 1, 1, 1, 1, 1)
ops.fix(3, 1, 1, 1, 1, 1, 1)
下面创建弹性梁单元，我使用了 combinations 这个函数用来从八个节点标签中取出任意两个顶点组成的组合，不重复！
from itertools import combinations
beam_link = list(combinations(nodeTags, 2))
print(beam_link)
for i, link in enumerate(beam_link):
eleTag = i + 1
coordI = ops.nodeCoord(link[0])
coordJ = ops.nodeCoord(link[1])
transfTag = beam_geom_transf(coordI, coordJ, transfType='Linear')
ops.element('elasticBeamColumn', eleTag, *link, 1., 1., 1., 1., 1., 1., transfTag, '-mass', 0.01)
[(0, 1),
(0, 2),
(0, 3),
(0, 4),
(0, 5),
(0, 6),
(0, 7),
(1, 2),
(1, 3),
(1, 4),
(1, 5),
(1, 6),
(1, 7),
(2, 3),
(2, 4),
(2, 5),
(2, 6),
(2, 7),
(3, 4),
(3, 5),
(3, 6),
(3, 7),
(4, 5),
(4, 6),
(4, 7),
(5, 6),
(5, 7),
(6, 7)]
使用我开发的 opstool 包来可视化看看，打开显示局部坐标系开关 show_local_crd=True， 关于 opstool 的包文档参见 https://opstool.readthedocs.io/en/latest/。
ModelData = opst.GetFEMdata(results_dir="opstool_output")
ModelData.get_model_data(save_file="ModelData.hdf5")
opsvis = opst.OpsVisPyvista(point_size=10, line_width=2,
colors_dict=None, theme="document",
color_map="Spectral", on_notebook=False,
results_dir="opstool_output")
opsvis.model_vis(input_file="ModelData.hdf5",
show_node_label=False, show_ele_label=False,
show_local_crd=True, label_size=8,
show_outline=False,
opacity=1.0,
save_fig=None)
Model data saved in opstool_output/ModelData.hdf5!
可以看出各梁单元的局部轴和我们预期的是一样的
竖向单元的局部z轴和整体x轴平行，局部y轴与整体y轴平行，
局部x轴与整体z轴平行！
非竖向单元，局部z轴位于局部x轴与整体z轴组成的平面内，且指向
位于同一半区，且水平单元局部z轴与整体z轴平行！
也可以可视化模态！
ModelData.get_eigen_data(mode_tag=5, solver="-genBandArpack",
save_file='EigenData.hdf5')
opsvis.eigen_vis(input_file='EigenData.hdf5',
mode_tags=[1, 5], subplots=True,
alpha=None, show_outline=False,
show_origin=False, opacity=1.0,
show_face_line=False, save_fig=None)
Eigen data saved in opstool_output/EigenData.hdf5 !
今天暂时到这里，再会！