用MATLAB手写空间桁架刚度法求解器：从原理到代码实现

简介：一套面向土木、机械与航空航天领域工程师及学生的MATLAB空间桁架计算源码包，基于结构力学方法实现空间桁架的静力分析，帮助用户理解节点坐标定义、杆件连接、材料属性赋值、荷载与约束处理，以及稀疏线性方程组的组装与求解。压缩包共9个文件，以7个.m脚本为主，另有2个.asv自动备份文件，整体仅4KB；主程序可一键调用全部子函数，完成单元刚度矩阵计算、总刚度矩阵装配、杆件内力与应力提取等任务，结构清晰，便于学习和二次改造。已有808人学习下载，适合MATLAB结构编程初学者参考。借助这套代码，读者可快速搭建自己的桁架分析框架，并可通过修改节点坐标、材料参数和边界条件，扩展到更复杂的工程实例。 第一次抱着“空间桁架”这个关键词找资料时，我看到的大半都是通用有限元软件的建模演示。我不否认图形界面有它的价值，可当我要连续比几十个杆件截面时，鼠标点选就成了最大的瓶颈。反正都要重复试算，不如把整个计算逻辑写成脚本。于是我放下软件，用 MATLAB 从零写了一个空间桁架刚度法求解器，整个过程比预想中顺得多。把它整理出来，希望能给正在做课程设计或想自建参数化计算工具的朋友一点参考。

“空间桁架”听起来高深，落到力学模型上其实就是：一系列只承受轴力的杆件，在三维空间里通过节点铰接成整体。它和平面桁架最大的区别，是每个节点多了一个自由度，杆件的方向不再受到二维平面的限制。如果你已经熟悉平面桁架的刚度法，那空间桁架不过是在原来的思路上加了一个维度。真正难的不是公式，而是把坐标、连接、约束、荷载组织成一套可靠的数据结构。下面我按自己实现的顺序，把整条链路拆开讲。

1. 从平面桁架到空间桁架：为什么我还是选择手动写一遍

1.1 空间桁架和平面桁架的计算差异

体育馆穹顶、大跨雨棚、通信塔架，这些结构都可以简化成空间桁架。从计算模型上看，平面桁架每个节点有两个平动自由度，单元刚度矩阵只和杆件与 x 轴的夹角有关；空间桁架每个节点有三个平动自由度，单元刚度矩阵要处理三个方向余弦。两者在逻辑上的关系，就像把一条直线方程推广到空间直线方程，方向从“一个角度”变成了“一个单位向量”。

很多人觉得手写空间桁架程序工作量吓人，其实恰好相反。空间杆单元是整个有限元里最简单的一种单元：每个单元只有两个节点，每个节点三个自由度，单元刚度矩阵本身甚至比平面梁单元还规整。真正麻烦的是后面组装全局矩阵时，需要反复和“哪个节点对应哪几个全局自由度”这种索引细节较劲。MATLAB 处理这种问题有天然优势：矩阵运算是语言底层能力，不需要引入额外库；调试时随便打印中间矩阵；画图也方便。用别的语言要写一大圈数据结构和遍历逻辑，在 MATLAB 里几段循环就能搞定。

1.2 手写刚度法程序到底能干什么

我给自己定的目标是：输入节点坐标、单元连接表、约束信息、外荷载和抗拉刚度 EA，程序输出节点位移、每根杆件的轴力、支座反力，并画出变形前后对比图。听起来像个简化版有限元软件，实际上核心代码加起来不到一百行。它的价值不在于取代 SAP2000、ANSYS 这类工具，而在于完全可控：想加温度荷载、初应变、自重，只需要在对应位置插入一段代码；想批量跑一百个截面方案，直接把截面参数放进循环就行。

如果你正处在“看得懂教材但不知道怎么落地”的阶段，这个例子比任何黑箱软件都更适合入门。因为你可以顺着代码一步步看到刚度矩阵是怎么组装的，荷载是怎么施加的，轴力又是怎么从位移结果里提取出来的。等这套逻辑跑通了，再回头看商业软件的内力云图，你会更有底气判断它到底对不对。

2. 计算前最关键的事：把坐标、连接表、约束和荷载写成四个变量

2.1 空间桁架模型的输入约定

我习惯先把结构和程序之间的“接口”定义清楚，再写任何算法。对空间桁架来说，最核心的输入就是四个变量：

