第七篇 4节点矩形平面单元分析-相同单元的均质材料(python,有限元)

引言

从本程序开始接触平面单元，程序种加入一些新的参数值，包括矩形单元的长aa和宽bb，单元厚度th
子程序geom_rect的作用是根据已知信息得出每个节点的坐标和对应的自由度序列号，该子程序在下一部分内容中也会广泛使用。

计算实例

由上可知，x方向分为两端，y方向分为两段，数据类型为1种
x方向单元长度为0.25，y方向单元长度为0.25，厚度为1.0
单元性质包含弹性模量E，泊松比V，分别为10.92和0.3
有8个约束点，0表示约束，1表示自由，例如对于节点1，只有xy方向有转角
节点9有垂直向下的力0.25
没有固定位移

代码

import numpy as np
import A
nxe=2
nye=2
nze=0
np_types=1
aa=0.25
bb=0.25
th=1.0
loaded_nodes=1
fixed_freedoms=0
nod=4
nprops=2
nr=8
nodof=4
ndim=2
nip=16
ndof=nod*nodof
element='quadrilateral'
nels=np.array([0])
nn=np.array([0])
A.mesh_size(element,nod,nels,nn,nxe,nye,nze)
nels=int(nels)
nn=int(nn)
#初始化定义数组
nf=np.ones((nodof,nn),dtype=np.int64)
g_coord=np.ones((ndim,nn))
g_num=np.ones((nod,nels))
g=np.ones((ndof,1),dtype=np.int64)
bm=np.zeros((3,1))
g_g=np.ones((ndof,nels))
coord=np.ones((nod,ndim))
km=np.zeros((ndof,ndof))
dtd=np.ones((ndof,ndof))
d2x=np.ones((ndof,1))
d2y=np.ones((ndof,1))
d2xy=np.ones((ndof,1))
num=np.ones((nod,1),dtype=np.int64)
x_coords=np.ones((nxe+1,1))
y_coords=np.ones((nxe+1,1))
points=np.ones((nip,ndim))
weights=np.ones((nip,1))
prop=np.ones((nprops,np_types))
etype=np.ones((nels,1),dtype=np.int64)
if np_types==1:
  etype[:,0]=1
else:
  etype[:,0]=(1,2,3,4,5)
if np_types==1:
  prop[:nprops,np_types-1]=(10.92,0.3)
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)
for i in range(1,nxe+2):
    x_coords[i-1]=(i-1)*aa
for i in range(1,nye+2):
    y_coords[i-1]=(i-1)*bb
#读取nr
if nr!=0:
  Dim_1=[1,2,3,4,6,7,8,9]
  nf_value=np.array([[0,0,0,0,1,0,1,1],[0,0,0,1,0,1,1,0],[0,1,1,0,1,0,0,0],[1,1,0,1,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):
    A.geom_rect(element,iel,x_coords,y_coords,coord,num,'x')
    A.num_to_g(num,nf,g)
    A.fkdiag(kdiag,g)
    g_num[:,iel-1]=num[:,0]
    g_coord=np.transpose(coord)
    g_g[:,iel-1]=g[:,0]
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,0]))
#单元刚度积分和安装
A.sample(element,points,weights)
for iel in range(1,nels+1):
    e=prop[0,etype[iel-1]-1]
    v=prop[1,etype[iel-1]-1]
    d=e*th**3/(12*(1-v**2))
    g[:,0]=g_g[:,iel-1]
    km[:]=0
    for i in range(1,nip+1):
        A.fmplat(d2x,d2y,d2xy,points,aa,bb,i)
        for k in range(1,ndof+1):
            dtd[k-1,:]=4*aa*bb*d*weights[i-1]*(d2x[k-1,0]*d2x[:,0]/(aa**4)+d2y[k-1,0]*d2y[:,0]/(bb**4)+   \
                (v*d2x[k-1,0]*d2y[:,0]+v*d2x[:,0]*d2y[k-1,0]+2.0*(1.0-v)*d2xy[k-1,0]*d2xy[:,0])/(aa**2*bb**2)) 
            dtd[:,k-1]=dtd[k-1,:]
            #print(km[15,15])
            #print(dtd[15,15])
        km=km+dtd
#call fsparv
    A.fsparv(kv,km,g,kdiag)
###读荷载和位移
if loaded_nodes!=0:
  Dim_2 = [9]
  load_value=np.array([[0.25],[0],[0],[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:
  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[:,0]=(1,3)
  sense[:,0]=(2,1)
  value[:,0]=(-0.001,-0.005)
####位移赋值
  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,0]=kv[kdiag[no-1,0]-1,0]*value
##荷载增量累加
A.sparin(kv,kdiag)
A.spabac(kv,loads,kdiag)
loads[neq]=0
print('节点 位移 x转角 y转角 xy转角')
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()
#取得单元中心点的力矩
nip=1
points=np.ones((nip,ndim))
A.sample(element,points,weights)
print('单元 x-力矩 y-力矩 xy-力矩')
for iel in range(1,nels+1):
    g[:,0]=g_g[:,iel-1]
    for i in range(1,nip+1):
        A.fmplat(d2x,d2y,d2xy,points,aa,bb,i)
        bm[:]=0
        for k in range(1,ndof+1):
            bm[0]=bm[0]+4.0*d*(d2x[k-1]/aa/aa+v*d2y[k-1]/bb/bb)*loads[g[k-1,0]-1,0]
            bm[1]=bm[1]+4.0*d*(v*d2x[k-1]/aa/aa+d2y[k-1]/bb/bb)*loads[g[k-1,0]-1,0]
            bm[2]=bm[2]+4.0*d*(1.0-v)*d2xy[k-1]/aa/bb*loads[g[k-1,0]-1,0]
    print(iel,end=' ')
    for m in range(1,4):
      print('{:9.4e}'.format(bm[m-1,0]),end=' ')
    print()

终端输出结果