1. 有限元分析基础与MATLAB实现概述
有限元分析（FEA）是工程仿真中最常用的数值计算方法之一，它通过将复杂结构离散化为有限数量的小单元来求解力学问题。在二维分析中，平面三角形单元因其网格划分灵活、算法成熟而广泛应用。MATLAB作为科学计算领域的标杆工具，凭借其矩阵运算优势和可视化能力，成为实现有限元分析的理想平台。
我在实际项目中发现，一个完整的MATLAB有限元程序通常包含以下核心模块：前处理（网格生成、边界条件定义）、求解器（刚度矩阵组装、方程求解）和后处理（结果可视化）。与商业软件相比，自主编程的优势在于可以完全掌控算法细节，比如这次我们要实现的INP文件读取功能，就能让程序兼容多种建模工具生成的网格数据。
初学者常遇到的第一个困惑是如何选择单元类型。平面三角形单元（如CST单元）虽然精度不如四边形单元，但在处理复杂几何形状时更具优势。我曾在一个悬臂梁分析项目中对比过，当模型存在曲线边界时，三角形单元的网格生成明显更简单，而计算精度通过加密网格也能满足工程需求。
2. INP文件解析与数据预处理
ABAQUS的INP文件格式因其结构清晰而被广泛借鉴。在我们的MATLAB实现中，自定义的INP文件包含六个关键部分：
* Config # 分析维度与单元类型 * Material # 材料参数 * Nodes # 节点坐标 * Elements # 单元连接关系 * Forces # 载荷信息 * Boundaries # 边界条件通过ReadData.m实现文件解析时，我采用了逐行读取+关键字识别的策略。例如处理节点数据时，先用fgetl获取行内容，当检测到* Nodes关键字后，下一行的数字表示节点总数，随后按节点编号, x坐标, y坐标格式读取数据。这里有个实用技巧：用strsplit处理逗号分隔的字符串，再用str2double转换为数值。
% 示例：读取节点行数据 line = '1, 0.0, 10.0'; data = strsplit(line, ','); node_id = str2double(data{1}); x_coord = str2double(data{2}); y_coord = str2double(data{3});在实际调试中发现，当处理大模型时，预分配数组内存能显著提升性能。比如已知节点数量为n时，先用nodes = zeros(n,2)分配空间，比动态扩展数组效率高出一个数量级。
3. 网格生成与可视化技巧
对于简单几何，可以完全用MATLAB实现网格生成。GenerateInpFromNumSeeds.m展示了一个等腰直角三角形的自动划分算法。核心思路是通过双层循环控制行列种子点，计算每个节点的坐标：
for row = 1:num_seeds+1 for col = 1:num_seeds+row x = side_len * (col-1); y = -side_len * (row-1); fprintf(f, "%d, %f, %f\n", n_dof, x, y); n_dof = n_dof + 1; end end可视化是验证网格质量的重要手段。推荐使用triplot函数显示拓扑关系：
tri = triangulation(elements(:,1:3), nodes); figure triplot(tri, 'LineWidth',1.5) axis equal title('网格拓扑检查')我曾遇到过一个典型问题：生成的网格中存在反转单元（内角大于180°），导致后续计算发散。通过添加面积正负检查避免了这种情况：
area = (x2*y3 - x3*y2) - (x1*y3 - x3*y1) + (x1*y2 - x2*y1); if area <= 0 error('单元%d存在反转，请检查节点顺序', elem_id) end4. 刚度矩阵组装与求解
平面应力问题的单元刚度矩阵计算涉及三个关键矩阵：
- 应变位移矩阵B：3×6矩阵，与节点坐标相关
- 弹性矩阵D：3×3矩阵，取决于材料参数
- 刚度矩阵Ke：6×6矩阵，通过积分得到
在ElementStiffness.m中，我采用解析积分法计算三角形单元刚度矩阵。核心代码片段：
% 计算雅可比行列式 Jacobi = [x1-x3, y1-y3; x2-x3, y2-y3]; detJ = det(Jacobi); % 构造B矩阵 B = [y2-y3, 0, y3-y1, 0, y1-y2, 0; 0, x3-x2, 0, x1-x3, 0, x2-x1; x3-x2, y2-y3, x1-x3, y3-y1, x2-x1, y1-y2] / detJ; % 计算单元刚度矩阵 Ke = 0.5 * thickness * abs(detJ) * B' * D * B;组装全局刚度矩阵时，采用稀疏存储能大幅降低内存占用。但在教学示例中，我建议初学者先用全矩阵理解组装过程：
K_global = zeros(2*n_nodes); for e = 1:n_elements nodes_e = elements(e,1:3); % 获取单元节点 dofs = [2*nodes_e-1, 2*nodes_e]; % 自由度索引 K_global(dofs,dofs) = K_global(dofs,dofs) + Ke; end处理边界条件时，罚函数法是简单有效的选择。通过给对角元素添加大数（如1e8）实现固定约束：
fixed_dof = boundaries(boundaries(:,2)==0, 1); K_global(2*fixed_dof-1, 2*fixed_dof-1) = K_global(2*fixed_dof-1, 2*fixed_dof-1) + 1e8;5. 结果可视化与验证
后处理阶段，通过patch函数绘制应力云图时，控制色彩映射范围很重要：
figure colormap('jet') for e = 1:n_elements % 获取单元节点位移后的坐标 patch('Faces', elements(e,1:3), 'Vertices', nodes_deformed,... 'FaceVertexCData', mises(e), 'FaceColor','interp') end colorbar clim([min(mises), max(mises)]) % 统一色标范围与ABAQUS结果对比时，我习惯用相对误差评估精度：
error_displacement = norm(u_matlab - u_abaqus) / norm(u_abaqus); fprintf('位移相对误差: %.4f%%\n', error_displacement*100)在最近的一个板件分析案例中，20节点模型的MATLAB计算结果与ABAQUS的Mises应力最大偏差仅0.2%，验证了程序的可靠性。不过要注意，三角形单元在弯曲问题中会出现剪切自锁现象，这时需要改用二阶单元或减缩积分。
通过这个全流程实现，不仅能深入理解有限元原理，还能培养解决实际工程问题的能力。建议读者尝试扩展程序功能，比如添加四边形单元支持或动态载荷分析，这对提升编程和力学理解都大有裨益。