• xyz：节点坐标矩阵，第一行表示节点1的坐标，三列分别是 x、y、z。
• elem：单元连接表，每一行表示一根杆件，例如[1 4]表示杆件从节点1连到节点4。
• clamped：被约束的自由度编号。自由度编号规则是：节点 i 的 x、y、z 分别编号为3*i-2、3*i-1、3*i。
• F：外荷载向量，长度等于总自由度数，单位与 EA 取一致。

这四个变量定义清楚后，后面的程序逻辑基本就不用动了。我见过很多初学者一上来就盯着单元刚度矩阵背公式，结果卡在模型输入上：节点坐标单位是毫米还是米没统一，导致结果差了一千倍；连接表方向写反，导致轴力符号全反。这些问题都不是算法问题，而是数据组织问题。

2.2 一个能随时复现的小算例

为了让后面每一步都有东西可验证，我构造了一个最简单的空间桁架：一个四面体，四个节点、六根杆件。三个底部节点完全固定，顶部节点作用一个向下的集中力。这个结构的自由度规模很小，适合拿来手算复核程序。

clear; close all; clc; % 节点坐标：节点1、2、3为底部三角，节点4为顶部 xyz = [ 0.0, 0.0, 0.0 3.0, 0.0, 0.0 1.5, 2.598076, 0.0 1.5, 0.866025, 2.5 ]; % 单元连接表：四面体的六条棱 elem = [ 1 2 1 3 2 3 1 4 2 4 3 4 ]; % 约束自由度：节点1、2、3的 x/y/z 全部固定 clamped = [1:3, 4:6, 7:9]; % 外荷载：节点4沿 z 负方向受 10000 N F = zeros(size(xyz,1)*3, 1); F(4*3) = -10000; % 抗拉刚度 EA，单位 N EA = 2e7;

节点编号本身没有对错之分，它只影响矩阵的稀疏程度和输出结果的查看顺序。但一个实用的建议是：先编约束节点，再编自由节点，这样clamped可以写成一整段连续编号，不容易漏。我这里把底部三个节点连续编号，就是出于这个考虑。

3. 单元刚度矩阵组装与整体方程组求解：核心代码

3.1 三维杆单元的整体刚度矩阵：不必绕行局部坐标

教材里推导三维杆单元时，通常先建立局部坐标系，再用坐标转换矩阵把它变到整体坐标系。但在程序里，我们其实可以直接得到整体坐标系下的单元刚度矩阵，因为杆件的方向余弦可以由两个端点的整体坐标直接算出来。

定义从节点 i 指向节点 j 的单位方向向量：

n = (xyz(j,:) - xyz(i,:)) / L

其中 L 是杆件长度。设这个单位向量为n = [nx, ny, nz]，则整体坐标系下的单元刚度矩阵可以写成：

Ke = EA/L * [ n'*n -n'*n; -n'*n n'*n ]

这里n'*n是一个 3×3 矩阵，它的元素是三个方向余弦的两两乘积。为什么可以直接这样写？因为桁架杆件只能承受轴力，而轴力的大小完全由杆件两端节点沿杆轴方向的相对位移决定。利用投影关系，这个矩阵正好把所有轴力贡献都包含进去了，不需要再额外地做一次旋转。这比先算局部矩阵再乘转换矩阵的方式少一次矩阵乘法，也少一个犯错的机会。

对应的 MATLAB 函数可以这样写：

function Ke = spaceBarKe(xyz, elem, e, EA) i = elem(e,1); j = elem(e,2); vec = xyz(j,:) - xyz(i,:); L = norm(vec); n = vec / L; Ke = EA / L * [n'*n -n'*n; -n'*n n'*n]; end

有一个细节需要注意：单元连接表的写入方向会影响方向向量n的正负。比如[1 4]和[4 1]两种写法，刚度矩阵本身完全一样，但后面提取轴力时的符号是相反的。这不是错误，只要提取轴力时保持同样的顺序即可。

3.2 组装整体刚度矩阵与求解位移

整体刚度矩阵的组装逻辑是固定的：对每个单元，算出局部 6×6 刚度矩阵，再把它叠加到全局矩阵对应的行和列上。对应关系就是“节点 i 的三个自由度”和“节点 j 的三个自由度”。

nd = size(xyz,1) * 3; K = zeros(nd, nd); for e = 1:size(elem,1) Ke = spaceBarKe(xyz, elem, e, EA); i = elem(e,1); j = elem(e,2); dofs = [3*i-2, 3*i-1, 3*i, 3*j-2, 3*j-1, 3*j]; K(dofs, dofs) = K(dofs, dofs) + Ke; end

