使用2节点梁或梁/杆单元分析弹塑性梁或框架(python,有限元)

第五篇 使用2节点梁或梁/杆单元分析弹塑性梁或框架

介绍

本程序属于第四篇的改进版，其中对构件能够承受的最大力矩进行了限制，随着结构上荷载的增加，形成塑性铰，结构会变得失稳，该程序目前仅限于荷载主导。
这是介绍的第一个非线性程序，假设材料为理想弹塑性，对于这种非线性问题，迭代的策略被使用，迭代法用于在给定的一组外加荷载下求节点位移和单元“作用”力，超过塑性极限的力矩将重新分配给其他仍保留力矩承载能力的铰头，当在公差范围内，构件节点处的力矩不超过其极限塑性值，且内部“作用”力与施加的外部荷载处于平衡状态时，迭代过程收敛。

解决这类问题的传统方法是，随着接头达到塑性极限，逐步修正全局刚度矩阵。修正是必要的，因为塑料接头被具有适当塑性力矩的铰接头代替，全局刚度矩阵只形成一次，通过迭代修正结构的受力荷载，直到收敛。
当问题涉及到常量左边矩阵和右侧向量求解，一般划分为两个阶段，乔列斯基分解sparin，和从前-从后迭代过程spabac。
本程序在需要读取全部的弹性性质和塑性力矩值，对于一维梁和二维骨架，仅仅需要一个塑性力矩，因此nprops在1维中为2，二维中为3，对于三维结构，三个塑性力矩值需要获取，分别为y和z轴的极限弯矩和横轴的极限扭矩，因此nprops=7,
载荷以增量incs的形式作用于节点，并且每个增量的大小被读取到矢量数据载荷dload中，荷载维持和塑性铰分析中同样的比例，因此节点荷载的相对大小读入node和val，其中node包含节点编号，val为每个自由度上的荷载权重。
在安装全局刚度矩阵和因式分解之后，程序进入荷载增量循环，通过iters来代表迭代次数，加上外加荷载增量之后重新分配荷载向量bdylds，使用spabac子程序解节点位移增量，然后通过子程序checon与前一次位移值比较，如果改变值小于容差tol，说明结构收敛。
在每次迭代中，遍历每个单元，并用节点位移和单元刚度矩阵的乘积计算其受力，子程序hinge将力添加到先前加载（保存在holdr中）已有的值中，并检查全部节点看是否超过塑性力矩值，如果超过塑性力矩值，则得到自平衡矢量react。在下图中，显示了一个典型的二维单元，其中一个特定的荷载增量使两个节点处的力矩值过其塑性极限，每个节点施加的修正向量力矩等于塑性力矩值的过盈量，然而，为了保持平衡，需要一对力矩，入下面第二个图的的局部坐标系所示。最后，如第三个图所示，在组装体荷载bdylds向量之前，将耦合转换为全局坐标，只有力矩超过塑性极限的构件才会增加体荷载bdylds。

如果在任何荷载阶段，该算法未能在规定的迭代上限内收敛，则表明结构失稳破坏，因为该算法在不违反塑性弯矩值的情况下不能满足平衡。

算例1

下图所示的第一个例子是受比例荷载作用的双开式门式刚架，每增加一次荷载后，得到完成收敛的加载因子λ(等于dload的累加值)，加载的节点位移和迭代次数值，为减少输出值量，给出加载节点(2、3、6)的位移。在图4.33中，A点的水平位移值被绘制到λ上，这与Horne(1971)给出的理论值λf = 1.375非常吻合。

代码1

import numpy as np
import A
import math
nels=7
np_types=3
loaded_nodes=3
fixed_freedoms=0
nod=2
nprops=3
nn=8
ndim=2
nr=4
limit=200
tol=0.0001
if ndim==2:
    nodof=3
elif ndim==3:
    nodof=6
ndof=nod*nodof
g_coord=np.array([[0.0,0.0,10.0,20.0,20.0,35.0,50.0,50.0],[0.0,15.0,15.0,15.0,0.0,15.0,15.0,0.0]])
g_num=np.array([[1,2,3,5,4,6,8],[2,3,4,4,6,7,7]])
#初始化定义数组
nf=np.ones((nodof,nn),dtype=np.int64)
g=np.ones((ndof,1),dtype=np.int64)
num=np.ones((nod,1),dtype=np.int64)
etype=np.ones((nels,1),dtype=np.int64)
eld=np.ones((ndof,1))
km=np.ones((ndof,ndof))
coord=np.ones((nod,ndim))
action=np.ones((ndof,1))
g_g=np.ones((ndof,nels))
prop=np.ones((nprops,np_types))
gamma=np.ones((nels,1))
holdr=np.zeros((ndof,nels))
react=np.zeros((ndof,1))
converged=np.array([True])
if np_types==1:
  etype[:,0]=1
else:
  etype[:,0]=(1,2,2,1,3,3,1)
if np_types==1:
  prop[:,0]=(1.0,1e4,1e4,1.0,1.0,1.0,1e8)
else:
  prop[0,:]=(1.0e10,1.0e10,1.0e10)
  prop[1,:]=(1.0e6,1.0e6,1.0e6)
  prop[2,:]=(20,50,80)
if ndim==3:
    gamma=np.array([[0],[0],[0],[0],[0],[0],[0],[0],[0],[0],[0],[0]])
#读取nr
if nr!=0:
  Dim_1=[1,5,8]
  nf_value=np.array([[0,0,0],[0,0,0],[0,0,0]])
  m=0
  for i in Dim_1:
      for j in range(1,</