简介面向固体力学有限元初学者的MATLAB源码包围绕单元矩阵组装这一核心环节展示从几何建模、网格离散、材料属性设置到全局刚度矩阵集成的完整流程适合力学、土木、机械方向学生与工程技术人员边读边练。压缩包共5个文件全部为.m脚本整体仅12KB轻量且便于逐行理解。脚本以多个一维、二维算例为线索常见几何对象包括矩形区域与圆柱体等读者可从局部单元矩阵出发逐步学习边界条件施加、载荷处理、线性方程组求解以及后处理绘图等关键步骤体会连续区域被离散为有限元后如何通过矩阵组装获得整体解。已有311人学习下载资源可作为课程作业、毕业设计或小型科研案例的改写模板帮助读者快速打通有限元理论与MATLAB编程实现之间的映射关系。1. 这个标题背后是一整套 MATLAB 固体有限元的起手式你搜到“FEM 有限元.zip”这类资源包时看到的几个关键词其实是同一件事的三个面向FEM 是方法MATLAB 是工具而“组装单元矩阵”才是整个求解流程里最容易被忽略、也最常出错的环节。网上流传的 matlab 有限元编程求解实例十有八九会在一维杆或二维平面问题上演一遍这个流程剖网格、写单刚、组总刚、加边界、求解、后处理。乍看每一步都有公式可真到自己动手时卡住的往往不是推导而是“这个数到底该加到刚度矩阵的哪一个位置”。这篇博文就沿着这条最朴素但完整的路线走一遍。目标是让一个只记得有限元基础概念的人也能用 MATLAB 把二维固体问题的求解闭环跑通从节点坐标和单元连接表出发依次完成单元刚度矩阵计算、全局组装、边界条件施加和应力后处理。适合两类人一类是刚学到刚度矩阵一章、想用代码把公式落到实处的学生另一类是工程中只想快速验证一个结构方案、没时间开大型有限元仿真软件的开发者。2. 固体有限元的数据结构与自由度布局组装前先想清楚编号2.1 弱形式与刚度矩阵的物理含义为什么“组装”是必经之路固体力学有限元的出发点是弹性体的平衡方程结合几何方程与材料本构关系后直接求解这类偏微分方程并不容易。有限元的做法是将求解区域划分成有限个单元在每个单元内假设位移的插值形式再通过虚位移原理得到一组代数方程。最终形成的就是经典的线性系统K 乘以所有节点的位移向量等于节点力向量。这里的 K 就是全局刚度矩阵它的来源是每个单元的单元刚度矩阵而把这些局部矩阵按自由度位置“叠加”进同一个全局矩阵的操作就是组装。组装不是单纯的数学技巧它背后是力平衡的物理要求。相邻两个单元共享的节点上来自两个单元的贡献必须相加而不是相互覆盖否则节点上的力就不平衡。换句话说全局刚度矩阵 K 的每一行代表的是某个自由度在所有相邻单元中的刚度贡献和。这也是为什么“组装单元矩阵”会成为这类 MATLAB 程序的关键词它直接决定了 K 是否描述了一个真实的结构整体。2.2 MATLAB 里用两个矩阵承载网格节点坐标与单元连接任何一个固体有限元程序不论一维还是三维都离不开两个最基本的数据结构节点坐标矩阵和单元连接矩阵。MATLAB 的惯例是让每个节点占一行行号就是节点的全局编号。以二维问题为例节点坐标矩阵是 n_nodes 行、2 列的数组单元连接矩阵的每一行存放一个单元的所有节点全局编号按逆时针或统一方向排列。这里用一个贯穿全文的网格例子一个 2×1 的矩形板由两个四节点平面单元Q4组成共 6 个节点。对应的 MATLAB 数据如下% nodes(i,:) 节点i的(x, y)坐标 nodes [ 0 0; 1 0; 2 0; 0 1; 1 1; 2 1 ]; % elems(i,:) 第i个单元的4个节点全局编号按逆时针给出 elems [ 1 2 5 4; 2 3 6 5 ];单元连接矩阵的节点顺序非常关键。第一个单元用节点 1、2、5、4形成左下角那个 1×1 正方形第二个单元用节点 2、3、6、5形成右下角那个正方形。如果写成 1 2 4 5得到的会是一个交叉扭曲的四边形Jacobian 行列式可能为负后面的积分计算就会出错。这类错误在 MATLAB 里通常不会报错只会表现为位移结果完全不合理。2.3 自由度编号二维固体每个节点有两个故事“组装”最容易被低估的一步是自由度编号。一维杆单元里每个节点只有一个位移自由度节点编号和自由度编号几乎可以混用这也让很多人在转向二维问题时吃了大亏二维固体的每个节点有 x 和 y 两个方向的位移所以节点 1 占全局自由度的第 1、2 列节点 2 占第 3、4 列依此类推。上面网格的自由度分布可以列出如下关系节点编号x 方向自由度y 方向自由度112234356478591061112当处理第一个单元节点 1、2、5、4时局部自由度 1 到 8 对应的全局自由度分别是 1、2、3、4、9、10、7、8。这里的顺序必须与单元刚度矩阵的行列顺序保持一致否则 K 矩阵装配出来之后载荷施加和位移提取都会对不上。节点位移排列方式没有统一规定但建议始终采用“按节点分组、节点内先 x 后 y”的方案。这样做的理由是便于阅读从解向量里挑出节点 i 的 x 位移只需取 u(2*i-1) 即可后处理和画云图时一步到位。下面是两种常见自由度布局的对比方案节点内自由度顺序全局向量的排列按节点分组[x, y][n1_x, n1_y, n2_x, n2_y, ...]按方向分组[所有 x, 所有 y][n1_x, n2_x, ..., n1_y, n2_y, ...]如果要在其他软件里交换数据按方向分组更方便某些求解器的内存布局但在 MATLAB 里按节点分组写起来最直白后续用 reshape 处理位移云图也更顺手。选定后不要中途切换否则组装函数和边界条件代码要全部重写。3. MATLAB 单元刚度矩阵与全局组装从一维杆热手到二维 Q4 实装3.1 一维杆单元先热手自由度映射与循环组装在跳进二维平面单元之前先用一维杆单元把“组装”的索引机制看明白。一根由三段节点、两个杆单元组成的结构每个单元的长度为 L弹性模量 E截面积 A。杆单元的单元刚度矩阵是 2×2 的% 一维杆单元组装示例单元自由度为节点编号本身 EA 1; L 1; elems [1 2; 2 3]; % 两个单元依次连接 K zeros(3, 3); for e 1:size(elems, 1) k EA / L * [1 -1; -1 1]; % 杆单元单刚 dofs elems(e, :); % 局部自由度映射到全局 K(dofs, dofs) K(dofs, dofs) k; end这段代码的要点在于 K(dofs, dofs) 用的是累加而不是赋值。如果写成K(dofs, dofs) k第二个单元会覆盖掉第一个单元在节点 2 上的贡献得到的总刚矩阵将是一个对角线矩阵物理上完全错误。一维杆问题里由于自由度编号恰好等于节点编号这个映射关系显得透明到几乎不存在但它正是“组装”的核心操作做一次局部坐标到全局坐标的索引查找然后叠加。3.2 二维 Q4 单元的单刚形函数、几何矩阵与高斯积分3.2.1 为什么用四节点等参单元二维连续体问题最常用的教学单元是四节点四边形单元 Q4。它的每个单元有 4 个节点、8 个自由度单元刚度矩阵是 8×8 的。相比三角形常应变单元 CSTQ4 能描述线性变化的应变场精度更高相比高阶六节点三角形或九节点四边形Q4 的公式推导和编码难度低得多。对初学者来说这是性价比最高的选择。Q4 单元在局部坐标系ξ, η下定义四个节点的局部坐标分别为 (-1,-1)、(1,-1)、(1,1)、(-1,1)。位移插值由双线性形函数完成而把局部坐标映射到物理坐标的桥梁是 Jacobian 矩阵。单元刚度矩阵的标准计算公式为k_e ∫ B^T D B t dΩ其中 D 是弹性矩阵B 是几何矩阵应变-位移矩阵t 是厚度积分数值上用高斯积分完成。3.2.2 可直接运行的 MATLAB 单刚函数下面这个函数接受单元四个节点的物理坐标输出 8×8 的单元刚度矩阵。这是整个有限元程序里复用率最高的一个模块function k q4_stiffness(nodes_e, E, nu, t, plane_type) % nodes_e : 4x2 矩阵第i行是第i个局部节点的(x, y)坐标 % plane_type: stress(平面应力) 或 strain(平面应变) if strcmp(plane_type, stress) D E/(1-nu^2) * [1 nu 0; nu 1 0; 0 0 (1-nu)/2]; else D E/((1nu)*(1-2*nu)) * [1-nu nu 0; nu 1-nu 0; 0 0 (1-2*nu)/2]; end gps [-1/sqrt(3), 1/sqrt(3)]; % 2点高斯积分坐标 k zeros(8,8); for i 1:2 for j 1:2 xi gps(i); eta gps(j); % 形函数对局部坐标的导数第1行是dN/dxi第2行是dN/deta dN 1/4 * [... -(1-eta), 1-eta, 1eta, -(1eta); -(1-xi), -(1xi), 1xi, 1-xi ]; J dN * nodes_e; % 2x2 Jacobian矩阵 invJ J \ eye(2); % 逆矩阵 dN_dx invJ(1,1)*dN(1,:) invJ(1,2)*dN(2,:); dN_dy invJ(2,1)*dN(1,:) invJ(2,2)*dN(2,:); B zeros(3,8); B(1,1:2:end) dN_dx; % eps_x 对应的形函数偏导 B(2,2:2:end) dN_dy; % eps_y B(3,1:2:end) dN_dy; % gamma_xy B(3,2:2:end) dN_dx; k k B * D * B * det(J) * t; % 2x2高斯积分权重积为1 end end end这段函数的数学逻辑是先在每个高斯点上计算形函数对自然坐标的导数通过 Jacobian 矩阵逆变换得到形函数对物理坐标的导数然后组装成 B 矩阵代入积分核心B*D*B*det(J)。高斯积分的权重在这里省略是因为二维二点高斯积分的四个点权重恰好都是 1乘积也是 1对 Q4 单元刚度矩阵的积分来说已经足够精确。D 矩阵的选取要按问题类型来。薄板结构受面内载荷时通常用平面应力公式里假设厚度方向应力为零而坝体、管道这类沿厚度方向尺寸很大的结构适合用平面应变模型。同一个 E 和 nu两种模型给出的 D 矩阵不同刚度结果的差异可能超过百分之十所以选错模型是常见的设计失误。如果手头问题既不是标准薄板也不是无限长柱体就应该考虑直接做三维分析而不是硬套二维假设。3.3 全局组装三种实现方式的取舍3.3.1 教科书三重循环写法拿到单刚函数后组装全局矩阵的方式决定了程序能撑多大的网格。下面这种三重循环写法最贴近教材也最容易对照公式理解function K assemble_stiffness(nodes, elems, E, nu, t, plane_type) n_dof 2 * size(nodes, 1); K sparse(n_dof, n_dof); for e 1:size(elems, 1) idx elems(e, :); % 单元节点全局编号 k q4_stiffness(nodes(idx, :), E, nu, t, plane_type); dofs [2*idx-1; 2*idx]; dofs dofs(:); % 关键列优先展开 K(dofs, dofs) K(dofs, dofs) k; end end注意这一行dofs dofs(:)。MATLAB 的(:)操作符按列展开矩阵所以 2×4 的矩阵会先取出第 1 列的 1、2再取第 3、4……得到的顺序是 1、2、3、4、9、10、7、8正好与单刚矩阵的排列节奏对齐。如果误写成dofs [2*idx-1, 2*idx]或用reshape按行展开自由度顺序会变成 1、3、9、7、2、4、10、8K 矩阵的行列含义将与单刚错位求解结果必然错误——而 MATLAB 不会给出任何警告。3.3.2 向量化组装用 sparse 三数组合并重复索引三重循环在自由度几万量级时会明显变慢问题出在稀疏矩阵的增量写入上。每执行一次K(dofs,dofs) K(dofs,dofs) k都要修改矩阵的存储结构开销远超预期。更可靠的做法是先把所有单刚的非零项收集到三个长数组中最后一次性调用 sparse 完成累加。function K assemble_stiffness_vec(nodes, elems, E, nu, t, plane_type) n_dof 2 * size(nodes, 1); n_ele size(elems, 1); n_edof 8; I zeros(n_ele * n_edof^2, 1); J zeros(n_ele * n_edof^2, 1); V zeros(n_ele * n_edof^2, 1); idx_start 0; for e 1:n_ele idx elems(e, :); k q4_stiffness(nodes(idx, :), E, nu, t, plane_type); dofs [2*idx-1; 2*idx]; dofs dofs(:); [ii, jj] ndgrid(dofs, dofs); block idx_start1 : idx_startn_edof^2; I(block) ii(:); J(block) jj(:); V(block) k(:); idx_start idx_start n_edof^2; end K sparse(I, J, V, n_dof, n_dof); endsparse(I, J, V) 这个构造方法有一个重要行为当多组 (I, J) 索引相同时V 中的对应值会被自动求和。这正是组装要做的事情——同一个全局自由度位置上的单元贡献相加。三重循环里的累加操作被这段代码移交给了 sparse 内部实现C 语言层面的哈希和排序累加比 MATLAB 脚本层循环快一到两个数量级。这个写法的适用范围很广从二维平面问题到三维六面体问题都可以套用只要把 n_edof 改成对应单元的自由度数即可。3.4 组装索引顺序的坑从行展开和列展开说起前面自由度展开已经出现了一个非常典型的踩坑场景。只因为换了一种展开方式K 矩阵就可能在没人注意的情况下被错误拼装。这类问题排查起来很费力因为位移结果不是完全离谱——它看起来像是某个方向刚度变软或变硬特别容易让人去怀疑材料参数而不是程序逻辑。我的经验是只要写完一个组装循环立刻用一个已知解析解的极简模型比如受拉杆或悬臂梁做验证不要直接上大网格。另外单元节点的局部排列方向也会影响组装结果。如果某个单元的节点是顺时针给出的det(J) 会是负值单刚的积分结果随之改变。上方 q4_stiffness 函数的实现会自动把这种现象体现出来k 矩阵的某些项变成负值最终总刚矩阵失去正定性。遇到这种情况把该单元的行向量反转一下比如 [1 2 5 4] 改成 [1 4 5 2]即可修复。4. 边界条件施加与后处理悬臂梁算例跑通验证4.1 边界条件的三类处理方法置零法、置大数法与自由度分离全局刚度矩阵 K 在未施加边界条件之前是奇异的因为结构存在刚体位移模式。求解前必须把足够的自由度约束住否则K\F会得到数值上无意义的结果。常见的处理方式有三种。置零法是教科书里最常见的思路把约束自由度对应的行列清零、对角线置 1同时把右端力向量对应项清零置大数法罚函数法则是在对角线加上一个很大的数让该自由度近似固定第三种是自由度分离法直接把自由度和约束自由度拆开求解。自由度分离法在 MATLAB 中实现最干净也最推荐使用。它不修改 K 矩阵本身只构造自由部分的子矩阵 K_ff 和载荷向量 F_f。固定端位移为零时核心逻辑只有两行u(free) K(free,free) \ F(free)。如果存在非零的已知位移公式变为F_f F(free) - K(free, fixed) * u_fixed即在回代时把已知位移换算成等效载荷。这种方法既避免了罚函数中系数选择的主观性也不会像置零法那样破坏稀疏矩阵的存储模式。约束自由度不足时K(free,free) 仍然奇异。MATLAB 的\运算在这种情形下会给出带警告的结果位移解里会混入刚体位移成分。常见于漏掉了某个节点的 x 方向自由度而节点本身又处于自由状态。检查方法是计算rank(K(free,free))若小于 free 自由度数量说明约束不够需要回到网格里逐个检查固定节点。4.2 完整算例两个 Q4 单元组成的悬臂梁受拉把前面所有的函数串起来跑一个可以用材料力学解析解验证的算例是检验整条链路是否正确的关键一步。沿用 2×1 的矩形网格左端两个节点完全固定右端两个节点各施加 0.5 的水平拉力总拉力为 1。材料参数取 E1、nu0.3、厚度 t1采用平面应力模型。根据杆件拉伸公式右端伸长的理论值为 2.0。% 主脚本两个 Q4 单元右端受拉 nodes [0 0; 1 0; 2 0; 0 1; 1 1; 2 1]; elems [1 2 5 4; 2 3 6 5]; E 1; nu 0.3; t 1; K assemble_stiffness_vec(nodes, elems, E, nu, t, stress); F zeros(12, 1); F(5) 0.5; % 节点3的x方向 F(11) 0.5; % 节点6的x方向 fixed_dofs [1 2 7 8]; % 节点1和节点4完全固定 free_dofs setdiff(1:12, fixed_dofs); u zeros(12, 1); u(free_dofs) K(free_dofs, free_dofs) \ F(free_dofs); fprintf(右端平均伸长 %.6f\n, mean(u([5 11]))); fprintf(理论伸长 %.6f\n, 2.0);同样一段代码完全适用于更复杂的网格。如果需要做悬臂梁弯曲问题只需把右端力改成 y 方向并且可以考虑使用 MATLAB 进行梁的有限元网格划分与计算中常见的梁单元或平面单元混合建模。平面单元模拟纯弯曲时Q4 单元由于只有双线性位移场弯曲变形会被“锁住”一部分这是单元本身的固有缺陷可以通过加密网格来缓解。注意以上参数全部采用自洽单位制长度、力、弹性模量的单位相互匹配。实际工程中如果混用毫米、牛顿和兆帕得到的位移值可能偏差几个数量级而程序本身不会报错。4.3 后处理位移云图与应力恢复求解完成后第一步是查看结构变形。将节点位移叠加到原始坐标上就得到变形后的形状。下面这段代码用 patch 函数画出受拉后的位移云图颜色表示 x 方向位移大小figure(Color, w); disp_x u(1:2:end); disp_y u(2:2:end); patch(Faces, elems, Vertices, nodes 50 * [disp_x disp_y], ... FaceColor, interp, CData, disp_x, EdgeColor, k); axis equal; colorbar; colormap(jet); xlabel(x); ylabel(y); title(x方向位移云图);这里的放大系数 50 只是显示效果用的视觉放大位移解本身很小不放大根本看不出变形趋势。实际后处理时建议把放大系数设为最大位移的 1% 到 5%使网格形状变化肉眼可辨。薄壁圆筒有限元等轴对称问题做后处理时常用 1/4 模型加对称边界条件也就是在对称面上约束法向位移这一设置同样可以从 fixed_dofs 中体现出来。应力恢复则要回到单元内部。位移解求出的是节点值而应变和应力在单元内是位置的函数。通常取单元中心点ξ0η0作为应力输出位置把 q4_stiffness 里计算 B 矩阵的代码抽出来代入该点形函数导数然后执行sigma D * B * u_e其中 u_e 是单元四个节点的位移向量。需要注意Q4 单元的应力在相邻单元之间是不连续的两个相邻单元在同一节点上的应力值会有差异这属于正常现象工程上常对节点做平均或使用超收敛点插值优化结果。4.4 从教学代码到能用代码保存与扩展当一个 MATLAB 脚本能跑通上述悬臂梁算例后再往正式应用方向扩展时有几个值得提前规划的点。建议把组装函数、单刚函数和边界条件处理封装为独立 .m 文件主脚本只负责描述网格、材料与载荷。这样更换网格尺寸或材料参数时不需要改动核心计算逻辑。对于多工况分析把 K 矩阵只组装一次再用不同的 F 向量分别求解能节省大量重复计算时间。另一个实际工程中常见的需求是处理非结构化网格。上面的两个 Q4 单元结构规则、排列整齐但实际网格往往来自 CAD 软件的自动剖分单元形状可能是任意四边形甚至包含少量退化单元。此时需要预先检查每个单元的 det(J) 值将非正的单元标记出来重新剖分或手动调整节点顺序。如果网格质量太差与其强行在 MATLAB 里修补不如回到网格生成软件中细化网格。5. 验证与提速把 MATLAB 有限元组装调顺的三个小实验5.1 刚体位移检验K 矩阵每行之和应当为零在施加任何边界条件之前全局刚度矩阵有一个非常强的数学性质如果结构整体发生刚体平移应力和应变都为零节点力也为零。对应到矩阵语言上K 矩阵每一行所有元素之和必须为零每一列同理。这个检验可以在组装完成后立刻执行max(abs(sum(K, 2))) % 若结果约等于0组装通过刚体检验理论上结果应为零实际中由于浮点舍入会在 1e-12 量级。如果出现明显的非零行说明该自由度对应的单元贡献不平衡多半是单元连接矩阵写错或自由度编号方案混乱。这个实验的价值在于它不依赖任何已知解纯粹从 K 的内部结构就能发现问题适合在每更换一套网格后运行一次。注意这个性质只在 K 矩阵组装完成后、边界条件处理前成立施加约束后行和自然不再为零。5.2 单刚自检对称性与特征值判断单元顺序单元刚度矩阵本身必然是实对称矩阵且是半正定的既有刚体位移对应的零特征值又有非负的变形特征值。对第一个单元跑一次自查能快速定位节点顺序问题k q4_stiffness(nodes(elems(1,:),:), E, nu, t, stress); max(max(abs(k - k))) % 应接近0 min(eig(k)) % 应不小于0如果第二行输出一个绝对值较大的负数几乎可以确定该单元的节点是以顺时针排列的导致 det(J) 为负。修复方法是把该行节点的顺序倒过来例如将[1 2 5 4]改为[1 4 5 2]。同样的检查也可以推广到所有单元在组装函数里添加一个 det(J) 的检查循环遇到非正值立即报错。这个前置检查的成本极低但能省下大量排查位移异常的时间。5.3 稀疏重排用 symamd 给自由度重新排序自由度编号顺序直接影响方程组求解的填充量。同样的网格如果节点编号沿短边连续排列刚度矩阵的带宽就小如果编号顺序杂乱LU 分解的中间因子会膨胀好几倍。对称近似最小度重排是一个几乎零成本的加速手段perm symamd(K(free_dofs, free_dofs)); u_free zeros(length(free_dofs), 1); u_free(perm) K(free_dofs(perm), free_dofs(perm)) \ F(free_dofs(perm)); u zeros(12, 1); u(free_dofs) u_free;重排后的解向量需要按 perm 回填到原自由度数序中否则位移结果会张冠李戴。这个技巧在几千自由度以下几乎看不出效果但网格到几万自由度、带宽较大的情况下求解时间可能有数量级差异。把上述三项检查整合成一个小函数check_assembly(nodes, elems, K)每次修改网格或材料参数后先跑一遍能挡住相当大一部分组装错误而不是等到位移云图完全变形后才回头查矩阵索引。本文还有配套的精品资源点击获取 MATLAB固体有限元编程：从单元刚度矩阵到全局组装实战
简介面向固体力学有限元初学者的MATLAB源码包围绕单元矩阵组装这一核心环节展示从几何建模、网格离散、材料属性设置到全局刚度矩阵集成的完整流程适合力学、土木、机械方向学生与工程技术人员边读边练。压缩包共5个文件全部为.m脚本整体仅12KB轻量且便于逐行理解。脚本以多个一维、二维算例为线索常见几何对象包括矩形区域与圆柱体等读者可从局部单元矩阵出发逐步学习边界条件施加、载荷处理、线性方程组求解以及后处理绘图等关键步骤体会连续区域被离散为有限元后如何通过矩阵组装获得整体解。已有311人学习下载资源可作为课程作业、毕业设计或小型科研案例的改写模板帮助读者快速打通有限元理论与MATLAB编程实现之间的映射关系。1. 这个标题背后是一整套 MATLAB 固体有限元的起手式你搜到“FEM 有限元.zip”这类资源包时看到的几个关键词其实是同一件事的三个面向FEM 是方法MATLAB 是工具而“组装单元矩阵”才是整个求解流程里最容易被忽略、也最常出错的环节。网上流传的 matlab 有限元编程求解实例十有八九会在一维杆或二维平面问题上演一遍这个流程剖网格、写单刚、组总刚、加边界、求解、后处理。乍看每一步都有公式可真到自己动手时卡住的往往不是推导而是“这个数到底该加到刚度矩阵的哪一个位置”。这篇博文就沿着这条最朴素但完整的路线走一遍。目标是让一个只记得有限元基础概念的人也能用 MATLAB 把二维固体问题的求解闭环跑通从节点坐标和单元连接表出发依次完成单元刚度矩阵计算、全局组装、边界条件施加和应力后处理。适合两类人一类是刚学到刚度矩阵一章、想用代码把公式落到实处的学生另一类是工程中只想快速验证一个结构方案、没时间开大型有限元仿真软件的开发者。2. 固体有限元的数据结构与自由度布局组装前先想清楚编号2.1 弱形式与刚度矩阵的物理含义为什么“组装”是必经之路固体力学有限元的出发点是弹性体的平衡方程结合几何方程与材料本构关系后直接求解这类偏微分方程并不容易。有限元的做法是将求解区域划分成有限个单元在每个单元内假设位移的插值形式再通过虚位移原理得到一组代数方程。最终形成的就是经典的线性系统K 乘以所有节点的位移向量等于节点力向量。这里的 K 就是全局刚度矩阵它的来源是每个单元的单元刚度矩阵而把这些局部矩阵按自由度位置“叠加”进同一个全局矩阵的操作就是组装。组装不是单纯的数学技巧它背后是力平衡的物理要求。相邻两个单元共享的节点上来自两个单元的贡献必须相加而不是相互覆盖否则节点上的力就不平衡。换句话说全局刚度矩阵 K 的每一行代表的是某个自由度在所有相邻单元中的刚度贡献和。这也是为什么“组装单元矩阵”会成为这类 MATLAB 程序的关键词它直接决定了 K 是否描述了一个真实的结构整体。2.2 MATLAB 里用两个矩阵承载网格节点坐标与单元连接任何一个固体有限元程序不论一维还是三维都离不开两个最基本的数据结构节点坐标矩阵和单元连接矩阵。MATLAB 的惯例是让每个节点占一行行号就是节点的全局编号。以二维问题为例节点坐标矩阵是 n_nodes 行、2 列的数组单元连接矩阵的每一行存放一个单元的所有节点全局编号按逆时针或统一方向排列。这里用一个贯穿全文的网格例子一个 2×1 的矩形板由两个四节点平面单元Q4组成共 6 个节点。对应的 MATLAB 数据如下% nodes(i,:) 节点i的(x, y)坐标 nodes [ 0 0; 1 0; 2 0; 0 1; 1 1; 2 1 ]; % elems(i,:) 第i个单元的4个节点全局编号按逆时针给出 elems [ 1 2 5 4; 2 3 6 5 ];单元连接矩阵的节点顺序非常关键。第一个单元用节点 1、2、5、4形成左下角那个 1×1 正方形第二个单元用节点 2、3、6、5形成右下角那个正方形。如果写成 1 2 4 5得到的会是一个交叉扭曲的四边形Jacobian 行列式可能为负后面的积分计算就会出错。这类错误在 MATLAB 里通常不会报错只会表现为位移结果完全不合理。2.3 自由度编号二维固体每个节点有两个故事“组装”最容易被低估的一步是自由度编号。一维杆单元里每个节点只有一个位移自由度节点编号和自由度编号几乎可以混用这也让很多人在转向二维问题时吃了大亏二维固体的每个节点有 x 和 y 两个方向的位移所以节点 1 占全局自由度的第 1、2 列节点 2 占第 3、4 列依此类推。上面网格的自由度分布可以列出如下关系节点编号x 方向自由度y 方向自由度112234356478591061112当处理第一个单元节点 1、2、5、4时局部自由度 1 到 8 对应的全局自由度分别是 1、2、3、4、9、10、7、8。这里的顺序必须与单元刚度矩阵的行列顺序保持一致否则 K 矩阵装配出来之后载荷施加和位移提取都会对不上。节点位移排列方式没有统一规定但建议始终采用“按节点分组、节点内先 x 后 y”的方案。这样做的理由是便于阅读从解向量里挑出节点 i 的 x 位移只需取 u(2*i-1) 即可后处理和画云图时一步到位。下面是两种常见自由度布局的对比方案节点内自由度顺序全局向量的排列按节点分组[x, y][n1_x, n1_y, n2_x, n2_y, ...]按方向分组[所有 x, 所有 y][n1_x, n2_x, ..., n1_y, n2_y, ...]如果要在其他软件里交换数据按方向分组更方便某些求解器的内存布局但在 MATLAB 里按节点分组写起来最直白后续用 reshape 处理位移云图也更顺手。选定后不要中途切换否则组装函数和边界条件代码要全部重写。3. MATLAB 单元刚度矩阵与全局组装从一维杆热手到二维 Q4 实装3.1 一维杆单元先热手自由度映射与循环组装在跳进二维平面单元之前先用一维杆单元把“组装”的索引机制看明白。一根由三段节点、两个杆单元组成的结构每个单元的长度为 L弹性模量 E截面积 A。杆单元的单元刚度矩阵是 2×2 的% 一维杆单元组装示例单元自由度为节点编号本身 EA 1; L 1; elems [1 2; 2 3]; % 两个单元依次连接 K zeros(3, 3); for e 1:size(elems, 1) k EA / L * [1 -1; -1 1]; % 杆单元单刚 dofs elems(e, :); % 局部自由度映射到全局 K(dofs, dofs) K(dofs, dofs) k; end这段代码的要点在于 K(dofs, dofs) 用的是累加而不是赋值。如果写成K(dofs, dofs) k第二个单元会覆盖掉第一个单元在节点 2 上的贡献得到的总刚矩阵将是一个对角线矩阵物理上完全错误。一维杆问题里由于自由度编号恰好等于节点编号这个映射关系显得透明到几乎不存在但它正是“组装”的核心操作做一次局部坐标到全局坐标的索引查找然后叠加。3.2 二维 Q4 单元的单刚形函数、几何矩阵与高斯积分3.2.1 为什么用四节点等参单元二维连续体问题最常用的教学单元是四节点四边形单元 Q4。它的每个单元有 4 个节点、8 个自由度单元刚度矩阵是 8×8 的。相比三角形常应变单元 CSTQ4 能描述线性变化的应变场精度更高相比高阶六节点三角形或九节点四边形Q4 的公式推导和编码难度低得多。对初学者来说这是性价比最高的选择。Q4 单元在局部坐标系ξ, η下定义四个节点的局部坐标分别为 (-1,-1)、(1,-1)、(1,1)、(-1,1)。位移插值由双线性形函数完成而把局部坐标映射到物理坐标的桥梁是 Jacobian 矩阵。单元刚度矩阵的标准计算公式为k_e ∫ B^T D B t dΩ其中 D 是弹性矩阵B 是几何矩阵应变-位移矩阵t 是厚度积分数值上用高斯积分完成。3.2.2 可直接运行的 MATLAB 单刚函数下面这个函数接受单元四个节点的物理坐标输出 8×8 的单元刚度矩阵。这是整个有限元程序里复用率最高的一个模块function k q4_stiffness(nodes_e, E, nu, t, plane_type) % nodes_e : 4x2 矩阵第i行是第i个局部节点的(x, y)坐标 % plane_type: stress(平面应力) 或 strain(平面应变) if strcmp(plane_type, stress) D E/(1-nu^2) * [1 nu 0; nu 1 0; 0 0 (1-nu)/2]; else D E/((1nu)*(1-2*nu)) * [1-nu nu 0; nu 1-nu 0; 0 0 (1-2*nu)/2]; end gps [-1/sqrt(3), 1/sqrt(3)]; % 2点高斯积分坐标 k zeros(8,8); for i 1:2 for j 1:2 xi gps(i); eta gps(j); % 形函数对局部坐标的导数第1行是dN/dxi第2行是dN/deta dN 1/4 * [... -(1-eta), 1-eta, 1eta, -(1eta); -(1-xi), -(1xi), 1xi, 1-xi ]; J dN * nodes_e; % 2x2 Jacobian矩阵 invJ J \ eye(2); % 逆矩阵 dN_dx invJ(1,1)*dN(1,:) invJ(1,2)*dN(2,:); dN_dy invJ(2,1)*dN(1,:) invJ(2,2)*dN(2,:); B zeros(3,8); B(1,1:2:end) dN_dx; % eps_x 对应的形函数偏导 B(2,2:2:end) dN_dy; % eps_y B(3,1:2:end) dN_dy; % gamma_xy B(3,2:2:end) dN_dx; k k B * D * B * det(J) * t; % 2x2高斯积分权重积为1 end end end这段函数的数学逻辑是先在每个高斯点上计算形函数对自然坐标的导数通过 Jacobian 矩阵逆变换得到形函数对物理坐标的导数然后组装成 B 矩阵代入积分核心B*D*B*det(J)。高斯积分的权重在这里省略是因为二维二点高斯积分的四个点权重恰好都是 1乘积也是 1对 Q4 单元刚度矩阵的积分来说已经足够精确。D 矩阵的选取要按问题类型来。薄板结构受面内载荷时通常用平面应力公式里假设厚度方向应力为零而坝体、管道这类沿厚度方向尺寸很大的结构适合用平面应变模型。同一个 E 和 nu两种模型给出的 D 矩阵不同刚度结果的差异可能超过百分之十所以选错模型是常见的设计失误。如果手头问题既不是标准薄板也不是无限长柱体就应该考虑直接做三维分析而不是硬套二维假设。3.3 全局组装三种实现方式的取舍3.3.1 教科书三重循环写法拿到单刚函数后组装全局矩阵的方式决定了程序能撑多大的网格。下面这种三重循环写法最贴近教材也最容易对照公式理解function K assemble_stiffness(nodes, elems, E, nu, t, plane_type) n_dof 2 * size(nodes, 1); K sparse(n_dof, n_dof); for e 1:size(elems, 1) idx elems(e, :); % 单元节点全局编号 k q4_stiffness(nodes(idx, :), E, nu, t, plane_type); dofs [2*idx-1; 2*idx]; dofs dofs(:); % 关键列优先展开 K(dofs, dofs) K(dofs, dofs) k; end end注意这一行dofs dofs(:)。MATLAB 的(:)操作符按列展开矩阵所以 2×4 的矩阵会先取出第 1 列的 1、2再取第 3、4……得到的顺序是 1、2、3、4、9、10、7、8正好与单刚矩阵的排列节奏对齐。如果误写成dofs [2*idx-1, 2*idx]或用reshape按行展开自由度顺序会变成 1、3、9、7、2、4、10、8K 矩阵的行列含义将与单刚错位求解结果必然错误——而 MATLAB 不会给出任何警告。3.3.2 向量化组装用 sparse 三数组合并重复索引三重循环在自由度几万量级时会明显变慢问题出在稀疏矩阵的增量写入上。每执行一次K(dofs,dofs) K(dofs,dofs) k都要修改矩阵的存储结构开销远超预期。更可靠的做法是先把所有单刚的非零项收集到三个长数组中最后一次性调用 sparse 完成累加。function K assemble_stiffness_vec(nodes, elems, E, nu, t, plane_type) n_dof 2 * size(nodes, 1); n_ele size(elems, 1); n_edof 8; I zeros(n_ele * n_edof^2, 1); J zeros(n_ele * n_edof^2, 1); V zeros(n_ele * n_edof^2, 1); idx_start 0; for e 1:n_ele idx elems(e, :); k q4_stiffness(nodes(idx, :), E, nu, t, plane_type); dofs [2*idx-1; 2*idx]; dofs dofs(:); [ii, jj] ndgrid(dofs, dofs); block idx_start1 : idx_startn_edof^2; I(block) ii(:); J(block) jj(:); V(block) k(:); idx_start idx_start n_edof^2; end K sparse(I, J, V, n_dof, n_dof); endsparse(I, J, V) 这个构造方法有一个重要行为当多组 (I, J) 索引相同时V 中的对应值会被自动求和。这正是组装要做的事情——同一个全局自由度位置上的单元贡献相加。三重循环里的累加操作被这段代码移交给了 sparse 内部实现C 语言层面的哈希和排序累加比 MATLAB 脚本层循环快一到两个数量级。这个写法的适用范围很广从二维平面问题到三维六面体问题都可以套用只要把 n_edof 改成对应单元的自由度数即可。3.4 组装索引顺序的坑从行展开和列展开说起前面自由度展开已经出现了一个非常典型的踩坑场景。只因为换了一种展开方式K 矩阵就可能在没人注意的情况下被错误拼装。这类问题排查起来很费力因为位移结果不是完全离谱——它看起来像是某个方向刚度变软或变硬特别容易让人去怀疑材料参数而不是程序逻辑。我的经验是只要写完一个组装循环立刻用一个已知解析解的极简模型比如受拉杆或悬臂梁做验证不要直接上大网格。另外单元节点的局部排列方向也会影响组装结果。如果某个单元的节点是顺时针给出的det(J) 会是负值单刚的积分结果随之改变。上方 q4_stiffness 函数的实现会自动把这种现象体现出来k 矩阵的某些项变成负值最终总刚矩阵失去正定性。遇到这种情况把该单元的行向量反转一下比如 [1 2 5 4] 改成 [1 4 5 2]即可修复。4. 边界条件施加与后处理悬臂梁算例跑通验证4.1 边界条件的三类处理方法置零法、置大数法与自由度分离全局刚度矩阵 K 在未施加边界条件之前是奇异的因为结构存在刚体位移模式。求解前必须把足够的自由度约束住否则K\F会得到数值上无意义的结果。常见的处理方式有三种。置零法是教科书里最常见的思路把约束自由度对应的行列清零、对角线置 1同时把右端力向量对应项清零置大数法罚函数法则是在对角线加上一个很大的数让该自由度近似固定第三种是自由度分离法直接把自由度和约束自由度拆开求解。自由度分离法在 MATLAB 中实现最干净也最推荐使用。它不修改 K 矩阵本身只构造自由部分的子矩阵 K_ff 和载荷向量 F_f。固定端位移为零时核心逻辑只有两行u(free) K(free,free) \ F(free)。如果存在非零的已知位移公式变为F_f F(free) - K(free, fixed) * u_fixed即在回代时把已知位移换算成等效载荷。这种方法既避免了罚函数中系数选择的主观性也不会像置零法那样破坏稀疏矩阵的存储模式。约束自由度不足时K(free,free) 仍然奇异。MATLAB 的\运算在这种情形下会给出带警告的结果位移解里会混入刚体位移成分。常见于漏掉了某个节点的 x 方向自由度而节点本身又处于自由状态。检查方法是计算rank(K(free,free))若小于 free 自由度数量说明约束不够需要回到网格里逐个检查固定节点。4.2 完整算例两个 Q4 单元组成的悬臂梁受拉把前面所有的函数串起来跑一个可以用材料力学解析解验证的算例是检验整条链路是否正确的关键一步。沿用 2×1 的矩形网格左端两个节点完全固定右端两个节点各施加 0.5 的水平拉力总拉力为 1。材料参数取 E1、nu0.3、厚度 t1采用平面应力模型。根据杆件拉伸公式右端伸长的理论值为 2.0。% 主脚本两个 Q4 单元右端受拉 nodes [0 0; 1 0; 2 0; 0 1; 1 1; 2 1]; elems [1 2 5 4; 2 3 6 5]; E 1; nu 0.3; t 1; K assemble_stiffness_vec(nodes, elems, E, nu, t, stress); F zeros(12, 1); F(5) 0.5; % 节点3的x方向 F(11) 0.5; % 节点6的x方向 fixed_dofs [1 2 7 8]; % 节点1和节点4完全固定 free_dofs setdiff(1:12, fixed_dofs); u zeros(12, 1); u(free_dofs) K(free_dofs, free_dofs) \ F(free_dofs); fprintf(右端平均伸长 %.6f\n, mean(u([5 11]))); fprintf(理论伸长 %.6f\n, 2.0);同样一段代码完全适用于更复杂的网格。如果需要做悬臂梁弯曲问题只需把右端力改成 y 方向并且可以考虑使用 MATLAB 进行梁的有限元网格划分与计算中常见的梁单元或平面单元混合建模。平面单元模拟纯弯曲时Q4 单元由于只有双线性位移场弯曲变形会被“锁住”一部分这是单元本身的固有缺陷可以通过加密网格来缓解。注意以上参数全部采用自洽单位制长度、力、弹性模量的单位相互匹配。实际工程中如果混用毫米、牛顿和兆帕得到的位移值可能偏差几个数量级而程序本身不会报错。4.3 后处理位移云图与应力恢复求解完成后第一步是查看结构变形。将节点位移叠加到原始坐标上就得到变形后的形状。下面这段代码用 patch 函数画出受拉后的位移云图颜色表示 x 方向位移大小figure(Color, w); disp_x u(1:2:end); disp_y u(2:2:end); patch(Faces, elems, Vertices, nodes 50 * [disp_x disp_y], ... FaceColor, interp, CData, disp_x, EdgeColor, k); axis equal; colorbar; colormap(jet); xlabel(x); ylabel(y); title(x方向位移云图);这里的放大系数 50 只是显示效果用的视觉放大位移解本身很小不放大根本看不出变形趋势。实际后处理时建议把放大系数设为最大位移的 1% 到 5%使网格形状变化肉眼可辨。薄壁圆筒有限元等轴对称问题做后处理时常用 1/4 模型加对称边界条件也就是在对称面上约束法向位移这一设置同样可以从 fixed_dofs 中体现出来。应力恢复则要回到单元内部。位移解求出的是节点值而应变和应力在单元内是位置的函数。通常取单元中心点ξ0η0作为应力输出位置把 q4_stiffness 里计算 B 矩阵的代码抽出来代入该点形函数导数然后执行sigma D * B * u_e其中 u_e 是单元四个节点的位移向量。需要注意Q4 单元的应力在相邻单元之间是不连续的两个相邻单元在同一节点上的应力值会有差异这属于正常现象工程上常对节点做平均或使用超收敛点插值优化结果。4.4 从教学代码到能用代码保存与扩展当一个 MATLAB 脚本能跑通上述悬臂梁算例后再往正式应用方向扩展时有几个值得提前规划的点。建议把组装函数、单刚函数和边界条件处理封装为独立 .m 文件主脚本只负责描述网格、材料与载荷。这样更换网格尺寸或材料参数时不需要改动核心计算逻辑。对于多工况分析把 K 矩阵只组装一次再用不同的 F 向量分别求解能节省大量重复计算时间。另一个实际工程中常见的需求是处理非结构化网格。上面的两个 Q4 单元结构规则、排列整齐但实际网格往往来自 CAD 软件的自动剖分单元形状可能是任意四边形甚至包含少量退化单元。此时需要预先检查每个单元的 det(J) 值将非正的单元标记出来重新剖分或手动调整节点顺序。如果网格质量太差与其强行在 MATLAB 里修补不如回到网格生成软件中细化网格。5. 验证与提速把 MATLAB 有限元组装调顺的三个小实验5.1 刚体位移检验K 矩阵每行之和应当为零在施加任何边界条件之前全局刚度矩阵有一个非常强的数学性质如果结构整体发生刚体平移应力和应变都为零节点力也为零。对应到矩阵语言上K 矩阵每一行所有元素之和必须为零每一列同理。这个检验可以在组装完成后立刻执行max(abs(sum(K, 2))) % 若结果约等于0组装通过刚体检验理论上结果应为零实际中由于浮点舍入会在 1e-12 量级。如果出现明显的非零行说明该自由度对应的单元贡献不平衡多半是单元连接矩阵写错或自由度编号方案混乱。这个实验的价值在于它不依赖任何已知解纯粹从 K 的内部结构就能发现问题适合在每更换一套网格后运行一次。注意这个性质只在 K 矩阵组装完成后、边界条件处理前成立施加约束后行和自然不再为零。5.2 单刚自检对称性与特征值判断单元顺序单元刚度矩阵本身必然是实对称矩阵且是半正定的既有刚体位移对应的零特征值又有非负的变形特征值。对第一个单元跑一次自查能快速定位节点顺序问题k q4_stiffness(nodes(elems(1,:),:), E, nu, t, stress); max(max(abs(k - k))) % 应接近0 min(eig(k)) % 应不小于0如果第二行输出一个绝对值较大的负数几乎可以确定该单元的节点是以顺时针排列的导致 det(J) 为负。修复方法是把该行节点的顺序倒过来例如将[1 2 5 4]改为[1 4 5 2]。同样的检查也可以推广到所有单元在组装函数里添加一个 det(J) 的检查循环遇到非正值立即报错。这个前置检查的成本极低但能省下大量排查位移异常的时间。5.3 稀疏重排用 symamd 给自由度重新排序自由度编号顺序直接影响方程组求解的填充量。同样的网格如果节点编号沿短边连续排列刚度矩阵的带宽就小如果编号顺序杂乱LU 分解的中间因子会膨胀好几倍。对称近似最小度重排是一个几乎零成本的加速手段perm symamd(K(free_dofs, free_dofs)); u_free zeros(length(free_dofs), 1); u_free(perm) K(free_dofs(perm), free_dofs(perm)) \ F(free_dofs(perm)); u zeros(12, 1); u(free_dofs) u_free;重排后的解向量需要按 perm 回填到原自由度数序中否则位移结果会张冠李戴。这个技巧在几千自由度以下几乎看不出效果但网格到几万自由度、带宽较大的情况下求解时间可能有数量级差异。把上述三项检查整合成一个小函数check_assembly(nodes, elems, K)每次修改网格或材料参数后先跑一遍能挡住相当大一部分组装错误而不是等到位移云图完全变形后才回头查矩阵索引。本文还有配套的精品资源点击获取