Go back

OpenSeesPy 快速施加重力荷载, 获取质量, 刚度矩阵

在最新的opstool版本中，实现了几个函数用来给 OpenSeesPy 模型快速施加重力荷载，获取质量、刚度以及阻尼系统矩阵。

在 OpenSees 中，我们有许多方式可以给模型施加质量，譬如：

创建节点时我们可以顺手指定：

node(nodeTag, *crds, '-mass', *mass)

可以专门指定节点质量（会覆盖掉创建节点时指定的质量）：

mass(nodeTag, *massValues)

可以在创建梁或桁架单元时顺手指定质量密度（单位长度，<'-mass', mass>, <'-cMass'>）：

element('elasticBeamColumn', eleTag, *eleNodes, Area, E_mod, G_mod, Jxx, Iy, Iz, transfTag, <'-mass', mass>, <'-cMass'>)

对于壳单元、平面以及实体单元等，我们还可以在创建nD材料时指定质量密度rho：

nDMaterial('ElasticIsotropic', matTag, E, nu, <rho>)

这些方式可谓五花八门，但是它们都可以被OpenSees处理，最终走进质量矩阵里，然后进行特征值分析、瞬态分析等。

朋友，如果您只是用来做模态分析和时程分析，那么您可以任意使用上述方式添加质量。 但是，如果您想把上述质量转化成重力荷载，恰巧模型又复杂，那么您的问题就来了，OpenSees 没有命令转化重力荷载，也没有命令可以提取出上述所有质量， 只有一个命令 nodeMass 只可以提取出 mass 创建的质量。

想必用 OpenSees 创建过复杂模型的朋友都有过这种痛苦，质量转化重力荷载非常繁琐。

但是，作为天选打工人，在下当然不会放过任何一个可以偷懒的机会，当然，独懒懒不如众懒懒，要快乐大家一起快乐。 因此，我实现了一个 create_gravity_load函数，能够自动将OpenSees模型中存在的所有质量转化为重力荷载。

以下是代码，别忘了更新opstool的版本，在新的版本里，还有线性屈曲分析、高斯点结果向节点插值等功能。

import openseespy.opensees as ops

import opstool as opst
import opstool.vis.pyvista as opsvis

一个简单的2D框架模型

模型定义

我们给节点1—6都显式施加了质量，也在创建梁单元时指定了质量线密度。 为了说明提取质量的鲁棒性，我们还故意用了两个多点约束，约束3和4号节点，5和6号节点。

ops.wipe()
ops.model("basic", "-ndm", 2, "-ndf", 3)

# Define the model nodes
ops.node(1, 0.0, 0.0)
ops.node(2, 2.0, 0.0)
ops.node(3, 0.0, 2.0)
ops.node(4, 0.0, 2.0)
ops.node(5, 2.0, 2.0)
ops.node(6, 2.0, 2.0)

# Constrain the nodes
ops.rigidLink("beam", 3, 4)  # Rigid link between nodes 3 and 4
ops.equalDOF(5, 6, 1, 2, 3)  # Equal DOF between nodes 5 and 6

ops.fix(1, 1, 1, 1)
ops.fix(2, 1, 1, 1)

# Assign the node masses
ops.mass(1, 1.0, 1.0, 0.0)
ops.mass(2, 1.0, 1.0, 0.0)
ops.mass(3, 1.0, 1.0, 0.0)
ops.mass(4, 1.0, 1.0, 0.0)
ops.mass(5, 10.0, 10.0, 1e-3)
ops.mass(6, 10.0, 10.0, 1e-3)

# Define the elements
ops.geomTransf("Linear", 1)

# element('elasticBeamColumn', eleTag, *eleNodes, Area, E_mod, Iz, transfTag, <'-mass', mass>)
Area, E_mod, Iz = 2.0, 2.0, 2.0
ops.element("elasticBeamColumn", 1, 1, 3, Area, E_mod, Iz, 1, "-mass", 2.0)
ops.element("elasticBeamColumn", 2, 2, 5, Area, E_mod, Iz, 1, "-mass", 2.0)
ops.element("elasticBeamColumn", 3, 4, 6, Area, E_mod, Iz, 1, "-mass", 2.0)