组装完成后，用setdiff(1:nd, clamped)得到自由自由度编号，把整体方程分成两部分：未知位移放在free，已知位移放在clamped。对固定支座来说，约束位移默认是零，所以可以直接把约束自由度对应的行和列删掉求解：

free = setdiff(1:nd, clamped); u = zeros(nd, 1); u(free) = K(free, free) \ F(free);

这里必须提一句：求解时不要用inv(K(free,free)) * F(free)。MATLAB 的反斜杠会自动根据矩阵特性选择 LU 分解、Cholesky 分解等合适算法，数值稳定性和速度都比显式求逆好。对于空间桁架这种中小规模问题，经验做法是直接用反斜杠，等模型节点数上到几万之后，再引入稀疏矩阵sparse和迭代求解器优化不迟。

如果你执行到这里发现矩阵奇异，大概率是约束不足。可以打印一下rcond(K(free,free))，这个值如果接近1e-16，说明存在机构位移或零刚度模式，需要回去检查clamped是否漏掉了某些节点。

4. 位移求出后的重头戏：轴力提取与支座反力校核

4.1 轴力公式与正负号约定

节点位移只是中间结果，工程上真正关心的是每根杆件的轴力。轴力提取本质上是一个反向过程：取出杆件两端节点的位移，计算它们沿杆轴方向的相对伸长，再乘上EA/L。

N = zeros(size(elem,1), 1); for e = 1:size(elem,1) i = elem(e,1); j = elem(e,2); vec = xyz(j,:) - xyz(i,:); L = norm(vec); n = vec / L; deltaL = (u(3*j-2:3*j) - u(3*i-2:3*i)) * n'; N(e) = EA / L * deltaL; end

这里面有一个很容易踩的坑：deltaL是“节点 j 相对节点 i 沿杆轴方向的位移投影差”，所以如果杆件受压，节点 j 会沿着指向节点 i 的方向移动，deltaL为负，轴力为负。我在代码里统一约定：正轴力代表拉力，负轴力代表压力。检查结果时，先看少数几根受力明确的杆件，确认符号方向符合直觉，再批量查看全部结果。

以我前面的四面体算例为例，设置EA = 2e7 N、顶部荷载 10000 N 时，得到的结果大致是：

• 节点4竖向位移约为-7.5e-4 m，即向下约 0.75 mm。
• 杆件 1-4、2-4、3-4 的轴力约为-4055 N，均为受压。
• 底部三根杆件 1-2、1-3、2-3 的轴力为 0。

底部杆件轴力为零经常让初学者觉得奇怪：明明整个结构在受力，为什么下面的杆子不受力？其实这是因为底部三个节点被完全约束后，位移全被限制为零，杆件两端没有位移差，自然不产生内力。这不是程序 bug，而是这道静定题目在固定支座模型下的真实结果。你若把这个模型放到通用有限元软件里算，也会得到同样的现象。

4.2 支座反力与整体平衡校核

位移解出来之后，支座反力可以非常简单地拿到：把整体刚度矩阵乘上完整位移向量，再减去外荷载向量。

R = K * u - F; Reaction = R(clamped);

这个式子的物理含义很直白：K*u是所有节点上的等效内力，减去已经施加的外荷载，剩下的就是在被约束节点上需要由支座提供的力。理论上，在自由自由度位置上，R应该等于零，这是验算程序的第一个信号。如果发现自由自由度上的残余力不是零，说明位移求解环节有问题，最常见的就是约束编号和外荷载位置错位。

第二个验算信号是整体平衡。把所有支座反力和外荷载分别合成合力，两者大小应该相等、方向相反。我实际调程序时经常在这里发现错误：某个集中力漏输在了错误节点上，或者某根杆件的连接表写成了[2 1]和[1 2]混用，看起来内力大小没问题，但某根杆的符号就反了。这些错误在矩阵层面很难一眼看出来，平衡校核却能在几秒内给出异常提示。

5. 把结果画出来：一张图能顶半页计算书

5.1 变形图与放大倍数

只给一串数字，很难判断模型整体变形模式是否合理。MATLAB 的三维绘图可以很好解决这个问题，核心思路就是把节点坐标加上位移乘以一个放大系数，再用线把杆件连起来。

scale = 300; % 位移放大倍数 def_xyz = xyz + scale * reshape(u, 3, [])'; figure; hold on; axis equal; grid on; for e = 1:size(elem,1) i = elem(e,1); j = elem(e,2); plot3(xyz([i j],1), xyz([i j],2), xyz([i j],3), 'b', 'LineWidth', 1.2); plot3(def_xyz([i j],1), def_xyz([i j],2), def_xyz([i j],3), 'r--', 'LineWidth', 1.5); end

