2节点杆-弹性铰接框架的分析(2维,3维)(python,有限元)

第二篇 2节点杆-弹性铰接框架的分析(2维,3维)

介绍

本文程序是在上文程序的基础上做的修改，把一维问题扩展到二维或者三维；一维问题中，节点数始终为比单元数大一，而这种情况在本文中不再适用；单元长度，即杆的长度不再直接读取，将根据坐标值计算；由于问题扩展到二维和三维，所以每个节点的自由度数量也变多。

计算实例1(2维)

单元个数=10，节点个数=6，2维，单元性质类型=1，为2.0e5；
约束点2个，节点1两个方向，节点2y方向；
节点6具有y方向荷载-10KN，无位移荷载。

代码1(2维)

import numpy as np
import A
nels=10
np_types=1
loaded_nodes=1
fixed_freedoms=0
nod=2
nprops=1
nn=6
ndim=2
nr=2
nodof=ndim
ndof=nod*nodof
g_coord=np.array([[0,4.0,4.0,8.0,8.0,12.0],[3.0,0.0,3.0,3.0,0.0,0.0]])
g_num=np.array([[1,1,3,3,3,2,2,5,4,5],[2,3,4,5,2,4,5,4,6,6]])
#初始化定义数组
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))
eld=np.ones((ndof,1))
km=np.ones((ndof,ndof))
coord=np.ones((nod,ndim),dtype=np.int64)
action=np.ones((ndof,1))
g_g=np.ones((ndof,nels))
prop=np.ones((nprops,np_types))
if np_types==1:
  etype[:,0]=1
else:
  etype[:,0]=(2,2,1,1)
if np_types==1:
  prop[:nprops,np_types-1]=2.0e5
else:
  prop[:nprops,np_types-2]=2000
  prop[:nprops,np_types-1]=1000
#读取nr
if nr!=0:
  Dim_1=[1,2]
  nf_value=np.array([[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]*nf.shape[1],1)))
kdiag=np.zeros((neq,1),dtype=np.int64)
loads=np.zeros((neq+1,1))
##累加单元去发现全局尺寸
for iel in range(1,nels+1):
  num[:,0]=g_num[:,iel-1]
  A.num_to_g(num,nf,g)
  g_g[:,iel-1]=g[:,0]
#call fkidiag
  A.fkdiag(kdiag,g)
for i in range(1,neq):
  kdiag[i]=kdiag[i]+kdiag[i-1]
kv=np.zeros((kdiag[neq-1,0],1),dtype=float)
print('一共有'+str(neq)+'等式和天际线存储个数为'+str(kdiag[neq-1]))
#全局刚度矩阵安装
for iel in range(1,nels+1):
    num=g_num[:,iel-1]
    coord=np.transpose(g_coord[:,num-1])
#call pin_jointed
    A.pin_jointed(km,prop[0,int(etype[iel-1])-1],coord)
##########
    g[:,0]=g_g[:,iel-1]
#call fsparv
    A.fsparv(kv,km,g,kdiag)
###读荷载和位移
if loaded_nodes!=0:
  Dim_2 = [6]
  load_value=np.array([[0.0],[-10.0]])
  m=0
  for i in Dim_2:
      for j in range(1,nodof+1):
        loads[nf[j-1,i-1]-1]=load_value[j-1,m]
      m=m+1
if fixed_freedoms!=0:
  Dim_3=[5]
  for i in Dim_3:
    node=np.ones((fixed_freedoms,1),dtype=np.int64)
    sense=np.ones((fixed_freedoms,1),dtype=np.int64)
    value=np.ones((fixed_freedoms,1))
    no=np.ones((fixed_freedoms,1),dtype=np.int64)
####位移赋值
    node[fixed_freedoms-len(Dim_3)]=i
    sense[fixed_freedoms-len(Dim_3)]=2
    value[fixed_freedoms-len(Dim_3)]=-0.0005
    for i in range(1,fixed_freedoms+1):
      no[i-1]=nf[sense[i-1]-1,node[i-1]-1]
    kv[kdiag[no-1]-1]=kv[kdiag[no-1]-1]+1e20
    loads[no-1]=kv[kdiag[no-1]-1]*value
#####方程求解
#call sparin  
A.sparin(kv,kdiag)
##call spabac
A.spabac(kv,loads,kdiag)
loads[neq]=0
print('节点   位移')
for k in range(1,nn+1):
  print(k,end=' ')
  for m in range(1,nodof+1):
    print('{:9.4e}'.format(loads[nf[m-1,k-1]-1,0]),end=' ')
  print( )
##重新取得单元和节点力
axial=0
print('单元  作用力')
for iel in range(1,nels+1):
#call rod_km
    num=g_num[:,iel-1]
    coord=np.transpose(g_coord[:,num-1])
#call pin_jointed
    A.pin_jointed(km,prop[0,int(etype[iel-1])-1],coord)
##########
    g[:,0]=g_g[:,iel-1]
    eld=loads[g-1,0]
    action=np.dot(km,eld)
    print(iel,end=' ')
    for i in range(1,ndof+1):
      print('{:9.4e}'.format(action[i-1,0]),end=' ')
    print( )
