代码收藏家 计算力学|采用python进行有限元模拟
从abaqus输出的inp文件中读取节点和单元信息
import meshio
mesh = meshio.read('Job-3.inp')
coords = mesh.points###coords即为各个节点的坐标
Edof = mesh.cells_dict['triangle']#Edof为三角形单元的节点号
1.单元刚度矩阵
def element_stiffness(n1,coords,E,v,t):
node1 = coords[n1[0]]
node2 = coords[n1[1]]
node3 = coords[n1[2]]
xi = node1[0]
yi = node1[1]
xj = node2[0]
yj = node2[1]
xm = node3[0]
ym = node3[1]
A=abs((xi*(yj-ym)+xj*(ym-yi)+xm*(yi-yj))*0.5)
bi=yj-ym
bj=ym-yi
bm=yi-yj
ci=-xi+xm
cj=-xm+xi
cm=-xi+xj
B=np.array([[bi,0,bj,0,bm,0],[0,ci,0,cj,0,cm],[ci,bi,cj,bj,cm,bm]])/(2*A)
BT=B.transpose()
D=np.array([[1,v,0],[v,1,0],[0,0,(1-v)/2]])*E/(1-v*v)
Ke=BT.dot(D).dot(B)*t*A
return Ke
2.总体刚度矩阵
def assemble_stiffness(coords,edof,E,v,t):
NoN = np.size(coords, 0)
NoE = np.size(edof, 0)
PD = np.size(coords, 1)
NPE = np.size(edof, 1)
K=np.zeros([NoN*PD,NoN*PD])
for i in range(len(edof)):
n1 = edof[i, 0:(NPE)]
Ke=element_stiffness(i,n1,coords, E,v,t)
for r in range(NPE):
for p in range(NPE):
K[2*n1[r] , 2*n1[p] ] = K[2*n1[r] – 0, 2*n1[p]-0] + Ke[2*r, 2*p]
K[2*n1[r] , 2*n1[p]+1] = K[2*n1[r] , 2*n1[p]+1] + Ke[2*r, 2*p+1]
K[2*n1[r] +1, 2*n1[p] ] = K[2*n1[r]+1, 2*n1[p] ] + Ke[2*r+1, 2*p]
K[2*n1[r] +1, 2*n1[p]+1] = K[2*n1[r]+1 , 2*n1[p]+1] + Ke[2*r+1, 2*p+1]
return K
3.施加载荷
在薄板的左右两端施加均布载荷，首先根据单元节点横坐标判断在左右两端，判断在左右两端的单元个数，然后在单元上设置近似的均布载荷。
def F(coords,NoN,PD,b,c):
G = np.zeros([NoN , PD])
U = np.zeros([NoN , PD])
# 施加左侧载荷
m=0
for i in range(NoN):
if int(coords[i,:][0]*10**2) == -(b/2)*10**2+1.:
m=m+1
p=c/m*689.5
k = 0
for i in range(NoN):
if int(coords[i][0]*10**2) == -(b/2)*10**2+1.:
print(i)
k = k + 1
G[i][0]= (-p)
# 施加右侧载荷
n=0
for i in range(NoN):
if int(coords[i,0:][0]*10**2) == (b/2)*10**2-1.:
n = n + 1
p = c/ n*689.5
for i in range(NoN):
if int(coords[i,0:][0]*10**2) == (b/2)*10**2-1.:
G[i][0] = p
return G
4.求解位移
U = np.zeros([2*NoN , PD-1])
T= np.linalg.inv(K)#T为K（总体刚度矩阵）逆矩阵
U=np.dot(T,GT)#U=[K]-1[G]
U=U.reshape(2*NoN,PD-1)#把位移转化为一维数组
5.计算应力应变
epsilonx = np.zeros(NoE)
epsilony = np.zeros(NoE)
epsilonxy = np.zeros(NoE)
A = np.zeros([3,3])
for i in range(NoE):
dot_id = Edof[i]
dot1 = np.array(coords[int(Edof[i, 0])])
dot2 = np.array(coords[int(Edof[i, 1])])
dot3 = np.array(coords[int(Edof[i, 2])])
beta1 = dot2[1] – dot3[1]
beta2 = dot3[1] – dot1[1]
beta3 = dot1[1] – dot2[1]
gamma1 = -dot2[0] + dot3[0]
gamma2 = -dot3[0] + dot1[0]
gamma3 = -dot1[0] + dot2[0]
B=np.array([[1, dot1[0], dot1[1]], [1, dot2[0], dot2[1]],[1, dot3[0], dot3[1]]]) / 2
A = np.linalg.det(B)
epsilonx[i – 1] = (beta1 * U[dot_id[0]] + beta2 * U[dot_id[1]] + beta3 * U[dot_id[2]]) / (2 * A)
epsilony[i – 1] = (gamma1 * U[NoN + dot_id[0]] +
gamma2 * U[NoN + dot_id[1]] + gamma3 * U[NoN + dot_id[2]]) / (2 * A)
epsilonxy[i – 1] = (beta1 * U[NoN + dot_id[0]] +
beta2 * U[NoN + dot_id[1]] + beta3 * U[NoN + dot_id[2]]
+ gamma1 * U[dot_id[0]] + gamma2 * U[dot_id[1]] + gamma3 * U[dot_id[2]]) / (2 *A)
sigmax = (E / (1 – v ** 2)) * (epsilonx + v * epsilony)
sigmay = (E / (1 – v ** 2)) * (epsilony + v * epsilonx)
sigmaxy = (E / (2 * (1 + v))) * epsilonxy
五、结果及分析
编程结果：
X方向位移（Abaqus）
y方向位移（Abaqus）
应力：
数值比较：最大x方向位移
|
编程 |
Abaqus |
|
1.79345518e-03 |
1.183e-03 |
数值比较：最大y方向位移
|
编程 |
Abaqus |
|
8.84718366e-04 |
1.838e-04 |
数值比较：最大x方向应力
|
编程 |
Abaqus |
|
2.82313695e+04 |
1.028e+04 |
分析：在用Abaqus求解时，将中间整个圆环边界位移设置为零，相当于减少了薄板的位移量，因此出现了编程位移数值大于Abaqus数值的问题。而在使用相同模型时，位移的大小决定了应变的大小，应力呈现与应变相同的变化趋势，所以编程的应力大于Abaqus的应力，但结果在同一个数量级，可以基本确定是以上原因。
作者：jackl的科研日常
