MFEAOOP V1.3

功能概览

  1. 支持荷载类型：
    1. 位移荷载
    2. 表面压强荷载
    3. 重力荷载
  2. 单元类型（共计支持33种单元类型）：
    1. 杆系单元：
      1. T2D2、T3D2
      2. B21、B23、B31、B33
    2. 平面单元：
      1. CPS4、CPS8、CPS3、CPS6、CPS4R、CPS8R、CPS4I、CPS4Bar、CPS4SR
      2. CPE4、CPE8、CPE3、CPE6、CPE4R、CPE8R、CPE4I、CPSEBar、CPE4SR
    3. 空间实体单元：
      1. C3D4、C3D10、C3D8、C3D20、C3D8R、C3D20R、C3D8I、C3D8Bar、C3D8SR
  3. 整体刚度矩阵组装：稀疏矩阵存储格式，使用三元组（Triplet）形式，即行索引数组、列索引数组和非零值数组
  4. 线性方程组求解技术：提供直接求解法和pcg迭代求解器选项，面对较大模型时可以大大提升刚度方程的求解效率
  5. 边界条件：使用罚函数技术
  6. 节点应力磨平处理:
    1. 对于单元形状不均匀的结构模型，对于共节点处的应力采用体积/面积加权磨平处理
    2. 对于单元形状均匀的结构模型，对于共节点处的应力采用直接绕节点平均磨平处理
  7. 支持多种材料属性的赋予
  8. 场输出:U、RF、UR、RFM、S
  9. 前处理：解析 Abaqus-inp 文件：
    1. Node节点坐标
    2. Element单元局部节点编号
    3. Nset节点集
    4. Elset单元集
    5. Surface单元面信息
  10. 配置文件：使用 json 文件作为程序的配置文件
    1. 材料信息
    2. 截面属性
    3. 边界条件
    4. 求解器
    5. 绘图控制
  11. 文件输出：
    1. 结果文件（.dat），记录各节点位移、支反力等
    2. vtk 文件，可导入 Paraview 中进行可视化显示，或者使用 Python 的 Pyvista 库进行可视化也都是可以的
    3. log 文件，记录计算各个过程的具体时间明细以及一些模型信息。

JSON配置文件

普通单一材料

平面单元需要指定厚度参数tt，空间实体单元无所谓，杆单元需要指定截面面积AA，梁单元需要指定截面类型，将会在下个小节详细介绍。

"materials" : [
        {
        "name" : "Material-1",
        "category" : "Elastic",
        "data" : {
            "E" : 1E6, // 弹性模量
            "v" : 0.3, // 泊松比
            "t" : 1, // 平面单元厚度
            "A" : 1 // 杆单元截面面积
            }
        }
    ]

梁单元截面

梁单元支持矩形截面、环形截面、圆形截面，需要指定几何参数。

"materials" : [
        {
        "name" : "Material-1",
        "category" : "Elastic",
        "data" : {
            "E" : 3e10, // 弹性模量
            "v" : 0.2, // 泊松比
            "Profile" : "Rectangular",  // 矩形截面
            "a" : 0.4, //宽度
            "b" : 0.4, // 高度
            "Profile" : "Pipe", //环形截面
            "D" : 0.2, //外半径
            "d" : 0.1, // 内半径
            "Profile" : "Circular", //圆形截面
            "r" : 0.1 // 半径
            }
        }
    ]

赋予多种材料属性

对于多种材料属性的施加，需要指定材料属性对应的单元集。

"materials" : [
    {
        "name" : "Material-1",
        "category" : "Elastic",
        "data" : {
            "E" : 1E6,
            "v" : 0.3,
            "t" : 1
        },
        "elSet" : "ZZ-Blank"
    },
    {
        "name" : "Material-2",
        "category" : "Elastic",
        "data" : {
            "E" : 1E7,
            "v" : 0.3,
            "t" : 1
        },
        "elSet" : "ZZA-Layer1"
    },
    {
        "name" : "Material-3",
        "category" : "Elastic",
        "data" : {
            "E" : 2E7,
            "v" : 0.3,
            "t" : 1
        },
        "elSet" : "ZZA-Layer2"
    }
],

以一个二维复合材料案例为例，前处理借助公众号：星辰北极星出品的POLARIS_MesoConcrete插件建立二维细观骨料模型，边界条件及网格划分如下图所示，采用CPS3单元：

二维细观骨料模型

Mises应力对比：

MFEAOOP

Abaqus

应力精度与Abaqus保持一致！

固支边界

type为constraint，需要指定不同的方向，nset为节点集（与inp文件一致）。比如想要在x方向和y方向均施加固定约束，可指定"direction": “xy”。

"bc" : [
    {
        "name" : "fix1",
        "type": "constraint",
        "direction": "xy",
        "nset" : "Set-XY"
    }]

位移边界

