题外话：因为微信的推荐机制变动，有可能大家不会第一时间看到我的文章，请大家给我的公众号标上⭐，以免错过好资源。
本期推文分享MFEAOOP V1.2的主要更新内容，对于一些功能细节可在以后推文中逐一介绍。
功能概览
支持荷载类型： 位移荷载 表面压强荷载(默认荷载方向垂直与单元表面，压为正，拉为负) 重力荷载(默认荷载方向竖直向下，即z轴负方向) 单元类型：C3D4、C3D10、C3D8、C3D20，后续的更多三维单元类型将逐步添加 内核求解：C3D20 单元增加了 14 点积分求解技术，保持精度的前提下，减少原有的高斯积分数量 稀疏矩阵存储格式：使用三元组（Triplet）形式，即行索引数组、列索引数组和非零值数组，分别对应矩阵中非零元素的位置和它们的值，大大提升了整体刚度矩阵的组装效率 线性方程组求解技术：提供直接求解法和pcg迭代求解器选项，面对较大模型时可以大大提升刚度方程的求解效率 边界条件：使用罚函数技术计算出节点位移后，一键解出节点的支反力 节点应力磨平处理: 对于单元形状不均匀的结构模型，对于共节点处的应力采用体积加权磨平处理 对于单元形状均匀的结构模型，对于共节点处的应力采用直接绕节点平均磨平处理 场输出: 节点位移U: Umag、U1、U2、U3 节点支反力RF: RFmag、RF1、RF2、RF3 节点应力S: S11、S22、S33、S12、S13、S23、Mises 前处理：支持一键解析 Abaqus-inp 文件 Node节点坐标 Element单元局部节点编号 Nset节点集 Elset单元集 Surface单元面信息 配置文件：使用 json 文件作为程序的配置文件，参数化控制: 材料信息 边界条件 求解器 绘图控制 文件输出： 结果文件（.dat），记录各节点的坐标、位移、支反力的值 vtk 文件，可导入 Paraview 中进行可视化显示，或者使用 Python 的 Pyvista 库进行可视化也都是可以的 log 文件，记录计算各个过程的具体时间明细以及一些模型信息。
JSON配置文件
现将配置文件的控制选项罗列如下:
{
"name": "C3D4", //Job name
"scale": 1, // Scale factor for the model
"inpFileName" : "./input/Abaqus/C3D4.inp", // Abaqus input file location
"mesh" : {
"etype" : "C3D4" // Element type
},
"materials" : [
{
"name" : "Material-1", // Material name
"category" : "Elastic", // Material category
"data" : {
"E" : 1E6, // Young's modulus
"v" : 0.3 // Poisson's ratio
}
}
],
"bc" : [
{
"name" : "fix1", // Boundary condition name
"type": "constraint", // Boundary condition type
"direction": "xyz", // Boundary condition direction
"nset" : "Set-FIX" // Node set name
},
{
"name" : "dis1", // Boundary condition name
"type": "displacement", // Boundary condition type
"direction": "xyz", // Boundary condition direction
"nset" : "Set-FIX", // Node set name
"value" : 2 // Displacement value
},
{
"name" : "pressure1", // Boundary condition name
"type": "pressure", // Boundary condition type
"surface" : "Surf-1", // Surface name
"value" : 100 // Pressure value
},
{
"name" : "Gravity1", // Boundary condition name
"type": "gravity", // Boundary condition type
"density" : 9.8, // Density value
"gravity" : 9.8 // Gravity value
}
],
"solve" : {
"calculateStress": false, // Whether to calculate stress
"solver" : "pcg", // Solver type(direct or pcg)
"tolerance": 1e-8, // Tolerance for the pcg solver
"maxIterations": 1000 // Maximum number of iterations for the pcg solver
},
"plot" : {
"fields" : ["Umag"], // Fields to plot(Umag U1 U2 U3 RFmag RF1 RF2 RF3 S11 S22 S33 S12 S13 S23 Mises)
"average" : "volume", // Averaging method(volume or simple)
"edgeColor" : "#000080", // Edge color
"faceColor" : "#77AC30", // Face color
"alpha" : 1, // Transparency
"colorMap" : "abaqus", // Color map
"discretize" : false, // Whether to discretize the color map
"discretizeNum" : 12 // Number of discretization intervals
}
}
精度对比
本次要测试的模型如下，一侧固定 XYZ 自由度方向的自由度，上表面承受表面压强荷载100，整个模型受到向下的重力荷载，密度和重力加速度均取9.8，以 C3D4、C3D10、C3D8、C3D20 单元为例，与 Abaqus 做对比验证 MFEA程序的精度。
注： 由于应力分量个数较多，一一对比时占据较多篇幅，故只取Mises应力对比分析。
C3D4单元
| 场信息 | MFEA | Abaqus |
|---|---|---|
| Umag | ||
| U1 | ||
| U2 | ||
| U3 | ||
| RFmag | ||
| RF1 | ||
| RF2 | ||
| RF3 | ||
| Mises |
C3D10单元
| 场信息 | MFEA | Abaqus |
|---|---|---|
| Umag | ||
| U1 | ||
| U2 | ||
| U3 | ||
| RFmag | ||
| RF1 | ||
| RF2 | ||
| RF3 | ||
| Mises |
C3D8单元
| 场信息 | MFEA | Abaqus |
|---|---|---|
| Umag | ||
| U1 | ||
| U2 | ||
| U3 | ||
| RFmag | ||
| RF1 | ||
| RF2 | ||
| RF3 | ||
| Mises |
C3D20单元
| 场信息 | MFEA | Abaqus |
|---|---|---|
| Umag | ||
| U1 | ||
| U2 | ||
| U3 | ||
| RFmag | ||
| RF1 | ||
| RF2 | ||
| RF3 | ||
| Mises |
有趣的小案例
如下图所示的“埃菲尔铁塔”模型，其材料属性为：弹性模量E=1.0E6，泊松比v=0.3，密度ρ=9.8，重力加速度取-9.8。现考虑自重作用下，塔结构的位移场变化，使用MFEA求解并对比Abaqus求解结果。（模型来自B站up主@MOYE_Zgc）
以上就是整理的MFEA和Abaqus的场信息对比，除了C3D8单元之外其余单元的精度均得到很好的验证。Abaqus对C3D8单元采用BAR修正，MFEA V1.2版尚未对单元的刚度矩阵进行特殊修正, 留在V1.3版本进行修改。
除此之外，MFEA V1.3将会配合《有限元基础编程百科全书》增加更多单元类型，将为读者用户展现更多有限元编程的乐趣。
有关程序更新的功能就介绍到这里，感谢您的阅读。整套程序已发布在知识星球中，后台回复：星球，即可加入，对源程序的疑问，可在星球内详细讨论。
参与更多互动交流，快快在下方留言区留下你的小脚印吧~
-End-
易木木响叮当
想陪你一起度过短暂且漫长的科研生活
