投稿
  1. 首页
  2. 随笔

矩形d性固体的平面或者轴对称应变分析(python,有限元)

质量控制计划 • 2022-11-10 • 随笔 • 阅读 16

矩形d性固体的平面或者轴对称应变分析(python,有限元),第1张

矩形d性固体的平面或者轴对称应变分析(python,有限元) 第八篇 矩形d性固体的平面或者轴对称应变分析(python,有限元) 引言

此部分程序主要用来解决小应变d性固体的平衡问题。和上一篇相似，当数据输入端提供单元类型和单元编号的放行，节点坐标和序列号会通过子程序产生，矩形结构的子程序为geom_rect，立方体为hexahearon_xz，通过子程序num_to_g由节点号num和每个节点的约束条件nf得到对应的自由度编号，从这部分开始，将提供有限元网格的划分图，包括未变形网格代码mesh，变形网格代码dismsh，位移向量图代码vecmsh。

计算实例1


2维平面单元，网格划分为三角形，3节点，节点自由度沿x轴计数
x方向网格条数为2，y方向为2，积分点为1，单元类型1种，值为d性模量E=1.06，泊松比v=0.3
节点x方向坐标值为0.0、0.5、1.0；y方向坐标值为0.0、-0.5、-1.0
约束节点为5个，分别为1、4节点x方向位移，8、9节点y方向位移，7几点两个方向都约束
荷载节点为3个，1、2、3节点沿y轴负方向值为0.25、0.50、0.25

代码1 
import numpy as np
import A
type_2d='plane'
element='triangle'
dir='x'
nxe=2;nye=2;nze=0
nst=3
np_types=1
loaded_nodes=3
fixed_freedoms=0
nod=3
nprops=2
nr=5
nodof=2
ndim=2
nip=1
nels=np.array([0])
nn=np.array([0])
A.mesh_size(element,nod,nels,nn,nxe,nye,nze)
ndof=nod*nodof
if type_2d=='axisymmetric':
    nst=4
nels=int(nels)
nn=int(nn)
#初始化定义数组
nf=np.ones((nodof,nn),dtype=np.int64)
points=np.ones((nip,ndim))
g=np.ones((ndof,1),dtype=np.int64)
g_coord=np.ones((ndim,nn))
fun=np.ones((nod,1))
coord=np.ones((nod,ndim))
jac=np.ones((ndim,ndim))
g_num=np.ones((nod,nels))
der=np.ones((ndim,nod))
deriv=np.ones((ndim,nod))
bee=np.zeros((nst,ndof))
km=np.ones((ndof,ndof))
eld=np.ones((ndof,1))
weights=np.ones((nip,1))
g_g=np.ones((ndof,nels))
prop=np.ones((nprops,np_types))
num=np.ones((nod,1),dtype=np.int64)
#x_coords=np.ones((nxe+1,1))
#y_coords=np.ones((nxe+1,1))
etype=np.ones((nels,1),dtype=np.int64)
gc=np.ones((ndim,1))
dee=np.zeros((nst,nst))
sigma=np.ones((nst,1))
km=np.zeros((ndof,ndof))
weights=np.ones((nip,1))
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.0e6,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)
x_coords=np.array([0.0,0.5,1.0])
y_coords=np.array([0.0,-0.5,-1.0])
#读取nr
if nr!=0:
    Dim_1=[1,4,7,8,9]
    nf_value=np.array([[0,0,0,1,1],[1,1,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,dir)
    A.num_to_g(num,nf,g)
    g_num[:,iel-1]=num[:,0]
    g_coord[:,num[:,0]-1]=np.transpose(coord[:,:])
    g_g[:,iel-1]=g[:,0]
    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,0]))
#单元刚度积分和安装
A.sample(element,points,weights)
for iel in range(1,nels+1):
    A.deemat(dee,prop[0,etype[iel-1]-1],prop[1,etype[iel-1]-1])
    num[:,0]=g_num[:,iel-1]
    coord[:,:]=np.transpose(g_coord[:,num[:,0]-1])
    g[:,0]=g_g[:,iel-1]
    km[:]=0
    for i in range(1,nip+1):
        A.shape_fun(fun,points,i)
        A.shape_der(der,points,i)
        jac=np.dot(der,coord)
        det=np.linalg.det(jac)
        A.invert(jac)
        deriv=np.dot(jac,der)
        A.beemat(bee,deriv)
        if type_2d=='axisymmetric':
            gc=np.dot(fun,coord)
            bee[3,0:ndof-1:2]=fun[:]/gc[0]
        km[:]=km[:]+np.dot(np.dot(np.transpose(bee),dee),bee)*det*weights[i-1]*gc[0]
