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

已于 2023-06-19 01:17:27 修改
阅读量1.4k 收藏 19
点赞数 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()

成田日上
关注
  • 1
    点赞
  • 踩
  • 19
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
专栏目录
poisson:使用有限元方法求解任意二维域上的泊松方程
06-12
任意二维域上的泊松 使用有限元方法求解具有 RHS f 和 Dirichlet 边界条件的任意二维域上的泊松方程。
有限元法计算二维圆柱绕流问题——Python代码实现
weixin_52046577的博客
05-11 1186
选取流函数Ψ为变量，对拉普拉斯方程进行求解（右边界为自然边界条件，其余边界为本质边界条件）；网格数据文件的生成暂时不在本文中详述。
python 求差分_用python实现简单的有限元方法（一）
weixin_39622710的博客
01-27 1041
华中师范大学 hahakity有限元算法（Finite Element Method，简称 FEM）是一种非常流行的求解偏微分方程的数值算法。有限元被广泛应用于结构受力分析、复杂边界的麦克斯韦方程求解以及热传导等问题。这一节介绍有限元方法的基本原理，以及如何用 Python 从头实现一个有限元算法，数值求解麦克斯韦方程。学习内容筑基：加权残差法 （Weighted Residual Method）...
python+feon有限元分析|单元类型
最新发布
T20151470的博客
03-18 382
2个节点，每个节点一个自由度（），局部坐标系的刚度矩阵和整体坐标系的刚度矩阵一致。假设弹簧单元的刚度为，则一次单元在局部坐标系和整体坐标系的刚度矩阵（2×2阶）为为局部坐标系中的单元力列阵为整体坐标系中的单元节点位移列阵。
meshpy:适用于Python的2D3D简单网格生成器接口（Triangle，TetGen，gmsh）
03-08
MeshPy：从Python生成简单网格 MeshPy为Python提供了高质量的三角形和四面体网格生成。 这种类型的网格主要用于有限元仿真代码中，但也有许多其他应用，从计算机图形学到机器人技术。 为了产生二维和三维网格，MeshPy提供的Python接口到三个广泛关注的网格生成器，由J. Shewchuk， 由杭司前二者都包含在包中稍微修改的版本。 包括用于后者的通用网格读取器，以及从Python脚本运行gmsh的简便方法。 有关Christophe 和Jean-Francois Remacle的的界面，请参见gmsh_interop < >。 除了C ++编译器和有效的Python安装，MeshPy没有任何依赖关系。 在安装meshpy之前，您可以使用以下命令安装pybind11： pip install pybind11 从0.91.2版开始，MeshPy还可以与Pyth
pyFEM:用python进行有限元分析
04-13
pyFEM 有限元方法。 目录 背景 有限元方法问题的实现。 安装 该项目使用和各种科学库。 建议安装 。 用法 您可以从Structure类创建模型。 from pyFEM import Structure structure = Structure () 贡献 随时潜入！ 或提交PR（拉动请求）。 pyFEM遵循《贡献者公约行为守则》。 执照 麻省理工学院
python有限元分析框架-《有限元分析基础教程》（曾攀）笔记一-二维杆单元有限元程序（基 .pdf
06-12
python有限元分析框架_《有限元分析基础教程》（曾攀）笔 记⼀-⼆维杆单元有限元程序（基。。。 曾攀⽼师的《有限元分析基础教程》第三章有⼆维杆单元的推导，并结合⼀个例题进⾏了解析解和基于Matlab的程序求解。但是我感觉书 中的MATLAB代码有点罗嗦，⽽且⼀些实现⽅法也⽐较⿇烦，⽐如已经知道了杆单元的起点和终点坐标，仍然需要另外给出单元局部坐标 与整体坐标的夹⾓，这完全没必要。于是我就⽤Python重构了这段程序，当然并不是把书中的MATLAB代码翻译成python(事实上完全可 以这么⼲，⽽且很快！)。⽐如我使⽤了⾯向对象的思想，把杆单元写成了⼀个类，这样思路⽐较清晰。 #! /usr/bin/python#-*- coding: utf-8 -*- importmathimportnumpy as np sqrt, cos, sin, pi=math.sqrt, math.cos, math.sin, math.pi"前处理"nodeNumber= 4KK= np.zeros((2*nodeNumber, 2*nodeNumber))"""边界条件，U表⽰节点的位移向量，如果
【python 数字孪生】之使用python进行有限元可视化
weixin_47510917的博客
07-24 1829
通过python来实现有限元计算或者有限元结果可视化
有限元python
NSidle的博客
05-06 562
self.C = 0self.S = 0EA = 1000[0,0,0,0],# print(T)# print(k)self.k = kreturn kP = -10EA = 1000try:
Python建立简单3维有限元框架
m0_65134578的博客
04-05 409
Python建立简单3维有限元框架
Python中的有限元方法：详细指南与代码实现，用于计算电磁学组建模电磁现象
qq_38334677的博客
09-22 681
有限元方法的基本思想是将一个连续的问题离散化，将其转化为在有限数量的点上求解的问题。虽然这只是一个简化的例子，但它为我们提供了一个框架，我们可以在此基础上添加更多的功能，如考虑更复杂的元素和基函数，或者考虑更复杂的物理现象。这些形状函数是这样定义的，它们在一个节点上的值为1，而在其他节点上的值为0。它的基本思想是将一个连续的问题离散化，将其转化为在有限数量的点上求解的问题。它们用于在每个元素上近似解。这两个类为我们提供了一个基础的框架，我们可以在此基础上添加更多的功能，如计算元素的刚度矩阵和载荷向量。
有限元求解二维调和方程——python
07-07
有限元求解二维调和方程——python
FEMstudies:二维线性有限元方法研究
05-24
有限元研究 2D线性有限元方法研究。
python简单的计算方法_用python实现简单的有限元方法（二）
weixin_39603357的博客
12-03 897
华中师范大学 hahakity第一部分介绍了加权残差法(Weighted Residual Method)并引出有限元算法中经常出现的伽辽金(Galerkin)法。简单说明了有限元与有限差分法的区别：有限差分是对微分算子做差分近似，有限元是对待求函数做基函数展开近似。这一节继续介绍伽辽金有限元算法。这里特意强调伽辽金有限元算法，是因为还有另一种基于泛函和变分的 Ritz 有限元算法。跳过 Ritz...
简单的有限元分析(python)
seventonight的博客
05-28 6267
第五十三篇 有限元法 有限元方法采用与有限差分完全不同的原理，即选择微分方程在解域的子区间上近似求解。这些子区域被称为“有限单元”，通常局限于简单的形状——三角形或四边形表示平面面积，四面体或六面体(“砖块”)表示体积。这些元素可能有弯曲的边(面)，但通常仅限于线性边和面。 以下图为例， 等效的有限元展示在下图。 该区域被划分为16个等大小的四边形单元(在这种情况下都是正方形)，有25个“节点”，相当于有限差分的“网格点”。在这种情况下，更简单的做法是将所有节点的温度视为“未知数”，即使是在边界上。一个典
python+feon有限元分析|求解实例
T20151470的博客
11-27 590
Quad2D11S ---- 平面应力四边形单元。Quad2D11 ---- 平面应变四边形单元。Plate3D11 ---Midline 板单元。Tri2d11S---- 平面应力三角形单元。Tri2D11 ---- 平面应变三角形单元。Tetra3D11 ---- 四面体单元。Spring1D11 - 一维弹簧单元。Spring2D11 - 二维弹簧单元。Spring3D11 - 三维弹簧单元。Brick3D11 --- 六面体单元。Link2D11 - 二维弹簧单元。
有限元-二维有限元编程（矩形区域、三角剖分）
栀彧
07-13 4274
有限元、矩形区域、三角剖分
利用python进行有限元分析（一）
weixin_51788042的博客
01-18 470
利用python进行有限元分析
『OPEN3D』1.2 mesh处理 python篇
NNNNNathan的博客
02-16 3448
介绍open3d对mesh的基本处理操作
用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、付费专栏及课程。

余额充值