题外话：因为微信的推荐机制变动，有可能大家不会第一时间看到我的文章，请大家给我的公众号标上⭐，以免错过好资源。
功能概览
支持荷载类型： 位移荷载 表面压强荷载 重力荷载 单元类型（共计支持33种单元类型）： T2D2、T3D2 B21、B23、B31、B33 CPS4、CPS8、CPS3、CPS6、CPS4R、CPS8R、CPS4I、CPS4Bar、CPS4SR CPE4、CPE8、CPE3、CPE6、CPE4R、CPE8R、CPE4I、CPSEBar、CPE4SR C3D4、C3D10、C3D8、C3D20、C3D8R、C3D20R、C3D8I、C3D8Bar、C3D8SR 整体刚度矩阵组装：稀疏矩阵存储格式，使用三元组（Triplet）形式，即行索引数组、列索引数组和非零值数组 线性方程组求解技术：提供直接求解法和pcg迭代求解器选项，面对较大模型时可以大大提升刚度方程的求解效率 边界条件：使用罚函数技术 节点应力磨平处理: 对于单元形状不均匀的结构模型，对于共节点处的应力采用体积/面积加权磨平处理 对于单元形状均匀的结构模型，对于共节点处的应力采用直接绕节点平均磨平处理 支持多种材料属性的赋予 场输出:U、RF、UR、RFM、S 前处理：解析 Abaqus-inp 文件： Node节点坐标 Element单元局部节点编号 Nset节点集 Elset单元集 Surface单元面信息 配置文件：使用 json 文件作为程序的配置文件 材料信息 截面属性 边界条件 求解器 绘图控制 文件输出： 结果文件（.dat），记录各节点位移、支反力等 vtk 文件，可导入 Paraview 中进行可视化显示，或者使用 Python 的 Pyvista 库进行可视化也都是可以的 log 文件，记录计算各个过程的具体时间明细以及一些模型信息。
JSON配置文件
普通单一材料
平面单元需要指定厚度参数，空间实体单元无所谓，杆单元需要指定截面面积，梁单元需要指定截面类型，将会在下个小节详细介绍。
"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应力对比：
应力精度与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 // 最大迭代次数
}
后处理
fields：场信息，需要指定要显示的场信息二维杆单元：U1, U2, Umag, RF1, RF2, RFmag, S11 三维杆单元：U1, U2, U3, Umag, RF1, RF2, RF3, RFmag, S11 二维梁单元：U1, U2, Umag, UR3, RF1, RF2, RFmag, RFM3 三维梁单元：U1, U2, U3, Umag, UR1, UR2, UR3, URmag, RF1, RF2, RF3, RFmag, RFM1, RFM2, RFM3, RFMmag 平面单元：U1, U2, Umag, RF1, RF2, RFmag, S11, S22, S12, Mises, model, deform 空间实体单元：U1, U2, U3, Umag, RF1, RF2, RF3, RFmag, S11, S22, S33, S12, S13, S23, Mises, model, deform average：平均化方式，可选：simple, volume，分别对应绕节点直接平均和体积/面积加权edgeColor：单元边界颜色faceColor：单元面颜色，当输出model, deform时，显示的单元面颜色alpha：单元面透明度colorMap：可选abaqus或者matlab支持的colormap，如jet、summer等discretize：若为true，则colormap分段显示；若为false，则colormap连续平滑显示discretizeNum：colormap分段数
"plot" : {
"fields" : ["Mises","Umag"],
"average" : "simple",
"edgeColor" : "#000080",
"faceColor" : "#77AC30",
"alpha" : 1,
"colorMap" : "abaqus",
"discretize" : true,
"discretizeNum" : 12
}
单元类型
共计支持33种单元类型，罗列如下：
杆系单元
T2D2：二维杆单元 T3D2：三维杆单元
梁单元
B21：二维Timoshenko梁单元（按照Abaqus帮助文档进行剪切刚度修正） B31：三维Timoshenko梁单元（按照Abaqus帮助文档进行剪切刚度修正） B23：二维Euler-Bernoulli 梁单元(cubic beam) B33：三维Euler-Bernoulli 梁单元(cubic beam)
平面单元
平面应力：
CPS3：线性三角形单元 CPS6：二次三角形单元 CPS4：线性四边形单元 CPS8：二次四边形单元 CPS8R：减缩积分单元 CPS4R：减缩积分单元，已添加沙漏刚度 CPS4I：非协调线性四边形单元，避免剪切自锁 CPS4SR：选择性减缩积分，避免体积自锁 CPS4Bar：对B矩阵加以Bar修正，避免体积自锁（与Abaqus帮助文档处理方法一致）
平面应变：
CPE3：线性三角形单元 CPE6：二次三角形单元 CPE4：线性四边形单元 CPE8：二次四边形单元 CPE8R：减缩积分单元 CPE4R：减缩积分单元，已添加沙漏刚度 CPE4I：非协调线性四边形单元，避免剪切自锁 CPE4SR：选择性减缩积分，避免体积自锁 CPE4Bar：对B矩阵加以Bar修正，避免体积自锁（与Abaqus帮助文档处理方法一致）
空间实体单元
C3D4：线性四面体单元 C3D8：线性六面体单元 C3D10：二次四面体单元 C3D20：二次六面体单元 C3D20：减缩积分单元 C3D8R：减缩积分单元，已添加沙漏刚度 C3D8I：非协调线性八面体单元，避免剪切自锁 C3D8SR：选择性减缩积分，避免体积自锁 C3D8Bar：对B矩阵加以Bar修正，避免体积自锁（与Abaqus帮助文档处理方法一致）
大模型测试
105万个C3D4单元，网格信息及边界条件设置：
计算结果对比：
计算过程信息：
********************* 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
由于添加的单元类型较多，精度验证需要占据大量篇幅，没有贴在本次推文中，不过大可放心，每种单元类型在添加时均和Abaqus做了验证分析，均验证通过！感兴趣的小伙伴可逐个验证。
后面会出推文介绍一些比较感兴趣的单元详解，比如三维Timoshenko梁单元是如何根据Abaqus帮助文档进行剪切刚度修正，也或是C3D8R是如何进行沙漏刚度施加，亦或是C3D8单元如何进行B矩阵的Bar修正等等等等，对程序内部的细节内容感兴趣的小伙伴也可以在留言区进行留言，木木会抽时间做详细介绍。
有关程序更新的功能就介绍到这里，感谢您的阅读。整套程序已发布在知识星球中，后台回复：星球，即可加入，对源程序的疑问，可在星球内详细讨论。
本次推文已上传至Web端，如果觉得公众号界面阅读不便，可点击下方“阅读原文”，即可跳转至木木的个人网页：
参与更多互动交流，快快在下方留言区留下你的小脚印吧~
-End-
易木木响叮当
想陪你一起度过短暂且漫长的科研生活