#call fsparv
    A.fsparv(kv,km,g,kdiag)
###读荷载和位移
if loaded_nodes!=0:
  Dim_2 = [1,2,3]
  load_value=np.array([[0.0,0.0,0.0],[-0.25,-0.50,-0.25]])
  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
if type_2d=='axisymmetric':
    print('节点 r-位移 z-位移')
else:
    print('节点 x-位移   y-位移')
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))
weights=np.ones((nip,1))
A.sample(element,points,weights)
print('积分点nip='+str(nip)+'应力为')
if type_2d=='axisymmetric':
    print('单元 r-坐标 z-坐标')
    print('sig_r sig_z tau_rz sig_t')
else:
    print('单元 x-坐标     y-坐标       sig_x       sig_y       tau_xy')
for iel in range(1,nels+1):
  A.deemat(dee,prop[0,etype[iel-1]-1],prop[1,etype[iel-1]-1])
  num[:,0]=g_num[:,iel-1]
  coord[:,:]=np.transpose(g_coord[:,num[:,0]-1])
  g[:,0]=g_g[:,iel-1]
  eld=loads[g-1,0]
  for i in range(1,nip+1):
    A.shape_fun(fun,points,i)
    A.shape_der(der,points,i)
    gc=np.dot(np.transpose(fun),coord)
    jac=np.dot(der,coord)
    A.invert(jac)
    deriv=np.dot(jac,der)
    A.beemat(bee,deriv)
    if type_2d=='axisymmetric':
        gc=np.dot(fun,coord)
        bee[3,0:ndof-1:2]=fun[:]/gc[0]
    sigma[:]=np.dot(dee[:],np.dot(bee[:],eld[:]))
    print(iel,end='  ')
    for i in range(1,ndim+1):
      print('{:9.4e}'.format(gc[0,i-1]),end='  ')
    for i in range(1,nst+1):
      print('{:9.4e}'.format(sigma[i-1,0]),end='  ')
    print( )
终端输出结果1

网格图1

1.未变形图

2.变形图

3.位移矢量图

计算实例2



2维平面单元，网格划分为三角形，15节点，节点自由度沿y轴计数
x方向网格条数为2，y方向为1，积分点为12，单元类型1种，值为d性模量E=1.05，泊松比v=0.2
节点x方向坐标值为0.0、1.0、6.0；y方向坐标值为0.0、-2.0
约束节点为17个，0为约束，1为自由，例如节点1x方向约束，y方向自由
荷载节点为5个，1、6、11、16、21节点沿y轴负方向值为0.0778、0.3556、0.1333、0.3556、0.0778

代码2 
import numpy as np
import A
type_2d='plane'
element='triangle'
dir='y'
nxe=2;nye=1;nze=0
nst=3
np_types=1
loaded_nodes=5
fixed_freedoms=0
nod=15
nprops=2
nr=17
nodof=2
ndim=2
nip=12
nels=np.array([0])
nn=np.array([0])
A.mesh_size(element,nod,nels,nn,nxe,nye,nze)
ndof=nod*nodof
if type_2d=='axisymmetric':
    nst=4
