简介：这是一套以MATLAB为平台的杆系结构有限元分析程序包，面向结构工程学习者与MATLAB初级用户，解决从结构建模、材料属性定义、边界条件施加到静力求解与结果可视化的完整分析需求，也适合作为有限元方法课程的配套练习代码。压缩包内共359个文件，主体为.m脚本与函数文件，并包含大量.asv自动备份文件及1个.doc说明文档；.m用于核心计算，.doc辅助阅读，整体仅402KB，轻量便于下载后直接运行和核对。目前已有871人学习或下载。程序内含结构读取、荷载读取等辅助脚本，并附带多个Exam示例数据文件，便于针对桁架、框架等不同杆系结构开展静力分析；说明文档则有助于理清函数调用关系、数据输入格式与运行步骤。结合建模、网格划分、求解与后处理流程，读者可在有限元理论学习之余，通过实际代码掌握MATLAB在结构分析中的典型用法。
1. 杆系结构有限元分析：一个能用MATLAB独立跑通的小型计算引擎
如果你正在做结构力学课程设计或毕业设计，手边需要一个能“讲得清原理、改得动代码”的有限元分析工具，而不是把模型丢进商业软件里当黑匣子用，那么这个用MATLAB开发的杆系结构分析程序会非常对胃口。它把二维桁架、刚架这类杆件结构从输入到输出的完整计算流程，压缩在几个清晰的.m文件里：你只需要定义节点坐标、单元连接关系、材料和截面参数，程序会自动组装整体刚度矩阵、处理边界条件、求解节点位移和杆件内力，最后给出变形后的可视化图形。对想搞懂有限元底层逻辑的开发者来说，它更像一份“会说话的源码教材”。
2. 程序架构与数据流：从节点坐标到内力结果
2.1 文件构成与职责划分：每个 .m 文件干什么
先不急着跑，打开源码包先把文件结构摸清楚。这个程序不是把所有逻辑塞进一个巨大脚本里，而是按有限元分析的标准流程拆成了几个模块，每个模块只干一件事。我在拆这类学术型项目时，第一件事永远是看main函数在哪里，然后用调试器在脑海里把流程走一遍。
通常你会看到这几个文件：
| 文件名 | 职责 |
|---|---|
main.m | 主入口，控制整体流程：定义数据 → 组装 → 求解 → 输出 |
TrussElement.m或elem_stiff.m | 计算单根杆件的单元刚度矩阵 |
assemble.m | 把单元刚度矩阵组装进全局刚度矩阵 |
applyBC.m | 施加边界约束（位移边界条件） |
solve.m | 求解线性方程组，得到节点位移 |
postprocess.m | 后处理：从位移反算杆件轴力、应力，绘制变形图 |
如果你拿到手的源码包结构略有差异，比如把所有函数都写在main.m里，也不要觉得奇怪——很多教学代码为了减少文件跳转，会用局部函数的方式实现。但无论是分散还是集中，搞清楚“谁在调用谁”是第一步。
从职责上看，这个程序把有限元分析分成了三个清晰的阶段：前处理（输入数据组织）、求解核心（刚度矩阵组装与线性求解）、后处理（结果提取与可视化）。这个分段思维特别值得学习，因为你后期想把程序从静力分析扩展到模态分析或非线性计算时，只需要替换“求解核心”这一层，前后处理都可以复用。
2.2 核心数据结构：节点表、单元表、约束表、载荷表
数据结构决定了一个有限元程序能处理多复杂的问题。这套程序的数据组织方式是典型的“表驱动”风格——所有几何信息、物理信息都以矩阵行的形式存储，每行代表一个对象。这种设计在MATLAB里尤其自然，因为矩阵运算本来就是MATLAB的主场。
节点表nodes用[节点编号, x坐标, y坐标]存储，比如nodes(1,:) = [1, 0, 0]表示节点 1 在坐标系原点。单元表elements用[单元编号, 起始节点号, 终止节点号, 弹性模量E, 截面积A]存储，这是杆系分析里最通用的定义格式——用两个节点号描述一根杆件的拓扑关系，用 E 和 A 描述它的截面力学属性。
约束表和载荷表的设计也很直白。约束表记录[节点号, 自由度方向(1代表x平动, 2代表y平动), 约束位移值]，对固定支座来说约束位移值就是 0；载荷表则是[节点号, 自由度方向, 载荷数值]。这种表结构的好处是直观：你在纸上画一个简支梁，在哪个节点加了多大力，直接对应到载荷表的一行，几乎不需要额外脑内翻译。
有一个细节值得注意：自由度编号的规则。每个节点有 2 个自由度，那么节点 i 的 x 方向自由度编号是2*i-1，y 方向是2*i。整套程序里所有矩阵组装、约束施加都依赖这个映射关系，一旦编号错误，后期排查会很痛苦。源码里如果这部分的注释写得足够清楚，说明作者想得很明白。
2.3 求解主流程：六个步骤的完整链路
我把主程序的执行流程梳理成了一条明确的链路，你对照源码看会非常有方向感：
- 第一步：定义节点、单元、约束、载荷数据。这些数据可以写在主程序里，也可以从外部数据文件读取。读取方式在第四章专门讲。
- 第二步：初始化全局刚度矩阵 K。总自由度数为
2 * 节点总数，矩阵大小就是2n × 2n，先全部置零。 - 第三步：循环每个单元，计算单元刚度矩阵 ke，再把 ke 按照自由度映射关系叠加进全局 K 中。这是整个程序的核心，第三章会展开讲。
- 第四步：根据约束表处理边界条件。常见做法有“划行划列法”和“置大数法”，这个程序具体用的哪种，源码一眼就能看出来。
- 第五步：求解线性方程组 K·u = F。节点位移 u 是未知向量，载荷 F 是按载荷表填充的已知向量。
- 第六步：后处理。从 u 中提取每个节点的位移，再回到每个单元计算杆件轴力、应力，最后调绘图函数把变形后的结构画出来。
3. 核心代码拆解：刚度矩阵从单元到整体的组装全过程
3.1 单元刚度矩阵：这条公式是整个程序的命根子
二维杆系单元的刚度矩阵公式在结构力学教材里很常见，它的计算逻辑就是：先用两节点坐标算出杆长和方向余弦，再代入弹性模量 E、截面积 A，用一条固定公式算出 4×4 的局部刚度矩阵。程序里通常把这段逻辑封装成一个独立函数，输入是两节点坐标、E 和 A，输出是 4×4 矩阵。下面是实际代码里最典型的实现：
function ke = elem_stiff(nodal_coords, E, A) % nodal_coords: [x1, y1; x2, y2] % 返回单元在全局坐标系下的刚度矩阵 dx = nodal_coords(2,1) - nodal_coords(1,1); dy = nodal_coords(2,2) - nodal_coords(1,2); L = sqrt(dx^2 + dy^2); c = dx / L; % 方向余弦 cos(theta) s = dy / L; % 方向余弦 sin(theta) % 局部坐标下的刚度矩阵(2D桁架单元) k_local = E * A / L * [ c^2, c*s, -c^2, -c*s; c*s, s^2, -c*s, -s^2; -c^2, -c*s, c^2, c*s; -c*s, -s^2, c*s, s^2]; ke = k_local; end这段代码的关键在于c和s的计算：dx是单元在 x 方向的投影长度，dy是 y 方向投影长度，L是单元的几何长度。方向余弦把单元从“局部坐标系”换算到“全局坐标系”，这样即使杆件是斜放的，刚度矩阵也能直接叠加到全局矩阵里。函数不长，但如果dx或dy出现符号错误，后面所有结果都会跟着错，这里值得反复核对。
3.2 全局组装：自由度编号与矩阵叠加
单个单元的刚度矩阵算出来后，接下来的活就是把一堆 4×4 的小矩阵拼成一个大矩阵。这个“拼”的过程在有限元里叫组装，核心依据是单元两端节点的全局自由度编号。我之前接手过一个类似程序，看到组装代码时差点被绕晕——它内部没有用任何稀疏矩阵技巧，就是一个纯存储的问题，但自由度对不齐就会出事。
% 假设 nodes 是 n x 2 的坐标矩阵, elements 是 m x 4 的单元表 ndof = 2 * size(nodes, 1); K = zeros(ndof, ndof); for e = 1:size(elements, 1) n1 = elements(e, 1); n2 = elements(e, 2); E = elements(e, 3); A = elements(e, 4); % 提取两个节点的坐标 coord = [nodes(n1,:); nodes(n2,:)]; ke = elem_stiff(coord, E, A); % 自由度映射: 节点 n1 对应 2*n1-1 和 2*n1 两个自由度 dof = [2*n1-1, 2*n1, 2*n2-1, 2*n2]; K(dof, dof) = K(dof, dof) + ke; end参数说明：ndof是总自由度数量，K是全局刚度矩阵。循环里先用单元的两个节点号定位坐标，算出单元刚度矩阵ke，再通过dof数组完成从局部到全局的索引映射。K(dof, dof) = K(dof, dof) + ke这行是精髓——刚度叠加是“加”的操作，因为杆件在节点处共享自由度。初学者最容易把=写成覆盖，导致后面组装的单元直接覆盖掉前面单元的刚度贡献，位移结果小到离谱。
3.3 边界条件处理与求解收敛性
边界条件的处理方式直接影响求解能否成功。如果所有约束都不施加，刚度矩阵是一个奇异矩阵，线性方程组的解不存在——这在工程直觉上很好理解：没有支座约束的结构一受力就会飞出去。程序里处理约束的方式很直观，这里展示最常用的“置换法”：
% 约束表: [节点号, 自由度方向(1或2), 位移值] % 将约束位移写入位移向量, 并把刚度矩阵对应行列强制为单位阵 for i = 1:size(constraints, 1) nid = constraints(i, 1); dof_dir = constraints(i, 2); dof_id = 2 * (nid - 1) + dof_dir; K(dof_id, :) = 0; K(:, dof_id) = 0; K(dof_id, dof_id) = 1; F(dof_id) = constraints(i, 3); % 约束位移值(通常为0) end % 用 MATLAB 内置求解器直接求解 u = K \ F;说明一下这一步在做什么：先把被约束的节点自由度那行那列全部清零，对角设置为 1，再把载荷向量里对应的位置改成约束位移值。这样方程K*u=F在这个自由度上会解出一个确定的已知位移。对于固定支座，约束位移值是 0，F对应位置也就是 0。
注意：如果程序里采用的是“置大数法”（把对角元素加一个极大值如
1e10），效果类似但数值稳定性不如置换法。我一般建议在代码里把置换法的注释写清楚，这样训练营也好向读者交代原理。
求解完成后，位移向量u就拿到了。接下来的后处理阶段，对于每一根单元，提取它的两个节点位移，代入几何方程算出轴应变，再用材料本构关系算出轴力。公式也简单：轴力 =E*A/L * (u_e2_x - u_e1_x) * c + (u_e2_y - u_e1_y) * s再乘以相关系数。不过要小心符号方向——拉力为正还是压力为正，不同教材习惯不同，程序输出的时候会有一个符号约定说明。
4. 运行与验证：从数据文件到可视化结果
4.1 环境准备与启动步骤
这个程序对MATLAB版本的要求不苛刻，从 R2016a 到最新版本我都跑过，核心代码几乎没有用到需要高版本才支持的语法特性。运行前你需要确认机器上装的是 64 位 MATLAB，然后把源码文件夹加入当前工作路径。在命令行里敲一句main就能看到流程结果——不要双击文件，main.m里如果调用子函数，双击运行有时候会因为路径问题找不到子文件。
运行后程序通常会在命令行打印几张表：节点位移表、单元内力表。如果程序带绘图功能，还会弹出一个 Figure 窗口，显示变形前后的结构对比图。
| 输出项 | 格式示例 | 含义 |
|---|---|---|
| 节点位移 | (0.0000, -0.0325) | 节点 x、y 方向的平移量，单位通常是米（毫米视输入而定） |
| 单元轴力 | N = 152.4 kN | 杆件轴力，正负号表示拉/压 |
| 单元应力 | σ = 152.4 MPa | 轴力除以截面积后得到的应力 |
值得注意的是，这个程序的物理量单位完全取决于输入：你用什么单位输入坐标、E、A 和载荷，结果就是什么单位体系。比如输入坐标用米、E 用 Pa、A 用平方米、载荷用牛顿，位移算出来是米，轴力是牛顿，应力是 Pa。很多同学跑完觉得结果跟教材对不上，十有八九是单位混用了。
4.2 输入数据怎么改：一个三杆桁架的完整示例
为了验证程序跑得是否正确，最稳妥的办法是拿一个手算过的简单结构试水。我习惯用三杆桁架做第一次测试：三个节点组成一个三角形，底部两个节点铰接，顶部受一个竖向载荷。这种结构手算一步能出结果，刚好拿来验程序。
% 简例: 三杆桁架 nodes = [ 1 0.0 0.0; % 左下节点 2 2.0 0.0; % 右下节点 3 1.0 1.5 % 顶部节点 ]; elements = [ 1 1 3 2.1e11 0.001; % 左斜杆 E=210GPa A=0.001m^2 2 2 3 2.1e11 0.001; % 右斜杆 3 1 2 2.1e11 0.001 % 底杆 ]; constraints = [ 1 1 0; % 节点1 x方向固定 1 2 0; % 节点1 y方向固定 2 2 0 % 节点2 y方向固定 ]; loads = [ 3 2 -100000 % 节点3施加向下的100kN ];改完这些矩阵，直接跑main，对比输出的位移和手算。对上了，说明程序的数据通路没问题；对不上，先进第五章的排查清单看一遍。这种“已知答案验证未知程序”的习惯值得保持，每次改完代码都拿着桁架简例跑一遍回归，比什么测试都好使。
4.3 结果可视化：把数字变成看得见的变形图
可视化是这个程序加分的地方——它能直接绘制原始结构和变形后的结构。绘图的实现逻辑有心的人可以学一下：先画原始结构，节点用连线画出杆件；然后把求解出的位移叠加到节点坐标上，再画一遍同样的拓扑结构。
% 后处理绘制: 实线为原始形状, 虚线为变形形状 hold on; for e = 1:size(elements, 1) n1 = elements(e, 1); n2 = elements(e, 2); % 原始结构 plot([nodes(n1,1), nodes(n2,1)], [nodes(n1,2), nodes(n2,2)], 'b-', 'LineWidth', 2); % 变形后结构(放大系数) scale = 50; % 位移放大系数, 让变形肉眼可见 coord_disp1 = nodes(n1,:) + scale * u(2*n1-1:2*n1)'; coord_disp2 = nodes(n2,:) + scale * u(2*n2-1:2*n2)'; plot([coord_disp1(1), coord_disp2(1)], [coord_disp1(2), coord_disp2(2)], 'r--'); end hold off;这里有个参数值得你注意：scale是缩放系数。因为结构位移通常只有毫米甚至微米量级，直接叠加到米为单位的坐标上看不出差别。放大 50 倍甚至 100 倍，变形形态才肉眼可见。程序默认值不一定贴合你的模型，跑出来发现变形图看不出变化，优先检查是不是放大系数太小。
5. 避坑与排查：从奇异矩阵到错误结果的检查清单
5.1 求解报错“矩阵接近奇异”或位移量级异常
现象：命令行提示Matrix is singular，或者位移解出来是Inf、NaN、1e16这种离谱值。
原因：结构约束不足，也就是平时说的存在“刚体位移”。比如你只约束了两个节点的 y 方向，结构在 x 方向还可以整车平移，刚度矩阵就是奇异的。这个最容易犯，尤其是平面结构忘了约束整体平动。
解决：先数一遍约束表里的自由度数量，一个平面结构至少要消除三个刚体位移——通常约束一个节点的 x 和 y，再约束另一个节点的 y。如果还不行，打开K矩阵用det(K)或rcond(K)看看是否接近零。
5.2 算出来的位移和手算差距很大但方向对
现象：数值在同一量级，但偏大或偏小百分之几十。
原因：多半是单元组装时写成了覆盖赋值。在K(dof, dof) = K(dof, dof) + ke这行里，如果把加法+漏掉了，组装第二个单元时第一个单元的贡献直接被清零，等效于结构只有最后一根杆件受力，位移当然会大得离谱。这类问题肉眼很难发现，因为代码语法完全没有错，只能一步步打断点看K矩阵的叠加过程。
解决：在循环里加一个临时检查点，第一次循环后nnz(K)非零元素数量应该增长；第二次循环后总刚度值应该在增加。要是非零数量不对，先确认自由度映射数组dof是不是算错了。
5.3 约束位置和载荷方向“看着对但结果反了”
现象：结构受力后应该下挠，结果变形图是向上鼓的。
原因：坐标轴定义的差异。有的程序用 y 轴向下为正，有的用 y 轴向上为正。教材里手算用的坐标系方向和程序的绘图方向不一致，看起来就是变形“反向”。另一个常见坑是载荷符号：你想施加向下的力，但在 y 轴向上的坐标系里应该是-100000，结果填成了正值。
解决：统一坐标轴约定，建议在代码头部的注释里写清“y 轴向上为正，重力载荷取负数”。检查载荷表里竖向力的符号，和教材例题逐一对上。值得注意的是，如果你把约束位移也写错成负值，可能会出现整体结构平移而不是变形，这个现象跟载荷符号错误要区分开。
5.4 单位换算造成数值悬殊
现象：手算得到的位移是 0.032mm，程序输出 3.2e-5，看起来差了好几个数量级。
原因：单位混合使用了。比如坐标用米，但材料弹性模量用了兆帕而不是帕，或者截面积用了平方厘米，整体数量级就乱套了。有限元程序里没有物理单位的概念，整个模型只是一个由数字组成的线性方程组，所以单位必须全程统一。
解决：建议统一使用一套国际单位制（米-牛顿-帕），并做一个简单的“单位换算表”和输入数据放一起。回看第四章里的三杆桁架例子，E 是2.1e11（帕），A 是0.001（平方米），载荷是100000（牛顿），输出就是标准的国际单位。改任何输入数据之前，先拿起计算器把数量级对齐。
6. 二次开发思路：从杆系走向通用有限元框架
这个程序最有价值的地方在于它的结构足够“薄”，每一层都可以单独替换。如果想把二维桁架扩展成二维刚架（梁单元），核心改动点在单元刚度矩阵：你需要把每个节点从 2 个自由度扩展为 3 个（加上转角自由度），单元矩阵从 4×4 变成 6×6。替换elem_stiff函数里的公式后，组装和求解流程几乎不用动——这得益于“数据驱动”的分层设计。
另一种扩展方向是载荷工况支持。现在程序只支持节点集中力，实际工程里还有均布载荷、温度载荷。常见做法是把均布载荷等效转换成节点力，转换公式在结构力学里有现成的“固端力”查表。你只需要在加载前写一个load_conversion.m函数，把均布力折算成节点等效载荷，填入载荷表即可，跟核心求解器完全解耦。
我给这个程序做“加螺栓”改造的时候，最大的感悟是：有限元程序的本质不是算法，而是数据结构。结构分析难点不在解方程，而在如何把物理世界里的约束、连接、载荷翻译成程序能理解的“表”。学会把问题抽象成表，你就掌握了这类程序的钥匙。
从那以后我每次拿到一个有限元源码包，都会强制自己先画一张数据流图，把每个矩阵的行列含义、每张表的字段含义标注清楚，再动代码。希望你这次拿到这个杆系分析程序，也能保持这个习惯——先把数据结构吃透，再谈运行和改造。希望帮到你。
本文还有配套的精品资源，点击获取
