MATLAB有限元编程实战：杆板组合结构与梯形薄壁板分析

简介：面向航空结构分析与有限元课程设计的MATLAB程序包，针对杆板（梯形板）薄壁结构静力求解，适合机械、航空航天、土木等专业本科生或研究生完成大作业、理解有限元编程实现。压缩包共5个文件，以m主程序源码、docx可运行备用程序与程序报告、txt输入数据、pdf输出结果说明组成，包体仅280KB，结构紧凑但内容闭环。已有859人学习下载，可作为独立完成同类题目的重要参考。程序涵盖有限元核心流程：几何建模、网格划分、单元与整体刚度矩阵组装、边界条件与载荷施加、线性方程组求解及后处理，报告中还包括误差分析与收敛性检验。通过该程序，读者可获得可运行的求解代码、完整报告模板、典型薄壁梯形板算例及输入输出数据，帮助快速复现结果并掌握将有限元理论转化为MATLAB实践的方法。 上学期有限元课的大作业里，有一道题是用 MATLAB 写程序，求解一个“杆板组合结构”，其中板还特意做成了梯形薄壁板。我拿到题的第一反应是头疼——前面上机题全是单一类型单元，这一下要把杆单元和板单元放到同一个整体刚度矩阵里组装，还要处理梯形边界，完全是新玩法。但做完之后回头看，这题出得是真值。它把单元刚度矩阵推导、坐标变换、整体组装、边界条件处理这些有限元最核心的步骤一口气全串了起来。

这篇文章完全按我实际写程序的顺序来整理，从“题目到底想考什么”到“MATLAB 代码怎么落地”，再到调试过程中撞过的墙，都会尽量讲清楚。目的是给你一份可以直接照着改、照着排错的参考，而不是那种贴了一大段程序却看不懂在干嘛的“答案”。

1. 先把这个作业拆开：到底要算一个什么结构

1.1 杆板组合结构是哪来的

所谓“杆板组合结构”，形象点说就是一块薄板加上若干根加强筋。梯形板主体承担面内载荷，比如拉伸、剪切，而杆件布置在边缘或者中间，用来分担轴向力、抑制局部变形。工程里这种构造很常见，汽车车门内板、飞机蒙皮加筋壁板，都是这个套路。大作业为了能把问题简化到课堂知识范围内，通常会把三维结构压成二维问题：板只考虑面内变形，厚度方向应力近似为零。

这个近似的前提就是“薄壁”。如果板的厚度远小于平面尺寸，且载荷只在板平面内作用，就可以当成平面应力问题处理。节点自由度变成两个：水平位移 ux 和垂直位移 uy。杆单元是平面杆单元，每个节点也是两个自由度，正好和板单元匹配，这是整个程序能组装起来的基石。

作业的典型形式一般是：结构一侧固定，另一侧施加分布力或集中力，要求计算节点位移场、单元应力，有的还会让校核某个危险点的等效应力。所以这道题真正考察的重点不是你背下来多少有限元公式，而是你知不知道“单元刚度矩阵”和“结构整体刚度矩阵”之间那层组装关系。

1.2 梯形板网格怎么想

梯形板的麻烦点不在于单元理论，而在于网格生成。它宽度沿高度线性变化，如果直接套用矩形单元，边界就会变成锯齿状，误差没法控制。

我当时跑了两个方案对比：一是全部划分成三角形单元，二是四边形等参元映射。方案一明显更划算。任意四边形都能切成三角形，3 节点三角形单元（CST）是恒定应变单元，公式简单、刚度矩阵好写，而且三角形边是直线，能严格贴合梯形斜边。方案二需要在自然坐标系里做雅可比变换，还得上高斯积分，编程量一下子大好几倍，大作业周期短，没必要在这个环节死磕。

我的做法是先把梯形“按行切开”。每一行的宽度随高度线性变化，上底窄、下底宽（反过来也行），每行再等分成若干段，形成一个一个的小四边形，然后把每个四边形沿对角线切成两个三角形。网格细一点之后，CST 的“恒定应变”缺陷会被摊薄，结果足够应付课程要求。

2. 两个单元的刚度矩阵，另一条主线

2.1 杆单元算起来最快

杆单元是一维单元，局部坐标系下刚度矩阵非常简单：

k_local = EA/L * [1 -1; -1 1]

其中 E 是弹性模量，A 是截面积，L 是杆长。但程序里的杆往往不水平，而是斜着布置的，所以必须把局部刚度矩阵变换到全局坐标系。

