简介面向需要开展有限元分析与四面体网格剖分的Matlab开发者、工程计算学习者这份源代码包提供了一套可直接运行的四面体剖分实现方案。包内共3个文件其中my_poufen.m为主程序涵盖读取几何数据、生成四面体网格、质量检查与可视化等完整环节NODE.txt记录节点三维坐标WN_NE.txt保存四面体单元连接信息两者共同构成后续刚度矩阵组装与求解的数据基础。整个压缩包大小仅2KB结构精简便于快速阅读和二次开发。目前已有601人学习下载适合希望从代码层面理解FEM离散化流程、掌握网格生成细节并优化求解效率的读者。通过研读源程序可以清晰看到节点编号、单元拓扑、边界条件处理如何衔接也能借鉴其数据结构设计思路迁移到更复杂的工程场景。整体注释简洁清晰变量命名规范是一份适合教学与自学的好模板。1. 四面体剖分从几何网格到有限元求解的第一道关在仿真项目里摸爬过几年的人都有这种体会模型画得再漂亮剖分这一步过不去后面全是空谈。四面体网格是三维有限元计算最通用的离散形式它能把任意复杂的三维实体切成一个个互不重叠的四面体单元让偏微分方程在每个小单元上获得近似解。这个标题里的“有限元分割”指的就是这一过程而Matlab实现四面体剖分的价值在于你不用离开Matlab环境就能完成从几何模型到可计算网格的全部准备工作。本文面向两类读者一类是做有限元编程求解实例、需要自己控制网格生成流程的研究者另一类是工程上做薄壁圆筒、梁结构或复杂铸件仿真被商业软件的黑盒网格逼疯的工程师。读完你可以用Matlab从零生成四面体网格理解剖分算法的选型逻辑并能定位那些让人挠头的网格质量问题的根源。这里不讨论用Patran或HyperMesh怎么点鼠标只讲用代码可控地完成这件事。2. 四面体剖分的算法选型Delaunay、八叉树与前沿推进的取舍2.1 为什么Delaunay是Matlab实现四面体剖分的主流选择四面体剖分在计算几何里已经有几十年的积累主流算法无非三支Delaunay四面体化、八叉树法Octree和前沿推进法Advancing Front Technique AFT。其中Delaunay是Matlab原生支持、上手最快的一条路原因在于它的数学性质足够“硬”。Delaunay四面体化满足空球准则任意一个四面体的外接球内部不包含其他顶点。这个准则保证了在非退化情况下剖分结果是唯一的并且能最大化所有四面体的最小二面角——通俗说它在“尽量不产生特别扁的单元”这个意义下是最优的。对有限元求解来说单元形状直接决定刚度矩阵的条件数和误差估计的可靠性一个接近退化薄片状的四面体足以让共轭梯度迭代次数翻几倍甚至在直接法里产生虚假的大位移。Matlab自带的delaunayTriangulation类底层基于CGAL的Delaunay实现支持三维点集的四面体化、约束边和约束面的嵌入、以及最近邻与插值查询。它不需要你额外安装任何工具箱语法也简单这是它在教学和原型验证里被用得最多的原因。2.2 八叉树和前沿推进法在什么场景下更适用八叉树法的思路是递归地把包围盒八等分直到每个叶子节点满足密度要求再把叶子节点里的点连接成单元。它的优点是剖分速度快、内存可控特别适合地形、城市建筑群这类大规模、表面不复杂的域。缺点也很明显表面单元与几何边界完全脱节要贴回原表面需要额外做投影和裁剪处理不好会在边界留下一层劣质单元。前沿推进法从边界三角形出发像“生长”一样逐层向内部插入节点、生成四面体。它能精确控制边界层网格的厚度和增长率因此在计算流体力学CFD的黏性边界层网格里几乎是不二之选。但AFT的鲁棒性依赖前沿面的拓扑管理实现复杂度远高于Delaunay而且当两个前沿面即将穿透时的冲突检测非常容易写错。算法核心思想表面保真边界层支持实现成本代表库Delaunay空球准则、逐点插入需约束恢复弱需各向异性扩展中CGAL、TetGen、Matlab自带Octree递归八等分包围盒差需投影修正弱低各开源网格器内部模块AFT从表面向内部推进好强高商业CFD网格器如果你做的是结构力学、传热、电磁场这类以体单元为主、边界层要求不高的仿真Matlab的Delaunay路线是成本最低的选择。要做边界层老老实实去调专业网格生成器别在Matlab里自己造AFT的轮子。2.3 从表面网格到体网格Matlab剖分的两种路线实际工程中几乎不会直接给出一堆散点让你剖而是先有CAD模型导出的表面网格通常是STL或OBJ再在内部填充四面体。针对这个流程Matlab下有两条技术路线。路线一是纯原生用triangulation读入表面用delaunayTriangulation做无约束四面体化再用几何谓词把位于表面之外的四面体剔除。这条路线代码量少、无第三方依赖适合几何不太复杂、对网格质量要求中等的场景也是后续三章的讲解主线。路线二是走第三方接口把表面网格导出为TetGen的.poly格式调用TetGen的可执行文件完成四面体化再读回Matlab分析。TetGen在约束Delaunay四面体化、网格质量优化和尺寸函数控制上比Matlab原生的Qhull路径成熟得多。代价是需要额外处理外部程序和文件格式而且TetGen的许可证对商用不友好你在项目里引入前要确认授权边界。3. 用MATLAB实现四面体剖分的最小可运行代码3.1 准备面网格输入从STL文件到三角剖分对象四面体剖分的输入通常是闭合的表面三角网格。第一步先把STL文件中重复的顶点合并因为STL用浮点数存储顶点同一个几何点在不同三角形里可能对应不同坐标值。以下代码完成STL读取与顶点去重。function [TR, Node, elem] read_stl_surface(filename) % filename: STL文件路径(ASCII或二进制) % 返回 TR: 表面三角剖分对象 % Node: 去重后的节点坐标, Nx3 % elem: 表面三角形节点索引, Mx3 % 使用Matlab自带的STL读取接口自动识别ASCII/二进制 stlData stlread(filename); tri stlData.ConnectivityList; % Mx3, 表面三角形 pt stlData.Points; % Px3, 含重复顶点 % 合并重复节点: 使用round将坐标对齐到容差网格 tol 1e-6; [Node, ~, ic] uniquetol(pt, tol, ByRows, true, ... OutputAllIndices, false); % 重新映射三角形顶点索引 elem ic(tri); TR triangulation(elem, Node); % 简单自检: 统计每条边出现的次数边界边应出现1次 edges edges(TR); fprintf(节点数: %d, 三角形数: %d\n, size(Node,1), size(elem,1)); end这里最关键的是uniquetol的参数tol是顶点合并的容差单位与模型几何单位一致。模型尺寸是米级1e-6通常够如果模型是微米级几何这个容差要把所有坐标先缩放到同一尺度再比较。直接用等号判断重复顶点几乎一定会漏因为不同三角形写入的坐标可能相差10的负16次方量级。3.2 Delaunay四面体化一条命令背后的合法性与限制拿到表面剖分对象后最直接的思路是对全部表面节点做Delaunay四面体化% 对表面节点直接做三维Delaunay DT delaunayTriangulation(Node); % 获取四面体连接关系 tet DT.ConnectivityList; fprintf(初始四面体数: %d\n, size(tet,1));这一步很快但有个致命问题无约束Delaunay四面体化只保证点集的凸包被剖分它完全不知道你原来的表面长什么样。对凹形体比如一个L形支架或带孔的法兰凸包区域里有大量四面体落在实体外部。第3.3节解决这个问题。另外如果点集存在退化情况多个点共球剖分结果不唯一后续单元质量评估会不稳定建议在剖分前对节点做微小随机扰动或规范化处理。delaunayTriangulation底层调用Qhull库默认启用合并共面点、处理精度误差的选项。你可以通过设置DT.Constraints来把表面三角形作为约束嵌入让剖分尽量保住表面三角形。但对复杂STL约束可能互相冲突Qhull会抛错或产生退化四面体这也是第5章要重点排查的内容。3.3 冻结外部单元pointLocation与重心坐标裁剪这一步是原生Delaunay路线里最容易被忽略、却决定成败的环节。思路如下剖分得到的每个四面体都有一个重心若是实体内部的四面体其重心必然在表面网格围成的区域内反之外部四面体的重心在区域外。用pointLocation查询每个重心落在哪个四面体里如果某个四面体的重心在它“自己”的内部但法向朝外需要结合面的朝向综合判断。更稳妥的做法是对每个表面三角形做射线法或使用inpolyhedron类工具判断重心是否在实体内部然后据此剔除外部单元。function [tet_in, valid] keep_inside_tets(Node, tet, TR) % Node: 节点坐标 % tet: 全部Delaunay四面体 % TR: 表面三角剖分对象 % tet_in: 保留的内部四面体 % valid: 逻辑向量标记原四面体是否通过检查 % 计算每个四面体的重心 P (Node(tet(:,1),:) Node(tet(:,2),:) ... Node(tet(:,3),:) Node(tet(:,4),:)) / 4; % 用表面三角剖分做点在多面体内部判定 % 原理: 以重心为端点朝任意方向发一条射线 % 统计与表面三角形交点数奇数为内部 valid false(size(tet,1),1); for i 1:size(tet,1) % 射线方向取(1,1,1)避开与三角形共面的退化情形 valid(i) point_in_polyhedron(TR, P(i,:), [1 1 1]); end tet_in tet(valid,:); fprintf(保留内部四面体: %d / %d\n, size(tet_in,1), size(tet,1)); endpoint_in_polyhedron的实现可以用Matlab File Exchange上的成熟函数或者自己写遍历表面三角形用Möller–Trumbore算法求射线与三角形交点统计与射线方向同向的交点个数。需要特别注意射线恰好通过三角形边或顶点的情况工程上常用的规避办法是对三个坐标轴分别发三条射线取多数表决结果。剔除外部单元后网格就真正“贴着”几何表面了但表面单元形状可能仍然很糟糕——因为内部节点的位置完全被Delaunay准则支配没有考虑表面三角形对单元形状的约束。这就是下一章的主题质量评估与密度控制。4. 四面体质量评估与网格密度控制别让一个坏单元毁掉一次求解4.1 用半径比和最小二面角量化单元质量四面体质量没有唯一标准但业界最常用的是半径比radius-edge ratio和最小二面角。半径比定义为四面体外接圆半径与最短边长的比值。正四面体的半径比最小约为0.612比值越大单元越扁。另一个直观指标是最小二面角六个二面角中最小者正四面体约70.5度有限元分析里通常要求不低于10度最好大于20度。function qual tet_quality(Node, tet) % 计算每个四面体的质量指标 % 返回结构体: radius_ratio, min_dihedral(度) n size(tet,1); qual.radius_ratio zeros(n,1); qual.min_dihedral zeros(n,1); for i 1:n p Node(tet(i,:), :); % 边向量 e1 p(2,:) - p(1,:); e2 p(3,:) - p(1,:); e3 p(4,:) - p(1,:); e4 p(3,:) - p(2,:); e5 p(4,:) - p(2,:); e6 p(4,:) - p(3,:); edges_len [norm(e1) norm(e2) norm(e3) ... norm(e4) norm(e5) norm(e6)]; hmin min(edges_len); % 体积: 用行列式符号可用于检测负体积 vol abs(det([e1; e2; e3])) / 6; % 外接球半径: 通过求解线性方程组得到球心 A 2 * [e1; e2; e3]; rhs [e1*e1; e2*e2; e3*e3]; center A \ rhs; R norm(center); qual.radius_ratio(i) R / hmin; end % 最小二面角: 用四面体每个面的法向量求夹角 % 详情略用r min(qual.radius_ratio)近似筛查 end参数经验值如下表你可以直接拿来做质量门禁半径比最小二面角评价建议动作 1.2 40°优秀直接用于求解1.2 ~ 1.820° ~ 40°可接受可求解误差略大1.8 ~ 2.510° ~ 20°勉强建议局部加密或平滑 2.5 10°较差必须修复否则求解器报错4.2 网格密度怎么控尺寸函数从简单到实用很多初学者的误区是“网格越密越准”实际上一味加密会让计算量爆炸而真正聪明的做法是让网格密度跟随解的梯度变化。剖分阶段我们能控制的是节点的空间分布。最简单的控制是全局尺寸在剖分前用meshgrid在包围盒内生成密度均匀的候选点再把表面节点和内部候选点合在一起做Delaunay。但这种方法对薄壁圆筒这类几何没什么用——壁厚方向一个单元都放不下。实用做法是引入尺寸函数size function每个空间位置定义该处允许的最大边长。有了尺寸函数就能在剖分前用泊松盘采样或蓝噪声采样生成非均匀分布的节点。以下代码演示如何生成一个基于到表面距离的尺寸函数并用它控制采样密度function P sample_with_size_function(TR, hmin, hmax, grow_rate) % TR: 表面三角剖分 % hmin: 表面附近最小边长 % hmax: 最远位置最大边长 % grow_rate: 从表面向外边长增长率(1.1 ~ 1.4) % 先取表面节点作为采样基础 bndPts TR.Points; % 计算表面节点到其余采样点的距离 % 用KDTree检索每个候选点距离表面的最近距离 kdtree KDTreeSearcher(bndPts); xmin min(bndPts(:,1)); xmax max(bndPts(:,1)); ymin min(bndPts(:,2)); ymax max(bndPts(:,2)); zmin min(bndPts(:,3)); zmax max(bndPts(:,3)); % 在包围盒内生成候选网格 [X,Y,Z] meshgrid(xmin:hmin:xmax, ymin:hmin:ymax, zmin:hmin:zmax); cand [X(:) Y(:) Z(:)]; % 逐点计算允许尺寸 [idx, dist] knnsearch(kdtree, cand); h_allow hmin (hmax - hmin) * (1 - exp(-dist / hmin / grow_rate)); % 按允许尺寸做空间均匀化采样: 简单方式是按比例丢弃候选点 keep rand(size(cand,1),1) (hmin ./ h_allow).^3; P cand(keep,:); % 把表面节点加回来保证表面几何被完整保留 P [P; bndPts]; end这段代码里采样密度与h_allow的三次方成反比——三维空间里边长减半单元数要翻8倍这个指数关系决定了加密的成本极高。所以尺寸函数的设置一定要考虑计算资源不要在远离关注区域的地方也加密。4.3 质量太差时的补救局部加密与单元翻转生成完网格后如果质量不达标最常见的补救手段是拉普拉斯平滑把每个内部节点移到其相邻节点坐标的平均位置。它不能修正拓扑错误但能把半径比从2.5以上大幅改善到1.5左右。代码如下function Node_new laplacian_smooth(Node, tet, iter) % Node: 节点坐标 % tet: 四面体连接关系 % iter: 平滑迭代次数 Node_new Node; % 建立节点到单元的邻接关系 for k 1:iter % 累加每个节点相邻节点的坐标和 node_sum zeros(size(Node_new)); node_cnt zeros(size(Node_new,1), 1); for i 1:size(tet,1) for j 1:4 v tet(i,j); node_sum(v,:) node_sum(v,:) ... sum(Node_new(tet(i,:),:), 1) - Node_new(v,:); node_cnt(v) node_cnt(v) 3; % 每个四面体有3个邻居 end end % 只移动内部节点表面节点保持固定 Node_new Node_new; interior setdiff(1:size(Node,1), unique(TR.ConnectivityList(:))); Node_new(interior,:) node_sum(interior,:) ./ node_cnt(interior); end end平滑时要特别注意移动节点后单元可能发生翻转体积变负。每次迭代后都应重算一次最小二面角如果质量指标变差就回退到上一轮坐标。另一个手段是边翻转edge flip——两个相邻四面体构成的凸六面体内可以交换共享面来改善最大二面角但实现复杂度较高工程上常用现成的网格优化库Matlab生态里可用的不多遇到极端质量问题时我一般直接导出到TetGen用它的-q选项优化。5. 避开Matlab四面体剖分的五个经典坑5.1 表面网格不封闭STL的缝隙和重叠面四面体剖分要求输入表面是一个封闭的二维流形即每条边恰好被两个三角形共享。STL从CAD导出时经常出现裂缝、重叠面和非流形边。先做诊断再剖分function diagnose_surface(TR) % 统计每条边被多少个三角形共享 edges TR.edges; % 每条边对应的两个节点 n_edges size(edges,1); % 构建边到三角形的映射 tri_nodes TR.ConnectivityList; edge_tag zeros(n_edges,1); for i 1:size(tri_nodes,1) t sort([tri_nodes(i,1) tri_nodes(i,2); tri_nodes(i,2) tri_nodes(i,3); tri_nodes(i,1) tri_nodes(i,3)], 2); for e 1:3 % 用unique行匹配来统计 end end % 实际实现可直接用triangulation的edgeAttachments [~, ~, att] edgeAttachments(TR); cnt cellfun(numel, att); fprintf(每条边被三角形共享次数分布: min%d, max%d\n, min(cnt), max(cnt)); if min(cnt) 2 error(存在边界边表面未闭合); end end如果存在边界边需要回到CAD软件补面或在Matlab里用patch的手工编辑修复。不要尝试用Delaunay剖分去“自动补洞”那只会生成一堆穿出表面的劣质单元。5.2 单位不统一几何尺寸与参数比例失调Qhull处理的是浮点坐标其内部容差基于坐标量级。如果你的模型坐标在1e-5量级或横跨1e-5到1e6的多个数量级精度误差会直接导致Qhull precision error。解决方法是剖分前把所有坐标平移到原点、缩放到最大边长1附近剖分后再把单元坐标缩放回去。这一条被很多人忽略却是单元质量差的最常见原因。5.3 Delaunay约束冲突与Qhull报错当你把TR.ConnectivityList作为Constraints属性传给delaunayTriangulation时约束面之间可能存在小角度相交或面与面几乎重合。Qhull会警告QH6214或直接出错。应对策略是按4.2节先做表面修复把过于接近的三角形合并或者放弃原生约束改用第3.3节的剔除方案。后者牺牲一点性能但鲁棒性好得多。5.4 判定函数写错重心法在凹形区域失效第3.3节用重心是否在内部来剔除外部四面体这个策略对凸形体完全正确对凹形体也基本正确但有一个陷阱如果重心恰好在表面附近由于浮点误差判定结果可能不稳定。做判定前把模型整体放大1000倍或对重心坐标加上一个远离表面的微扰能显著降低误判率。更可靠的做法是统计表面三角形法向所有法向朝外时内部点发出的任意射线与表面交点是奇数个工程上用三条不同射线同时判定。5.5 负体积单元方向约定不一致Delaunay剖分生成的四面体通常满足右手定则节点的顺序使体积为正但在做局部加密或平滑操作后某些节点的移动可能使单元翻转。负体积单元在有限元装配时会导致刚度矩阵出现负的特征值求解器直接报错。每次修改网格后都要检查带符号体积function is_negative check_negative_volume(Node, tet) % 计算带符号体积正为右手系负为翻转 v zeros(size(tet,1),1); for i 1:size(tet,1) p Node(tet(i,:), :); v(i) det([p(2,:)-p(1,:); p(3,:)-p(1,:); p(4,:)-p(1,:)]) / 6; end is_negative v 0; end这条检查也要加入你的网格流水线末端和单元质量评估一起作为输出门禁。很多人的网格在视觉上毫无问题一进求解器就报“Element distorted”根因就是这里。6. 从剖分到求解把四面体网格交给有限元装配器6.1 导出为求解器可读格式写出TetGen/Abaqus风格文件网格生成不是终点真正要跑的是装配和求解。Matlab的assembleFEMatrices能直接吃网格对象但更通用的做法是导出为文本格式供外部求解器使用。下面是导出为TetGen.node和.ele格式的代码function export_tetgen(filename_base, Node, tet) % 导出为TetGen格式 fid_node fopen([filename_base .node], w); fprintf(fid_node, %d 3 0 0\n, size(Node,1)); fprintf(fid_node, %d %g %g %g\n, ... [(0:size(Node,1)-1) Node]); fclose(fid_node); fid_ele fopen([filename_base .ele], w); fprintf(fid_ele, %d 4 0\n, size(tet,1)); fprintf(fid_ele, %d %d %d %d %d\n, ... [(0:size(tet,1)-1) tet-1]); fclose(fid_ele); end注意这里的节点编号从0开始而Matlab里索引从1开始所以输出时要减1。四面体四个节点的顺序要和求解器的约定一致大多数有限元求解器要求单元顶点按右手系排列顺序不对体积就是负的——前面检查负体积的代码在这里复用。6.2 节点编号重排带宽优化与Cuthill-McKee稀疏线性方程组的求解效率高度依赖刚度矩阵的非零元分布。默认剖分编号的矩阵带宽往往很大用Matlab的symrcm做一次节点重排能显著压缩带宽% 构建单元连接图的邻接矩阵 A sparse(size(Node,1), size(Node,1)); for i 1:size(tet,1) t tet(i,:); A(t, t) 1; end % 对称Cuthill-McKee重排 perm symrcm(A); new_node Node(perm, :); new_tet zeros(size(tet)); for i 1:numel(perm) new_tet(tet perm(i)) i; end重排不影响数值结果但能大幅降低直接法分解的填充量和迭代法的收敛时间。对10万节点以上的网格这一步带来的加速可能达到2到5倍值得固化成剖面流程的一个标准环节。6.3 自适应细化闭环用一次粗糙解驱动下一次剖分四面体剖分的终极形态是自适应先跑一次粗糙网格的计算得到误差指示子比如相邻单元应力梯度过大处在误差大区域插入节点重新剖分再求解。这个闭环在Matlab里的实现流程是求解→计算每个四面体的误差指示子→提取高误差四面体的中心点作为新增采样点→合并进节点集重新做Delaunay→回到第4章做质量检查。整个过程可以用一个while循环驱动直到最大误差低于阈值。这里有一个经验值可以参考每次只加密误差最大的前20%单元迭代3到5次比一次性全量加密获得同样精度要省一半以上的节点数。做的时候记得把尺寸函数也同步更新否则只加节点不更新密度约束剖分器会把新增节点“平均”掉白白浪费计算量。本文还有配套的精品资源点击获取