可视化一下，jupyter_backend="static"只是为了生成静态图形（懒得截图），实际中不要使用。

fig = opsvis.plot_model(show_node_numbering=True, show_ele_numbering=True)
fig.show()

直接转化重力荷载

timeSeries和pattern需要事先定义，然后使用create_gravity_load来一键创建重力荷载。

内部会使用 load命令来创建节点荷载。

ops.timeSeries("Constant", 1)  # Define a constant time series
ops.pattern("Plain", 1, 1)  # Define a load pattern
node_loads = opst.pre.create_gravity_load(direction="Y", factor=-9.81, exclude_nodes=None)  # Create gravity loads in the pattern 1

检查一下返回的荷载信息：

for node, load in node_loads.items():
    print(f"Applying load at node {node}: {load}")

Applying load at node 1: [0.0, -29.43, 0.0]
Applying load at node 2: [0.0, -29.43, 0.0]
Applying load at node 3: [0.0, -29.43, 0.0]
Applying load at node 4: [0.0, -29.43, 0.0]
Applying load at node 5: [0.0, -117.72, 0.0]
Applying load at node 6: [0.0, -117.72, 0.0]

获取质量

实际上，create_gravity_load内部先使用get_node_mass来提取节点质量，然后再转化荷载。

get_node_mass 帮助我们提取当前模型中存在的所有质量（会等效为节点质量）：

node_mass = opst.pre.get_node_mass()  # get the node mass dict

打印出来返回的节点质量：

print("Node Masses:")
for node, mass in node_mass.items():
    print(f"Node {node}: {mass}")

Node Masses:
Node 1: [3.0, 3.0, 0.0]
Node 2: [3.0, 3.0, 0.0]
Node 3: [3.0, 3.0, 0.0]
Node 4: [3.0, 3.0, 0.0]
Node 5: [12.0, 12.0, 0.001]
Node 6: [12.0, 12.0, 0.001]

让我们检查下：

─── ⋆⋅☆⋅⋆ ── ✍️

▶️ 首先，一部分质量来源于直接指定的节点质量

ops.mass(1, 1.0, 1.0, 0.0)
ops.mass(2, 1.0, 1.0, 0.0)
ops.mass(3, 1.0, 1.0, 0.0)
ops.mass(4, 1.0, 1.0, 0.0)
ops.mass(5, 10.0, 10.0, 1e-3)
ops.mass(6, 10.0, 10.0, 1e-3)

▶️ 其次来源于梁单元的质量线密度，每个单元的：

长度 = 2, 线密度 = 2.0

每个端节点获得的质量:

2 * 2.0 / 2 = 2.0

加起来，和获取的质量是一样的，完美！

🤖 因此，您可以在创建模型时随意施加质量，怎么方便怎么来，opstool 将会为那扫清障碍！

获取系统矩阵

实际上，节点质量是从模型的质量矩阵推算出来的，至于怎么推算的，可以去看一下 opstool 源代码。

新版，我也实现了获取系统矩阵的函数。

在获取系统矩阵之前，我需要强调，不同的constraints命令会改变系统矩阵的维度，system命令必须用"FullGeneral"才可以完整的矩阵，numberer命令会改变矩阵的拓扑（节点出现的行列位置），因此获取的系统矩阵自然是不同的。 至于怎么样的影响，可以翻翻以前的公众号推文。

选择如下的参数：

constraints_args = ["Penalty", 1e3, 1e3]  # This affects the dimensions of the matrix
system_args = ["FullGeneral"]  # Full matrix system must be used
numberer_args = ["RCM"]  # It affects the topology of the matrix and the order of dof.

打印质量矩阵

M = opst.pre.get_mck("m", constraints_args=constraints_args, system_args=system_args, numberer_args=numberer_args)
print("Mass Matrix (M):")
print(M.to_numpy().shape)

Mass Matrix (M):
(18, 18)

看一下 节点号-自由度号 的排列，可以看到在 Penalty 约束下，所有节点的所有自由的都被保留下来了。

