有限元算法依据常见的有限元法教材,简单复现悬臂梁在重力作用下的形变(为了变形更明显,重力大小扩大了10倍),还没来得及写注释。【卧槽快跑,没注释!】
节点是随机函数撒的点,完全没有优化;
meshpy库中的Delauny优化算法计算得到三角单元;
pygame实现图形绘制,图形如下(文字是自己后来写上去的):文章来源:https://www.toymoban.com/news/detail-741883.html
文章来源地址https://www.toymoban.com/news/detail-741883.html
import numpy as np
import copy
import pygame,sys
from pygame.locals import *
import meshpy.triangle as triangle
import json
class Node:
def __init__(self):
self.id = -1
self.type = -1
self.coord = -1
class Material:
def __init__(self):
self.id = -1
self.E = -1
self.h = -1
self.miu = -1
self.mass = -1
self.thickness = -1
class Element:
def __init__(self):
self.id = -1
self.type = -1
self.nodes = []
self.material = -1
self.C = -1
self.S = -1
self.k = -1
self.runType = 0
self.A = -1
self.G = -1
def get_k(self,nodes):
material = self.material
if material == -1:
return
if self.type == 0:
E = material.E
miu = material.miu
else :
E = material.E/(1-material.miu**2)
miu = material.miu/(1-material.miu)
D = np.mat([[1,miu,0],[miu,1,0],[0,0,(1-miu)/2]])*E/(1-miu**2)
A = np.mat([[1,nodes[self.nodes[0]].coord[0],nodes[self.nodes[0]].coord[1]],
[1,nodes[self.nodes[1]].coord[0],nodes[self.nodes[1]].coord[1]],
[1,nodes[self.nodes[2]].coord[0],nodes[self.nodes[2]].coord[1]],])
A = 0.5*np.abs(np.linalg.det(A))
node_abc = [[],[],[]]
B = np.zeros((3,6))
for i in range(3):
if i == 0:
j = 1
m = 2
elif i == 1:
j = 2
m = 0
else:
j = 0
m = 1
a = nodes[self.nodes[j]].coord[0]*nodes[self.nodes[m]].coord[1] - nodes[self.nodes[m]].coord[0]*nodes[self.nodes[j]].coord[1]
b = nodes[self.nodes[j]].coord[1] - nodes[self.nodes[m]].coord[1]
c = -nodes[self.nodes[j]].coord[0] + nodes[self.nodes[m]].coord[0]
node_abc[i] = [a,b,c]
Bi = np.mat([[b,0],[0,c],[c,b]])/2/A
B[:,i*2:i*2+2] = Bi
self.k = B.T*D*B*material.h*A
self.G = np.mat([0,1,0,1,0,1]).T*-1/3*material.mass*material.h*A
return self.k
class Main:
def __init__(self):
self.nodes = []
self.elements = []
self.structureHeight = -1
self.structureWidth = -1
self.structureThickness = -1
self.material = Material()
self.material.id = 1
self.material.E = 100000000 #材料弹性模量
self.material.mass = 100 #材料密度
self.material.miu = 0.2 #材料泊松比
self.material.h = 0.2 #忘了是啥了
self.meshSize = -1
self.runType = 0
self.g = -9.806*10 #为了使形变更明显,将重力大小扩大10倍
self.X = -1
def run(self):
self.structureHeight = 0.5
self.structureWidth = 4
self.structureThickness = 1
self.meshSize = [3,3]
node_row = int(np.ceil(self.structureHeight/self.meshSize[1]))+1
node_col = int(np.ceil(self.structureWidth / self.meshSize[0]))+1
node_count = 0
np.set_printoptions(precision=2)
mesh_info = triangle.MeshInfo()
points = [[0, 0], [0, self.structureHeight], [self.structureWidth, self.structureHeight], [self.structureWidth, 0]]
random_points = np.random.rand(300, 2)
random_points = np.mat(random_points)
random_points[:, 0] = random_points[:, 0] * self.structureWidth
random_points[:, 1] = random_points[:, 1] * self.structureHeight
points.extend(random_points.tolist())
# 将节点存储到json,可供第三方软件读取节点信息
# with open('data_random12172.json', 'w') as json_file:
# json_file.write(json.dumps(points))
mesh_info.set_points(points)
facets = [[0,1],[1,2],[2,3],[3,0]]
mesh_info.set_facets(facets)
mesh = triangle.build(mesh_info)
node_count = 0
for i, p in enumerate(mesh.points):
node = Node()
node.id = node_count
node.type = 1
node.coord = p
self.nodes.append(node)
node_count = node_count + 1
element_count = 0
K = np.zeros((node_count*2,node_count*2))
F = np.zeros((2 * node_count, 1))
for i,element_item in enumerate(mesh.elements):
element = Element()
element.id = element_count
element.nodes = element_item
element.material = self.material
k = element.get_k(self.nodes)
for j in range(3):
F[2*element_item[j]:2*element_item[j]+2,:] = F[2*element_item[j]:2*element_item[j]+2,:]+element.G[2*j:2*j+2,:]*self.g
for m in range(3):
K[2*element_item[j]:2*element_item[j]+2,2*element_item[m]:2*element_item[m]+2] = \
K[2 * element_item[j]:2 * element_item[j] + 2, 2 * element_item[m]:2 * element_item[m] + 2] + k[2*j:2*j+2,2*m:2*m+2]
self.elements.append(element)
element_count = element_count + 1
for i in range(node_count):
if np.abs(self.nodes[i].coord[0] - 0) < 1e-5:
K[:,2*i:2*i+2] = 0
K[2 * i:2 * i + 2,:] = 0
K[2 * i:2 * i + 2,2 * i:2 * i + 2] = np.eye(2)
F[2*i:2*i+2,:] = 0
self.X = np.mat(K).I*np.mat(F)
a = 0
if __name__ == '__main__':
main = Main()
main.run()
pygame.init()
geo_scale = [300,300]
geo_offset_mixed = [55,50]
geo_offset_origin = [55,300]
geo_offset_onload = [55,550]
screen = pygame.display.set_mode((1400, 800))
screen.fill(pygame.Color(255, 255, 255))
line_color_orogin = (50, 50, 255)
line_color_onload = (50, 50, 255)
flag = 0
while True:
for event in pygame.event.get():
if event.type == QUIT:
pygame.quit()
sys.exit()
my_font = pygame.font.SysFont('arial', 16)
### 无荷载图形与带荷载图形重叠绘制
pygame.draw
for i in range(len(main.nodes)):
if flag == 1:
break
pointpos_origin = (np.multiply(np.mat(main.nodes[i].coord),np.mat(geo_scale))+np.mat(geo_offset_mixed)).tolist()[0]
pointpos_onload = (np.multiply(np.mat(main.nodes[i].coord+main.X[2*i:2*i+2,:].T),np.mat(geo_scale))+np.mat(geo_offset_mixed)).tolist()[0]
pygame.draw.circle(screen, [255, 0, 0], [int(pointpos_origin[0]), int(pointpos_origin[1])], 3, 2)
pygame.draw.circle(screen, [50, 200, 50], [int(pointpos_onload[0]), int(pointpos_onload[1])], 3, 2)
for i in range(len(main.elements)):
if flag == 1:
break
node1 = (np.multiply(np.mat(main.nodes[main.elements[i].nodes[0]].coord),np.mat(geo_scale))
+np.mat(geo_offset_mixed)).tolist()[0]
node2 = (np.multiply(np.mat(main.nodes[main.elements[i].nodes[1]].coord),np.mat(geo_scale))
+np.mat(geo_offset_mixed)).tolist()[0]
node3 = (np.multiply(np.mat(main.nodes[main.elements[i].nodes[2]].coord),np.mat(geo_scale))
+np.mat(geo_offset_mixed)).tolist()[0]
pygame.draw.lines(screen, line_color_orogin, False, [node1, node2], 1)
pygame.draw.lines(screen, line_color_orogin, False, [node1, node3], 1)
pygame.draw.lines(screen, line_color_orogin, False, [node2, node3], 1)
node1 = (np.multiply(np.mat(main.nodes[main.elements[i].nodes[0]].coord)+main.X[2*main.elements[i].nodes[0]:2*main.elements[i].nodes[0]+2,:].T, np.mat(geo_scale))
+ np.mat(geo_offset_mixed)).tolist()[0]
node2 = (np.multiply(np.mat(main.nodes[main.elements[i].nodes[1]].coord)+main.X[2*main.elements[i].nodes[1]:2*main.elements[i].nodes[1]+2,:].T, np.mat(geo_scale))
+ np.mat(geo_offset_mixed)).tolist()[0]
node3 = (np.multiply(np.mat(main.nodes[main.elements[i].nodes[2]].coord)+main.X[2*main.elements[i].nodes[2]:2*main.elements[i].nodes[2]+2,:].T, np.mat(geo_scale))
+ np.mat(geo_offset_mixed)).tolist()[0]
pygame.draw.lines(screen, line_color_orogin, False, [node1, node2], 1)
pygame.draw.lines(screen, line_color_orogin, False, [node1, node3], 1)
pygame.draw.lines(screen, line_color_orogin, False, [node2, node3], 1)
### 无荷载图形与带荷载图形分开绘制
for i in range(len(main.nodes)):
if flag == 1:
break
pointpos_origin = (np.multiply(np.mat(main.nodes[i].coord),np.mat(geo_scale))+np.mat(geo_offset_origin)).tolist()[0]
pointpos_onload = (np.multiply(np.mat(main.nodes[i].coord+main.X[2*i:2*i+2,:].T),np.mat(geo_scale))+np.mat(geo_offset_onload)).tolist()[0]
pygame.draw.circle(screen, [255, 0, 0], [int(pointpos_origin[0]), int(pointpos_origin[1])], 3, 2)
pygame.draw.circle(screen, [50, 200, 50], [int(pointpos_onload[0]), int(pointpos_onload[1])], 3, 2)
for i in range(len(main.elements)):
if flag == 1:
break
node1 = (np.multiply(np.mat(main.nodes[main.elements[i].nodes[0]].coord),np.mat(geo_scale))
+np.mat(geo_offset_origin)).tolist()[0]
node2 = (np.multiply(np.mat(main.nodes[main.elements[i].nodes[1]].coord),np.mat(geo_scale))
+np.mat(geo_offset_origin)).tolist()[0]
node3 = (np.multiply(np.mat(main.nodes[main.elements[i].nodes[2]].coord),np.mat(geo_scale))
+np.mat(geo_offset_origin)).tolist()[0]
pygame.draw.lines(screen, line_color_orogin, False, [node1, node2], 1)
pygame.draw.lines(screen, line_color_orogin, False, [node1, node3], 1)
pygame.draw.lines(screen, line_color_orogin, False, [node2, node3], 1)
node1 = (np.multiply(np.mat(main.nodes[main.elements[i].nodes[0]].coord)+main.X[2*main.elements[i].nodes[0]:2*main.elements[i].nodes[0]+2,:].T, np.mat(geo_scale))
+ np.mat(geo_offset_onload)).tolist()[0]
node2 = (np.multiply(np.mat(main.nodes[main.elements[i].nodes[1]].coord)+main.X[2*main.elements[i].nodes[1]:2*main.elements[i].nodes[1]+2,:].T, np.mat(geo_scale))
+ np.mat(geo_offset_onload)).tolist()[0]
node3 = (np.multiply(np.mat(main.nodes[main.elements[i].nodes[2]].coord)+main.X[2*main.elements[i].nodes[2]:2*main.elements[i].nodes[2]+2,:].T, np.mat(geo_scale))
+ np.mat(geo_offset_onload)).tolist()[0]
pygame.draw.lines(screen, line_color_orogin, False, [node1, node2], 1)
pygame.draw.lines(screen, line_color_orogin, False, [node1, node3], 1)
pygame.draw.lines(screen, line_color_orogin, False, [node2, node3], 1)
flag = 1
pygame.display.update()
到了这里,关于python实现简单的二维有限元计算的文章就介绍完了。如果您还想了解更多内容,请在右上角搜索TOY模板网以前的文章或继续浏览下面的相关文章,希望大家以后多多支持TOY模板网!