简介这套基于MATLAB的有限元程序设计算例包覆盖从杆系到空间块体的典型单元实现适合正在学习有限元理论、需要结合代码验证矩阵推导的本科生与工程技术人员。压缩包内共56个文件其中41个m脚本对应各算例主程序与子函数6个doc文档为算例说明另有少量dat数据、txt文本及asv自动备份文件整体仅333KB便于快速下载与本地运行。具体算例包括三梁平面框架Beam2D2Node、四杆桁架Bar2D2Node、基于三角形、四边形单元的矩形薄板、基于四面体、六面体单元的空间块体分析等也包含斜支座处理等典型例题通过说明文档与MATLAB代码联动可直观对比不同单元类型的构造差异与计算结果。目前已有1249人学习下载适合配合有限元课程或自学使用是理解节点编号、单元刚度矩阵与后处理流程的实用参考。1. 用 MATLAB 写有限元程序六面体网格为什么比想象中更好上手用 MATLAB 写有限元程序的人第一眼看到 quad2d4node 这个函数名很容易产生两种反应一个认为它是二维四边形四节点单元和六面体没有关系另一个直接按字面猜测它是某种六面体单元的变种。实际在网上下载的带算例源码包里quad2d4node 往往只是作者早期的平面单元函数它的形函数结构一旦看明白扩展成三维六面体八节点单元其实就是多乘一项的问题。本文围绕六面体 8 节点等参元把形函数、Jacobian 矩阵、高斯积分、MATLAB 刚度矩阵实现和悬臂梁算例串成一条能直接跑通的路线最后给出负 Jacobian 排查与分片验证的方法。适合刚读完有限元教材、想在 MATLAB 里亲手跑通第一个实体单元分析的工程师和学生。2. 六面体有限元基础8节点等参元的形函数与Jacobian矩阵2.1 从quad2d4node到hex8二维四边形到六面体的维度扩展quad2d4node 的本质是二维四节点四边形单元的形函数实现四个节点的局部坐标是 (-1,-1)、(1,-1)、(1,1)、(-1,1)每个节点对应一个双线性形函数。六面体 8 节点单元与它的亲缘关系非常直接把 xi、eta 两个方向扩展到 xi、eta、zeta 三个方向每个节点的形函数变成三个方向双线性因子的乘积即 N_i (1/8)(1xi_i·xi)(1eta_i·eta)(1zeta_i·zeta)。这个单元就是常说的 hex8也有人叫它三线性六面体单元。从 quad2d4node 迁移到 hex8差异集中在三处。一是形函数多乘一项 (1zeta_i·zeta)因此形函数对局部坐标的导数矩阵 dNdxi 也增加了一行二是高斯积分点从二维平面的 2x2 变成三维立方体里的 2x2x2三是应变位移矩阵 B 从 3x8 变成 6x24对应三个线应变 eps_x、eps_y、eps_z 和三个剪应变 gamma_xy、gamma_yz、gamma_xz。理解这三处差异之后单元刚度矩阵的程序框架不需要做任何结构性修改循环、组装、求解逻辑全部沿用。这也是为什么很多开源算例包里既有 quad2d4node又有 hex8 版本两个函数放在一起对比代码骨架完全一致。2.2 高斯积分从一维2点积分到三维2x2x2刚度矩阵中的积分项 K_e ∫ B^T D B dV 没有显式原函数程序里用高斯数值积分在单元内部逐点累加。一维两点高斯规则用两个积分点加权求和对不超过三次的多项式可以精确积分三维情况由三个方向的张量积得到2x2x2 表示每个方向取两个点共八个积分点。下面的表是常用的一维高斯积分点设置三维权重直接取三个方向权重的乘积。积分点数 n积分点位置 xi权重 w1022-1/sqrt(3), 1/sqrt(3)1, 13-sqrt(3/5), 0, sqrt(3/5)5/9, 8/9, 5/9例如在 MATLAB 里实现一维两点积分代码结构是gp [-1/sqrt(3), 1/sqrt(3)]; % 2点高斯积分点 gw [1, 1]; % 对应权重 W 0; for i 1:2 W W gw(i) * f(gp(i)); % 累加 f 在各积分点的加权值 end代码里的 f 是任意被积函数gp 和 gw 数量必须一致换三点点只需把 gp 和 gw 替换成表中第三行。对六面体单元嵌套三层循环分别遍历 xi、eta、zeta 方向的积分点每个积分点的总权重就是三个方向权重的乘积。2.3 Jacobian矩阵与单元坐标变换形函数给出的是局部坐标下的插值关系而单元刚度矩阵需要全局坐标下的导数这个转换依赖 Jacobian 矩阵 J。J 的每一行是局部坐标对全局坐标的映射具体计算用 dNdxi 乘以单元节点坐标矩阵 node_e。J 的行列式 det(J) 是局部体积到全局体积的缩放因子是所有高斯积分点处体积分的基础。det(J) 如果不是正值单元刚度矩阵就会出现不正定问题求解出的位移完全不可信。最典型的成因是单元节点顺序不符合逆时针约定或者单元发生了内凹、扭转。类似薄壁圆筒有限元这类厚度方向尺寸很小的模型径向只有一层单元时 det(J) 会非常小计算出的应力常常对网格畸变特别敏感。所以在写单元函数时我一般会在刚度矩阵循环里同时输出每个积分点的 detJ先定位负值再继续求解。3. MATLAB实现六面体刚度矩阵从quad2d4node到hex83.1 形函数与B矩阵四节点循环改写为八节点先写形函数函数这个函数同时返回 N 和 dNdxi后面的 hex8 单元刚度矩阵直接调用它。节点编号顺序必须固定前四个节点是 zeta-1 面上的四个角点后四个是 zeta1 面上的对应角点每个面内按逆时针排列。function [N, dNdxi] shape_hex8(xi, eta, zeta) % 8节点六面体单元形函数 % xi, eta, zeta: 局部坐标取值 [-1, 1] % N: 1x8 形函数值 % dNdxi: 3x8 形函数对局部坐标的导数 xi_i [-1 1 1 -1 -1 1 1 -1]; eta_i [-1 -1 1 1 -1 -1 1 1]; zeta_i [-1 -1 -1 -1 1 1 1 1]; for i 1:8 N(i) (1/8) * (1 xi_i(i)*xi) * (1 eta_i(i)*eta) * (1 zeta_i(i)*zeta); dNdxi(1,i) (1/8) * xi_i(i) * (1 eta_i(i)*eta) * (1 zeta_i(i)*zeta); dNdxi(2,i) (1/8) * eta_i(i) * (1 xi_i(i)*xi) * (1 zeta_i(i)*zeta); dNdxi(3,i) (1/8) * zeta_i(i) * (1 xi_i(i)*xi) * (1 eta_i(i)*eta); end end这里的 xi_i、eta_i、zeta_i 是八个节点的局部坐标值与节点编号一一对应。不想手写三个数组的话也可以把 2x2x2 角点坐标放进 8x3 矩阵再用循环读出但上面的写法更直观也方便对照教材公式排错。3.2 单元刚度矩阵的MATLAB代码与高斯积分循环有了形函数再写单元刚度矩阵。材料矩阵 D 用各向同性线弹性假设ngp 控制每个方向的积分点数默认 2 对应 2x2x2 完整积分。整个函数返回 24x24 的单元刚度矩阵因为每个节点三个平动自由度。function Ke hex8_stiffness(node_e, E, nu, ngp) % 返回 24x24 六面体单元刚度矩阵 % node_e: 8x3 单元节点全局坐标 % E, nu: 弹性模量与泊松比 % ngp: 每个方向高斯积分点数量通常取2 C E/((1nu)*(1-2*nu)) * ... [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, gw] gauss_points(ngp); Ke zeros(24, 24); for i 1:ngp for j 1:ngp for k 1:ngp [N, dNdxi] shape_hex8(gp(i), gp(j), gp(k)); J dNdxi * node_e; % 3x3 Jacobian 矩阵 invJ inv(J); dNdx invJ * dNdxi; % 全局坐标下的形函数导数3x8 B zeros(6, 24); for a 1:8 B(:, 3*a-2:3*a) ... [dNdx(1,a) 0 0; 0 dNdx(2,a) 0; 0 0 dNdx(3,a); dNdx(2,a) dNdx(1,a) 0; 0 dNdx(3,a) dNdx(2,a); dNdx(3,a) 0 dNdx(1,a)]; end detJ det(J); Ke Ke gw(i)*gw(j)*gw(k) * (B*C*B) * detJ; end end end end function [gp, gw] gauss_points(ngp) % 一维高斯积分点和权重 if ngp 1 gp 0; gw 2; elseif ngp 2 gp [-1/sqrt(3), 1/sqrt(3)]; gw [1, 1]; else gp [-sqrt(3/5), 0, sqrt(3/5)]; gw [5/9, 8/9, 5/9]; end end代码里 B 矩阵的每一列对应一个节点的三个自由度排列顺序是 ux、uy、uz。dNdx invJ * dNdxi 这一步容易写错注意是转置而不是直接左乘矩阵乘法的维度对应关系是 (3x3)^T * (3x8)结果仍是 3x8。detJ 必须大于零因此上面的代码没有做符号判断实际用的时候建议在每个积分点处加一个 if detJ 0 的警告。注意单元节点顺序一旦颠倒detJ 变负求解出的刚度矩阵不正定。出现这种情况时先检查 node_e 的行顺序而不是怀疑公式。3.3 整体刚度矩阵组装与节点编号规则单元刚度矩阵算出来后要按自由度编号放入整体矩阵。MATLAB 里用稀疏矩阵存储可以显著降低内存占用尤其是网格规模超过几千个节点时。整体刚度矩阵 K 的组装逻辑是遍历所有单元把单元自由度编号映射到全局自由度编号。% node: 节点坐标矩阵size nnode x 3 % elem: 单元连接矩阵size nelem x 8 ndof 3 * size(node, 1); K sparse(ndof, ndof); for e 1:size(elem, 1) idx elem(e, :); % 当前单元的8个节点编号 edof zeros(24, 1); for a 1:8 edof(3*a-2:3*a) 3*idx(a)-2 : 3*idx(a); % 节点a的三个自由度 end Ke hex8_stiffness(node(idx, :), E, nu, ngp); K(edof, edof) K(edof, edof) Ke; endedof 的生成是组装最需要小心的部分。3*idx(a)-2 是该节点第一个自由度的全局编号例如节点编号 1 对应自由度 1、2、3节点编号 2 对应 4、5、6依此类推。K(edof, edof) 必须用累加而不是赋值因为一个节点通常被多个单元共享。4. 算例MATLAB悬臂梁六面体网格划分与求解4.1 悬臂梁参数表与节点坐标生成用悬臂梁做六面体有限元编程求解实例是验证程序正确性最直观的做法。梁沿 x 方向长度为 1 米截面为 0.1m x 0.1m 的正方形左端固定右端施加向下的总力 1000N。为了兼顾计算速度和网格形态单元数取 16x4x4即长度方向 16 个单元高度和宽度方向各 4 个单元。参数取值说明梁长 Lx1.0 mx 方向截面 Ly x Lz0.1 x 0.1 my 与 z 方向单元数 nx x ny x nz16 x 4 x 4长度 x 高度 x 宽度节点总数17 x 5 x 5 425与单元数对应自由度数12753 x 节点总数弹性模量 E200 GPa各向同性材料泊松比 nu0.3各向同性材料端部载荷-1000 N沿 y 方向均布在自由端面节点网格生成函数的关键是节点编号按 i 方向最快、j 次之、k 最慢的方式推进这样单元连接矩阵的规律性最强也最容易检查和排错。nx 16; ny 4; nz 4; Lx 1.0; Ly 0.1; Lz 0.1; dx Lx/nx; dy Ly/ny; dz Lz/nz; % 生成节点坐标 nnode (nx1)*(ny1)*(nz1); node zeros(nnode, 3); n 0; for k 0:nz for j 0:ny for i 0:nx n n 1; node(n, :) [i*dx, j*dy, k*dz]; end end end % 生成单元连接矩阵 nelem nx*ny*nz; elem zeros(nelem, 8); nxy (nx1)*(ny1); e 0; for k 0:nz-1 for j 0:ny-1 for i 0:nx-1 n0 k*nxy j*(nx1) i 1; e e 1; elem(e, :) [n0, n01, n01(nx1), n0(nx1), ... n0nxy, n0nxy1, n0nxy1(nx1), n0nxy(nx1)]; end end end单元连接矩阵里的 n0 是当前单元的起始节点括号里的偏移量按照“先 x 后 y 再 z”的顺序推导。比如 n01 是同一层右侧节点n0(nx1) 是后一排节点n0nxy 是上一层对应节点。这个顺序和 shape_hex8 里的节点顺序保持一致如果调整遍历顺序必须同步修改形函数的坐标数组。4.2 边界条件处理与线性方程组求解边界条件分两步找出左端面 x0 的节点约束这些节点的全部自由度找出右端面 xLx 的节点把 1000N 均分到每个节点上。MATLAB 的 find 可以按坐标筛选节点再通过自由度编号映射到全局方程组。ndof 3 * size(node, 1); K sparse(ndof, ndof); for e 1:size(elem, 1) idx elem(e, :); edof zeros(24, 1); for a 1:8 edof(3*a-2:3*a) 3*idx(a)-2 : 3*idx(a); end Ke hex8_stiffness(node(idx, :), 200e9, 0.3, 2); K(edof, edof) K(edof, edof) Ke; end % 固定端x0 的所有节点六个自由度全部约束 fixed_nodes find(node(:,1) 1e-10); fixed_dofs []; for i 1:length(fixed_nodes) fixed_dofs [fixed_dofs, ... 3*fixed_nodes(i)-2, 3*fixed_nodes(i)-1, 3*fixed_nodes(i)]; end % 载荷x1.0 自由端面的节点沿 y 方向均匀分配 load_nodes find(node(:,1) 1.0 - 1e-10); n_load length(load_nodes); F sparse(ndof, 1); for i 1:n_load F(3*load_nodes(i)-1) -1000.0 / n_load; end % 求解自由位移 free_dofs setdiff(1:ndof, fixed_dofs); d zeros(ndof, 1); d(free_dofs) K(free_dofs, free_dofs) \ F(free_dofs);fprintf 输出最大位移时可以直接查看 d 中对应自由度的值。固定端全部自由度约束后整体矩阵可逆反斜杠求解会走直接法对 1275 个自由度来说计算时间在一秒以内。这里没有使用缩减自由度的方法因为先组装完整 K 再取子块代码更容易理解也方便后续做约束方程扩展。4.3 位移云图绘制与结果合理性检查求解结束后把位移分量从 d 里拆出来再按节点坐标绘制变形后的散点云图。位移数量级很小直接绘制几乎看不出变形所以需要乘一个放大系数。Ux d(1:3:end); Uy d(2:3:end); Uz d(3:3:end); umag sqrt(Ux.^2 Uy.^2 Uz.^2); % 变形放大到梁长的 10% 左右 scale 0.1 / max(abs(Uy)); node_d node scale * [Ux, Uy, Uz]; figure; scatter3(node_d(:,1), node_d(:,2), node_d(:,3), 36, umag, filled); axis equal; colorbar; xlabel(x); ylabel(y); zlabel(z);理论上梁端挠度可以用欧拉梁公式估算I Ly^3 * Lz / 12 8.33e-6 m^4δ P·L^3 / (3·E·I)代入数值得到约 0.2mm。六面体实体单元在 16x4x4 网格下得到的端部位移会略大于这个值因为实体单元的剪切变形和端部载荷的等效方式都会让挠度偏大。只要误差在百分之几到十几之间程序逻辑就是对的。5. 验证与进阶负Jacobian排查、单单元测试与降阶积分5.1 负Jacobian的快速检测网格单元数较多时逐单元打印 detJ 太啰嗦。常见做法是在组装循环里加一个中心积分点检测单元中心点处 detJ 最容易反映整体畸变程度。下面的代码可以在组装前单独跑一遍[~, dNdxi0] shape_hex8(0, 0, 0); % 单元中心点 for e 1:size(elem,1) Jc dNdxi0 * node(elem(e,:), :); if det(Jc) 1e-12 fprintf(单元 %d 负Jacobian: %e\n, e, det(Jc)); end enddetJ 出现负值最常见的原因是单元节点顺序错误其次是单元内凹、长宽比过大。网格长宽比超过 10 时即使 detJ 为正应力精度也会明显下降尤其是薄壁圆筒有限元这类厚度方向尺寸远小于其他方向的模型。5.2 单单元压缩测试与分片验证整体算例跑通之前先做一个单单元压缩测试。用坐标 [0,1]^3 的单个八节点立方体左端面 x0 全约束右端面 x1 的节点直接赋位移 ux1e-3然后求解支反力。理想情况下应力 σx E·ε 200e9×1e-3 2e8 Pa支反力除以面积应等于 2e8。这个测试只需要一个单元几分钟就能定位公式错误。5.3 降阶积分、优化工具箱与批处理程序跑通后把 hex8_stiffness 的 ngp 改为 1 就是降阶积分。降阶积分能缓解弯曲问题中的剪切锁闭但会引入零能模式实体结构只用中心点积分时刚度偏柔。我一般对弯曲占优的梁板问题用 1 点积分体积变形占优的实体问题保持 2x2x2。如果要做尺寸优化把整个求解流程封装成函数 solve_cantilever(ny, nz)目标函数返回最大位移再用优化工具箱的 fmincon 自动搜索截面尺寸多个网格方案对比时用 parfor 替换 for 循环即可。本文还有配套的精品资源点击获取
MATLAB六面体有限元编程：从形函数到刚度矩阵与悬臂梁算例
简介这套基于MATLAB的有限元程序设计算例包覆盖从杆系到空间块体的典型单元实现适合正在学习有限元理论、需要结合代码验证矩阵推导的本科生与工程技术人员。压缩包内共56个文件其中41个m脚本对应各算例主程序与子函数6个doc文档为算例说明另有少量dat数据、txt文本及asv自动备份文件整体仅333KB便于快速下载与本地运行。具体算例包括三梁平面框架Beam2D2Node、四杆桁架Bar2D2Node、基于三角形、四边形单元的矩形薄板、基于四面体、六面体单元的空间块体分析等也包含斜支座处理等典型例题通过说明文档与MATLAB代码联动可直观对比不同单元类型的构造差异与计算结果。目前已有1249人学习下载适合配合有限元课程或自学使用是理解节点编号、单元刚度矩阵与后处理流程的实用参考。1. 用 MATLAB 写有限元程序六面体网格为什么比想象中更好上手用 MATLAB 写有限元程序的人第一眼看到 quad2d4node 这个函数名很容易产生两种反应一个认为它是二维四边形四节点单元和六面体没有关系另一个直接按字面猜测它是某种六面体单元的变种。实际在网上下载的带算例源码包里quad2d4node 往往只是作者早期的平面单元函数它的形函数结构一旦看明白扩展成三维六面体八节点单元其实就是多乘一项的问题。本文围绕六面体 8 节点等参元把形函数、Jacobian 矩阵、高斯积分、MATLAB 刚度矩阵实现和悬臂梁算例串成一条能直接跑通的路线最后给出负 Jacobian 排查与分片验证的方法。适合刚读完有限元教材、想在 MATLAB 里亲手跑通第一个实体单元分析的工程师和学生。2. 六面体有限元基础8节点等参元的形函数与Jacobian矩阵2.1 从quad2d4node到hex8二维四边形到六面体的维度扩展quad2d4node 的本质是二维四节点四边形单元的形函数实现四个节点的局部坐标是 (-1,-1)、(1,-1)、(1,1)、(-1,1)每个节点对应一个双线性形函数。六面体 8 节点单元与它的亲缘关系非常直接把 xi、eta 两个方向扩展到 xi、eta、zeta 三个方向每个节点的形函数变成三个方向双线性因子的乘积即 N_i (1/8)(1xi_i·xi)(1eta_i·eta)(1zeta_i·zeta)。这个单元就是常说的 hex8也有人叫它三线性六面体单元。从 quad2d4node 迁移到 hex8差异集中在三处。一是形函数多乘一项 (1zeta_i·zeta)因此形函数对局部坐标的导数矩阵 dNdxi 也增加了一行二是高斯积分点从二维平面的 2x2 变成三维立方体里的 2x2x2三是应变位移矩阵 B 从 3x8 变成 6x24对应三个线应变 eps_x、eps_y、eps_z 和三个剪应变 gamma_xy、gamma_yz、gamma_xz。理解这三处差异之后单元刚度矩阵的程序框架不需要做任何结构性修改循环、组装、求解逻辑全部沿用。这也是为什么很多开源算例包里既有 quad2d4node又有 hex8 版本两个函数放在一起对比代码骨架完全一致。2.2 高斯积分从一维2点积分到三维2x2x2刚度矩阵中的积分项 K_e ∫ B^T D B dV 没有显式原函数程序里用高斯数值积分在单元内部逐点累加。一维两点高斯规则用两个积分点加权求和对不超过三次的多项式可以精确积分三维情况由三个方向的张量积得到2x2x2 表示每个方向取两个点共八个积分点。下面的表是常用的一维高斯积分点设置三维权重直接取三个方向权重的乘积。积分点数 n积分点位置 xi权重 w1022-1/sqrt(3), 1/sqrt(3)1, 13-sqrt(3/5), 0, sqrt(3/5)5/9, 8/9, 5/9例如在 MATLAB 里实现一维两点积分代码结构是gp [-1/sqrt(3), 1/sqrt(3)]; % 2点高斯积分点 gw [1, 1]; % 对应权重 W 0; for i 1:2 W W gw(i) * f(gp(i)); % 累加 f 在各积分点的加权值 end代码里的 f 是任意被积函数gp 和 gw 数量必须一致换三点点只需把 gp 和 gw 替换成表中第三行。对六面体单元嵌套三层循环分别遍历 xi、eta、zeta 方向的积分点每个积分点的总权重就是三个方向权重的乘积。2.3 Jacobian矩阵与单元坐标变换形函数给出的是局部坐标下的插值关系而单元刚度矩阵需要全局坐标下的导数这个转换依赖 Jacobian 矩阵 J。J 的每一行是局部坐标对全局坐标的映射具体计算用 dNdxi 乘以单元节点坐标矩阵 node_e。J 的行列式 det(J) 是局部体积到全局体积的缩放因子是所有高斯积分点处体积分的基础。det(J) 如果不是正值单元刚度矩阵就会出现不正定问题求解出的位移完全不可信。最典型的成因是单元节点顺序不符合逆时针约定或者单元发生了内凹、扭转。类似薄壁圆筒有限元这类厚度方向尺寸很小的模型径向只有一层单元时 det(J) 会非常小计算出的应力常常对网格畸变特别敏感。所以在写单元函数时我一般会在刚度矩阵循环里同时输出每个积分点的 detJ先定位负值再继续求解。3. MATLAB实现六面体刚度矩阵从quad2d4node到hex83.1 形函数与B矩阵四节点循环改写为八节点先写形函数函数这个函数同时返回 N 和 dNdxi后面的 hex8 单元刚度矩阵直接调用它。节点编号顺序必须固定前四个节点是 zeta-1 面上的四个角点后四个是 zeta1 面上的对应角点每个面内按逆时针排列。function [N, dNdxi] shape_hex8(xi, eta, zeta) % 8节点六面体单元形函数 % xi, eta, zeta: 局部坐标取值 [-1, 1] % N: 1x8 形函数值 % dNdxi: 3x8 形函数对局部坐标的导数 xi_i [-1 1 1 -1 -1 1 1 -1]; eta_i [-1 -1 1 1 -1 -1 1 1]; zeta_i [-1 -1 -1 -1 1 1 1 1]; for i 1:8 N(i) (1/8) * (1 xi_i(i)*xi) * (1 eta_i(i)*eta) * (1 zeta_i(i)*zeta); dNdxi(1,i) (1/8) * xi_i(i) * (1 eta_i(i)*eta) * (1 zeta_i(i)*zeta); dNdxi(2,i) (1/8) * eta_i(i) * (1 xi_i(i)*xi) * (1 zeta_i(i)*zeta); dNdxi(3,i) (1/8) * zeta_i(i) * (1 xi_i(i)*xi) * (1 eta_i(i)*eta); end end这里的 xi_i、eta_i、zeta_i 是八个节点的局部坐标值与节点编号一一对应。不想手写三个数组的话也可以把 2x2x2 角点坐标放进 8x3 矩阵再用循环读出但上面的写法更直观也方便对照教材公式排错。3.2 单元刚度矩阵的MATLAB代码与高斯积分循环有了形函数再写单元刚度矩阵。材料矩阵 D 用各向同性线弹性假设ngp 控制每个方向的积分点数默认 2 对应 2x2x2 完整积分。整个函数返回 24x24 的单元刚度矩阵因为每个节点三个平动自由度。function Ke hex8_stiffness(node_e, E, nu, ngp) % 返回 24x24 六面体单元刚度矩阵 % node_e: 8x3 单元节点全局坐标 % E, nu: 弹性模量与泊松比 % ngp: 每个方向高斯积分点数量通常取2 C E/((1nu)*(1-2*nu)) * ... [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, gw] gauss_points(ngp); Ke zeros(24, 24); for i 1:ngp for j 1:ngp for k 1:ngp [N, dNdxi] shape_hex8(gp(i), gp(j), gp(k)); J dNdxi * node_e; % 3x3 Jacobian 矩阵 invJ inv(J); dNdx invJ * dNdxi; % 全局坐标下的形函数导数3x8 B zeros(6, 24); for a 1:8 B(:, 3*a-2:3*a) ... [dNdx(1,a) 0 0; 0 dNdx(2,a) 0; 0 0 dNdx(3,a); dNdx(2,a) dNdx(1,a) 0; 0 dNdx(3,a) dNdx(2,a); dNdx(3,a) 0 dNdx(1,a)]; end detJ det(J); Ke Ke gw(i)*gw(j)*gw(k) * (B*C*B) * detJ; end end end end function [gp, gw] gauss_points(ngp) % 一维高斯积分点和权重 if ngp 1 gp 0; gw 2; elseif ngp 2 gp [-1/sqrt(3), 1/sqrt(3)]; gw [1, 1]; else gp [-sqrt(3/5), 0, sqrt(3/5)]; gw [5/9, 8/9, 5/9]; end end代码里 B 矩阵的每一列对应一个节点的三个自由度排列顺序是 ux、uy、uz。dNdx invJ * dNdxi 这一步容易写错注意是转置而不是直接左乘矩阵乘法的维度对应关系是 (3x3)^T * (3x8)结果仍是 3x8。detJ 必须大于零因此上面的代码没有做符号判断实际用的时候建议在每个积分点处加一个 if detJ 0 的警告。注意单元节点顺序一旦颠倒detJ 变负求解出的刚度矩阵不正定。出现这种情况时先检查 node_e 的行顺序而不是怀疑公式。3.3 整体刚度矩阵组装与节点编号规则单元刚度矩阵算出来后要按自由度编号放入整体矩阵。MATLAB 里用稀疏矩阵存储可以显著降低内存占用尤其是网格规模超过几千个节点时。整体刚度矩阵 K 的组装逻辑是遍历所有单元把单元自由度编号映射到全局自由度编号。% node: 节点坐标矩阵size nnode x 3 % elem: 单元连接矩阵size nelem x 8 ndof 3 * size(node, 1); K sparse(ndof, ndof); for e 1:size(elem, 1) idx elem(e, :); % 当前单元的8个节点编号 edof zeros(24, 1); for a 1:8 edof(3*a-2:3*a) 3*idx(a)-2 : 3*idx(a); % 节点a的三个自由度 end Ke hex8_stiffness(node(idx, :), E, nu, ngp); K(edof, edof) K(edof, edof) Ke; endedof 的生成是组装最需要小心的部分。3*idx(a)-2 是该节点第一个自由度的全局编号例如节点编号 1 对应自由度 1、2、3节点编号 2 对应 4、5、6依此类推。K(edof, edof) 必须用累加而不是赋值因为一个节点通常被多个单元共享。4. 算例MATLAB悬臂梁六面体网格划分与求解4.1 悬臂梁参数表与节点坐标生成用悬臂梁做六面体有限元编程求解实例是验证程序正确性最直观的做法。梁沿 x 方向长度为 1 米截面为 0.1m x 0.1m 的正方形左端固定右端施加向下的总力 1000N。为了兼顾计算速度和网格形态单元数取 16x4x4即长度方向 16 个单元高度和宽度方向各 4 个单元。参数取值说明梁长 Lx1.0 mx 方向截面 Ly x Lz0.1 x 0.1 my 与 z 方向单元数 nx x ny x nz16 x 4 x 4长度 x 高度 x 宽度节点总数17 x 5 x 5 425与单元数对应自由度数12753 x 节点总数弹性模量 E200 GPa各向同性材料泊松比 nu0.3各向同性材料端部载荷-1000 N沿 y 方向均布在自由端面节点网格生成函数的关键是节点编号按 i 方向最快、j 次之、k 最慢的方式推进这样单元连接矩阵的规律性最强也最容易检查和排错。nx 16; ny 4; nz 4; Lx 1.0; Ly 0.1; Lz 0.1; dx Lx/nx; dy Ly/ny; dz Lz/nz; % 生成节点坐标 nnode (nx1)*(ny1)*(nz1); node zeros(nnode, 3); n 0; for k 0:nz for j 0:ny for i 0:nx n n 1; node(n, :) [i*dx, j*dy, k*dz]; end end end % 生成单元连接矩阵 nelem nx*ny*nz; elem zeros(nelem, 8); nxy (nx1)*(ny1); e 0; for k 0:nz-1 for j 0:ny-1 for i 0:nx-1 n0 k*nxy j*(nx1) i 1; e e 1; elem(e, :) [n0, n01, n01(nx1), n0(nx1), ... n0nxy, n0nxy1, n0nxy1(nx1), n0nxy(nx1)]; end end end单元连接矩阵里的 n0 是当前单元的起始节点括号里的偏移量按照“先 x 后 y 再 z”的顺序推导。比如 n01 是同一层右侧节点n0(nx1) 是后一排节点n0nxy 是上一层对应节点。这个顺序和 shape_hex8 里的节点顺序保持一致如果调整遍历顺序必须同步修改形函数的坐标数组。4.2 边界条件处理与线性方程组求解边界条件分两步找出左端面 x0 的节点约束这些节点的全部自由度找出右端面 xLx 的节点把 1000N 均分到每个节点上。MATLAB 的 find 可以按坐标筛选节点再通过自由度编号映射到全局方程组。ndof 3 * size(node, 1); K sparse(ndof, ndof); for e 1:size(elem, 1) idx elem(e, :); edof zeros(24, 1); for a 1:8 edof(3*a-2:3*a) 3*idx(a)-2 : 3*idx(a); end Ke hex8_stiffness(node(idx, :), 200e9, 0.3, 2); K(edof, edof) K(edof, edof) Ke; end % 固定端x0 的所有节点六个自由度全部约束 fixed_nodes find(node(:,1) 1e-10); fixed_dofs []; for i 1:length(fixed_nodes) fixed_dofs [fixed_dofs, ... 3*fixed_nodes(i)-2, 3*fixed_nodes(i)-1, 3*fixed_nodes(i)]; end % 载荷x1.0 自由端面的节点沿 y 方向均匀分配 load_nodes find(node(:,1) 1.0 - 1e-10); n_load length(load_nodes); F sparse(ndof, 1); for i 1:n_load F(3*load_nodes(i)-1) -1000.0 / n_load; end % 求解自由位移 free_dofs setdiff(1:ndof, fixed_dofs); d zeros(ndof, 1); d(free_dofs) K(free_dofs, free_dofs) \ F(free_dofs);fprintf 输出最大位移时可以直接查看 d 中对应自由度的值。固定端全部自由度约束后整体矩阵可逆反斜杠求解会走直接法对 1275 个自由度来说计算时间在一秒以内。这里没有使用缩减自由度的方法因为先组装完整 K 再取子块代码更容易理解也方便后续做约束方程扩展。4.3 位移云图绘制与结果合理性检查求解结束后把位移分量从 d 里拆出来再按节点坐标绘制变形后的散点云图。位移数量级很小直接绘制几乎看不出变形所以需要乘一个放大系数。Ux d(1:3:end); Uy d(2:3:end); Uz d(3:3:end); umag sqrt(Ux.^2 Uy.^2 Uz.^2); % 变形放大到梁长的 10% 左右 scale 0.1 / max(abs(Uy)); node_d node scale * [Ux, Uy, Uz]; figure; scatter3(node_d(:,1), node_d(:,2), node_d(:,3), 36, umag, filled); axis equal; colorbar; xlabel(x); ylabel(y); zlabel(z);理论上梁端挠度可以用欧拉梁公式估算I Ly^3 * Lz / 12 8.33e-6 m^4δ P·L^3 / (3·E·I)代入数值得到约 0.2mm。六面体实体单元在 16x4x4 网格下得到的端部位移会略大于这个值因为实体单元的剪切变形和端部载荷的等效方式都会让挠度偏大。只要误差在百分之几到十几之间程序逻辑就是对的。5. 验证与进阶负Jacobian排查、单单元测试与降阶积分5.1 负Jacobian的快速检测网格单元数较多时逐单元打印 detJ 太啰嗦。常见做法是在组装循环里加一个中心积分点检测单元中心点处 detJ 最容易反映整体畸变程度。下面的代码可以在组装前单独跑一遍[~, dNdxi0] shape_hex8(0, 0, 0); % 单元中心点 for e 1:size(elem,1) Jc dNdxi0 * node(elem(e,:), :); if det(Jc) 1e-12 fprintf(单元 %d 负Jacobian: %e\n, e, det(Jc)); end enddetJ 出现负值最常见的原因是单元节点顺序错误其次是单元内凹、长宽比过大。网格长宽比超过 10 时即使 detJ 为正应力精度也会明显下降尤其是薄壁圆筒有限元这类厚度方向尺寸远小于其他方向的模型。5.2 单单元压缩测试与分片验证整体算例跑通之前先做一个单单元压缩测试。用坐标 [0,1]^3 的单个八节点立方体左端面 x0 全约束右端面 x1 的节点直接赋位移 ux1e-3然后求解支反力。理想情况下应力 σx E·ε 200e9×1e-3 2e8 Pa支反力除以面积应等于 2e8。这个测试只需要一个单元几分钟就能定位公式错误。5.3 降阶积分、优化工具箱与批处理程序跑通后把 hex8_stiffness 的 ngp 改为 1 就是降阶积分。降阶积分能缓解弯曲问题中的剪切锁闭但会引入零能模式实体结构只用中心点积分时刚度偏柔。我一般对弯曲占优的梁板问题用 1 点积分体积变形占优的实体问题保持 2x2x2。如果要做尺寸优化把整个求解流程封装成函数 solve_cantilever(ny, nz)目标函数返回最大位移再用优化工具箱的 fmincon 自动搜索截面尺寸多个网格方案对比时用 parfor 替换 for 循环即可。本文还有配套的精品资源点击获取
