MATLAB三角形单元有限元编程：从刚度矩阵到悬臂梁分析

简介一套用于悬臂梁应力应变分析的MATLAB有限元程序基于三角形单元求解结构力学中的位移、应力与应变分布。面向学习有限元方法FEM的工科学生、科研人员或需要快速搭建二维弹性力学算例的MATLAB开发者适合用于理解从网格离散、刚度矩阵组装到结果后处理的完整流程。压缩包内共4个文件全部为后缀为.m的MATLAB源程序分别实现三角形单元刚度矩阵计算、单元应力求解、整体刚度矩阵组装以及主程序调用整体仅3KB代码精简便于逐行研读与二次修改。自发布以来已有510人学习下载适合作为有限元编程入门或课程作业参考。通过运行该程序读者可以直观看到悬臂梁在荷载作用下的变形结果掌握如何利用MATLAB内置求解器完成线性方程组求解并进一步将节点位移转换为单元应力与应变为后续扩展到更复杂结构的数值分析打下基础。1. 悬臂梁分析里的三角形单元MATLAB 怎么从零拼出一套有限元程序拿到这个新建文件夹.zip的时候里面的文件结构其实很清晰LinearTriangleElementStiffness.m、LinearTriangleElementStresses.m、LinearTriangleAssemble.m加上一个Untitled.m主程序正好构成一套完整的线性三角形单元CST平面应力分析流程。很多人第一次接触 MATLAB 有限元遇到的第一个坎不是理论而是「刚度矩阵算出来了怎么放回全局矩阵」「自由端节点力怎么加进去」「后处理出来的应力为什么和手算对不上」。这套代码恰好把这三个问题都覆盖了。适用对象很清楚正在学有限元课程、需要做平面应力算例验证或者想把自己的理论推导变成可运行代码的工程师。悬臂梁是最适合用来验证三角形单元精度的算例——它既有弯曲变形又有剪切变形固定端约束和一端自由边界条件直接、理论解现成拿它来检验程序里每一步实现是否正确比任何复杂结构都有效。下面从单元刚度矩阵开始把这套小程序逐段拆开。2. 三角形单元刚度矩阵的实现LinearTriangleElementStiffness 的数学内核2.1 常应变三角形单元的刚度矩阵推导这套程序用的是三节点线性三角形单元也就是 CSTConstant Strain Triangle单元。平面应力问题的核心公式是单元刚度矩阵[ \mathbf{k} t A \mathbf{B}^T \mathbf{D} \mathbf{B} ]其中 (t) 是单元厚度(A) 是三角形面积(\mathbf{B}) 是应变-位移矩阵(\mathbf{D}) 是弹性矩阵。对平面应力状态弹性矩阵只依赖弹性模量 (E) 和泊松比 (\nu)[ \mathbf{D} \frac{E}{1-\nu^2} \begin{bmatrix} 1 \nu 0 \ \nu 1 0 \ 0 0 \frac{1-\nu}{2} \end{bmatrix} ]三个节点的坐标是 ((x_1,y_1))、((x_2,y_2))、((x_3,y_3))每个节点两个自由度所以单元刚度矩阵是 (6 \times 6)。(\mathbf{B}) 矩阵由三个形函数对 (x)、(y) 求导得到常数项全部由节点坐标差构成。三角形面积可以用行列式计算(2A x_1(y_2-y_3) x_2(y_3-y_1) x_3(y_1-y_2))如果按逆时针排列节点面积值恒为正这也是程序里要求节点按逆时针输入的隐含约定。LinearTriangleElementStiffness.m函数的完整签名一般是function k LinearTriangleElementStiffness(E, NU, t, node1_coord, node2_coord, node3_coord)这里node1_coord、node2_coord、node3_coord是 (1 \times 2) 的行向量内部先计算面积、组装 (\mathbf{B})、组装 (\mathbf{D})最后做三次矩阵乘法得到 6×6 的k。函数返回值直接交给组装环节使用。2.2 B 矩阵和 D 矩阵的具体组装方式(\mathbf{B}) 矩阵的结构是每个节点对应一个 (3 \times 2) 子块。对节点 i子块为[ \mathbf{B}_i \frac{1}{2A} \begin{bmatrix} b_i 0 \ 0 c_i \ c_i b_i \end{bmatrix} ]其中 (b_i y_j - y_k)(c_i x_k - x_j)这里 (j)、(k) 是另外两个节点的编号按逆时针循环取。MATLAB 实现时可以直接用坐标差来计算不需要先求形函数再求导因为线性单元的形函数导数就是常数。实际代码里常见的写法是x1 node1_coord(1); y1 node1_coord(2); x2 node2_coord(1); y2 node2_coord(2); x3 node3_coord(1); y3 node3_coord(2); A 0.5 * abs(x1*(y2-y3) x2*(y3-y1) x3*(y1-y2)); beta1 y2 - y3; gamma1 x3 - x2; beta2 y3 - y1; gamma2 x1 - x3; beta3 y1 - y2; gamma3 x2 - x1; B 1/(2*A) * [beta1 0 beta2 0 beta3 0; 0 gamma1 0 gamma2 0 gamma3; gamma1 beta1 gamma2 beta2 gamma3 beta3];这段代码里的abs是为了防止节点输入顺序不一致导致面积为负。注意这里A用的是绝对值但beta、gamma的计算依赖于节点顺序——如果顺序颠倒应变矩阵的部分分量会变号。因此最稳妥的做法是一开始就约定逆时针输入并在注释里写明。程序包里如果没写我一般会建议加一行注释说明这个约定避免后续对结果产生困惑。2.3 单元刚度矩阵组装时需要厘清的自由度编号每个单元有 6 个自由度排列顺序是[u1 v1 u2 v2 u3 v3]对应三个节点各自的 (x) 方向位移和 (y) 方向位移。这个顺序必须和LinearTriangleAssemble.m里的编号规则保持一致。常见做法是给每个节点一个全局自由度编号节点 i 的 (x) 方向自由度(2i-1)节点 i 的 (y) 方向自由度(2i)举例来说如果悬臂梁被划分为 200 个三角形单元节点总数为 (N)那么全局刚度矩阵 (\mathbf{K}) 的尺寸就是 (2N \times 2N)。单元刚度矩阵 k 中的第 1、2 行/列对应节点 1 的 (x)、(y) 自由度第 3、4 行/列对应节点 2第 5、6 行/列对应节点 3。组装时把这些局部自由度映射到全局编号逐个填进大矩阵对应位置。我见过不少人在这个环节出错——单元矩阵算对了但组装时局部自由度对错了全局位置结果整个刚度矩阵不对称或者解出来的位移完全没意义。一个简单的自检方法是组装完成后检查 K 是否对称非对称就说明编号映射有问题。3. LinearTriangleAssemble 组装循环与固定端约束的处理3.1 全局刚度矩阵的组装策略全局刚度矩阵 (\mathbf{K}) 是一个稀疏的带状矩阵。200 个三角形单元的网格规模不算大但理论上说直接用一个 (2N \times 2N) 的稠密矩阵存储当节点数超过 500 时内存占用就会明显上升。程序包里的LinearTriangleAssemble.m一般写成循环遍历所有单元、逐一累加的形式function K LinearTriangleAssemble(K, k, node1, node2, node3) % K : 当前全局刚度矩阵尺寸为 2N x 2N % k : 单元刚度矩阵尺寸为 6 x 6 % node1/2/3 : 单元的三个全局节点编号 dof1 [2*node1-1, 2*node1]; dof2 [2*node2-1, 2*node2]; dof3 [2*node3-1, 2*node3]; dofs [dof1, dof2, dof3]; K(dofs, dofs) K(dofs, dofs) k; end这段代码的逻辑是把单元的 6 个自由度编号拼成一个向量dofs然后利用 MATLAB 的矩阵索引一次性完成 6×6 分块的累加。理解这个函数的关键在于dofs的构造方式每个节点对应两个连续编号所有单元共用同一套编号规则才能保证位移连续性。如果有两个相邻单元共享一条边、两个节点那么它们的单元刚度矩阵贡献会累加到全局矩阵的相同位置这个过程叫“叠加”。主程序Untitled.m里调用这段组装的思路通常是K zeros(2*N, 2*N); for e 1:numel(elements) node1 elements(e, 1); node2 elements(e, 2); node3 elements(e, 3); k LinearTriangleElementStiffness(E, NU, t, nodes(node1,:), nodes(node2,:), nodes(node3,:)); K LinearTriangleAssemble(K, k, node1, node2, node3); end3.2 固定端约束的实现方式悬臂梁左端固定意味着所有固定端节点的u和v自由度位移都为零。在 MATLAB 里常见的做法是“置一法”把全局刚度矩阵中对应固定自由度的行和列全部清零对角线置 1荷载向量中对应位置也置 0。这样做的好处是保留矩阵的稀疏结构避免重新编号带来的麻烦。写出来大致是这样fixed_dofs []; % 收集固定端节点的所有自由度编号 fixed_nodes 1:6; % 示例固定端包含节点1到6 for i fixed_nodes fixed_dofs [fixed_dofs, 2*i-1, 2*i]; end K(fixed_dofs, :) 0; K(:, fixed_dofs) 0; K(fixed_dofs, fixed_dofs) eye(length(fixed_dofs)); F(fixed_dofs) 0;这个做法的底层逻辑是固定自由度的位移已知为零把刚度矩阵中这些自由度与其他自由度的耦合项全部切断只保留对角线上的 1使得求解方程 ( \mathbf{K} \mathbf{U} \mathbf{F} ) 中该自由度的方程变为 (1 \times u 0)自然解出零位移。耦合信息丢失不影响其他自由度的计算结果因为物理上固定端的位移本来就为零不需要参与求解。需要注意的是固定端节点集合要根据实际网格编号来定。如果把左端视为 y 轴上的节点编号一般在 1 到某个数之间具体编号取决于网格生成时的节点排序方式。我会建议先在Untitled.m里加一行disp(nodes(1:20, :))检查节点坐标和编号的对应关系再确定固定节点的范围这样比肉眼数网格可靠得多。3.3 自由端荷载向量的构造自由端的节点力需要根据等效节点荷载的原则施加。如果你是集中力作用在自由端中点那么这个力就直接加在对应节点的 y 自由度上。如果荷载是均布力就需要按静力等效原则分配到多个节点上。程序包里的场景是悬臂梁自由端受集中载荷荷载向量 F 的构造要留心方向沿着 y 负方向施加在 MATLAB 中对应负值。常见的错误是忘记在荷载向量里把载荷除以节点数——如果你把总荷载 P 平均分到若干个自由端节点每个节点分到的力是 P 除以节点数而不是每个节点都写 P。下表整理了不同荷载情况下的处理方式供参考荷载类型处理方法操作示例自由端集中力直接加到节点的 y 自由度F(2*free_node, 1) -P;均布荷载按长度分配到多个节点F(2*node_list, 1) -q*dx/2;端部节点减半自重荷载转换为体积力加到每个单元节点需要单元级等效节点力再组装到 F节点力矩转化为等效节点力偶需要额外处理CST 单元不支持转角自由度对于悬臂梁来说最常见的组合是固定端零位移 自由端集中力或均布力。如果你的题目是均布荷载我一般会建议不要偷懒把所有力加在一个节点上而是按相邻节点的间距分配到多个节点——虽然总力一样但节点分布不同会影响局部应力峰值。4. 200 个三角形单元的网格生成与求解流程4.1 结构化网格生成与节点编号要对悬臂梁划分三角形单元最直接的方式不是用delaunay随机生成而是先构造四边形网格再对每个四边形切分成两个三角形。这样做的好处是节点分布规则、单元编号可控后处理提取某个截面上的应力也方便。假设悬臂梁长 L、高 H沿长度方向划分 nx 段沿高度方向划分 ny 段nx 20; ny 10; % 长度和高度方向分段数 L 4; H 0.4; % 几何尺寸单位按工程需要 x_nodes linspace(0, L, nx1); y_nodes linspace(0, H, ny1); nodes zeros((nx1)*(ny1), 2); idx 0; for j 1:ny1 for i 1:nx1 idx idx 1; nodes(idx, :) [x_nodes(i), y_nodes(j)]; end end这样生成的节点矩阵nodes每一行是一个节点的坐标编号先沿 x 方向再沿 y 方向递增。下面把每个矩形单元切成两个三角形。一个矩形有四个顶点(i,j)、(i1,j)、(i,j1)、(i1,j1)对应的节点编号分别是p1、p2、p3、p4。切成两个三角形时公共边可以沿 p1-p4 或者沿 p2-p3前者单元形状更稳定推荐使用。elements []; for j 1:ny for i 1:nx p1 (j-1)*(nx1) i; p2 p1 1; p3 p1 (nx1); p4 p3 1; elements [elements; p1, p2, p4; p1, p4, p3]; end end每个矩形切成两个三角形后单元总数正好是 (2 \times nx \times ny)当 nx20、ny10 时就是 400 个三角形。如果题目明确说划分为 200 个三角形单元那可以将 nx10、ny10对应 100 个矩形单元。这个映射关系要注意区分200 个三角形单元与 200 个节点不是同一个概念。三角形单元排列顺序影响的是单元矩阵的组装顺序不影响最终结果。但一个问题值得注意如果对角线方向不统一相邻单元的公共边方向可能不一致导致应力结果在某些单元交界处出现轻微不连续这是 CST 单元的固有特性不是程序的 bug。想改善这个现象可以加密网格或者在几何形状允许时统一对角线方向。4.2 用 mldivide 求解位移场全局刚度矩阵 K 和荷载向量 F 都组装好、边界条件处理完之后求解节点位移就是一个线性代数问题。用 MATLAB 的\运算符是最省事的写法U K \ F;这个运算符会根据矩阵的特性自动选择合适的求解器。当 K 是稀疏对称正定矩阵时MATLAB 通常走 Cholesky 分解路径如果 K 因为是稠密存储而变得规模较大求解会明显变慢。我的习惯是在组装前先声明K sparse(2*N, 2*N)把 K 定义为稀疏矩阵后续的累加操作不变但内存和求解速度都会改善。对 200 个单元量级的问题这个优化感受不明显但当网格超过 1000 个单元时差别会非常显著。U的排列顺序和自由度的编号规则一致U(2*i-1)是节点 i 的 x 方向位移U(2*i)是节点 i 的 y 方向位移。提取某个节点的位移时u_free U(2*free_node - 1); % 自由端节点的水平位移 v_free U(2*free_node); % 自由端节点的竖向位移求解完成后一个值得做的检查是观察固定端节点的位移是否严格为零——由于浮点运算和置一法处理这个值通常在 (10^{-16}) 量级可以视为零。如果发现固定端位移很大那一定是边界条件处理有误。4.3 求解前后各阶段的调试要点有限元程序最常见的出错点集中在三个阶段网格生成、刚度矩阵组装、边界条件处理。网格阶段要检查的是单元内有没有面积为零的退化三角形一个简单方法是求每个单元的最小边长小于某个阈值就报错。刚度矩阵组装阶段要检查 K 是否对称、对角线元素是否为正。边界条件阶段要检查固定节点的自由度在 U 中是否解出接近零的值。我建议在Untitled.m中加一个中间验证步骤% 检查全局刚度矩阵的对称性 err_sym norm(K - K, fro); if err_sym 1e-10 warning(刚度矩阵不对称误差: %e, err_sym); end % 检查自由度编号是否越界 NN size(nodes, 1); if max(max(elements)) NN || min(min(elements)) 1 error(单元节点编号超出范围); end这段检查不要删——程序包本身可能没有这个功能但在你自己的调试过程中这两个检查几乎是免费的保险。K 不对称几乎可以断定是LinearTriangleAssemble里自由度映射写错了。5. 应力应变后处理与精度验证的实用技巧5.1 单元应力应变恢复的基本流程节点位移解出后每个单元的应力和应变可以用单元刚度矩阵类似的方式恢复。应变由 (\boldsymbol{\varepsilon} \mathbf{B} \mathbf{u}^{(e)}) 计算其中 (\mathbf{u}^{(e)}) 是单元的 6 个节点位移分量应力由 (\boldsymbol{\sigma} \mathbf{D} \mathbf{B} \mathbf{u}^{(e)}) 得到。程序包里的LinearTriangleElementStresses.m做的事情就是这两步function [sigma, epsilon] LinearTriangleElementStresses(E, NU, node1_coord, node2_coord, node3_coord, u) % u : 单元节点位移向量 [u1 v1 u2 v2 u3 v3]6x1 % 先构造 B 和 D再计算应变和应力 epsilon B * u; sigma D * epsilon; end对应变结果要有一个心理预期CST 单元在每个单元内应力和应变是常数所以画云图的时候每个单元内部的颜色是一块色块不会像四边形或高阶单元那样平滑过渡。自由端附近如果有集中力局部区域会出现应力集中的假象因为单元在该处应力值高、相邻单元过渡陡峭。要减少这种误差可以在自由端附近适当加密网格或者对结果做节点平均后再绘制。5.2 与悬臂梁理论解的对比验证验证程序正确性的核心是理论解对比。悬臂梁受端部集中力 P 作用时梁理论给出固定端最大弯曲应力为[ \sigma_{\text{max}} \frac{6PL}{b h^2} ]关键点在于程序里使用的单位必须一致。如果长度用米、力用牛顿弹性模量用帕斯卡Pa得出的应力和位移就是国际单位制。一个常见的坑是弹性模量使用了 GPa 而长度用了毫米两者没有统一换算导致结果差了几个数量级。自由端挠度的理论公式是[ v_{\text{free}} \frac{PL^3}{3EI} ]其中 (I b h^3 / 12)(b) 是梁的宽度(h) 是高度。将程序求出的v_free与理论值对比误差通常在 5%~15% 之间具体取决于三角形单元网格的划分方式。CST 单元是常应变单元对弯曲问题存在“剪切锁死”现象导致计算出的挠度偏小、结构偏刚。网格越密误差越小。如果你发现误差达到 30% 以上优先怀疑的是单元数量不够或边界条件处理有误而不是理论公式用错了。我习惯用下面的方式做快速验证% 理论解 I b*h^3/12; v_theory P*L^3/(3*E*I); sigma_theory 6*P*L/(b*h^2); % 数值解 v_num max(abs(U(2:2:end))); % 所有节点的竖向位移最大值通常出现在自由端 sigma_top max(abs(sigma(:, 1))); % 取 x 方向应力的最大值 % 输出对比 fprintf(自由端挠度: 理论 %e, 数值 %e, 误差 %f%%\n, v_theory, v_num, abs(v_num-v_theory)/v_theory*100); fprintf(最大应力: 理论 %e, 数值 %e, 误差 %f%%\n, sigma_theory, sigma_top, abs(sigma_top-sigma_theory)/sigma_theory*100);这里v_num取所有节点竖向位移的最大值物理上对应自由端的最大挠度sigma_top取所有单元 x 方向应力的最大值理论上应该出现在固定端附近的单元。如果最大应力出现在自由端单元说明荷载施加区域存在严重的应力集中需要细化网格或改用面积等效方式加载。5.3 应力云图绘制与结果导出的收尾操作最后一步是把位移和应力画出来。MATLAB 里用patch或trisurf都可以。我推荐patch因为它可以直接按单元填充颜色并且能显示节点坐标和单元连接关系figure; patch(Faces, elements, Vertices, nodes(:,1:2), FaceVertexCData, stress_avg, ... FaceColor, flat, EdgeColor, none); colorbar; colormap(jet); axis equal; xlabel(x (m)); ylabel(y (m)); title(x方向应力分布 \sigma_x);stress_avg是每个单元的一个标量值比如 x 方向应力FaceVertexCData按单元填充颜色。如果想让云图更平滑可以对节点做应力平均把共享该节点的所有单元应力取平均再赋给节点。两种方式各有优劣单元应力云图更能反映 CST 单元的离散特性节点平均云图更接近视觉直觉。导出结果时fprintf配合格式化控制符就够了。如果想导出完整数据供后处理用writematrix或xlswrite会更快results [(1:size(elements,1)), sigma]; writematrix(results, element_stress.csv); headers {ElementID, sigma_x, sigma_y, tau_xy}; writecell(headers, element_stress_headers.csv);至此从单元刚度矩阵到全局组装、从边界条件到应力恢复的完整流程已经打通。往后在这个基础上换材料参数、换荷载工况或者把三角形单元换成四边形单元都可以在现有框架内扩展。本文还有配套的精品资源点击获取