python实现简单的二维有限元计算-Toy模板网

这篇具有很好参考价值的文章主要介绍了python实现简单的二维有限元计算。希望对大家有所帮助。如果存在错误或未考虑完全的地方，请大家不吝赐教，您也可以点击"举报违法"按钮提交疑问。

有限元算法依据常见的有限元法教材，简单复现悬臂梁在重力作用下的形变（为了变形更明显，重力大小扩大了10倍），还没来得及写注释。【卧槽快跑，没注释！】

节点是随机函数撒的点，完全没有优化；

meshpy库中的Delauny优化算法计算得到三角单元；

pygame实现图形绘制，图形如下（文字是自己后来写上去的）：

python有限元编程,曾经笔记,python,悬臂梁,有限元,结构力学文章来源地址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模板网！