设杆的两个端点分别是 P1(x1,y1)、P2(x2,y2)，那么杆长 L、方向余弦 c、s 为：

L = sqrt((x2-x1)^2 + (y2-y1)^2); c = (x2-x1) / L; s = (y2-y1) / L;

变换后，平面内任意方向杆单元的全局刚度矩阵是：

k_global = EA/L * [c*c c*s -c*c -c*s; c*s s*s -c*s -s*s; -c*c -c*s c*c c*s; -c*s -s*s c*s s*s];

这个 4×4 矩阵对应两个节点的顺序是 [u1x; u1y; u2x; u2y]，正好和后面板单元的自由度顺序一致。这里最容易犯错的是 c 和 s 的符号，建议代码里直接由节点坐标计算，千万不要手填。

2.2 板单元就用 CST 三角形

3 节点三角形单元的每个节点有两个自由度，单元自由度总共 6 个，对应位移向量：

u_e = [u1x; u1y; u2x; u2y; u3x; u3y]

单元内位移场用形函数插值，形函数是坐标的线性函数。因此应变矩阵 B 是常数矩阵，这就是“常应变三角形”名字的由来。先算三角形面积：

detJ = (x2-x1)*(y3-y1) - (x3-x1)*(y2-y1); A = 0.5 * abs(detJ);

B 矩阵可以直接用节点坐标差写出来。定义每边的差分：

b1 = y2 - y3; c1 = x3 - x2; b2 = y3 - y1; c2 = x1 - x3; b3 = y1 - y2; c3 = x2 - x1;

那么：

B = 1/(2*A) * [b1 0 b2 0 b3 0; 0 c1 0 c2 0 c3; c1 b1 c2 b2 c3 b3];

平面应力问题下，弹性矩阵 D 为：

D = E/(1-nu^2) * [1 nu 0; nu 1 0; 0 0 (1-nu)/2];

如果题目给的是薄壁板，默认就是平面应力，千万不要随手用平面应变的 D 矩阵，两者差挺多。单元刚度矩阵是：

ke = t * A * B' * D * B;

这里 t 是板厚。由于 B、D 都是常数矩阵，积分直接退化成一次乘法。这就是为什么大作业推荐 CST——你不需要写任何数值积分函数，就能把单元刚度矩阵算出来。

2.3 总刚矩阵的组装就一两行循环

有限元的灵魂不在单个单元，而是在“整体组装”。组装前先规定自由度全局编号，我的习惯是：

节点 i 的水平自由度 -> 2*i - 1 节点 i 的垂直自由度 -> 2*i

假设某个单元的节点编号是 n1、n2、n3，那么它的单元自由度索引是：

edof = [2*n1-1 2*n1 2*n2-1 2*n2 2*n3-1 2*n3];

整体刚度矩阵初始化为全零，然后循环单元，把每个单元的 ke 按 edof 索引累加进去：

K(edof, edof) = K(edof, edof) + ke;

这段代码看起来简单，但信息量很大。为什么是加不是赋值？因为一个节点往往被多个单元共享，不同单元的贡献必须叠加。为什么只操作 K(edof, edof)？这是 Rabinowitz 组装法的核心，本质就是把单元自由度位置映射到全局自由度位置。

杆单元也是同样的逻辑，自由度索引同样是[2*n1-1 2*n1 2*n2-1 2*n2]，所以杆和板在组装层面毫无冲突。

3. MATLAB 程序怎么落地

3.1 网格生成函数：把梯形切开

我先说几何参数：上底半宽 0.2 m，下底半宽 0.4 m，高度 0.3 m，厚度 0.005 m。网格规模取 nx=4（水平方向 4 段）、ny=3（竖直方向 3 段），先用小网格验证逻辑，后面再加密。

% 梯形板网格生成 a_top = 0.2; % 上底半宽，单位 m b_bot = 0.4; % 下底半宽，单位 m H = 0.3; % 板高度 nx = 4; ny = 3; % 水平分段数、竖直分段数 coord = zeros((nx+1)*(ny+1), 2); conn = zeros(nx*ny*2, 3); % 每个四边形切两个三角形 knode = 0; % 生成节点坐标 for i = 1:ny+1 t = (i-1)/ny; y = t * H; halfW = a_top + (b_bot - a_top) * t; % 该行半宽线性插值 xs = linspace(-halfW, halfW, nx+1); for j = 1:nx+1 knode = knode + 1; coord(knode, :) = [xs(j), y]; end end % 生成单元连接矩阵，每个四边形切成两个三角形 kelem = 0; for i = 1:ny for j = 1:nx % 四边形四个顶点编号（逆时针） p1 = (i-1)*(nx+1) + j; p2 = p1 + 1; p3 = p1 + (nx+1) + 1; p4 = p1 + (nx+1); % 三角形1：p1 p2 p3 kelem = kelem + 1; conn(kelem, :) = [p1 p2 p3]; % 三角形2：p1 p3 p4 kelem = kelem + 1; conn(kelem, :) = [p1 p3 p4]; end end

