一阶四面体单元,共有4个节点,每个节点有ux, uy, uz 3个自由度,共有12个自由度。一阶四面体单元的位移函数u(x,y,z), v(x,y,z) 和w(x,y,z)均为线性函数,故其单元应变场和单元应力场皆为常量。一阶四面体单元的单元刚度矩阵Ke的公式推导可参考《有限元方法基础教程(第5版)》相关的章节。
一阶四面体单元的单元刚度矩阵求解的python代码:
#一阶四面体实体单元
#共4个节点,每个节点有且仅有 UX,UY,UZ三个自由度import numpy as np
from numpy import dot #矩阵乘法
from numpy.linalg import det # 求行列式
class Tetra3D(object):
def __init__(self,nodes, E, nu):
self.nodes = nodes
self.E = E
self.nu = nu
self.sixV = self.six_V()
self.B = self.calc_B()
self.D = self.calc_D()
self.Ke = self.calc_Ke()
def six_V(self): # 四面体体积的6倍
m = np.ones((4,4),dtype = np.float)
for i in range(4):
m[i,1:] = self.nodes[i]
return det(m) def calc_B(self):
R = (0,1,2,3) # 替代range(4)
#Alpha = [ (-1)**i * det(np.delete(self.nodes,i,axis=0)) for i in R] #用不上 a = np.ones((4,3))
for i in R:
a[i,1:] = self.nodes[i,(1,2)]
Belta = [ (-1)**(i+1) *det(np.delete(a,i,axis=0)) for i in R] a = np.ones((4,3))
for i in R:
a[i,1:] = self.nodes[i,(0,2)]
Gama = [ (-1)**i *det(np.delete(a,i,axis=0)) for i in R]
a = np.ones((4,3))
for i in R:
a[i,1:] = self.nodes[i,(0,1)]
Delta = [ (-1)**(i+1) *det(np.delete(a,i,axis=0)) for i in R]
B = np.array(((Belta[0], 0, 0, Belta[1], 0, 0, Belta[2], 0, 0, Belta[3], 0, 0),
( 0, Gama[0],0, 0, Gama[1],0, 0, Gama[2],0, 0, Gama[3],0),
( 0, 0, Delta[0],0, 0, Delta[1],0, 0,Delta[2],0,0,Delta[3]),
(Gama[0],Belta[0],0,Gama[1],Belta[1],0,Gama[2],Belta[2],0,Gama[3],Belta[3],0,),
(0,Delta[0],Gama[0],0,Delta[1],Gama[1],0,Delta[2],Gama[2],0,Delta[3],Gama[3]),
(Delta[0],0,Belta[0],Delta[1],0,Belta[1],Delta[2],0,Belta[2],Delta[3],0,Belta[3])),
dtype = np.float64)
B *= 1.0/self.sixV return B
def calc_D(self):
nu = self.nu
a = 1-nu
b = (1-2*nu)/2.0
D = np.array(((a, nu, nu, 0, 0, 0),
(nu,a, nu, 0, 0, 0),
(nu,nu, a, 0, 0, 0),
(0, 0, 0, b, 0, 0),
(0, 0, 0, 0, b, 0),
(0, 0, 0, 0, 0, b)),
dtype = np.float64)
D *= self.E/(1+nu)/(2*b)
return D def calc_Ke(self):
Ke = self.sixV/6.0 * dot(self.B.T, dot(self.D,self.B))
return Ke def calc_Strain(self,disp):#disp大小12x1,单元位移矩阵
return dot(self.B,disp)
def calc_Stress(self,disp):
Strain = self.calc_Strain(disp)
return dot(self.D,Strain)
if __name__ =="__main__":
# 节点坐标 x1,y1; x2,y2; x3,y3; x4,y4 ,须逆时针排 #import time
#t0 = time.time()
#t1 = time.time() #print(f"耗时{t1-t0}") nodes= np.array([[0.5, 0.5,0],[0.25, 0.25,0.4],[0, 0,0],[0.5, 0,0]])
#nodes= np.array([[0, 0,0],[0.5, 0,0],[0.5, 0.5,0],[0.25, 0.25,0.4]])
e = Tetra3D(nodes, E= 1.0E6, nu = 0.25)
print(e.sixV)
print(e.B[0,0])
下面开始求解一个简单的有限元模型。如下图,该模型共有8个节点,5个一阶四面体单元。左边的四个节点约束住全部自由度,右边的4个节点施加y向的外力。
有限元模型的创建和求解的python代码:
from numpy import array, mat,zeros, double, integer,float64,sqrt,ravel
from numpy.linalg import solve #,det
from random import random
from tetra3d import Tetra3D
#from readFromFemap_Tetra import readNeuFile
import time
time1 = time.time()
#单位体系 N,m, kg
#Nodes info.: x,y,z. #Node Id must start from 1,the step is 1.
#读入节点坐标和单元信息
#NODE, ELEM = readNeuFile("3d.neu")
#NODE = array(NODE,dtype=double)
#ELEM = array(ELEM,dtype=integer)
dofs_per_node = 3 #x,y,z
NODE= array([[0,0,0],
[0.025,0,0],
[0,0.5,0],
[0.025,0.5,0],
[0.,0.,0.25],
[0.025,0,0.25],
[0,0.5,0.25],
[0.025,0.5,0.25]],
dtype=float64)
#Node ID 从0开始,步长1为1自增
node_qty = NODE.shape[0] #节点数量
#MAT ID,TYPE ID,4 Nodes' ID
ELEM = array([[1,1,3,5,0,1],
[1,1,0,3,2,6],
[1,1,5,4,6,0],
[1,1,5,6,7,3],
[1,1,0,5,3,6]],
dtype=integer)
#单元ID 从0开始,步长1为1自增
elem_qty = ELEM.shape[0] #单元数量
#Boundary conditions
#节点自由度 dof CONSTRAINtrain
fixed_nodeID = [0,1,4,5]
CONSTRAIN = [(i,j,0.0) for i in fixed_nodeID for j in range(dofs_per_node)] #约束x,y,z
#外力 Force矩阵
F = mat(zeros((node_qty * dofs_per_node,1)), dtype=double)
F[dofs_per_node*2 + 1]= 3.125 # Fy on node 2,第1个自由度上(每节点3自由度时)
F[dofs_per_node*3 + 1]= 6.25 # Fy on node 3,第1个自由度上(每节点3自由度时)
F[dofs_per_node*7 + 1]= 3.125 # Fy on node 7,第1个自由度上(每节点3自由度时)
F[dofs_per_node*6 + 1]= 6.25 # Fy on node 6,第1个自由度上(每节点3自由度时)
print(F)
#材料属性
E= 210e6 # 弹性模量
nu = 0.30 #泊松比
time2 = time.time()
print(f"有限元模型输入完成!耗时{time2-time1}")
def assembleK(NODE,ELEM,CONSTRAIN):
elems = [] #用于存储单元对象
K= mat(zeros((node_qty*dofs_per_node ,node_qty *dofs_per_node)), dtype=float64)#初始化总刚度矩阵
#遍历单元
for ie in range(elem_qty):
nodesID = ELEM[ie,2:6] #6 =2+4, 单元的4个节点的ID:
i,j,m,n = nodesID
nodes = array([NODE[i],NODE[j],NODE[m],NODE[n]])
elem = Tetra3D(nodes,E, nu)#生成单元
elems.append(elem)
#总刚度矩阵组装(更新4x4个区域)
M = 3
for N1,I in enumerate(nodesID):
for N2,J in enumerate(nodesID):
K[dofs_per_node*I:dofs_per_node*I+M,dofs_per_node*J:dofs_per_node*J+M] += elem.Ke[M*N1:M*(N1+1), M*N2:M*(N2+1)]
#将边界条件(力,位移约束)更新到刚度矩阵;更新外力矩阵
BigNum = 1.0e150#大数法
#cr = CONSTRAIN.shape[0]
#a = (K-K.T)
#print(K==K.T)
#遍历约束节点
for nodeID,dofID,nodeDisp in CONSTRAIN:
jj = nodeID* dofs_per_node + dofID # 约束在总刚度矩阵的行号
jj = int(jj)#CONSTRAIN 是浮点型数组,所有jj的结果是浮点型。做索引需是整数
K[jj, jj] *= BigNum #总刚度矩阵对角线上的元素乘以大数
F[jj] = K[jj, jj]*nodeDisp #载荷也用更新后的K[jj,jj]乘以指定位移
return K, elems
# SOLVE
# 求解 位移矩阵Deformation Matrix DISP
#K,elems = assembleK(NODE,ELEM,CONSTRAIN)
K,elems= assembleK(NODE,ELEM,CONSTRAIN)
time3 = time.time()
print(f"总刚度矩阵组装完成!耗时{time3-time2}")
DISP= solve(K, F)#位移矩阵
#将位移在边界条件节点上的值用输入的约束值修正
for nodeID,dofID,nodeDisp in CONSTRAIN:
jj = nodeID* dofs_per_node + dofID
jj=int(jj)
DISP[jj] = nodeDisp
time4 = time.time()
print(f"位移矩阵求解完成!耗时{time4-time3}")
print("位移矩阵:")
print(DISP)
#提高可读性:现在每个节点Ux,Uy,Uz 显示在同一行
Delta = DISP.reshape((-1,3))
X= NODE[:,0]
Y= NODE[:,1]
Z= NODE[:,2]
Delta_X = array(Delta[:,0])
Delta_Y = array(Delta[:,1])
Delta_Z = array(Delta[:,2])
#节点总位移
DISP_total = sqrt(Delta_X**2 + Delta_Y**2,Delta_Z**2)
Stress_elems = zeros((elem_qty,6),dtype=float64)
#遍历单元,求各个单元上各个节点上的应力和应变
for ie in range(elem_qty):
i,j,m,n = ELEM[ie,2:2+4] #单元的4个节点的ID:
print(DISP[dofs_per_node*i+2,0])
DISPe = array([[DISP[dofs_per_node*i,0]],
[DISP[dofs_per_node*i+1,0]],
[DISP[dofs_per_node*i+2,0]],
[DISP[dofs_per_node*i+2,0]],
[DISP[dofs_per_node*j+1,0]],
[DISP[dofs_per_node*j+2,0]],
[DISP[dofs_per_node*m,0]],
[DISP[dofs_per_node*m+1,0]],
[DISP[dofs_per_node*m+2,0]],
[DISP[dofs_per_node*n,0]],
[DISP[dofs_per_node*n+1,0]],
[DISP[dofs_per_node*n+2,0]]])
elem = elems[ie] #当前计算的单元
#elem.calc_Strain(DISPe)
Stress = elem.calc_Stress(DISPe)
Stress_elems[ie] = ravel(Stress) #按C语言的顺序展平
求得位移矩阵(24x24 reshape为8x3)为:
y向最大变形为:
应力矩阵(5x6,5为单元个数)为:
最大的y向正应力为:
本文分享自 Python可视化编程机器学习OpenCV 微信公众号,前往查看
如有侵权,请联系 cloudcommunity@tencent.com 删除。
本文参与 腾讯云自媒体同步曝光计划 ,欢迎热爱写作的你一起参与!