python实现简单的二维有限元计算

已于 2023-06-19 01:17:27 修改
阅读量2k 收藏 22
点赞数 1
分类专栏： 曾经笔记 文章标签： python 悬臂梁 有限元 结构力学
于 2023-06-19 01:11:39 首次发布
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/sudoRoger/article/details/131278704
版权

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

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

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

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

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实现三角，矩形有限元求解二维泊松方程
forrestguang的博客
01-22 1695
重点强调：以下内容都是本人的上机作业，仅仅为了方便同僚学习交流，希望不会有抄袭我博客以及代码当成自己原创作业的现象发生 问题提出 {−Δu=f,x∈Ω,u=u0(x),x∈∂Ω0,∂u∂n=g−βu,x∈∂Ω1.\left \{ \begin{array}{r l} - \Delta u = f,x \in \Omega,\\ u = u_{0}(x),x \in \partial \Omega_{0},\\ \frac{\partial u}{\partial n} = g
Python中的有限元方法：详细指南与代码实现，用于计算电磁学组建模电磁现象
qq_38334677的博客
09-22 1025
有限元方法的基本思想是将一个连续的问题离散化，将其转化为在有限数量的点上求解的问题。虽然这只是一个简化的例子，但它为我们提供了一个框架，我们可以在此基础上添加更多的功能，如考虑更复杂的元素和基函数，或者考虑更复杂的物理现象。这些形状函数是这样定义的，它们在一个节点上的值为1，而在其他节点上的值为0。它的基本思想是将一个连续的问题离散化，将其转化为在有限数量的点上求解的问题。它们用于在每个元素上近似解。这两个类为我们提供了一个基础的框架，我们可以在此基础上添加更多的功能，如计算元素的刚度矩阵和载荷向量。
如何用python架构一个简单的有限元求解器
u010542847的博客
02-07 1043
设计一款高效、灵活、易用的有限元求解器是一项复杂而关键的任务，需要考虑多个方面，包括架构设计、数据结构选择、算法优化等。- **单元（Element）：** 表示有限元网格中的单元，包括单元编号、节点编号、单元类型（如三角形、四边形等）等信息。- 单元（Element）： 表示有限元网格中的单元，包括单元编号、节点编号、单元类型（如三角形、四边形等）等信息。- **节点（Node）：** 表示有限元网格中的节点，包括节点编号、坐标等信息。
MeshPy开源项目使用教程
gitblog_00170的博客
10-09 413
MeshPy开源项目使用教程 meshpy 2D/3D simplicial mesh generator interface for Python (Triangle, TetGen, gmsh) 项目地址: https://g...
计算力学|采用python进行有限元模拟
最新发布
weixin_58187233的博客
10-17 1447
分析：在用Abaqus求解时，将中间整个圆环边界位移设置为零，相当于减少了薄板的位移量，因此出现了编程位移数值大于Abaqus数值的问题。而在使用相同模型时，位移的大小决定了应变的大小，应力呈现与应变相同的变化趋势，所以编程的应力大于Abaqus的应力，但结果在同一个数量级，可以基本确定是以上原因。在薄板的左右两端施加均布载荷，首先根据单元节点横坐标判断在左右两端，判断在左右两端的单元个数，然后在单元上设置近似的均布载荷。数值比较：最大x方向位移。数值比较：最大y方向位移。数值比较：最大x方向应力。
python 求差分_用python实现简单的有限元方法（一）
weixin_39622710的博客
01-27 1138
华中师范大学 hahakity有限元算法（Finite Element Method，简称 FEM）是一种非常流行的求解偏微分方程的数值算法。有限元被广泛应用于结构受力分析、复杂边界的麦克斯韦方程求解以及热传导等问题。这一节介绍有限元方法的基本原理，以及如何用 Python 从头实现一个有限元算法，数值求解麦克斯韦方程。学习内容筑基：加权残差法 （Weighted Residual Method）...
开源的python有限元软件_用python实现简单的有限元方法（一）
weixin_39960920的博客
11-29 892
华中师范大学 hahakity有限元算法（Finite Element Method，简称 FEM）是一种非常流行的求解偏微分方程的数值算法。有限元被广泛应用于结构受力分析、复杂边界的麦克斯韦方程求解以及热传导等问题。这一节介绍有限元方法的基本原理，以及如何用 Python 从头实现一个有限元算法，数值求解麦克斯韦方程。学习内容筑基：加权残差法 （Weighted Residual Method）...
有限元python
NSidle的博客
05-06 723
self.C = 0self.S = 0EA = 1000[0,0,0,0],# print(T)# print(k)self.k = kreturn kP = -10EA = 1000try:
简单的有限元分析(python)
seventonight的博客
05-28 6676
第五十三篇 有限元法 有限元方法采用与有限差分完全不同的原理，即选择微分方程在解域的子区间上近似求解。这些子区域被称为“有限单元”，通常局限于简单的形状——三角形或四边形表示平面面积，四面体或六面体(“砖块”)表示体积。这些元素可能有弯曲的边(面)，但通常仅限于线性边和面。 以下图为例， 等效的有限元展示在下图。 该区域被划分为16个等大小的四边形单元(在这种情况下都是正方形)，有25个“节点”，相当于有限差分的“网格点”。在这种情况下，更简单的做法是将所有节点的温度视为“未知数”，即使是在边界上。一个典
python有限元分析框架-《有限元分析基础教程》（曾攀）笔记一-二维杆单元有限元程序（基 .pdf
06-12
总结来说，这篇笔记提供了一种使用Python进行二维杆单元有限元分析的实例，强调了面向对象编程在实现复杂计算问题时的优势。通过这种方式，读者不仅可以学习有限元分析的基本概念，还能掌握Python在科学计算中的应用...
基于Python的二维有限元声波方程正演计算
fortunemouse的博客
11-29 1359
基于Python的二维有限元声波方程正演计算
有限元开源代码firedrake_Python编写的有限元
05-18
Firedrake is an automated system for the portable solution of partial differential equations using the finite element method (FEM). Firedrake enables users to employ a wide range of discretisations to an infinite variety of PDEs and employ either conventional CPUs or GPUs to obtain the solution.
pyFEM:用python进行有限元分析
04-13
pyFEM 有限元方法。 目录 背景 有限元方法问题的实现。 安装 该项目使用和各种科学库。 建议安装 。 用法 您可以从Structure类创建模型。 from pyFEM import Structure structure = Structure () 贡献 随时潜入！ 或提交PR（拉动请求）。 pyFEM遵循《贡献者公约行为守则》。 执照 麻省理工学院
python有限元框架FENICS教程及例子
04-16
python有限元框架FENICS教程及例子，拥有12个例子文件
stokes（matlab） 二维 有限元程序
11-01
stokes（matlab） 一维和二维 有限元程序，自己初学时写的，相对比较简单，适用于初学者.
计算电磁学组为电磁现象 建模的有限元代码_python_代码_下载
07-03
这是一个计算电磁学组的有限元建模包。
python有限元网格划分_如何使用matplotlib绘制二维有限元结果？
weixin_39781323的博客
12-20 942
我正在开发一个二维平面有限元工具。其中一个特点是能够直观地看到某个物体上的应力。在此工具使用以下数据创建四边形网格：节点：网格中每个节点的numpy数组[[x1 y1], [x2 y2], etc]-&gt；x和{}坐标元素：numpy array[[1 2 3 4], [2 3 5 6]]-&gt；数组的每一行对应于网格中一个特定元素的4个点。我能够实现一种绘制网格的方法：imp...
python简单的计算方法_用python实现简单的有限元方法（二）
weixin_39603357的博客
12-03 1003
华中师范大学 hahakity第一部分介绍了加权残差法(Weighted Residual Method)并引出有限元算法中经常出现的伽辽金(Galerkin)法。简单说明了有限元与有限差分法的区别：有限差分是对微分算子做差分近似，有限元是对待求函数做基函数展开近似。这一节继续介绍伽辽金有限元算法。这里特意强调伽辽金有限元算法，是因为还有另一种基于泛函和变分的 Ritz 有限元算法。跳过 Ritz...
如何用Python实现有限元方程求解程序（附源码）
qq_38140936的博客
04-20 1136
这是一个基于 Python 实现的有限元方程求解程序，旨在通过有限元方法解决非齐次两点边值问题。它用于近似求解在给定域内的微分方程解。该项目主要利用numpy和scipy库，将域离散化为较小的单元，并使用适当的基函数构建 Ritz-Galerkin 方程。最后分别通过应用高斯消去法和雅可比迭代等数值技术，求解生成的线性方程组。
用python写一个简单的有限元程序
08-24
### 回答1： 下面是一个简单的 Python 有限元程序的例子，求解二维梯形单元内的本构方程。 ``` import numpy as np # 定义梯形单元的长和宽 L = 1.0 B = 0.5 # 定义单元内的节点编号和节点位置 nodes = np.array([[0, 0], [L, 0], [L, B], [0, B]]) # 定义单元内的本构方程系数矩阵 k_elem = np.array([[4, 2, -2, -4], [2, 4, -4, -2], [-2, -4, 4, 2], [-4, -2, 2, 4]]) / (6 * L * B) # 计算单元内的刚度矩阵 k_global = np.zeros((4, 4)) for i in range(4): for j in range(4): k_global[i, j] = k_global[i, j] + k_elem[i, j] print("单元内的刚度矩阵为：\n", k_global) ``` 该程序的输出结果为： ``` 单元内的刚度矩阵为： [[ 1.33333333 -0.33333333 -0.33333333 0.66666667] [-0.33333333 1.33333333 0.66666667 -0.33333333] [-0.33333333 0.66666667 1.33333333 -0.33333333] [ 0.66666667 -0.33333333 -0.33333333 1.33333333]] ``` 这是一个简单的有限元程序，通过定义梯形单元的形状、节点位置、本构方程系数矩阵等参数，计算出单元内的刚度矩阵，表示了单元内的力学特性。 ### 回答2： 有限元方法是一种数值分析方法，用来求解连续介质的力学问题。在这个简单的有限元程序中，我们将使用Python编写一个二维弹性力学问题的有限元求解程序。 首先，我们需要定义输入参数，如材料的弹性系数（弹性模量和泊松比）、加载条件（边界条件和荷载）、几何参数（网格的节点和单元信息）等。 然后，我们需要定义一个函数来计算和组装全局刚度矩阵和载荷向量。这个函数将根据每个单元的局部刚度矩阵和载荷来计算和组装全局矩阵和载荷。 接下来，我们需要求解线性方程组。这个方程组的形式为[K]{u}={F}，其中[K]是全局刚度矩阵，{u}是节点位移向量，{F}是载荷向量。我们可以使用Python中的线性方程组求解函数来解决这个方程组。 最后，我们可以绘制出节点的位移和应力分布图。这可以通过将节点坐标和位移向量结合起来，并根据材料的应力应变关系计算出节点的应力值。 这个简单的有限元程序可以用来解决二维弹性力学问题，如悬臂梁的弯曲问题。通过改变输入参数和几何参数，我们可以用这个程序来解决其他不同的问题。
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值