本篇讨论一个轴对称或者反轴对称的固体遭受一个非对称的简谐波荷载的情况
将引入三个新的参数，lth表示表示承受的对应简谐波荷载的振幅值，iflag值为1表示轴对称，-1表示反向轴对称，chi表示平面所表示的二维图形在轴对称固体中的位置，用度数表示
一个通用的轴对称偏导矩阵B为
径向划分为1，竖向网格划分为5
轴对称，简谐波荷载为1，所选平面为主视图所在面0度，单元只有一种类型
性质为弹性模量E1.0e5，泊松比0.3
径向坐标0.0，0.5；竖向坐标10.0、8.0、6.0、4.0、2.0、0.0
约束节点有13个，0为约束，1为自由，例如1节点，径向约束，竖向自由
荷载节点1个，3节点径向力0.3183，为1/pai
import numpy as np
import A
import math
type_2d='plane'
element='quadrilateral'
dir=='x'
nxe=1;nye=5;nze=0
nst=6
np_types=1
loaded_nodes=1
fixed_freedoms=0
nod=8
nprops=2
nr=13
nodof=3
ndim=2
nip=4
lth=1
iflag=1
chi=0.0
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))
radius=np.ones((1,1))
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.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])
y_coords=np.array([10.0,8.0,6.0,4.0,2.0,0.0])
#读取nr
if nr!=0:
Dim_1=[1,4,6,9,11,14,16,19,21,24,26,27,28]
nf_value=np.array([[1,1,1,1,1,1,1,1,1,1,0,0,0],[0,0,0,0,0,0,0,0,0,0,0,0,0],[1,1,1,1,1,1,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))
pi=math.acos(-1)
chi=chi*pi/180
ca=math.cos(chi)
sa=math.sin(chi)
##累加单元去发现全局尺寸
for iel in range(1,nels+1):
A.geom_rect(element,iel,x_coords,y_coords,coord,num,'r')
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]))
#单元刚度积分和安装
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.bmat_nonaxi(bee,radius,coord,deriv,fun,iflag,lth)
km[:]=km[:]+np.dot(np.dot(np.transpose(bee),dee),bee)*det*weights[i-1]*radius[0,0]
#call fsparv
A.fsparv(kv,km,g,kdiag)
###读荷载和位移
if loaded_nodes!=0:
Dim_2 = [3]
load_value=np.array([[0.3183],[0.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('节点 r-位移 z-位移 t-位移')
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)+'应力为')
print('单元 r-坐标 z-坐标 sig_r sig_z sig_t tau_rz tau_zt tau_tr')
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.bmat_nonaxi(bee,radius,coord,deriv,fun,iflag,lth)
bee[0:4,:]=bee[0:4,:]*ca
bee[4:6,:]=bee[4:6,:]*sa
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( )
终端输出结果