nels=int(nels)
nn=int(nn)
#初始化定义数组
nf=np.ones((nodof,nn),dtype=np.int64)
points=np.ones((nip,ndim))
g=np.ones((ndof,1),dtype=np.int64)
g_coord=np.ones((ndim,nn))
fun=np.ones((nod,1))
coord=np.ones((nod,ndim))
jac=np.ones((ndim,ndim))
g_num=np.ones((nod,nels))
der=np.ones((ndim,nod))
deriv=np.ones((ndim,nod))
bee=np.zeros((nst,ndof))
km=np.ones((ndof,ndof))
eld=np.ones((ndof,1))
weights=np.ones((nip,1))
g_g=np.ones((ndof,nels))
prop=np.ones((nprops,np_types))
num=np.ones((nod,1),dtype=np.int64)
#x_coords=np.ones((nxe+1,1))
#y_coords=np.ones((nxe+1,1))
etype=np.ones((nels,1),dtype=np.int64)
gc=np.ones((ndim,1))
dee=np.zeros((nst,nst))
sigma=np.ones((nst,1))
km=np.zeros((ndof,ndof))
weights=np.ones((nip,1))
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.0e5,0.2)
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)
x_coords=np.array([0.0,1.0,6.0])
y_coords=np.array([0.0,-2.0])
#读取nr
if nr!=0:
    Dim_1=[1,2,3,4,5,10,15,20,25,30,35,40,41,42,43,44,45]
    nf_value=np.array([[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],[1,1,1,1,0,0,0,0,0,0,0,0,1,1,1,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)
loads=np.zeros((neq+1,1))
##累加单元去发现全局尺寸
for iel in range(1,nels+1):
    A.geom_rect(element,iel,x_coords,y_coords,coord,num,dir)
    A.num_to_g(num,nf,g)
    g_num[:,iel-1]=num[:,0]
    g_coord[:,num[:,0]-1]=np.transpose(coord[:,:])
    g_g[:,iel-1]=g[:,0]
    A.fkdiag(kdiag,g)
A.mesh(g_coord,g_num)
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):
    A.deemat(dee,prop[0,etype[iel-1]-1],prop[1,etype[iel-1]-1])
    num[:,0]=g_num[:,iel-1]
    coord[:,:]=np.transpose(g_coord[:,num[:,0]-1])
    g[:,0]=g_g[:,iel-1]
    km[:]=0
    for i in range(1,nip+1):
        A.shape_fun(fun,points,i)
        A.shape_der(der,points,i)
        jac=np.dot(der,coord)
        det=np.linalg.det(jac)
        A.invert(jac)
        deriv=np.dot(jac,der)
        A.beemat(bee,deriv)
        if type_2d=='axisymmetric':
            gc=np.dot(fun,coord)
            bee[3,0:ndof-1:2]=fun[:]/gc[0]
        km[:]=km[:]+np.dot(np.dot(np.transpose(bee),dee),bee)*det*weights[i-1]*gc[0]
#call fsparv
    A.fsparv(kv,km,g,kdiag)
###读荷载和位移
if loaded_nodes!=0:
  Dim_2 = [1,6,11,16,21]
  load_value=np.array([[0.0,0.0,0.0,0.0,0.0],[-0.0778,-0.3556,-0.1333,-0.3556,-0.0778]])
  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
if type_2d=='axisymmetric':
    print('节点 r-位移 z-位移')
else:
    print('节点 x-位移   y-位移')
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))
weights=np.ones((nip,1))
A.sample(element,points,weights)
print('积分点nip='+str(nip)+'应力为')
if type_2d=='axisymmetric':
    print('单元 r-坐标 z-坐标')
    print('sig_r sig_z tau_rz sig_t')
else:
    print('单元 x-坐标     y-坐标       sig_x       sig_y       tau_xy')
for iel in range(1,nels+1):
  A.deemat(dee,prop[0,etype[iel-1]-1],prop[1,etype[iel-1]-1])
  num[:,0]=g_num[:,iel-1]
  coord[:,:]=np.transpose(g_coord[:,num[:,0]-1])
  g[:,0]=g_g[:,iel-1]
  eld=loads[g-1,0]
  for i in range(1,nip+1):
    A.shape_fun(fun,points,i)
    A.shape_der(der,points,i)
    gc=np.dot(np.transpose(fun),coord)
    jac=np.dot(der,coord)
    A.invert(jac)
    deriv=np.dot(jac,der)
    A.beemat(bee,deriv)
    if type_2d=='axisymmetric':
        gc=np.dot(fun,coord)
        bee[3,0:ndof-1:2]=fun[:]/gc[0]
    sigma[:]=np.dot(dee[:],np.dot(bee[:],eld[:]))
    print(iel,end='  ')
    for i in range(1,ndim+1):
      print('{:9.4e}'.format(gc[0,i-1]),end='  ')
    for i in range(1,nst+1):
      print('{:9.4e}'.format(sigma[i-1,0]),end='  ')
    print( )
A.dismsh(g_coord,g_num,loads,nf)
A.vecmsh(loads,nf,0.1,0.05,g_coord,g_num)
终端输出结果2


网格划分图

1.未变形图

2.变形图

3.位移矢量图

计算实例3