type为displacement，需要指定不同的方向，nset为节点集（与inp文件一致），value为位移值。比如想要在x方向和y方向均施加位移，可指定"direction": “xy”。

"bc" : [
    {
        "name" : "dis1", 
        "type": "displacement", 
        "direction": "x", 
        "nset" : "Set-U", 
        "value" : 1 
    }]

表面压强边界

type为pressure，需要指定表面名字（与inp文件一致），以及压强值。（不支持杆单元、三维梁单元）

"bc" : [
    {
        "name" : "pressure1", 
        "type": "pressure", 
        "surface" : "Surf-1", 
        "value" : 1000
    }]

重力荷载

type为gravity，需要指定重力加速度，以及重力方向。(不支持杆单元)

"bc" : [
    {
        "name" : "Gravity1", 
        "type": "gravity", 
        "density" : 9.8, 
        "gravity" : {
            "x": 0.0,
            "y": 0.0,
            "z": -9.8
        } 
    }]

求解器

"solve" : {
        "calculateStress": true, // 应力计算开关
        "solver" : "direct", // 求解器类型，可选：direct、pcg
        "tolerance": 1e-8, // 容差
        "maxIterations": 1000 // 最大迭代次数
    }

后处理

  1. fields：场信息，需要指定要显示的场信息
    1. 二维杆单元：U1, U2, Umag, RF1, RF2, RFmag, S11
    2. 三维杆单元：U1, U2, U3, Umag, RF1, RF2, RF3, RFmag, S11
    3. 二维梁单元：U1, U2, Umag, UR3, RF1, RF2, RFmag, RFM3
    4. 三维梁单元：U1, U2, U3, Umag, UR1, UR2, UR3, URmag, RF1, RF2, RF3, RFmag, RFM1, RFM2, RFM3, RFMmag
    5. 平面单元：U1, U2, Umag, RF1, RF2, RFmag, S11, S22, S12, Mises, model, deform
    6. 空间实体单元：U1, U2, U3, Umag, RF1, RF2, RF3, RFmag, S11, S22, S33, S12, S13, S23, Mises, model, deform
  2. average：平均化方式，可选：simple, volume，分别对应绕节点直接平均和体积/面积加权
  3. edgeColor：单元边界颜色
  4. faceColor：单元面颜色，当输出model, deform时，显示的单元面颜色
  5. alpha：单元面透明度
  6. colorMap：可选abaqus或者matlab支持的colormap，如jet、summer等
  7. discretize：若为true，则colormap分段显示；若为false，则colormap连续平滑显示
  8. discretizeNum：colormap分段数
"plot" : {
        "fields" : ["Mises","Umag"],
        "average" : "simple",
        "edgeColor" : "#000080",
        "faceColor" : "#77AC30",
        "alpha" : 1,
        "colorMap" : "abaqus",
        "discretize" : true,
        "discretizeNum" : 12
    }

单元类型

共计支持33种单元类型，罗列如下：

杆系单元

  1. T2D2：二维杆单元
  2. T3D2：三维杆单元

梁单元

  1. B21：二维Timoshenko梁单元（按照Abaqus帮助文档进行剪切刚度修正）
  2. B31：三维Timoshenko梁单元（按照Abaqus帮助文档进行剪切刚度修正）
  3. B23：二维Euler-Bernoulli 梁单元(cubic beam)
  4. B33：三维Euler-Bernoulli 梁单元(cubic beam)

平面单元

平面应力：

  1. CPS3：线性三角形单元
  2. CPS6：二次三角形单元
  3. CPS4：线性四边形单元
  4. CPS8：二次四边形单元
  5. CPS8R：减缩积分单元
  6. CPS4R：减缩积分单元，已添加沙漏刚度
  7. CPS4I：非协调线性四边形单元，避免剪切自锁
  8. CPS4SR：选择性减缩积分，避免体积自锁
  9. CPS4Bar：对B矩阵加以Bar修正，避免体积自锁（与Abaqus帮助文档处理方法一致）

平面应变：

  1. CPE3：线性三角形单元
  2. CPE6：二次三角形单元
  3. CPE4：线性四边形单元
  4. CPE8：二次四边形单元
  5. CPE8R：减缩积分单元
  6. CPE4R：减缩积分单元，已添加沙漏刚度
  7. CPE4I：非协调线性四边形单元，避免剪切自锁
  8. CPE4SR：选择性减缩积分，避免体积自锁
  9. CPE4Bar：对B矩阵加以Bar修正，避免体积自锁（与Abaqus帮助文档处理方法一致）

空间实体单元

  1. C3D4：线性四面体单元
  2. C3D8：线性六面体单元
  3. C3D10：二次四面体单元
  4. C3D20：二次六面体单元
  5. C3D20：减缩积分单元
  6. C3D8R：减缩积分单元，已添加沙漏刚度
  7. C3D8I：非协调线性八面体单元，避免剪切自锁
  8. C3D8SR：选择性减缩积分，避免体积自锁
  9. C3D8Bar：对B矩阵加以Bar修正，避免体积自锁（与Abaqus帮助文档处理方法一致）

