2节点弹性梁的稳定性分析(弹性基础)(python,有限元)

此程序计算遭受轴向压缩荷载的基础屈曲荷载和相关联梁的模态形状，本程序与第三篇相似但包括弹性选项，子程序beam_gm产生单元几何矩阵gm，组装几何矩阵储存在向量gv中，梁刚度矩阵由km和mm组装为kv矩阵。
子程序stability使用简单迭代去得到最小特征值eval，对应最低的屈曲荷载，和对应的特征向量evec。其他相关联的标量tol为收敛容差，iters为迭代次数，limit为迭代上限。

第一个例子为单位长度梁，一端固定一端铰接。

其中单元个数=4，单元性质种类为1种，性质类型也只有1种EI，其值为1.0
单元长度分别为0.25，0.25，0.25，0.25
约束点2个，节点1x方向的约束,节点5两个方向都约束；
迭代最高次数为100，收敛容差为1.0e-5。

import numpy as np
import A
nels=4
np_types=1
loaded_nodes=1
fixed_freedoms=0
nod=2
nprops=1
nn=nels+1
nr=2
nodof=2
ndof=nod*nodof
limit=100
tol=1e-5
ell=np.array([[0.25],[0.25],[0.25],[0.25]])
#初始化定义数组
iters=np.array([0])
eval1=np.array([0])
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)
km=np.ones((ndof,ndof))
gm=np.ones((ndof,ndof))
mm=np.zeros((ndof,ndof))
g_g=np.ones((ndof,nels))
prop=np.ones((nprops,np_types))
if np_types==1:
  etype[:,0]=1
else:
  etype[:,0]=(1,2,3,4,5)
if np_types==1:
  prop[:nprops,np_types-1]=(1.0)
else:
  prop[0,:]=(1.92e4,1.92e4,1.92e4,1.92e4,1.92e4)
  prop[1,:]=(0.2,0.6,1.0,1.4,1.8)
#读取nr
if nr!=0:
  Dim_1=[1,5]
  nf_value=np.array([[0,0],[1,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)
evec=np.zeros((neq+1,1))
##累加单元去发现全局尺寸
for iel in range(1,nels+1):
  num[:nod,0]=(iel,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)
gv=np.zeros((kdiag[neq-1,0],1),dtype=float)
print('一共有'+str(neq)+'等式和天际线存储个数为'+str(kdiag[neq-1]))
#全局刚度矩阵安装
for iel in range(1,nels+1):
#call pin_jointed
  A.beam_km(km,prop[0,etype[iel-1,0]-1],ell[iel-1,0])
  g[:,0]=g_g[:,iel-1]
  if nprops>1:
    A.beam_mm(mm,prop[1,etype[iel-1,0]-1],ell[iel-1,0])
  A.beam_gm(gm,ell[iel-1])
#call fsparv
  A.fsparv(kv,km+mm,g,kdiag)
  A.fsparv(gv,gm,g,kdiag)
###解特征值问题
A.stability(kv,gv,kdiag,tol,limit,iters,evec,eval1)
print('屈曲荷载='+str(eval1[0]))
evec[neq]=0
print('屈曲模态')
print('节点 位移 转角')
for k in range(1,nn+1):
  print(k,end=' ')
  for m in range(1,nodof+1):
    print('{:9.4e}'.format(evec[nf[m-1,k-1]-1,0]),end=' ')
  print( )
print('收敛在'+str(iters[0])+'次迭代')

第二个例子为一个弹性基础的简单支持梁
其中单元个数=4，单元性质为1种，性质类型也只有2种，EI=1.0,k=800.0
单元长度分别为0.25，0.25，0.25，0.25
约束点2个，节点1x方向的约束,节点5x方向的约束；
迭代最高次数为100，收敛容差为1.0e-5。

import numpy as np
import A
nels=4
np_types=1
loaded_nodes=1
fixed_freedoms=0
nod=2
nprops=2
nn=nels+1
nr=2
nodof=2
ndof=nod*nodof
limit=100
tol=1e-5
ell=np.array([[0.25],[0.25],[0.25],[0.25]])
#初始化定义数组
iters=np.array([0])
eval1=np.array([0])
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)
km=np.ones((ndof,ndof))
gm=np.ones((ndof,ndof))
mm=np.zeros((ndof,ndof))
g_g=np.ones((ndof,nels))
prop=np.ones((nprops,np_types))
if np_types==1:
  etype[:,0]=1
else:
  etype[:,0]=(1,2,3,4,5)
if np_types==1:
  prop[:nprops,np_types-1]=(1.0,800.0)
else:
  prop[0,:]=(1.92e4,1.92e4,1.92e4,1.92e4,1.92e4)
  prop[1,:]=(0.2,0.6,1.0,1.4,1.8)
#读取nr
if nr!=0:
  Dim_1=[1,5]
  nf_value=np.array([[0,0],[1,1]])
  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)
evec=np.zeros((neq+1,1))
##累加单元去发现全局尺寸
for iel in range(1,nels+1):
  num[:nod,0]=(iel,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)
gv=np.zeros((kdiag[neq-1,0],1),dtype=float)
print('一共有'+str(neq)+'等式和天际线存储个数为'+str(kdiag[neq-1]))
#全局刚度矩阵安装
for iel in range(1,nels+1):
#call pin_jointed
  A.beam_km(km,prop[0,etype[iel-1,0]-1],ell[iel-1,0])
  g[:,0]=g_g[:,iel-1]
  if nprops>1:
    A.beam_mm(mm,prop[1,etype[iel-1,0]-1],ell[iel-1,0])
  A.beam_gm(gm,ell[iel-1])
#call fsparv
  A.fsparv(kv,km+mm,g,kdiag)
  A.fsparv(gv,gm,g,kdiag)
###解特征值问题
A.stability(kv,gv,kdiag,tol,limit,iters,evec,eval1)
print('屈曲荷载='+str(eval1[0]))
evec[neq]=0
print('屈曲模态')
print('节点 位移 转角')
for k in range(1,nn+1):
  print(k,end=' ')
  for m in range(1,nodof+1):
    print('{:9.4e}'.format(evec[nf[m-1,k-1]-1,0]),end=' ')
  print( )
print('收敛在'+str(iters[0])+'次迭代')