2维平面单元，网格划分四边形，4节点，节点自由度沿y轴计数
x方向网格条数为3，y方向为2，积分点为4，单元类型1种，值为d性模量E=1.06，泊松比v=0.3
节点x方向坐标值为0.0、10.0、20.0、30.0；y方向坐标值为0.0、-5.0、-10.0
约束节点为8个，0为约束，1为自由，例如节点1x方向约束，y方向自由
固定节点位移2个，1、4节点沿y轴负方向值为1.0e-5,1.0e-5

代码3 
import numpy as np
import A
type_2d='plane'
element='quadrilateral'
dir='y'
nxe=3;nye=2;nze=0
nst=3
np_types=1
loaded_nodes=0
fixed_freedoms=2
nod=4
nprops=2
nr=8
nodof=2
ndim=2
nip=4
nels=np.array([0])
nn=np.array([0])
A.mesh_size(element,nod,nels,nn,nxe,nye,nze)
ndof=nod*nodof
if type_2d=='axisymmetric':
    nst=4
nels=int(nels)
nn=int(nn)
#初始化定义数组
nf=np.ones((nodof,nn),dtype=np.int64)
points=np.ones((nip,ndim))
g=np.ones((ndof,1),dtype=np.int64)
g_coord=np.ones((ndim,nn))
fun=np.ones((nod,1))
coord=np.ones((nod,ndim))
jac=np.ones((ndim,ndim))
g_num=np.ones((nod,nels))
der=np.ones((ndim,nod))
deriv=np.ones((ndim,nod))
bee=np.zeros((nst,ndof))
km=np.ones((ndof,ndof))
eld=np.ones((ndof,1))
weights=np.ones((nip,1))
g_g=np.ones((ndof,nels))
prop=np.ones((nprops,np_types))
num=np.ones((nod,1),dtype=np.int64)
#x_coords=np.ones((nxe+1,1))
#y_coords=np.ones((nxe+1,1))
etype=np.ones((nels,1),dtype=np.int64)
gc=np.ones((ndim,1))
dee=np.zeros((nst,nst))
sigma=np.ones((nst,1))
km=np.zeros((ndof,ndof))
weights=np.ones((nip,1))
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.0e6,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)
x_coords=np.array([0.0,10.0,20.0,30.0])
y_coords=np.array([0.0,-5.0,-10.0])
#读取nr
if nr!=0:
    Dim_1=[1,2,3,6,9,10,11,12]
    nf_value=np.array([[0,0,0,0,0,0,0,0],[1,1,0,0,0,1,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)
loads=np.zeros((neq+1,1))
##累加单元去发现全局尺寸
for iel in range(1,nels+1):
    A.geom_rect(element,iel,x_coords,y_coords,coord,num,dir)
    A.num_to_g(num,nf,g)
    g_num[:,iel-1]=num[:,0]
    g_coord[:,num[:,0]-1]=np.transpose(coord[:,:])
    g_g[:,iel-1]=g[:,0]
    A.fkdiag(kdiag,g)
A.mesh(g_coord,g_num)
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):
    A.deemat(dee,prop[0,etype[iel-1]-1],prop[1,etype[iel-1]-1])
    num[:,0]=g_num[:,iel-1]
    coord[:,:]=np.transpose(g_coord[:,num[:,0]-1])
    g[:,0]=g_g[:,iel-1]
    km[:]=0
    for i in range(1,nip+1):
        A.shape_fun(fun,points,i)
        A.shape_der(der,points,i)
        jac=np.dot(der,coord)
        det=np.linalg.det(jac)
        A.invert(jac)
        deriv=np.dot(jac,der)
        A.beemat(bee,deriv)
        if type_2d=='axisymmetric':
            gc=np.dot(fun,coord)
            bee[3,0:ndof-1:2]=fun[:]/gc[0]
        km[:]=km[:]+np.dot(np.dot(np.transpose(bee),dee),bee)*det*weights[i-1]*gc[0]
#call fsparv
    A.fsparv(kv,km,g,kdiag)
###读荷载和位移
if loaded_nodes!=0:
  Dim_2 = [1,2,3]
  load_value=np.array([[0.0,0.0,0.0],[-0.25,-0.50,-0.25]])
  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,4)
  sense[:,0]=(2,2)
  value[:,0]=(-1.0e-5,-1.0e-5)
####位移赋值
  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
