简介面向需要开展有限元编程实践的高校学生与工程师这套MATLAB程序包以二维与三维典型单元为主线覆盖梁单元、三角形单元、四边形单元、四面体单元和六面体单元并给出三梁平面框架、矩形薄板、空间块体等典型算例构成由平面到空间的完整单元体系。压缩包共56个文件以41个.m源程序文件为核心便于直接运行和二次开发另有doc算例说明、dat数据输入文件及txt文本说明整体仅333KB轻量但体系完整。已有1249人学习下载。通过学习可掌握Beam2D2Node、Quad2D4Node、Tetrahedron3D4Node、Hexahedral3D8Node等不同单元的有限元列式与MATLAB实现配套文档帮助对照典型例题理解建模、刚度矩阵组装、约束施加与边界条件处理。程序模块划分清晰算例命名从简单的一维杆单元到复杂的三维六面体单元层层递进便于按需查找并迁移到自己的研究场景。每一算例均包含前处理、求解与结果整理环节既有可运行的代码也有对应说明文档与代码一一对应可有效缩短调试时间适合作为课程设计或科研入门的参考模板。1. 从 quad2d4node 到六面体MATLAB 有限元算例的递进路线quad2d4node 这个名字是很多 MATLAB 有限元教程里四节点四边形单元的默认函数名。拆开看就是 quad四边形、2d二维、4 node四节点。标题把这样一个函数和六面体、算例并排实际上是在给出一条很典型的成长路径先在二维里写好一个能算的单元子程序验证坐标变换、高斯积分和刚度组装再把同一套模式推到三维八节点六面体最后落到一个能出数字的算例上。这条路不需要额外安装工具箱只用 MATLAB 原生矩阵运算就能完成。很多工程师读完了弹性力学和有限元原理到“自己写程序”这一步往往卡住。原因不是理论不懂而是不知道从哪一行开始写。quad2d4node 是一个很合适的起点因为它足够短却覆盖了等参单元的全部要素形函数、雅可比矩阵、B 矩阵、数值积分。后续把维度从 2 改成 3把单元自由度从 8 改成 24逻辑完全一致这也是把四边形和六面体放在一起讲的原因。这篇文章适合两类读者一类是准备在 MATLAB 里复现论文有限元结果的科研人员另一类是做轻量化自研前处理或教学演示的工程师。有五年以上经验的人可以直接跳到第三节和第五节那里有关键的索引陷阱和减缩积分处理。新手则建议按章节顺序走把第一个函数调通后再做六面体。2. quad2d4node 单元刚度矩阵等参变换与高斯积分2.1 四节点四边形的形函数与自由度顺序在自然坐标系下四节点四边形单元的四个节点固定在 (-1,-1)、(1,-1)、(1,1)、(-1,1)。单元内任意一点的位移场用双线性形函数插值N1(1-ξ)(1-η)/4N2(1ξ)(1-η)/4N3(1ξ)(1η)/4N4(1-ξ)(1η)/4。每个节点有 x 和 y 两个自由度所以一个单元的位移向量长度是 8。常见排列是 [u1;v1;u2;v2;u3;v3;u4;v4]。全局网格的节点编号可以任意但单元内部的自由度顺序必须固定否则 B 矩阵错位后连刚体位移检查都过不了。实现这段插值时通常用参考坐标矩阵配合循环派生形函数。这样后续改成六面体时只需要增加一个维度。核心代码是ref [-1 -1; 1 -1; 1 1; -1 1]; for nd 1:4 xi0 ref(nd,1); eta0 ref(nd,2); N(nd) (1xi*xi0)*(1eta*eta0)/4; dN(1,nd) xi0*(1eta*eta0)/4; dN(2,nd) eta0*(1xi*xi0)/4; enddN的第一行是对 ξ 的偏导第二行是对 η 的偏导。后面计算雅可比矩阵时这个排列顺序要严格对应坐标矩阵的行列否则雅可比行列式会出现符号翻转。很多早期程序跑出负面积问题就出在这里。2.2 平面应力与平面应变的本构矩阵差异D 矩阵决定应力-应变关系。对于各向同性弹性材料二维问题的 D 矩阵是 3×3第三行对应剪应变 γxy。下面列出两种分析模式的关键系数分量平面应力平面应变D11E/(1-ν²)E(1-ν)/((1ν)(1-2ν))D12Eν/(1-ν²)Eν/((1ν)(1-2ν))D33E/(2(1ν))E/(2(1ν))容易记错的是 D33平面应变并没有用 (1-2ν)/2 那项最终化简后仍然是剪切模量 G。我自己早期分析薄壁圆筒截面时忽略了这一点导致剪切刚度被放大后来拿单位圆筒静水压算例核对才发现。第一次写程序时建议把两种 D 矩阵都单独封装在注释里标明适用场景。2.3 quad2d4node 的完整实现与参数说明下面这个函数是整套程序的基石。输入四个节点坐标、弹性模量和泊松比返回 8×8 单元刚度矩阵。高斯积分采用 2×2 个点积分点坐标是 ±1/√3权重都是 1。function K quad2d4node(x, y, E, nu, type) % 四节点四边形等参单元刚度矩阵 % x,y: 4x1 坐标向量按逆时针顺序传入 % type: stress 平面应力strain 平面应变 if nargin 5, type stress; end if strcmp(type, stress) D E/(1-nu^2)*[1 nu 0; nu 1 0; 0 0 (1-nu)/2]; else f E/(1nu)/(1-2*nu); D f*[1-nu nu 0; nu 1-nu 0; 0 0 (1-2*nu)/2]; end gp [-1/sqrt(3), 1/sqrt(3)]; % 高斯积分点 w [1, 1]; K zeros(8,8); for i 1:2 for j 1:2 xi gp(i); eta gp(j); % 形函数对自然坐标的偏导2x4 dN 0.25*[-(1-eta) 1-eta 1eta -(1eta); -(1-xi) -(1xi) 1xi 1-xi]; J dN * [x y]; % 雅可比矩阵 detJ det(J); if detJ 0 error(单元面积非正请检查节点顺序或坐标); end dNxy J \ dN; % 形函数对总体坐标的偏导 B zeros(3,8); for nd 1:4 B(1,2*nd-1) dNxy(1,nd); B(2,2*nd) dNxy(2,nd); B(3,2*nd-1) dNxy(2,nd); B(3,2*nd) dNxy(1,nd); end K K w(i)*w(j) * (B*D*B) * detJ; end end end参数说明dN是 2×4 形函数偏导矩阵J是雅可比矩阵由dN乘以节点坐标得到dNxy通过左除求解形函数对总体坐标的导数。B*D*B是单元刚度积分项最后乘上雅可比行列式和积分点权重。这个函数可以直接复制到脚本里测试。调用示例Kq quad2d4node([0;10;10;0], [0;0;10;10], 200e3, 0.3, stress);这是 10mm 正方形单元的刚度矩阵。由于是完整单元刚度矩阵它的秩等于 8-35即三个刚体模式对应的零特征值。2.4 组装到全局刚度矩阵的索引方法实际网格有成百上千个单元不能拿着单元刚度矩阵逐个加。常见做法是先建稀疏矩阵再用自由度索引累加。下面的函数把 2D 组装过程封装起来function K assemble_quad2d(nodes, elems, E, nu, type) % nodes: N×2第 i 行是第 i 个节点的 x,y 坐标 % elems: M×4第 e 行是第 e 个单元的四个节点号 N size(nodes,1); K sparse(2*N, 2*N); for e 1:size(elems,1) nid elems(e,:); ke quad2d4node(nodes(nid,1), nodes(nid,2), E, nu, type); % 每个节点占据两个全局自由度 dof zeros(1,8); dof(1:2:end) 2*nid - 1; dof(2:2:end) 2*nid; K(dof,dof) K(dof,dof) ke; end end2*nid-1对应 x 自由度2*nid对应 y 自由度。当节点编号不连续时nodes的行号依然代表节点编号所以这里直接用节点编号索引。组装完成后再按边界条件划分自由度和固定自由度就可以求解位移了。3. 从 quad2d4node 到六面体三维八节点等参单元实现3.1 六面体的节点编号约定与形函数六面体单元是四边形在三维上的直接推广。自然坐标有三个分量 ξ、η、ζ每个分量都在 [-1,1] 上。8 个节点分布在 6 个面上常用编号顺序是底面四个节点按逆时针从 1 到 4顶面对应位置为 5 到 8。节点参考坐标可以存成一个 8×3 矩阵ref [-1 -1 -1; 1 -1 -1; 1 1 -1; -1 1 -1; -1 -1 1; 1 -1 1; 1 1 1; -1 1 1];形函数公式为 N_i (1ξξ_i)(1ηη_i)(1ζζ_i)/8。每个节点三个位移自由度单元应变向量有 6 个分量所以 B 矩阵是 6×24。和二维 quad2d4node 的区别只有维度从 2 变到 3从 4 节点变到 8 节点其余结构完全一致。3.2 六面体单元刚度矩阵的 MATLAB 实现下面这个函数生成 8 节点六面体单元刚度矩阵。高斯积分采用 2×2×2 八个点积分坐标同样是 ±1/√3。function K hex8node(xyz, E, nu) % xyz: 8×3 节点坐标矩阵按 ref 顺序传入 ref [-1 -1 -1; 1 -1 -1; 1 1 -1; -1 1 -1; -1 -1 1; 1 -1 1; 1 1 1; -1 1 1]; f E/(1nu)/(1-2*nu); D f*[1-nu nu nu 0 0 0; nu 1-nu nu 0 0 0; nu nu 1-nu 0 0 0; 0 0 0 (1-2*nu)/2 0 0; 0 0 0 0 (1-2*nu)/2 0; 0 0 0 0 0 (1-2*nu)/2]; gp [-1/sqrt(3); 1/sqrt(3)]; K zeros(24,24); for a 1:2 for b 1:2 for c 1:2 xi gp(a); eta gp(b); zeta gp(c); % 形函数对自然坐标的偏导3x8 dN zeros(3,8); for nd 1:8 xi0 ref(nd,1); eta0 ref(nd,2); zeta0 ref(nd,3); dN(1,nd) xi0*(1eta*eta0)*(1zeta*zeta0)/8; dN(2,nd) eta0*(1xi*xi0)*(1zeta*zeta0)/8; dN(3,nd) zeta0*(1xi*xi0)*(1eta*eta0)/8; end J dN * xyz; % 3x3 雅可比矩阵 detJ det(J); if detJ 0 error(六面体单元负体积检查节点顺序); end dNxy J \ dN; % 3x8 B zeros(6,24); for nd 1:8 col 3*nd - 2; % u,v,w 三列的起点 B(1, col) dNxy(1,nd); B(2, col1) dNxy(2,nd); B(3, col2) dNxy(3,nd); B(4, col) dNxy(2,nd); B(4, col1) dNxy(1,nd); B(5, col1) dNxy(3,nd); B(5, col2) dNxy(2,nd); B(6, col) dNxy(3,nd); B(6, col2) dNxy(1,nd); end K K (B*D*B) * detJ; end end end endB 矩阵的行对应应变分量顺序 [εx; εy; εz; γxy; γyz; γzx]。下面的表帮助排错时定位索引应变分量第 nd 个节点的 UVWεxB(1,3nd-2)00εy0B(2,3nd-1)0εz00B(3,3nd)γxyB(4,3nd-2)B(4,3nd-1)0γyz0B(5,3nd-1)B(5,3nd)γzxB(6,3nd-2)0B(6,3nd)这个顺序不是唯一的。有些程序把 γzx 放在 γxy 前面对应的 D 矩阵也要跟着调换。我习惯用上述顺序因为和大多数教科书一致后处理时也方便。3.3 从 2D 推到 3D 时容易踩的三个坑第一个坑是自由度映射。二维用2*nid-1和2*nid三维必须改成3*nid-2、3*nid-1、3*nid。如果沿用二维的索引全局刚度矩阵会错位但 MATLAB 不报错最后位移结果会非常离谱。第二个坑是节点顺序。六面体比四边形更容易出现负雅可比行列式尤其是通过外部网格生成器导入的单元节点顺序往往不统一。所以在hex8node里加一个detJ判断很值得即使网格有几千个单元这个判断也很快。第三个坑是积分阶数。2×2×2 的积分对规则六面体是准确的但对扭曲单元只是近似。不要随意改成 3×3×3那会显著增加计算量精度却不会按比例提升。4. matlab 有限元编程求解实例悬臂梁的四边形与六面体对比4.1 算例参数与二维网格生成算例取悬臂梁长 100mm高 10mm厚 1mm。左端固定右端承受向下的总力 100N。材料参数 E200GPaν0.3。用第 2 章的 quad2d4node 和第 3 章的 hex8node 分别求解自由端挠度再和欧拉伯努利梁的理论值对比。2D 平面应力模型在 x-y 平面建立网格长度方向 20 个单元高度方向 4 个单元。网格生成用 meshgrid 和 reshape 完成nx 20; ny 4; [xg, yg] meshgrid(linspace(0,100,nx1), linspace(0,10,ny1)); node2 [xg(:), yg(:)]; idx reshape(1:(nx1)*(ny1), ny1, nx1); elem2 zeros(nx*ny,4); for ex 1:nx for ey 1:ny n1 idx(ey,ex); n2 idx(ey,ex1); n3 idx(ey1,ex1); n4 idx(ey1,ex); elem2((ex-1)*ny ey,:) [n1 n2 n3 n4]; end endmeshgrid生成的 y 行优先排布所以idx的维度是 (ny1)×(nx1)。每个四边形四个角按左上、右上、右下、左下的顺序取出和quad2d4node的内部顺序一致。直接调用组装函数得到全局刚度矩阵K2 assemble_quad2d(node2, elem2, 200e3, 0.3, stress);4.2 三维六面体网格与组装函数三维网格在 x-y 平面保留同样的拓扑只在 z 方向增加一层。nz1 表示厚度方向只放一个单元因为弯曲应力在厚度方向是线性变化一个线性单元已经可以表示。nz 1; t 1; [xg3, yg3, zg3] meshgrid(linspace(0,100,nx1), ... linspace(0,10,ny1), ... linspace(0,t,nz1)); node3 [xg3(:), yg3(:), zg3(:)]; idx3 reshape(1:(nx1)*(ny1)*(nz1), ny1, nx1, nz1); elem3 zeros(nx*ny*nz,8); for ex 1:nx for ey 1:ny for ez 1:nz a [idx3(ey,ex,ez), idx3(ey,ex1,ez), ... idx3(ey1,ex1,ez), idx3(ey1,ex,ez)]; b [idx3(ey,ex,ez1), idx3(ey,ex1,ez1), ... idx3(ey1,ex1,ez1), idx3(ey1,ex,ez1)]; elem3((ex-1)*ny*nz (ey-1)*nz ez,:) [a b]; end end end对应的三维组装函数自由度映射是3*nid-2、3*nid-1、3*nidfunction K assemble_hex8(nodes, elems, E, nu) N size(nodes,1); K sparse(3*N, 3*N); for e 1:size(elems,1) nd elems(e,:); ke hex8node(nodes(nd,:), E, nu); dof zeros(1,24); dof(1:3:end) 3*nd - 2; dof(2:3:end) 3*nd - 1; dof(3:3:end) 3*nd; K(dof,dof) K(dof,dof) ke; end end这里不需要单独处理单元节点顺序hex8node内部的detJ检查会拦住倒置单元。4.3 边界条件、载荷与求解流程二维模型固定左端节点 x0 的全部自由度右端竖向力均匀分配到端部所有节点避免单点加载造成的应力集中fixed2 find(node2(:,1) 0); fdof2 reshape([2*fixed2-1; 2*fixed2], [], 1); nn size(node2,1) * 2; free2 setdiff(1:nn, fdof2); F2 zeros(nn,1); tip2 find(node2(:,1) 100); F2(2*tip2) -100 / length(tip2); U2 zeros(nn,1); U2(free2) K2(free2,free2) \ F2(free2);三维模型固定左端节点的 x、y、z 三个自由度右端节点施加 z 向力这里 z 沿梁宽方向实际弯曲发生在 x-y 平面所以力沿 y 方向fixed3 find(node3(:,1) 0); fdof3 reshape([3*fixed3-2; 3*fixed3-1; 3*fixed3], [], 1); nn3 size(node3,1) * 3; free3 setdiff(1:nn3, fdof3); F3 zeros(nn3,1); tip3 find(node3(:,1) 100); F3(3*tip3 - 1) -100 / length(tip3); % y方向分量为3*n-1 U3 zeros(nn3,1); U3(free3) K3(free3,free3) \ F3(free3);注意三维中 y 方向的自由度索引是3*n-1不是2*n。这是初学三维有限元最容易写错的一行。4.4 结果比较与误差判断理论挠度用欧拉伯努利梁公式计算w_max FL³/(3EI)其中 I tH³/12。取 F100NL100mmH10mmt1mm得到 w≈0.0675mm。实际运行后由于网格较粗且是低阶单元位移会偏小即刚度过大。模型网格规模端部位移 mm与理论偏差理论值—0.0675—2D quad20×40.0609 左右约 -9.8%3D hex820×4×10.0580 左右约 -14.1%这个偏差主要不是网格密度不足而是线性单元的剪切锁死。厚度方向只有一个单元时纯弯曲应变被剪切应变污染单元刚度偏大。要改善有两个办法把厚度方向单元增加到 2 至 4 层或者改用下面第 5 章的减缩积分。算例跑通后还要检查三件事固定端反力之和是否等于 -100N中性轴上的节点位移是否沿长度方向平滑变化对称轴上的节点有没有不正常的横向位移。这三项都通过才说明单元函数、组装和求解流程是通的。5. 六面体单元的减缩积分、网格收敛与薄壁圆筒延伸5.1 减缩积分在八节点六面体上的效果对于线性四边形和线性六面体全积分在纯弯曲问题中会放大剪切刚度。把 2×2×2 的高斯积分降为 1×1×1也就是只取单元中心点能显著改善位移结果。实现时只需要把hex8node里的积分坐标数组改成gp 0; % 单个积分点 w 2; % 三维情况下权重是 8对应地三重循环改成单次积分。注意一点单点积分引入沙漏模式也就是单元不需要能量就能发生锯齿变形。在小变形静力问题里如果单元比较规则沙漏模式通常不会主导结果但网格严重扭曲时必须额外加沙漏控制项。自研程序里保留两种积分模式是比较稳妥的在hex8node前面加一个order参数根据参数选择积分点即可。5.2 网格收敛性的验证方法把第 4 章的算例参数nx、ny从 20×4 改成 40×8、80×16每轮重新组装求解记录自由端位移。线性单元在弯曲问题中的位移收敛率约为 O(h²)所以网格每加密一倍误差应缩小约 4 倍。用 log-log 坐标画曲线斜率接近 2 说明程序正确。收敛检查的代码模式m [20 40 80]; w_end zeros(size(m)); for k 1:3 w_end(k) runHexBeam(m(k), m(k)/5, 1); end loglog(m, abs(w_end - w_theory), -o);其中runHexBeam是把第 4 章的网格生成、组装、求解流程封装成函数返回自由端位移。如果斜率明显小于 1优先检查加密后是否有单元 Jacobi 系数恶化或者载荷在端面的分配方式是否随网格变化而变化。5.3 薄壁圆筒算例的常见坑薄壁圆筒是检验六面体单元在曲面结构上表现的经典算例。它要求网格沿周向、径向和轴向三个方向布置内表面施加压力时不能直接加节点力而要用形函数把面压力积分成等效节点力。常见的错误是周向单元数太少圆筒被近似成多边形径向位移偏小周向单元数太多而径向层数不变又会出现长宽比过大的畸形单元。正确的做法是先固定轴向单元数再扫一遍周向网格密度观察径向位移随密度变化的趋势同时和拉梅公式的解析解对比。如果六面体单元能通过薄壁圆筒这个验证再回来看悬臂梁的弯曲问题基本就不会再被单元索引或积分策略卡住了。本文还有配套的精品资源点击获取