大模型测试

105万个C3D4单元，网格信息及边界条件设置：

模型设置

计算结果对比：

MFEAOOP位移云图

Abaqus位移云图

计算过程信息：

*********************  MFEAOOP Version 1.3  ********************
About: 基于Matlab平台的小型有限元通用分析程序，采用面向对象的编程风格
Theory: 《有限元基础编程百科全书》
Author: 易公子
欢迎关注微信公众号/B站: 易木木响叮当
Email: yimumumfea@163.com
****************************************************************
Submitted: 2024-10-11 00:20:59
Analysis Input File                                 15.1943 s
----Number of nodes:                                197340
----Number of elements:                             1051568
----Element type:                                   C3D4
----Number of dofs:                                 592020
----Constraint:                                     0.2698 s
----Pressure:                                       0.0696 s
----Pressure:                                       0.0598 s
----Pressure:                                       0.0581 s
----Pressure:                                       0.0502 s
Assembly global stifness matrix:                    165.0002 s
Using PCG solver for large systems...
Solving Systems of Linear Equations:                106.6774 s
Stress calculation skipped based on JSON settings!
Write result File                                        1.7861 s
Write vtk File                                           4.3930 s
Plot                                                     2.2932 s
Completed: 2024-10-11 00:25:57                      297.9422 s
****************************************************************

项目文件结构图

V1-3
├─ class
│  ├─ FEAData.m
│  ├─ FEASolver.m
│  └─ writeFile.m
├─ elements
│  ├─ 2D
│  │  ├─ CPE3.m
│  │  ├─ CPE4.m
│  │  ├─ CPE4Bar.m
│  │  ├─ CPE4I.m
│  │  ├─ CPE4R.m
│  │  ├─ CPE4SR.m
│  │  ├─ CPE6.m
│  │  ├─ CPE8.m
│  │  ├─ CPE8R.m
│  │  ├─ CPS3.m
│  │  ├─ CPS4.m
│  │  ├─ CPS4Bar.m
│  │  ├─ CPS4I.m
│  │  ├─ CPS4R.m
│  │  ├─ CPS4SR.m
│  │  ├─ CPS6.m
│  │  ├─ CPS8.m
│  │  └─ CPS8R.m
│  ├─ 3D
│  │  ├─ C3D10.m
│  │  ├─ C3D20.m
│  │  ├─ C3D20R.m
│  │  ├─ C3D4.m
│  │  ├─ C3D8.m
│  │  ├─ C3D8Bar.m
│  │  ├─ C3D8I.m
│  │  ├─ C3D8R.m
│  │  └─ C3D8SR.m
│  ├─ beam
│  │  ├─ B21.m
│  │  ├─ B23.m
│  │  ├─ B31.m
│  │  └─ B33.m
│  └─ truss
│     ├─ T2D2.m
│     └─ T3D2.m
├─ input
│  ├─ Abaqus
│  │  ├─ aiFeiEr.inp
│  │  ├─ B21.inp
│  │  ├─ B23.inp
│  │  ├─ B31.inp
│  │  ├─ B33.inp
│  │  ├─ C3D10.inp
│  │  ├─ C3D20.inp
│  │  ├─ C3D4.inp
│  │  ├─ C3D4test.inp
│  │  ├─ C3D8.inp
│  │  ├─ CPS3.inp
│  │  ├─ CPS4.inp
│  │  ├─ CPS6.inp
│  │  ├─ CPS8.inp
│  │  ├─ hg4.inp
│  │  ├─ hg8.inp
│  │  ├─ T1D3.inp
│  │  ├─ T2D2.inp
│  │  ├─ T3D2.inp
│  │  └─ T3D2test.inp
│  ├─ beam2D.json
│  ├─ beam3D.json
│  ├─ plane.json
│  ├─ solid.json
│  ├─ truss2D.json
│  └─ truss3D.json
├─ main.m
├─ output
├─ README.md
└─ utitls
   ├─ Bmatrix2D.m
   ├─ Bmatrix3D.m
   ├─ Dmatrix2D.m
   ├─ Dmatrix3D.m
   ├─ gauss1D.m
   ├─ gauss2D.m
   ├─ gauss3D.m
   ├─ hammer2D.m
   ├─ hammer3D.m
   ├─ Jacobian.m
   ├─ shapeFunction1D.m
   ├─ shapeFunction2D.m
   ├─ shapeFunction3D.m
   └─ writeLog.m
有限元
#MFEA
MFEAOOP V1.3
http://yimumu.top/2024/10/08/MFEAOOP-V1.3/
作者
易公子
发布于
2024年10月8日
更新于
2024年10月11日
许可协议