MFEAOOP V1.4版本震撼来袭！更丰富的单元库！

题外话：因为微信的推荐机制变动，有可能大家不会第一时间看到我的文章，请大家给我的公众号标上⭐，以免错过好资源。

1.4版本更新内容

1. 增加单元类型：
1. 轴对称单元: CAX3、CAX4、CAX6、CAX8、CAX8R，均支持表面压强荷载和重力荷载
2. 考虑剪切变形的 Euler-Bernoulli 平面梁单元B23M
2. 增加89种colormap配色方案：matlab内置+Matplotlib官网+自定义配色方案
3. 增加梁单元内力(轴力、剪力、弯矩)计算输出
1. 二维梁：Fx Fy RM3
2. 三维梁：Fx Fy Fz RM1 RM2 RM3
4. 所有单元类型，除了梁单元，其余均采用mex文件加速，组装整体刚度矩阵的时间可以提升一倍多，梁单元因为截面类型不同，所需参数不同，不利于程序的统一实现，因此不采用mex文件加速
5. 后处理部分增加节点标号显示："nodeLabel" : true
6. 增加MPC多点约束功能（对应于Abaqus的TIE功能，即指定节点的自由度被约束）

MFEAOOP V1.4所有内容

1. 支持荷载类型：
1. 位移荷载
2. 表面压强荷载
3. 重力荷载
2. 单元类型（共计支持39种单元类型）：
1. C3D4、C3D10、C3D8、C3D20、C3D8R、C3D20R、C3D8I、C3D8Bar、C3D8SR
2. CAX3、CAX4、CAX6、CAX8、CAX8R
3. CPS4、CPS8、CPS3、CPS6、CPS4R、CPS8R、CPS4I、CPS4Bar、CPS4SR
4. CPE4、CPE8、CPE3、CPE6、CPE4R、CPE8R、CPE4I、CPSEBar、CPE4SR
5. T2D2、T3D2
6. B21、B23、B31、B33、B23M
7. 杆系单元：
8. 平面单元：
9. 轴对称单元：
10. 空间实体单元：
3. 整体刚度矩阵组装：稀疏矩阵存储格式，使用三元组（Triplet）形式，即行索引数组、列索引数组和非零值数组。使用mex文件加速，组装整体刚度矩阵的时间可以提升一倍多。
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 文件，记录计算各个过程的具体时间明细以及一些模型信息。

MPC多点约束

MPC约束节点集的格式如下：

   {
   "name": "fix2",
   "type": "mpc",
   "couplings": [
         {
            "master_nodes": "Set-M",   // 约束节点集合1      
            "slave_nodes": "Set-S",       // 约束节点集合2
            "dof": [1, 2]  // 约束的自由度（说明1自由度和2自由度均被约束）
         },
         {
            "master_nodes": "Set-p11",         
            "slave_nodes": "Set-p12",      
            "dof": [1]                       
         }
   ]
}

上面的约束条件表示：约束节点集合Set-M和节点集合Set-S的1、2自由度进行TIE连接（变形过过程中节点集Set-M和Set-S的1、2自由度值保持一致），约束节点集合Set-p11和节点集合Set-p12的1自由度进行TIE连接。

接下来以一个平面模型为例，展示程序的MPC功能。

程序中将11-12号节点以及4-5号节点的1、2自由度进行TIE连接，位移云图中，这两对节点值是完全一样的。

该部分程序实现方式的理论参考：罚函数处理边界条件

mex文件加速

在MATLAB中，MEX文件是一种利用C、C++或Fortran编写的外部函数，可以在MATLAB中直接调用。这种文件具有显著的加速效果，尤其适合需要大量计算或重复运算的程序。通过MEX文件，将计算密集型的代码段编译为机器码，MATLAB可以直接调用这些代码，而无需解释执行，从而获得接近原生C/C++程序的执行速度。

将程序中的所有单元刚度矩阵计算函数一起转换为mex文件，可编制以下脚本：