放大系数的选择没有统一标准。我的习惯是让最大位移在图上显示为结构跨度的十分之一到五分之一左右，这样既能看清变形趋势，又不至于变形夸张到和原结构重叠。像上面那个四面体，最大位移约 0.75 mm，跨度 3 m，放大 300 倍后大约 0.225 m，在图上已经有可辨识的偏离。如果你用通用有限元软件画后处理图，会发现软件给的默认变形放大系数也约在这个量级。

axis equal这一行容易被忽略，但对三维桁架特别重要。不写它的话，MATLAB 会自动拉伸坐标轴，一个正方体看起来会变成扁盒子，变形方向判断会受到严重干扰。加上axis equal后，三个轴的比例一致，才能直观看出结构真实的空间形状。

5.2 用线宽和颜色表达受力大小

线宽在三维桁架渲染里是一个很好用的视觉变量。把所有杆件按轴力绝对值归一化，再映射到线宽，受力大的杆件一眼就能分辨出来：

maxN = max(abs(N)); for e = 1:size(elem,1) lw = 1 + 3 * abs(N(e)) / maxN; plot3(xyz([elem(e,1) elem(e,2)],1), ... xyz([elem(e,1) elem(e,2)],2), ... xyz([elem(e,1) elem(e,2)],3), ... 'Color', [0 0.45 0.75], 'LineWidth', lw); end

如果你想进一步区分拉压，可以把拉杆和压杆分成两组，分别用不同颜色画出。这样结构内部哪根杆受拉、哪根杆受压、哪些杆处于零杆状态，一眼就能看明白，非常适合放进计算书或者方案汇报里。

5.3 用实时脚本做参数化试算

我后面越来越依赖 MATLAB 实时脚本，因为可以把“改变荷载数值”改成拖动滑块操作。在实时编辑器里插入一个数值滑块，然后让绘图代码引用这个滑块的值，点一下就能看到变形和内力变化。对于方案比选阶段，这种交互比每次改代码再运行快很多，尤其在和甲方沟通时，现场拖一个滑块演示不同荷载下的受力响应，比单纯放几张静态图更有说服力。

6. 空间桁架程序调试中的一张问题排查表

6.1 我在实际调试里遇到最多的几类问题

手写数值程序，一次跑通是运气，多数情况需要反复排查。下面这张表是我这几年调试空间桁架程序时反复用到的诊断清单，基本覆盖了最常见的问题来源：

最容易被忽视的问题是单位。很多人写程序时一会用米、一会用毫米，结果位移和轴力完全对不上。我的习惯是写死一套单位制：几何尺寸用米，力用牛顿，弹性模量用帕，面积用平方米。这样 EA 的单位就是牛顿，位移的单位就是米。所有输入数据在进入程序前统一换算，程序内部不做任何单位转换，能省掉大半调试时间。

6.2 程序稳定之后再往哪个方向扩展

空间桁架刚度法程序一旦跑通，后续扩展是顺水推舟的事。最自然的下一步是加入自重荷载：把每根杆件的质量平均分配到两端节点，形成等效节点力叠加到 F 向量上。再进一步，可以加入温度效应，把温度引起的初应变加入轴力提取公式，变成deltaL - alpha * deltaT * L。

我目前最常做的是把整个求解器包成一个函数，输入参数里多一个截面面积数组，然后用循环批量计算不同截面组合下的结构响应。这样在做杆件优选时，程序会自动输出每组截面的最大轴力、最大位移和总用钢量。对实际项目来说，这个能力比单纯算一次内力有用得多。

另外要提醒一句：这个程序基于线弹性小变形假设，只适用于节点位移远小于构件长度、且杆件不发生整体失稳的情况。如果结构位移明显变大，或者你正在分析受压杆件失稳临界荷载，就需要引入几何刚度矩阵，甚至直接上非线性求解器。MATLAB 手写程序的价值在于把线性问题的每个细节看得清清楚楚，遇到非线性情况时，你也能知道自己卡在哪个环节。

最后再分享一个小习惯。我在每个算例脚本的第一行都会用注释写明单位、EA 来源、荷载取值依据，并把模型草图对应的坐标写在旁边。三个月后再打开脚本，你绝不会记得当初为什么这么建模型，但注释会帮你把当时的设计意图完整还原。空间桁架这类问题的代码节奏一旦跑通，后面扩展到平面框架、空间刚架，都只是换单元类型的事，真正沉淀下来的，是你对刚度法全过程的理解。

本文还有配套的精品资源，点击获取