“RCM” 的编号处理器为了减少带宽，并没有按照节点号建立的顺序来排列自由度。

print(M.coords["nodeTagsDofs-1"])
<xarray.DataArray 'nodeTagsDofs-1' (nodeTagsDofs-1: 18)> Size: 288B
    array(['2-0', '2-1', '2-2', '5-3', '5-4', '5-5', '6-6', '6-7', '6-8', '4-9',
           '4-10', '4-11', '3-12', '3-13', '3-14', '1-15', '1-16', '1-17'],
          dtype='<U4')
    Coordinates:
      * nodeTagsDofs-1  (nodeTagsDofs-1) <U4 288B '2-0' '2-1' ... '1-16' '1-17'

刚度矩阵

K = opst.pre.get_mck("k", constraints_args=constraints_args, system_args=system_args, numberer_args=numberer_args)
print("Stiffness Matrix (K):")
print(K.to_numpy().shape)

Stiffness Matrix (K):
(18, 18)

K.plot.imshow()

再换一套参数：

constraints_args = ["Transformation"]  # This affects the dimensions of the matrix
system_args = ["FullGeneral"]  # Full matrix system must be used
numberer_args = ["Plain"]  # It affects the topology of the matrix and the order of dof.
K = opst.pre.get_mck("k", constraints_args=constraints_args, system_args=system_args, numberer_args=numberer_args)
print("Stiffness Matrix (K):")
print(K.to_numpy())
    Stiffness Matrix (K):
    [[ 8.  0.  6. -2.  0.  0.]
     [ 0.  8.  6.  0. -6.  6.]
     [ 6.  6. 16.  0. -6.  4.]
     [-2.  0.  0.  8.  0.  6.]
     [ 0. -6. -6.  0.  8. -6.]
     [ 0.  6.  4.  6. -6. 16.]]
M = opst.pre.get_mck("M", constraints_args=constraints_args, system_args=system_args, numberer_args=numberer_args)
print("Mass Matrix (M):")
print(M)
    Mass Matrix (M):
    <xarray.DataArray (nodeTagsDofs-1: 6, nodeTagsDofs-2: 6)> Size: 288B
    array([[6.0e+00, 0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00],
           [0.0e+00, 6.0e+00, 0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00],
           [0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00],
           [0.0e+00, 0.0e+00, 0.0e+00, 2.4e+01, 0.0e+00, 0.0e+00],
           [0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00, 2.4e+01, 0.0e+00],
           [0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00, 0.0e+00, 2.0e-03]])
    Coordinates:
      * nodeTagsDofs-1  (nodeTagsDofs-1) <U3 72B '4-0' '4-1' '4-2' '6-3' '6-4' '6-5'
      * nodeTagsDofs-2  (nodeTagsDofs-2) <U3 72B '4-0' '4-1' '4-2' '6-3' '6-4' '6-5'

Coordinates 是 “nodeTag-dof” 的列表，可以看出在 “Transformation” 约束下，只保留了 4 号 和 6号节点，这是因为：

1和2号节点被 fix 约束了，3和4号节点和5、6节点是多点约束，Transformation 约束会消去它们。

通过质量的值，我们还可以发现，消去之前，OpenSees 内部把质量也给搬运叠加了。

如果试一下其他的参数，您会发现：

  • constraints(‘Plain’) 系统矩阵维度会减少（消去约束自由度）；
  • constraints(‘Lagrange’, alphaS, alphaM) 会增加矩阵维度，增加的个数为自由度约束数量；
  • constraints(‘Penalty’, alphaS, alphaM) 矩阵维度不变。
  • constraints(‘Transformation’) 系统矩阵维度会减少（消去约束自由度）；
  • numberer(‘Plain’) 按节点号顺序排列自由度
  • numberer(‘RCM’) 节点号在矩阵中出现的位置不确定
  • numberer(‘AMD’) 节点号在矩阵中出现的位置不确定
Share this post on X Share this post via Telegram Share this post on Mastodon
Tags:
Next Post
我的诗词小集