clc; clear
% 梁单元不进行mex处理，因为截面参数多变，不利于统一封装
% 设置基本配置
cfg = coder.config('mex');
cfg.EnableOpenMP = true;    % 启用多线程支持

% 定义单元类型列表和相应的节点数量
elementTypes = {
    'T2D2', 'T3D2', ...
    'CPS4', 'CPS8', 'CPS3', 'CPS6', 'CPS4R', 'CPS8R', 'CPS4I', 'CPS4Bar', 'CPS4SR', ...
    'CPE4', 'CPE8', 'CPE3', 'CPE6', 'CPE4R', 'CPE8R', 'CPE4I', 'CPE4Bar', 'CPE4SR', ...
    'C3D4', 'C3D10', 'C3D8', 'C3D20', 'C3D8R', 'C3D20R', 'C3D8I', 'C3D8Bar', 'C3D8SR'
};

% 定义默认输入参数结构体示例
 % 示例材料属性

% 遍历每种单元类型并生成 MEX 文件
for i = 1:numel(elementTypes)
    elementType = elementTypes{i};
    
    % 根据不同的单元类型设置节点数和维度
    switch elementType
        case {'T2D2'}
            materialProps = struct('E', 210e9, 'A', 0.3);
            nodeCoordinates = zeros(2, 2); % 2D 桁架单元，2 个节点，每节点 2 个坐标
        case {'T3D2'}
            materialProps = struct('E', 210e9, 'A', 0.3);
            nodeCoordinates = zeros(2, 3); % 3D 桁架单元，2 个节点，每节点 3 个坐标
        case {'CPS4', 'CPS4R', 'CPS4I', 'CPS4Bar', 'CPS4SR', ...
              'CPE4', 'CPE4R', 'CPE4I', 'CPE4Bar', 'CPE4SR'}
            materialProps = struct('E', 210e9,'v', 0.3,'t', 0.3);
            nodeCoordinates = zeros(4, 2); % 2D 四节点平面单元，4 个节点，每节点 2 个坐标
        case {'CPS8', 'CPS8R'}
            materialProps = struct('E', 210e9,'v', 0.3,'t', 0.3);
            nodeCoordinates = zeros(8, 2); % 2D 八节点平面单元，8 个节点，每节点 2 个坐标
        case {'CPS3', 'CPE3'}
            materialProps = struct('E', 210e9,'v', 0.3,'t', 0.3);
            nodeCoordinates = zeros(3, 2); % 2D 三节点平面单元，3 个节点，每节点 2 个坐标
        case {'CPS6', 'CPE6'}
            materialProps = struct('E', 210e9,'v', 0.3,'t', 0.3);
            nodeCoordinates = zeros(6, 2); % 2D 六节点平面单元，6 个节点，每节点 2 个坐标
        case {'CPE8', 'CPE8R'}
            materialProps = struct('E', 210e9,'v', 0.3,'t', 0.3);
            nodeCoordinates = zeros(8, 2); % 2D 八节点平面单元，8 个节点，每节点 2 个坐标
        case 'C3D4'
            materialProps = struct('E', 210e9, 'v', 0.3);
            nodeCoordinates = zeros(4, 3); % 3D 四节点四面体单元，4 个节点，每节点 3 个坐标
        case 'C3D10'
            materialProps = struct('E', 210e9, 'v', 0.3);
            nodeCoordinates = zeros(10, 3); % 3D 十节点四面体单元，10 个节点，每节点 3 个坐标
        case {'C3D8', 'C3D8R', 'C3D8I', 'C3D8Bar', 'C3D8SR'}
            materialProps = struct('E', 210e9, 'v', 0.3);
            nodeCoordinates = zeros(8, 3); % 3D 八节点六面体单元，8 个节点，每节点 3 个坐标
        case {'C3D20', 'C3D20R'}
            materialProps = struct('E', 210e9, 'v', 0.3);
            nodeCoordinates = zeros(20, 3); % 3D 二十节点六面体单元，20 个节点，每节点 3 个坐标
        otherwise
            warning('Unknown element type: %s. Skipping...', elementType);
            continue;
    end
    
    try
         eval(['codegen -config cfg ' elementType ' -args {materialProps, nodeCoordinates};']);
        fprintf('Successfully generated MEX file for %s.\n', elementType);
    catch ME
        fprintf('Failed to generate MEX file for %s: %s\n', elementType, ME.message);
    end