if type_2d=='axisymmetric':
    print('节点 r-位移 z-位移')
else:
    print('节点 x-位移   y-位移')
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))
weights=np.ones((nip,1))
A.sample(element,points,weights)
print('积分点nip='+str(nip)+'应力为')
if type_2d=='axisymmetric':
    print('单元 r-坐标 z-坐标')
    print('sig_r sig_z tau_rz sig_t')
else:
    print('单元 x-坐标     y-坐标       sig_x       sig_y       tau_xy')
for iel in range(1,nels+1):
  A.deemat(dee,prop[0,etype[iel-1]-1],prop[1,etype[iel-1]-1])
  num[:,0]=g_num[:,iel-1]
  coord[:,:]=np.transpose(g_coord[:,num[:,0]-1])
  g[:,0]=g_g[:,iel-1]
  eld=loads[g-1,0]
  for i in range(1,nip+1):
    A.shape_fun(fun,points,i)
    A.shape_der(der,points,i)
    gc=np.dot(np.transpose(fun),coord)
    jac=np.dot(der,coord)
    A.invert(jac)
    deriv=np.dot(jac,der)
    A.beemat(bee,deriv)
    if type_2d=='axisymmetric':
        gc=np.dot(fun,coord)
        bee[3,0:ndof-1:2]=fun[:]/gc[0]
    sigma[:]=np.dot(dee[:],np.dot(bee[:],eld[:]))
    print(iel,end='  ')
    for i in range(1,ndim+1):
      print('{:9.4e}'.format(gc[0,i-1]),end='  ')
    for i in range(1,nst+1):
      print('{:9.4e}'.format(sigma[i-1,0]),end='  ')
    print( )
A.dismsh(g_coord,g_num,loads,nf)
A.vecmsh(loads,nf,0.1,0.05,g_coord,g_num)
终端输出结果3

网格划分图3

1.未变形图

2.变形图

3.位移矢量图

计算实例4



2维平面单元，网格划分四边形，8节点，节点自由度沿x轴计数
x方向网格条数为2，y方向为3，积分点为4，单元类型1种，值为d性模量E=1.06，泊松比v=0.3
节点x方向坐标值为0.0、3.0、6.0；y方向坐标值为0.0、-3.0、-6.0、-9.0
约束节点为17个，0为约束，1为自由，例如节点1x方向约束，y方向自由
荷载节点3个，1、2、3节点沿y轴负方向值为0.5、2.0、0.5

代码4 
import numpy as np
import A
type_2d='plane'
element='quadrilateral'
dir='x'
nxe=2;nye=3;nze=0
nst=3
np_types=1
loaded_nodes=3
fixed_freedoms=0
nod=8
nprops=2
nr=17
nodof=2
ndim=2
nip=4
nels=np.array([0])
nn=np.array([0])
A.mesh_size(element,nod,nels,nn,nxe,nye,nze)
ndof=nod*nodof
if type_2d=='axisymmetric':
    nst=4