这段代码里最容易被忽略的是halfW的线性插值。每一行的半宽都不同，底边是 0.4，顶边是 0.2，中间行是两者的线性混合。如果这里写死成同一个宽度，那就变成矩形网格，梯形斜边就丢了。

三角形节点顺序必须是逆时针，否则面积 A 算出来是负的，B 矩阵符号反转，刚度矩阵跟着全错。这段代码里 p1、p2、p3、p4 是逆时针取的，所以没问题。

3.2 组装、边界条件和求解

单元刚度函数按前面理论写，然后执行组装。

E = 210e9; % 钢弹性模量，Pa nu = 0.3; % 泊松比 t = 0.005; % 板厚，m A_bar = 1e-4; % 杆截面积，m^2 nNodes = size(coord, 1); K = zeros(2*nNodes, 2*nNodes); % 组装板单元 for e = 1:size(conn, 1) n1 = conn(e,1); n2 = conn(e,2); n3 = conn(e,3); ec = coord([n1 n2 n3], :); % 单元节点坐标 ke = triStiffness(E, nu, t, ec); edof = [2*n1-1 2*n1 2*n2-1 2*n2 2*n3-1 2*n3]; K(edof, edof) = K(edof, edof) + ke; end % 组装杆单元（例如上底边缘的加强筋） % 假设上底左节点编号 nodeA、右节点编号 nodeB nodeA = 1; nodeB = nx+1; Lbar = norm(coord(nodeB,:) - coord(nodeA,:)); kbar = barStiffness(E, A_bar, coord(nodeA,:), coord(nodeB,:)); edofBar = [2*nodeA-1 2*nodeA 2*nodeB-1 2*nodeB]; K(edofBar, edofBar) = K(edofBar, edofBar) + kbar;

边界条件处理我用的是“划行划列置 1 法”，简单直观。比如左侧一列节点固定，节点编号是第 1 列所有节点：

fixedNodes = 1 : nx+1 : nNodes; % 左侧列节点编号 fixedDofs = []; for n = fixedNodes fixedDofs = [fixedDofs, 2*n-1, 2*n]; end F = zeros(2*nNodes, 1); % 施加载荷：例如下底边受向右总力 1000 N，等效到节点 % 下底边节点：节点编号为第 ny+1 行所有节点 bottomNodes = (nNodes - nx) : nNodes; % 这行索引需按生成规则修正 % 均布载荷等效：两端取半，中间取整，这里先给集中力示例 F(2*bottomNodes(end) - 1) = 1000; % 示意，实际应等效分配 % 处理固定边界 K(fixedDofs, :) = 0; K(:, fixedDofs) = 0; K(fixedDofs, fixedDofs) = eye(length(fixedDofs)); F(fixedDofs) = 0; % 求解 U = K \ F;

注意，上面代码里 bottomNodes 那行注释没有展开写全，实际作业里要按你自己的节点编号顺序重新推导。我的建议是把载荷等效单独写一个函数，统一处理，别在求解段里临时手算。

求解用K \ F就够了，自由度规模不会太大。如果网格加密到几千自由度，建议把 K 改成稀疏矩阵，使用前面思路预分配spalloc(2*nNodes, 2*nNodes, 稀疏带宽估算)，不然全零矩阵配几千自由度会明显拖慢速度。

3.3 后处理看云图，别只盯着数据

算完位移，接下来要算单元应力，这是大作业评分重点之一。

stress = zeros(size(conn,1), 3); % 每行 [Sx Sy Txy] for e = 1:size(conn,1) edof = [2*conn(e,1)-1 2*conn(e,1) 2*conn(e,2)-1 2*conn(e,2) 2*conn(e,3)-1 2*conn(e,3)]; ue = U(edof); ec = coord(conn(e,:), :); B = computeB(ec); D = computeD(E, nu); stress(e,:) = (D * B * ue)'; end

CST 单元应力在单元内部是常数，云图看起来是一块一块的色块，这是正常现象。想看位移云图，用 triangulation 和 trisurf 最省事：

