2节点梁/杆单元的弹性刚接框架(二维、三维)（python,有限元）

第四篇 2节点梁/杆单元的弹性刚接框架(二维、三维)

介绍

本程序使用的梁-杆单元刚度矩阵是通过叠加之前描述的梁和杆刚度矩阵形成的，可以分析二维或三维的框架结构，其单元刚度矩阵由子程序rigid_jointed生成。

计算实例1

第一个算例如下图所示，为一个受均布荷载和点荷载作用的二维刚节点框架，在2D中，单元有6个自由度，如下图2所示，每个节点有x和y上两个自由度和一个旋转自由度(按此顺序)，自由节点编号包含在向量g，因此,单元1的编号自由度为[0 0 1 2 3 4)，单元2的编号自由度为(2 3 4 5 6 7)。输入数据喝之前类似，第一行数据是单元的的数量(nel)、节点的数量(nn)、维数(ndim)、性质nprops的数量和性质类型个数np_types的数量。
所需的材料属性的个数和维度有关，因此nprops表示每种属性类型所需的材料属性的个数，在一个2D框架问题中，需要两个材质属性(EA和EI)，因此nprops = 2，然后将每种类型的材料性质值读入二维数组prop，材料属性数据后面跟着etype向量(如果需要的话)，为每个单元赋予对应的材料性质，之后有单元顺序对应的坐标值g_coord和每个单元对应的节点序号g_num，约束自由度节点号和等效节点荷载。

单元个数=6，节点个数=6；维度=2；单元性质类型=2，性质个数=2，为4.0e9的ea和6.0e4的ei;1.0e9的ea和2.0e4的EI，前者为第一类，后者为第二类；
每个单元的性质种类为1，1，1，2，2，2；
每个节点坐标值为(0.0,0.0),(6.0,0.0),(6.0,-4.0),(12.0,0.0),(12.0,-5.0),(14.0,0.0),
按顺序每个单元对应节点号为(1,2),(2,4),(4,6),(3,2),(3,4),(5,4);
约束点3个，节点1x,y方向的约束,节点3和5位移和旋转完全约束；
节点1具有y方向荷载-60和力矩-60，节点2具有y方向荷载-180和力矩-80，节点4具有y方向荷载-140和力矩133.33，节点1具有y方向荷载-20和力矩-6.67；
没有位移和转角。

代码1

import numpy as np
import A
import math
nels=6
np_types=2
loaded_nodes=4
fixed_freedoms=0
nod=2
nprops=2
nn=6
ndim=2
nr=3
if ndim==2:
    nodof=3
elif ndim==3:
    nodof=6
ndof=nod*nodof
g_coord=np.array([[0,6.0,6.0,12.0,12.0,14.0],[0.0,0.0,-4.0,0.0,-5.0,0.0]])
g_num=np.array([[1,2,4,3,3,5],[2,4,6,2,4,4]])
#初始化定义数组
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))
if np_types==1:
  etype[:,0]=1
else:
  etype[:,0]=(1,1,1,2,2,2)
if np_types==1:
  prop[:,0]=(5.0e9,6.0e4)
  prop[:,1]=(1.0e9,2.0e4)
else:
  prop[:,0]=(5.0e9,6.0e4)
  prop[:,1]=(1.0e9,2.0e4)
if ndim==3:
    gamma=np.array([[0],[0],[90]])
#读取nr
if nr!=0:
  Dim_1=[1,3,5]
  nf_value=np.array([[0,0,0],[0,0,0],[1,0,0]])
  m=0
  for i in Dim_1:
      for j in range(1,nodof+1):
        nf[j-1,i-1]=nf_value[j-1,m]
      m=m+1  
#form nf
A.formnf(nf)
neq=int(max(nf.reshape(nf.shape[0]