end

从以上程序可以看到在转换过程中，需要明确给定形参的初始值，在设计刚度矩阵时，我均采用这样的形参，以CPS4单元为例：

function k=CPS4(prop, elemNodeCoordinate)

    DOF = 2;
    elenode = length(elemNodeCoordinate);
    k = zeros(DOF*elenode,DOF*elenode);

    t = prop.t;
    stressType = 'planeStress';
    D = Dmatrix2D(prop, stressType);

    ngp = 4; 
    [weights,locations]=gauss2D(ngp);

    for n=1:size(weights,1)  
        xi_Gauss=locations(n,1);
        eta_Gauss=locations(n,2);

        [B, Jacob] = Bmatrix2D(elenode, elemNodeCoordinate, DOF, xi_Gauss, eta_Gauss);
        k = k+B.'*D*t*B*weights(n)*det(Jacob);
    end

end

其中prop是从json读取的材料参数，是个可变参数，单元类型不同，用到的参数也不同，平面单元中用到的是弹性模量、泊松比、厚度，所以我在设计：materialProps = struct('E', 210e9,'v', 0.3,'t', 0.3)。

轴对称单元

为丰富程序的单元库，本次更新特意新添加了五种轴对称单元：CAX3、CAX4、CAX6、CAX8、CAX8R，具体的理论部分我会找时间更新至《有限元基础编程百科全书》，以CAX8单元为例，工况如下：

大家可能注意到MFEAOOP的应力场有些怪怪的，即x=0的轴上为空白，这是为什么呢？程序刚运行出这样的结果的时候我也会很惊讶，后来仔细思考，发现这是合理的，因为我在节点应力计算时，我直接将节点的坐标代入应力计算公式，而不是高斯积分点，所以会导致x=0处的节点应力出现奇异的情况。

这个时候可以更改应力计算方案，将高斯积分点代入应力计算公式，然后进行外插即可，但是这样以来，程序需要根据不同的单元类型进行选择合适的形函数，多少有些繁琐。

我采用的方法时，将x=0的节点x坐标设为0.0001，做一个小小的修正，即可以避免数值奇异的情况，效果如下：

colormap

为增加后处理的灵活性，我将matplotlib的colormap添加到MFEAOOP中，供用户使用，效果如下：

bug修复

1. 梁单元在约束自由度时，比如要约束全部自由度，需要xyzUR1UR2UR3

程序心得

1. 编码格式要统一
2. 尽量使用英文报错、警告信息
3. 程序出现bug，不要慌，根据报错信息，查询错误源头，如果不会高级调试就使用最简单的方式：打印变量
4. 如果不想再看程序了，那就关闭程序，隔几天后再来编写你的程序，这样你的思路会清晰很多
5. 不追求华丽花哨的代码，追求简洁、易读、易维护的代码

由于添加的单元类型较多，精度验证需要占据大量篇幅，没有贴在本次推文中，不过大可放心，每种单元类型在添加时均和Abaqus做了验证分析，均验证通过！感兴趣的小伙伴可逐个验证。

程序目前先不继续更新功能了，接下来就是补充完善理论部分，及时更新《有限元基础编程百科全书》，程序的思维导图也在绘制中，不久将全部公布，感谢大家关注。

有关程序更新的功能就介绍到这里，感谢您的阅读。整套程序已发布在知识星球中，后台回复：星球，即可加入，对源程序的疑问，可在星球内详细讨论。

本次推文已上传至Web端，如果觉得公众号界面阅读不便，可点击下方“阅读原文”，即可跳转至木木的个人网页：