T = triangulation(conn, coord(:,1), coord(:,2)); Ux = U(1:2:end); Uy = U(2:2:end); subplot(1,2,1); trisurf(T, Ux); title('Ux'); subplot(1,2,2); trisurf(T, Uy); title('Uy');

提交前强烈建议做一步“残差验证”：把算出来的 U 代回去，计算F_check = K * U，检查在非约束自由度上 F_check 是否等于施加的 F。如果这里都对不上，说明组装或者边界条件一定有问题，不要着急往下推进。

4. 最容易翻车的几个地方

4.1 刚阵奇异，多半是约束没给够

Matrix is singular to working precision是有限元作业里最高频的报错。原因九成是结构存在刚体位移，约束不足，K 矩阵不可逆。

我第一次遇到这个问题时，找了半天错在哪儿，最后发现是固定自由度的编号集合算错了。节点编号不是从 1 开始的自然增长顺序吗？我写了个 fixedNodes = 1:nx+1，结果那是最上边一行，不是左边，实际左侧节点应该是1:nx+1:end这个模式。这种错误不仔细看根本发现不了。

排查手段很简单：算一下 K 的零特征值数量。完整未约束结构应该至少有 3 个零特征值（两个平动、一个转动），施加足够约束后，如果仍然有接近 0 的特征值，就说明约束方向没约束上。

4.2 单位制不统一，结果全是笑话

这个错最隐蔽，因为它不报错，只是结果完全不可信。大作业常用单位容易混：几何尺寸用 mm，弹性模量用 MPa，最后力却用了 N，算出来位移数量级怎么都不对。

我自己的教训是把 E 写成 2.1e5 MPa，几何却全用的 m，然后板厚 0.005 m 又没换算，结果位移凭空大了三个数量级。最后统一成国际单位制（米、牛、帕），把所有的量纲列一遍：

在程序开头写清楚所有常量的单位，算完之后再在绘图或输出时统一换算，这样就不用动不动回改数据。

4.3 单元退化，detJ 不能是负的

CST 单元理论上不会像高阶单元那样因为畸变而严重降低精度，但如果三角形太扁，面积很小，刚度矩阵条件数会变大，结果还是可能出问题。

更常见的是网格生成时节点顺序写反，面积算成负值，B 矩阵整体变号。检测方法是在组装前统一计算每个单元的面积并打印最小值：

areaAll = zeros(size(conn,1),1); for e = 1:size(conn,1) ec = coord(conn(e,:), :); detJ = (ec(2,1)-ec(1,1))*(ec(3,2)-ec(1,2)) - ... (ec(3,1)-ec(1,1))*(ec(2,2)-ec(1,2)); areaAll(e) = 0.5 * detJ; end if any(areaAll < 0) error('存在负面积单元，节点顺序有误'); end

如果网格太粗导致位移云图出现明显的块状跳跃，先把 nx、ny 加密，比如 8×6、12×8，同时看关心的位移值是否趋于稳定。如果加密后数值还在很大范围内波动，就要回去查单元刚度矩阵推导，别指望靠网格把错误掩盖掉。

4.4 载荷分配，别把力重复加

载荷加载是个不高深但特别容易错的环节。如果题目给的是均布力，不能直接选一个节点把总量怼上去，必须按等效原则分配。

以底边受均布拉力为例，把底边分成 nx 段，总力为 F_total。内部节点承担的等效拉力是相邻两段上均布力之和的一半，而底边两端节点只承担半段，所以是内点的一半。更小的网格意味着更接近真实分布，这个等效不是近似，而是有限元里“一致节点载荷”的最基本处理。

我见过不少同学把底边每个节点都加上 F_total，最后结构受力是实际值的 nx+1 倍以上。提交前最简单的自检是：把所有载荷分量求和，看总合力是否等于已知外力，如果不等，多半就是重复加载或遗漏半段载荷的问题。

做这个作业时我最大的体会是：不要一上来就加密网格，先用 4×3 的小网格把位移和关键应力手算一遍，哪怕只用“纯拉伸”这种粗糙解析解对照，确认思路没跑偏，再逐步加密看收敛性。后面我又加了这个思路：把每个单元的面积、连接节点坐标、局部刚度矩阵对称性逐个打印出来，虽然土，但出问题找起来特别快。有限元程序调试顺序永远是“先小模型、再大模型，先单独单元、再组合结构”，这套习惯，到我现在做更复杂的仿真时仍然受用。

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