#call glob_to_axial
    A.glob_to_axial(axial,action,coord)
    

终端输出结果如下

计算实例2(3维)

单元个数=4，节点个数=5，3维，单元性质类型=1，为5.0e5；
约束点4个，节点1,2,3,4都为3个方向的约束；
节点5具有x方向荷载10KN，y方向荷载-20KN,z方向荷载30.0KN；
节点5具有y方向位移为-0.0005。

代码2(3维)

import numpy as np
import A
nels=4
np_types=1
loaded_nodes=1
fixed_freedoms=1
nod=2
nprops=1
nn=5
ndim=3
nr=4
nodof=ndim
ndof=nod*nodof
g_coord=np.array([[0,1.25,3.5,4,2],[0,3,2,1,1.5],[0,0,0,0,3]])
g_num=np.array([[1,2,3,4],[5,5,5,5]])
#初始化定义数组
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))
eld=np.ones((ndof,1))
km=np.ones((ndof,ndof))
coord=np.ones((nod,ndim),dtype=np.int64)
action=np.ones((ndof,1))
g_g=np.ones((ndof,nels))
prop=np.ones((nprops,np_types))
if np_types==1:
  etype[:,0]=1
else:
  etype[:,0]=(2,2,1,1)
if np_types==1:
  prop[:nprops,np_types-1]=5e5
else:
  prop[:nprops,np_types-2]=2000
  prop[:nprops,np_types-1]=1000
#读取nr
if nr!=0:
  Dim_1=[1,2,3,4]
  nf_value=np.array([[0,0,0,0],[0,0,0,0],[0,0,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]*nf.shape[1],1)))
kdiag=np.zeros((neq,1),dtype=np.int64)
loads=np.zeros((neq+1,1))
##累加单元去发现全局尺寸
for iel in range(1,nels+1):
  num[:,0]=g_num[:,iel-1]
  A.num_to_g(num,nf,g)
  g_g[:,iel-1]=g[:,0]
#call fkidiag
  A.fkdiag(kdiag,g)
for i in range(1,neq):
  kdiag[i]=kdiag[i]+kdiag[i-1]
kv=np.zeros((kdiag[neq-1,0],1),dtype=float)
print('一共有'+str(neq)+'等式和天际线存储个数为'+str(kdiag[neq-1]))
#全局刚度矩阵安装
for iel in range(1,nels+1):
    num=g_num[:,iel-1]
    coord=np.transpose(g_coord[:,num-1])
#call pin_jointed
    A.pin_jointed(km,prop[0,int(etype[iel-1])-1],coord)
##########
    g[:,0]=g_g[:,iel-1]
#call fsparv
    A.fsparv(kv,km,g,kdiag)
###读荷载和位移
if loaded_nodes!=0:
  Dim_2 = [5]
  load_value=np.array([[20],[-20],[30]])
  m=0
  for i in Dim_2:
      for j in range(1,nodof+1):
        loads[nf[j-1,i-1]-1]=load_value[j-1,m]
      m=m+1
if fixed_freedoms!=0:
  Dim_3=[5]
  for i in Dim_3:
    node=np.ones((fixed_freedoms,1),dtype=np.int64)
    sense=np.ones((fixed_freedoms,1),dtype=np.int64)
    value=np.ones((fixed_freedoms,1))
    no=np.ones((fixed_freedoms,1),dtype=np.int64)
####位移赋值
    node[fixed_freedoms-len(Dim_3)]=i
    sense[fixed_freedoms-len(Dim_3)]=2
    value[fixed_freedoms-len(Dim_3)]=-0.0005
    for i in range(1,fixed_freedoms+1):
      no[i-1]=nf[sense[i-1]-1,node[i-1]-1]
    kv[kdiag[no-1]-1]=kv[kdiag[no-1]-1]+1e20
    loads[no-1]=kv[kdiag[no-1]-1]*value
#####方程求解
#call sparin  
A.sparin(kv,kdiag)
##call spabac
A.spabac(kv,loads,kdiag)
loads[neq]=0
print('节点   位移')
for k in range(1,nn+1):
  print(k,end=' ')
  for m in range(1,nodof+1):
    print('{:9.4e}'.format(loads[nf[m-1,k-1]-1,0]),end=' ')
  print( )
##重新取得单元和节点力
axial=0
print('单元  作用力')
for iel in range(1,nels+1):
#call rod_km
    num=g_num[:,iel-1]
    coord=np.transpose(g_coord[:,num-1])
#call pin_jointed
    A.pin_jointed(km,prop[0,int(etype[iel-1])-1],coord)
##########
    g[:,0]=g_g[:,iel-1]
    eld=loads[g-1,0]
    action=np.dot(km,eld)
    print(iel,end=' ')
    for i in range(1,ndof+1):
      print('{:9.4e}'.format(action[i-1,0]),end=' ')
    print( )
#call glob_to_axial
    A.glob_to_axial(axial,action,coord)
    

终端输出结果

(再次申明，本部分编码为本人初学时编成，如有代码复杂化或过程错误欢迎指正)