nels=int(nels)
nn=int(nn)
#初始化定义数组
nf=np.ones((nodof,nn),dtype=np.int64)
points=np.ones((nip,ndim))
g=np.ones((ndof,1),dtype=np.int64)
g_coord=np.ones((ndim,nn))
fun=np.ones((nod,1))
coord=np.ones((nod,ndim))
jac=np.ones((ndim,ndim))
g_num=np.ones((nod,nels))
der=np.ones((ndim,nod))
deriv=np.ones((ndim,nod))
bee=np.zeros((nst,ndof))
km=np.ones((ndof,ndof))
eld=np.ones((ndof,1))
weights=np.ones((nip,1))
g_g=np.ones((ndof,nels))
prop=np.ones((nprops,np_types))
num=np.ones((nod,1),dtype=np.int64)
#x_coords=np.ones((nxe+1,1))
#y_coords=np.ones((nxe+1,1))
etype=np.ones((nels,1),dtype=np.int64)
gc=np.ones((ndim,1))
dee=np.zeros((nst,nst))
sigma=np.ones((nst,1))
km=np.zeros((ndof,ndof))
weights=np.ones((nip,1))
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.0e6,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)
x_coords=np.array([0.0,3.0,6.0])
y_coords=np.array([0.0,-3.0,-6.0,-9.0])
#读取nr
if nr!=0:
    Dim_1=[1,6,9,14,17,22,25,26,27,28,5,8,13,16,21,24,29]
    nf_value=np.array([[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],[1,1,1,1,1,1,0,0,0,0,1,1,1,1,1,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)
loads=np.zeros((neq+1,1))
##累加单元去发现全局尺寸
for iel in range(1,nels+1):
    A.geom_rect(element,iel,x_coords,y_coords,coord,num,dir)
    A.num_to_g(num,nf,g)
    g_num[:,iel-1]=num[:,0]
    g_coord[:,num[:,0]-1]=np.transpose(coord[:,:])
    g_g[:,iel-1]=g[:,0]
    A.fkdiag(kdiag,g)
A.mesh(g_coord,g_num)
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):
    A.deemat(dee,prop[0,etype[iel-1]-1],prop[1,etype[iel-1]-1])
    num[:,0]=g_num[:,iel-1]
    coord[:,:]=np.transpose(g_coord[:,num[:,0]-1])
    g[:,0]=g_g[:,iel-1]
    km[:]=0
    for i in range(1,nip+1):
        A.shape_fun(fun,points,i)
        A.shape_der(der,points,i)
        jac=np.dot(der,coord)
        det=np.linalg.det(jac)
        A.invert(jac)
        deriv=np.dot(jac,der)
        A.beemat(bee,deriv)
        if type_2d=='axisymmetric':
            gc=np.dot(fun,coord)
            bee[3,0:ndof-1:2]=fun[:]/gc[0]
        km[:]=km[:]+np.dot(np.dot(np.transpose(bee),dee),bee)*det*weights[i-1]*gc[0]
#call fsparv
    A.fsparv(kv,km,g,kdiag)
###读荷载和位移
if loaded_nodes!=0:
  Dim_2 = [1,2,3]
  load_value=np.array([[0.0,0.0,0.0],[-0.5,-2.0,-0.5]])
  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,4)
  sense[:,0]=(2,2)
  value[:,0]=(-1.0e-5,-1.0e-5)
####位移赋值
  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
if type_2d=='axisymmetric':
    print('节点 r-位移 z-位移')
else:
    print('节点 x-位移   y-位移')
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))
weights=np.ones((nip,1))
A.sample(element,points,weights)
print('积分点nip='+str(nip)+'应力为')
if type_2d=='axisymmetric':
    print('单元 r-坐标 z-坐标')
    print('sig_r sig_z tau_rz sig_t')
else:
    print('单元 x-坐标     y-坐标       sig_x       sig_y       tau_xy')
for iel in range(1,nels+1):
  A.deemat(dee,prop[0,etype[iel-1]-1],prop[1,etype[iel-1]-1])
  num[:,0]=g_num[:,iel-1]
  coord[:,:]=np.transpose(g_coord[:,num[:,0]-1])
  g[:,0]=g_g[:,iel-1]
  eld=loads[g-1,0]
  for i in range(1,nip+1):
    A.shape_fun(fun,points,i)
    A.shape_der(der,points,i)
    gc=np.dot(np.transpose(fun),coord)
    jac=np.dot(der,coord)
    A.invert(jac)
    deriv=np.dot(jac,der)
    A.beemat(bee,deriv)
    if type_2d=='axisymmetric':
        gc=np.dot(fun,coord)
        bee[3,0:ndof-1:2]=fun[:]/gc[0]
    sigma[:]=np.dot(dee[:],np.dot(bee[:],eld[:]))
    print(iel,end='  ')
    for i in range(1,ndim+1):
      print('{:9.4e}'.format(gc[0,i-1]),end='  ')
    for i in range(1,nst+1):
      print('{:9.4e}'.format(sigma[i-1,0]),end='  ')
    print( )
A.dismsh(g_coord,g_num,loads,nf)
A.vecmsh(loads,nf,0.1,0.05,g_coord,g_num)
终端输出结果4

网格划分图4

1.未变形图

2.变形图

3.位移矢量图

欢迎分享，转载请注明来源：内存溢出

原文地址: http://outofmemory.cn/zaji/4830396.html

赞 (0)
打赏 微信扫一扫 微信扫一扫 支付宝扫一扫 支付宝扫一扫
质量控制计划 质量控制计划 一级用户组
0 0
上一篇 2022-11-10
下一篇 2022-11-10

发表评论

登录后才能评论

评论列表（0条）

保存