尧图网站建设 
简介一套基于MATLAB的平面桁架有限元分析代码包面向结构工程、力学及相关专业的本科与硕士阶段教研学习可用于平面桁架建模、单元刚度矩阵组装、整体求解与结果可视化等典型有限元课程设计环节。资源支持MATLAB 2014/2019a/2021a等版本压缩包共包含4个文件以两个.m脚本为主提供完整的有限元分析实现配有一份.txt说明文档便于理解程序结构和参数含义另含一张PNG运行结果图可直观对照输出。整套资源仅50KB轻量紧凑下载后即可快速运行。目前已有225人学习适合课程设计、毕业设计或有限元入门研习使用。脚本与结果图结合既能帮助读者掌握平面桁架有限元程序的组织方式也可作为二次开发或教学演示的基础模板。1. 平面桁架有限元分析值得用MATLAB手写一遍拿到一个名为 PlaneTrussfem.m 的脚本第一反应往往是改一改节点坐标、跑通看到位移图就收工。但真到了毕业设计或工程复核需要自定义支座、增加多个载荷工况时脚本稍加改动就出各种奇怪结果。平面桁架虽然单元简单却是理解有限元方法最小但完整的载体每个节点两个自由度单元刚度矩阵不过 4×4却能完整走通“离散化→单元刚度→整体组装→约束施加→求解位移→恢复内力”全流程对教学和快速验证都极有价值。这套基于 MATLAB 的平面桁架有限元程序适合两类人本科或硕士做结构有限元课程设计需要一个能反复拆解、逐行验证的算例工程师想快速验证任意拓扑平面桁架的静力响应又不愿为一个十几根杆的小问题打开 ANSYS Workbench。MATLAB 脚本最大的优势是可以批量生成拓扑、修改材料和载荷并嵌入后续优化循环。下面按我拆这类代码的习惯把单元刚度、组装、边界条件和后处理逐层展开。2. 桁架单元刚度矩阵与方向余弦的MATLAB实现2.1 单元自由度与局部坐标系下的刚度矩阵平面桁架单元是典型二力杆只传递轴向力不需要考虑弯矩和剪力。每个节点有两个位移分量分别是 x 方向和 y 方向因此单元共有 4 个自由度。为了在整体结构里统一组装必须约定自由度顺序。我习惯把第 n 个节点的 x 向位移编号为 2n-1y 向位移编号为 2n。这样单元两端 i、j 的位移向量顺序就是 (ui, vi, uj, vj)自由度编号对应关系如下节点号x 方向自由度y 方向自由度112234n2n-12n在只沿杆轴方向的局部坐标系中单元刚度矩阵非常简单只有轴向一项k_local (E * A / L) * [1 -1; -1 1];这个矩阵描述的是杆端两个轴向位移与两个轴向力之间的关系。如果直接把 k_local 展开成 4×4 再往全局矩阵里放结果必然错误因为局部坐标轴与整体坐标轴通常不一致。杆在整体坐标中既有水平投影又有竖直投影必须先把局部刚度变换到整体坐标系。2.2 方向余弦与整体坐标下的 4×4 刚度矩阵设杆件从 i 节点指向 j 节点整体坐标差为 dx xj - xidy yj - yi杆长 L sqrt(dx^2 dy^2)方向余弦为 cx dx / Lcy dy / L。整体坐标系下单元刚度矩阵的解析式为k_e (EA/L) * [ cxcx cxcy -cxcx -cxcy; cxcy cycy -cxcy -cycy; -cxcx -cxcy cxcx cxcy; -cxcy -cycy cxcy cycy ]这个式子本质上就是 k_e T * k_local * T 的展开结果其中 T 是坐标变换矩阵。用 MATLAB 封装成函数可以让主程序循环更干净function ke truss2d_stiffness(E, A, xi, yi, xj, yj) % 计算平面桁架单元在整体坐标下的刚度矩阵 % 输入 % E - 弹性模量单位 Pa % A - 截面积单位 m^2 % xi, yi - i 节点整体坐标单位 m % xj, yj - j 节点整体坐标单位 m % 输出 % ke - 4x4 单元刚度矩阵自由度顺序为 (ui, vi, uj, vj) L sqrt((xj - xi)^2 (yj - yi)^2); % 杆长 cx (xj - xi) / L; % 整体坐标下的方向余弦 cy (yj - yi) / L; ke (E * A / L) * [cx*cx cx*cy -cx*cx -cx*cy; ... cx*cy cy*cy -cx*cy -cy*cy; ... -cx*cx -cx*cy cx*cx cx*cy; ... -cx*cy -cy*cy cx*cy cy*cy]; end这里 E 和 A 作为参数传入是为了后续做材料灵敏度分析或拓扑优化时只需在循环中修改数值不必改动单元函数。还要提醒一点L 不能为 0也就是两个节点不能重合否则会出现除零。实际模型里两根杆交叉但不共享节点或者节点坐标误输入为同一位置都会触发这个问题。2.3 全局刚度矩阵的组装与对号入座总刚矩阵的大小是 2N×2NN 为节点总数。组装时先初始化全零矩阵再逐单元取出两端节点编号计算单元刚度矩阵最后按自由度编号叠加到总刚的对应位置K zeros(2 * size(nodes, 1)); % 全局刚度矩阵初始化 for e 1:size(elements, 1) i elements(e, 1); j elements(e, 2); ke truss2d_stiffness(E, A, nodes(i,1), nodes(i,2), nodes(j,1), nodes(j,2)); dof [2*i-1 2*i 2*j-1 2*j]; % 单元四个自由度在全局中的编号 K(dof, dof) K(dof, dof) ke; % 叠加到全局矩阵 enddof 数组的排列顺序必须与 ke 中的 (ui, vi, uj, vj) 严格一致这是组装中最容易出错、也最不容易被察觉的地方。如果 dof 写成 [2i 2i-1 2j 2j-1]位移顺序发生变化K 矩阵就会错位。调试时可以用最简单的两节点单杆模型验证此时总刚 K 只有 4×4 非零展开后应当与 ke 完全一样。组装完成后K 一定是对称矩阵。用 issymmetric(K) 检查一下能快速发现写入顺序错误。对于节点数较多的模型我更倾向于在组装前先生成一个 N×2 的自由度编号表后续求解和后处理都基于这张表取值避免到处出现 2*i-1 这种魔法表达式。3. 桁架边界条件处理与静力凝聚求解的两种写法3.1 固定自由度与自由自由度的划分桁架结构在没有施加边界条件时刚度矩阵 K 是奇异的因为结构存在刚体位移模式。平面刚体位移包括两个平动和一个转动所以至少需要约束 3 个独立的自由度而且约束方向不能共线否则仍然可能绕汇交点转动。常见做法是用两个铰支座加一个滚动支座形成稳定的简支支承。MATLAB 里用 setdiff 划分自由度和固定自由度非常简洁fixed_dofs [1 2 4]; % 示例节点1固定节点2约束y向 free_dofs setdiff(1:size(K,1), fixed_dofs); % 剩余自由度fixed_dofs 的选取要检查是否完全消除了刚体位移。判断方法是对 K 做 rank 计算如果秩为 2N 减去约束数说明约束数量足够。但更稳妥的方式是用 null(K, r) 查看零空间基零空间对应的是允许的刚体或机构位移模式如果零空间里存在非零位移模式说明约束不足。3.2 划行划列法直接求解自由子矩阵最常用也最稳定的边界条件处理方式是划行划列。把 K 和载荷向量 F 按固定、自由两类自由度重新分块只对自由部分求解Kff K(free_dofs, free_dofs); Ff F(free_dofs); U zeros(size(F)); U(free_dofs) Kff \ Ff; % 只求自由自由度位移这个方法数值特性好求解规模也小。对于约束位移不为零的支座沉降情况需要增加一项修正Uc zeros(size(F)); Uc(fixed_dofs) ... % 已知支座位移比如0.01 F_eff F(free_dofs) - K(free_dofs, fixed_dofs) * Uc(fixed_dofs); U(free_dofs) Kff \ F_eff; U(fixed_dofs) Uc(fixed_dofs);这一步常被忽略。如果直接照搬零位移的代码去算支座沉降得到的位移场是错的。工程上做基础不均匀沉降分析时固定自由度上的已知位移必须按上面方式进入等效载荷。3.3 惩罚法适合快速检查和自动生成模型惩罚法是最容易写出来的边界条件处理方法思路是在固定自由度对应的对角元上添加一个大数让该自由度位移趋近于零penalty 1e12 * max(K(:)); for d fixed_dofs K(d, d) K(d, d) penalty; F(d) penalty * 0; % 等效于约束该自由度 end U K \ F;如果约束自由度有已知非零位移 uc要把 F(d) 改成 penalty * uc否则会强制位移为零。惩罚法对于程序自动生成约束集的场景很方便不需要反复维护自由度索引但缺点是 penalty 的选取直接影响求解精度和条件数。取太小则约束不彻底取太大则数值误差会掩盖真实位移。两种方法的对比见下表方法实现难度支持非零约束数值稳定性适用场景划行划列法中需要额外处理好常规静力求解、精确分析惩罚法低修改 F(d) 即可受惩罚数影响快速检查、拓扑自动生成我通常在主程序里用注释开关切换这两种写法。先用惩罚法跑通流程再切换到划行划列法做最终计算并比较两种结果。如果两种求解得到的最大位移差超过 1e-6 量级说明载荷或约束定义可能存在矛盾。4. 桁架轴力应力后处理与算例验证流程4.1 轴力和应力的计算公式平面桁架杆件的内力恢复必须在后处理中计算不能直接从节点位移读数判断。每个单元两端的位移向量 d_e [ui vi uj vj]沿杆轴方向的轴向变形为两个端点位移在杆轴方向的投影差deltaL (uj - ui) * cx (vj - vi) * cy轴力 N (EA/L) * deltaL应力 sigma N / A。在 MATLAB 中实现如下for e 1:size(elements, 1) i elements(e, 1); j elements(e, 2); L norm(nodes(j, :) - nodes(i, :)); cx (nodes(j,1) - nodes(i,1)) / L; cy (nodes(j,2) - nodes(i,2)) / L; dei [U(2*i-1); U(2*i); U(2*j-1); U(2*j)]; deltaL (dei(3) - dei(1)) * cx (dei(4) - dei(2)) * cy; axial_force(e) E * A / L * deltaL; axial_stress(e) axial_force(e) / A; end轴向变形 deltaL 是两个端点在杆轴方向的相对位移反映的是杆件真实的伸长或缩短。轴力为正表示受拉为负表示受压。这样计算出来的轴力会自动满足杆端平衡因为单元刚度矩阵本身就建立在轴向变形基础上。4.2 三杆三角形桁架算例用一个经典的三杆三角形桁架验证整个求解流程。节点1 (0,0)、节点2 (4,0) 为固定铰支座节点3 (2,2) 受竖直向下 20kN 载荷。E210GPaA0.01m^2。这个算例的优势在于可以手算校核快速判断程序是否正确。nodes [0 0; 4 0; 2 2]; elements [1 2; 1 3; 2 3]; E 210e9; A 0.01; F zeros(6, 1); F(6) -20000; % 节点3 y方向载荷 fixed_dofs [1 2 3 4]; % 节点1、2全部固定 free_dofs setdiff(1:6, fixed_dofs); K assemble_truss(nodes, elements, E, A); U zeros(6,1); U(free_dofs) K(free_dofs, free_dofs) \ F(free_dofs);其中 assemble_truss 是第 2.3 节组装代码封装后的函数内容就是循环叠加单位刚度矩阵这里不再重复。运行后得到节点3的竖向位移约为 -2.69e-5 m水平位移为 0符合对称结构的定性判断。杆件轴力结果如下单元i 节点j 节点轴力/kN判断1120.00零力杆213-14.14压杆323-14.14压杆水平下弦杆为零力杆是因为结构对称且没有水平载荷竖直载荷完全由两根斜杆承担。两根斜杆各承担 14.14kN竖直分量之和正好是 20kN与载荷平衡。如果手动算出的节点位移与代码结果不一致先检查 cx、cy 的符号特别是斜杆从节点1到节点3和从节点2到节点3方向不同cy 符号相反。4.3 支座反力复验求解位移后可以用 K*U 恢复节点力进而得到支座反力R F - K * U;R 在自由自由度上的分量应接近 0在固定自由度上的分量就是支座反力。对三杆算例节点1竖直反力约 10kN节点2竖直反力约 10kN水平支座反力均为 0。如果反力总和不等于外载荷说明组装或边界条件处理有误。这一检查对于任意规模的桁架模型都适用而且实现成本极低。5. 桁架刚度矩阵奇异排查与位移验证的三个技巧再复杂的平面桁架程序最终问题大多集中在矩阵奇异和位移结果异常上。分享三个我经常用到的排查技巧可以直接加到主程序尾部。第一个技巧是检查刚度矩阵的秩和零空间。对未约束的 K 计算 null(K, r)能得到刚体位移和机构位移模式。如果零空间基向量数大于 3说明结构存在机构也就是某些杆件布置方式让结构产生了多余的自由移动能力。这种情况在静定桁架中尤其常见比如四杆矩形框架不设斜杆时就是一个典型的可动机构。如果只想看秩可以用 rank(K)对于完全约束后的 Kff其秩应等于 free_dofs 的长度。第二个技巧是用最小特征值判断近奇异。约束后的 Kff 理论上正定但如果条件数过大最小特征值会趋近于零。实际代码里我会写min_eig eigs(K(free_dofs, free_dofs), 1, smallestabs); if min_eig 1e-10 * max(K(:)) warning(刚度矩阵近奇异请检查约束和杆件拓扑); end这里用 eigs 求一个最小特征值速度快也够用。如果触发警告优先检查是否存在孤立节点比如某个节点没有和任何单元相连或者两个节点之间缺少斜杆形成局部机构。用 spy(K) 看非零分布也能直观发现孤立的自由度行。第三个技巧是用应变能守恒验证位移结果。线性静力问题的应变能等于外力功的一半表达式为 0.5 * U * K * U应与 0.5 * F * U 一致。在 MATLAB 里比较这两个值如果相对误差大于 1e-8说明求解或自由度映射有问题。这个方法对局部错误特别敏感因为任何一行的位移或力弄错都会直接反映到能量上。把这三个检查写成一个独立函数在每次求解后自动调用能让程序从“跑通”变成“可信”。平面桁架虽然简单但这类线性静力程序做好之后稍加扩展就能处理二维刚架、空间桁架甚至更多类型的杆系结构。本文还有配套的精品资源点击获取 MATLAB平面桁架有限元分析：从单元刚度到轴力后处理
简介一套基于MATLAB的平面桁架有限元分析代码包面向结构工程、力学及相关专业的本科与硕士阶段教研学习可用于平面桁架建模、单元刚度矩阵组装、整体求解与结果可视化等典型有限元课程设计环节。资源支持MATLAB 2014/2019a/2021a等版本压缩包共包含4个文件以两个.m脚本为主提供完整的有限元分析实现配有一份.txt说明文档便于理解程序结构和参数含义另含一张PNG运行结果图可直观对照输出。整套资源仅50KB轻量紧凑下载后即可快速运行。目前已有225人学习适合课程设计、毕业设计或有限元入门研习使用。脚本与结果图结合既能帮助读者掌握平面桁架有限元程序的组织方式也可作为二次开发或教学演示的基础模板。1. 平面桁架有限元分析值得用MATLAB手写一遍拿到一个名为 PlaneTrussfem.m 的脚本第一反应往往是改一改节点坐标、跑通看到位移图就收工。但真到了毕业设计或工程复核需要自定义支座、增加多个载荷工况时脚本稍加改动就出各种奇怪结果。平面桁架虽然单元简单却是理解有限元方法最小但完整的载体每个节点两个自由度单元刚度矩阵不过 4×4却能完整走通“离散化→单元刚度→整体组装→约束施加→求解位移→恢复内力”全流程对教学和快速验证都极有价值。这套基于 MATLAB 的平面桁架有限元程序适合两类人本科或硕士做结构有限元课程设计需要一个能反复拆解、逐行验证的算例工程师想快速验证任意拓扑平面桁架的静力响应又不愿为一个十几根杆的小问题打开 ANSYS Workbench。MATLAB 脚本最大的优势是可以批量生成拓扑、修改材料和载荷并嵌入后续优化循环。下面按我拆这类代码的习惯把单元刚度、组装、边界条件和后处理逐层展开。2. 桁架单元刚度矩阵与方向余弦的MATLAB实现2.1 单元自由度与局部坐标系下的刚度矩阵平面桁架单元是典型二力杆只传递轴向力不需要考虑弯矩和剪力。每个节点有两个位移分量分别是 x 方向和 y 方向因此单元共有 4 个自由度。为了在整体结构里统一组装必须约定自由度顺序。我习惯把第 n 个节点的 x 向位移编号为 2n-1y 向位移编号为 2n。这样单元两端 i、j 的位移向量顺序就是 (ui, vi, uj, vj)自由度编号对应关系如下节点号x 方向自由度y 方向自由度112234n2n-12n在只沿杆轴方向的局部坐标系中单元刚度矩阵非常简单只有轴向一项k_local (E * A / L) * [1 -1; -1 1];这个矩阵描述的是杆端两个轴向位移与两个轴向力之间的关系。如果直接把 k_local 展开成 4×4 再往全局矩阵里放结果必然错误因为局部坐标轴与整体坐标轴通常不一致。杆在整体坐标中既有水平投影又有竖直投影必须先把局部刚度变换到整体坐标系。2.2 方向余弦与整体坐标下的 4×4 刚度矩阵设杆件从 i 节点指向 j 节点整体坐标差为 dx xj - xidy yj - yi杆长 L sqrt(dx^2 dy^2)方向余弦为 cx dx / Lcy dy / L。整体坐标系下单元刚度矩阵的解析式为k_e (EA/L) * [ cxcx cxcy -cxcx -cxcy; cxcy cycy -cxcy -cycy; -cxcx -cxcy cxcx cxcy; -cxcy -cycy cxcy cycy ]这个式子本质上就是 k_e T * k_local * T 的展开结果其中 T 是坐标变换矩阵。用 MATLAB 封装成函数可以让主程序循环更干净function ke truss2d_stiffness(E, A, xi, yi, xj, yj) % 计算平面桁架单元在整体坐标下的刚度矩阵 % 输入 % E - 弹性模量单位 Pa % A - 截面积单位 m^2 % xi, yi - i 节点整体坐标单位 m % xj, yj - j 节点整体坐标单位 m % 输出 % ke - 4x4 单元刚度矩阵自由度顺序为 (ui, vi, uj, vj) L sqrt((xj - xi)^2 (yj - yi)^2); % 杆长 cx (xj - xi) / L; % 整体坐标下的方向余弦 cy (yj - yi) / L; ke (E * A / L) * [cx*cx cx*cy -cx*cx -cx*cy; ... cx*cy cy*cy -cx*cy -cy*cy; ... -cx*cx -cx*cy cx*cx cx*cy; ... -cx*cy -cy*cy cx*cy cy*cy]; end这里 E 和 A 作为参数传入是为了后续做材料灵敏度分析或拓扑优化时只需在循环中修改数值不必改动单元函数。还要提醒一点L 不能为 0也就是两个节点不能重合否则会出现除零。实际模型里两根杆交叉但不共享节点或者节点坐标误输入为同一位置都会触发这个问题。2.3 全局刚度矩阵的组装与对号入座总刚矩阵的大小是 2N×2NN 为节点总数。组装时先初始化全零矩阵再逐单元取出两端节点编号计算单元刚度矩阵最后按自由度编号叠加到总刚的对应位置K zeros(2 * size(nodes, 1)); % 全局刚度矩阵初始化 for e 1:size(elements, 1) i elements(e, 1); j elements(e, 2); ke truss2d_stiffness(E, A, nodes(i,1), nodes(i,2), nodes(j,1), nodes(j,2)); dof [2*i-1 2*i 2*j-1 2*j]; % 单元四个自由度在全局中的编号 K(dof, dof) K(dof, dof) ke; % 叠加到全局矩阵 enddof 数组的排列顺序必须与 ke 中的 (ui, vi, uj, vj) 严格一致这是组装中最容易出错、也最不容易被察觉的地方。如果 dof 写成 [2i 2i-1 2j 2j-1]位移顺序发生变化K 矩阵就会错位。调试时可以用最简单的两节点单杆模型验证此时总刚 K 只有 4×4 非零展开后应当与 ke 完全一样。组装完成后K 一定是对称矩阵。用 issymmetric(K) 检查一下能快速发现写入顺序错误。对于节点数较多的模型我更倾向于在组装前先生成一个 N×2 的自由度编号表后续求解和后处理都基于这张表取值避免到处出现 2*i-1 这种魔法表达式。3. 桁架边界条件处理与静力凝聚求解的两种写法3.1 固定自由度与自由自由度的划分桁架结构在没有施加边界条件时刚度矩阵 K 是奇异的因为结构存在刚体位移模式。平面刚体位移包括两个平动和一个转动所以至少需要约束 3 个独立的自由度而且约束方向不能共线否则仍然可能绕汇交点转动。常见做法是用两个铰支座加一个滚动支座形成稳定的简支支承。MATLAB 里用 setdiff 划分自由度和固定自由度非常简洁fixed_dofs [1 2 4]; % 示例节点1固定节点2约束y向 free_dofs setdiff(1:size(K,1), fixed_dofs); % 剩余自由度fixed_dofs 的选取要检查是否完全消除了刚体位移。判断方法是对 K 做 rank 计算如果秩为 2N 减去约束数说明约束数量足够。但更稳妥的方式是用 null(K, r) 查看零空间基零空间对应的是允许的刚体或机构位移模式如果零空间里存在非零位移模式说明约束不足。3.2 划行划列法直接求解自由子矩阵最常用也最稳定的边界条件处理方式是划行划列。把 K 和载荷向量 F 按固定、自由两类自由度重新分块只对自由部分求解Kff K(free_dofs, free_dofs); Ff F(free_dofs); U zeros(size(F)); U(free_dofs) Kff \ Ff; % 只求自由自由度位移这个方法数值特性好求解规模也小。对于约束位移不为零的支座沉降情况需要增加一项修正Uc zeros(size(F)); Uc(fixed_dofs) ... % 已知支座位移比如0.01 F_eff F(free_dofs) - K(free_dofs, fixed_dofs) * Uc(fixed_dofs); U(free_dofs) Kff \ F_eff; U(fixed_dofs) Uc(fixed_dofs);这一步常被忽略。如果直接照搬零位移的代码去算支座沉降得到的位移场是错的。工程上做基础不均匀沉降分析时固定自由度上的已知位移必须按上面方式进入等效载荷。3.3 惩罚法适合快速检查和自动生成模型惩罚法是最容易写出来的边界条件处理方法思路是在固定自由度对应的对角元上添加一个大数让该自由度位移趋近于零penalty 1e12 * max(K(:)); for d fixed_dofs K(d, d) K(d, d) penalty; F(d) penalty * 0; % 等效于约束该自由度 end U K \ F;如果约束自由度有已知非零位移 uc要把 F(d) 改成 penalty * uc否则会强制位移为零。惩罚法对于程序自动生成约束集的场景很方便不需要反复维护自由度索引但缺点是 penalty 的选取直接影响求解精度和条件数。取太小则约束不彻底取太大则数值误差会掩盖真实位移。两种方法的对比见下表方法实现难度支持非零约束数值稳定性适用场景划行划列法中需要额外处理好常规静力求解、精确分析惩罚法低修改 F(d) 即可受惩罚数影响快速检查、拓扑自动生成我通常在主程序里用注释开关切换这两种写法。先用惩罚法跑通流程再切换到划行划列法做最终计算并比较两种结果。如果两种求解得到的最大位移差超过 1e-6 量级说明载荷或约束定义可能存在矛盾。4. 桁架轴力应力后处理与算例验证流程4.1 轴力和应力的计算公式平面桁架杆件的内力恢复必须在后处理中计算不能直接从节点位移读数判断。每个单元两端的位移向量 d_e [ui vi uj vj]沿杆轴方向的轴向变形为两个端点位移在杆轴方向的投影差deltaL (uj - ui) * cx (vj - vi) * cy轴力 N (EA/L) * deltaL应力 sigma N / A。在 MATLAB 中实现如下for e 1:size(elements, 1) i elements(e, 1); j elements(e, 2); L norm(nodes(j, :) - nodes(i, :)); cx (nodes(j,1) - nodes(i,1)) / L; cy (nodes(j,2) - nodes(i,2)) / L; dei [U(2*i-1); U(2*i); U(2*j-1); U(2*j)]; deltaL (dei(3) - dei(1)) * cx (dei(4) - dei(2)) * cy; axial_force(e) E * A / L * deltaL; axial_stress(e) axial_force(e) / A; end轴向变形 deltaL 是两个端点在杆轴方向的相对位移反映的是杆件真实的伸长或缩短。轴力为正表示受拉为负表示受压。这样计算出来的轴力会自动满足杆端平衡因为单元刚度矩阵本身就建立在轴向变形基础上。4.2 三杆三角形桁架算例用一个经典的三杆三角形桁架验证整个求解流程。节点1 (0,0)、节点2 (4,0) 为固定铰支座节点3 (2,2) 受竖直向下 20kN 载荷。E210GPaA0.01m^2。这个算例的优势在于可以手算校核快速判断程序是否正确。nodes [0 0; 4 0; 2 2]; elements [1 2; 1 3; 2 3]; E 210e9; A 0.01; F zeros(6, 1); F(6) -20000; % 节点3 y方向载荷 fixed_dofs [1 2 3 4]; % 节点1、2全部固定 free_dofs setdiff(1:6, fixed_dofs); K assemble_truss(nodes, elements, E, A); U zeros(6,1); U(free_dofs) K(free_dofs, free_dofs) \ F(free_dofs);其中 assemble_truss 是第 2.3 节组装代码封装后的函数内容就是循环叠加单位刚度矩阵这里不再重复。运行后得到节点3的竖向位移约为 -2.69e-5 m水平位移为 0符合对称结构的定性判断。杆件轴力结果如下单元i 节点j 节点轴力/kN判断1120.00零力杆213-14.14压杆323-14.14压杆水平下弦杆为零力杆是因为结构对称且没有水平载荷竖直载荷完全由两根斜杆承担。两根斜杆各承担 14.14kN竖直分量之和正好是 20kN与载荷平衡。如果手动算出的节点位移与代码结果不一致先检查 cx、cy 的符号特别是斜杆从节点1到节点3和从节点2到节点3方向不同cy 符号相反。4.3 支座反力复验求解位移后可以用 K*U 恢复节点力进而得到支座反力R F - K * U;R 在自由自由度上的分量应接近 0在固定自由度上的分量就是支座反力。对三杆算例节点1竖直反力约 10kN节点2竖直反力约 10kN水平支座反力均为 0。如果反力总和不等于外载荷说明组装或边界条件处理有误。这一检查对于任意规模的桁架模型都适用而且实现成本极低。5. 桁架刚度矩阵奇异排查与位移验证的三个技巧再复杂的平面桁架程序最终问题大多集中在矩阵奇异和位移结果异常上。分享三个我经常用到的排查技巧可以直接加到主程序尾部。第一个技巧是检查刚度矩阵的秩和零空间。对未约束的 K 计算 null(K, r)能得到刚体位移和机构位移模式。如果零空间基向量数大于 3说明结构存在机构也就是某些杆件布置方式让结构产生了多余的自由移动能力。这种情况在静定桁架中尤其常见比如四杆矩形框架不设斜杆时就是一个典型的可动机构。如果只想看秩可以用 rank(K)对于完全约束后的 Kff其秩应等于 free_dofs 的长度。第二个技巧是用最小特征值判断近奇异。约束后的 Kff 理论上正定但如果条件数过大最小特征值会趋近于零。实际代码里我会写min_eig eigs(K(free_dofs, free_dofs), 1, smallestabs); if min_eig 1e-10 * max(K(:)) warning(刚度矩阵近奇异请检查约束和杆件拓扑); end这里用 eigs 求一个最小特征值速度快也够用。如果触发警告优先检查是否存在孤立节点比如某个节点没有和任何单元相连或者两个节点之间缺少斜杆形成局部机构。用 spy(K) 看非零分布也能直观发现孤立的自由度行。第三个技巧是用应变能守恒验证位移结果。线性静力问题的应变能等于外力功的一半表达式为 0.5 * U * K * U应与 0.5 * F * U 一致。在 MATLAB 里比较这两个值如果相对误差大于 1e-8说明求解或自由度映射有问题。这个方法对局部错误特别敏感因为任何一行的位移或力弄错都会直接反映到能量上。把这三个检查写成一个独立函数在每次求解后自动调用能让程序从“跑通”变成“可信”。平面桁架虽然简单但这类线性静力程序做好之后稍加扩展就能处理二维刚架、空间桁架甚至更多类型的杆系结构。本文还有配套的精品资源点击获取