trullekrul各向异性网格自适应拓扑优化与有限元装配

简介：这套开源 MATLAB/Octave 脚本面向从事拓扑优化、各向异性网格自适应与有限元分析的研究生、科研人员和工程开发者，围绕二维/三维最佳导热与线性弹性最小柔度问题提供可运行的实现参考。包内共 5 个文件，以 3 个 .m 脚本为主体，另含 .md 与 .rst 说明文档，压缩包约 173KB，体积轻量却保留了完整的算例调用方式与参数入口，便于直接改参复现实验。已有 324 人学习，说明其在相关课程与课题中具备一定参考价值。读者可借助示例脚本理解各向异性网格自适应如何与拓扑优化耦合，掌握导热、最小柔度两类目标的参数组织与结果复现思路，并依据说明文档快速定位环境要求与许可信息，适合作为有限元与优化课程的辅助材料或科研原型起点。

1. 拓扑优化出图全是锯齿之后：我为什么把 trullekrul 装进求解链

第一次用密度法跑悬臂梁，迭代 80 步收敛得挺漂亮，导出图一看，结构边界像被狗啃过——棋盘格、灰度过渡带、边界一圈锯齿。换加密网格，拓扑变了；再换回来，又变一个样。问题不在优化器，在网格：各向同性均匀网格既喂不饱边界的高曲率，又在平坦区浪费自由度。trullekrul 这套 MATLAB/Octave 脚本正好对着这个痛点——它把各向异性网格自适应、拓扑优化、有限元装配三件事放在同一个仓库里，用度量张量驱动网格在边界拉长加密、在内部粗化，并且整套脚本是开源的，改一行就能看到一行的影响，没有黑匣子。适合手上有 MATLAB 或 Octave、想自己掌控求解链而不是调黑盒 API 的结构优化从业者。

2. 环境与冒烟测试：MATLAB/Octave 双栈跑通的三步闭环

拿到一个脚本包，最忌讳上来就run example。脚本类资源没有install、没有工具箱签名校验，出错信息往往是一句Undefined function，你连它到底缺哪一层都看不出来。这套代码同时标了 MATLAB 和 Octave，意味着作者刻意避开了专有工具箱依赖，但也意味着你在两边会踩到不同的坑：Octave 的稀疏求解器后端、图形后端、路径大小写规则都和 MATLAB 不一样。所以我把第一步拆成三件事：自检解释器、挂对路径、跑通一个最小闭环。

2.1 先自检，不要先跑算例

拓扑优化的墙钟时间八成花在稀疏求解那一行反斜杠上。如果这一步在你的环境里退化（比如稀疏矩阵被稀疏化失败，走了稠密求解），后面所有"为什么不收敛"的排查都是白费。下面这个自检脚本我一般放在tools/下，跑算例前先过一遍。

% tools/check_env.m —— 三项自检，全过再谈跑算例 function check_env() % 1) 解释器识别：两个环境的部分内置函数不同名，先分清自己在哪边 isOctave = exist('OCTAVE_VERSION', 'builtin') ~= 0; if isOctave fprintf('Interpreter : Octave %s\n', OCTAVE_VERSION); else fprintf('Interpreter : MATLAB %s\n', version); end % 2) 稀疏求解：拓扑优化的主耗时项，必须确认没退化成稠密 n = 2000; A = spdiags([-ones(n,1) 4*ones(n,1) -ones(n,1)], -1:1, n, n); b = ones(n,1); tic; x = A \ b; tSolve = toc; assert(norm(A*x - b) < 1e-10, 'sparse backslash failed'); fprintf('Sparse solve: %.3f s, nnz(A) = %d\n', tSolve, nnz(A)); if tSolve > 1 warning('稀疏求解偏慢，检查是否误用 full() 或矩阵病态'); end % 3) 图形后端：无显示环境下 patch 会直接抛错，但只影响出图不影响求解 try h = figure('visible', 'off'); patch('Faces', [1 2 3], 'Vertices', [0 0; 1 0; 0 1], 'FaceColor', 'flat'); close(h); fprintf('Graphics : ok\n'); catch err fprintf('Graphics : %s（只影响可视化，求解不受影响）\n', err.message); end end

代码里的三个判断各有指向：exist('OCTAVE_VERSION','builtin')是跨平台识别解释器最稳的写法，比判断matlabroot可靠；稀疏求解那一段构造的是标准三对角刚度矩阵，规模 2000×2000，正常应该在毫秒到几十毫秒量级，如果超过 1 秒，八成是你所在的环境里A\b走了稠密分解，或者是包里的装配函数返回了full矩阵——这一点在后面的装配章节会再撞见；图形部分故意用try/catch包住，是因为服务器上跑批量算例时 headless 环境报错会直接中断整个脚本，而出图本来就不是求解链的必要环节。

2.2 路径挂载：genpath 的坑比你想的多

脚本包没有安装流程，靠addpath生效。常见的做法是一行addpath(genpath(root))完事，但这一行在两类场景下会翻车：一是包里带.git、.svn目录，genpath会把它们一起挂进搜索路径，路径变量膨胀，启动明显变慢；二是包内不同子目录存在同名函数（比如mesh/和adapt/各有一个edge_length），MATLAB 按路径顺序取第一个命中，谁在前谁生效，出问题时表现是"函数存在但行为不对"，非常难查。

% run_example.m —— 所有算例的统一入口 root = fileparts(mfilename('fullpath')); % 只挂一层源码目录，避免 genpath 把 .git / 备份目录也收进来 srcDirs = {'mesh', 'fem', 'adapt', 'topopt', 'io', 'tools'}; for k = 1:numel(srcDirs) d = fullfile(root, 'src', srcDirs{k}); if isfolder(d), addpath(d); end end addpath(fullfile(root, 'examples')); check_env(); % 同名函数排查：命中多于一个就打印出来 cand = which('edge_length', '-all'); if numel(cand) > 1 fprintf('重复定义 edge_length：\n'); fprintf(' %s\n', cand{:}); end

which -all这一行值得养成习惯。我一般会在排查阶段对包内所有函数名跑一遍重复检测，把命中数大于 1 的列出来，一次性把遮蔽问题清掉，比等到结果诡异时再回头找要省事得多。另外注意isfolder在较早版本的 MATLAB 里叫isdir，如果你手上的解释器版本偏老，改成exist(d, 'dir') == 7更保险。

2.3 最小闭环：网格、装配、求解三段式

真正验证"这套脚本能用"的标准，不是它能出图，而是它能完成一次从几何到柔度数值的完整链路。我给自己的验收脚本是一个矩形域上的弹性静力问题，规模小到秒级出结果，但走的函数路径和大算例完全一致。

% examples/smoke_test.m —— 最小闭环，验证三段式是否打通 Lx = 2; Ly = 1; nx = 80; ny = 40; % 1) 网格：结构化四边形，节点坐标 Nx2，单元连接 Nx4 [node, elem] = rect_mesh(Lx, Ly, nx, ny); % 2) 材料与本构：归一化参数，E 取 1 方便看相对量 mat.E0 = 1.0; mat.nu = 0.3; mat.Emin = 1e-9; [D, Ke] = fem_matrices(mat, elem); % D: 3x3 本构, Ke: 8x8xNel 单元刚度 % 3) 装配：输出必须保持稀疏 K = fem_assemble(node, elem, Ke); assert(issparse(K), '装配返回了稠密矩阵，大算例会直接爆内存'); % 4) 边界条件：按坐标找节点，别缓存索引 fixed = find(abs(node(:,1)) < 1e-9); % 左端固支 tip = find(abs(node(:,1) - Lx) < 1e-9); % 右端加载 F = sparse(2*size(node,1), 1); F(2*tip-1) = 1 / numel(tip); % 均布剪力，避免单点应力奇异 % 5) 求解并统计柔度 U = solve_dirichlet(K, F, fixed); c = full(F' * U); fprintf('dof = %d, compliance = %.6e\n', size(K,1), c);

这段代码里的函数名是按功能归类的通用命名，你手上的包里可能叫别的名字，对照目录找语义最近的即可——网格生成在mesh/下，本构与单元刚度在fem/下，装配和约束施加也在fem/下，这套职责划分在同类脚本里基本是共识。

几个参数值得说清。nx=80, ny=40是个刻意的选择：80×40 的四边形网格是 3280 个节点、6560 个自由度，稀疏直接求解在一秒内能出结果，同时单元数又足够让后面的自适应算法有东西可做，太小了看不出度量张量的效果。mat.Emin = 1e-9是 SIMP 里的空洞刚度下限，这里先设上是为了让冒烟测试和后面的优化算例共用同一份材料结构，但它本身就是个坑——Emin 太小会让刚度矩阵条件数飙升，稀疏求解器的数值精度会掉，后面避坑章节会专门讲。右端加载用1/numel(tip)做均布而不是在单点施加集中力，是因为单点集中力在 Q1 单元上会产生局部应力奇异，柔度值会随网格加密单调发散，你用它做网格无关性验证永远收敛不了。

三段式跑通的标准是：柔度数值有限、不为负、量级在 1e-2 到 1e2 之间（归一化材料参数下），并且把nx翻倍后柔度变化在个位数百分比。如果这三条里任何一条不满足，先别往下走自适应，回去查装配的稀疏性和约束施加方式。

3. 各向异性网格自适应：从 Hessian 到度量张量的落地步骤

各向同性自适应只回答一个问题：这里误差大不大。各向异性自适应多回答两个：在哪个方向大、大多少倍。这两问的答案就是度量张量。理解这一点，后面所有参数才有落脚处：为什么单元可以被拉成 10:1 的长条而不影响精度，为什么长宽比不能无限拉，为什么单元数预算是靠行列式控制的。

3.1 各向异性优于各向同性的前提条件

结论先说：只有当解的二阶导在两个主方向上量级差得足够远时，各向异性才有收益。典型场景是薄壁结构、带边界层的场、剪切主导的位移场——这些地方解沿壁面方向几乎线性，沿法向剧烈变化，各向同性加密会在壁面方向白白堆单元。反过来，如果解的 Hessian 两个特征值量级接近（比如各向同性的点载荷附近），强行各向异性只会产生一堆瘦长单元，最小角掉到 10 度以下，刚度矩阵条件数爆炸，求解精度反而下降。

所以开工前我会先做一件事：在初始网格上解一次，把每个节点的 Hessian 特征值比打出来看分布。

% 诊断：Hessian 各向异性比的分布，决定要不要开各向异性 H = recover_hessian(node, elem, U); % 每个节点一个 2x2 对称阵 ratio = zeros(size(node,1), 1); for i = 1:size(node,1) ev = eig(0.5*(H(:,:,i) + H(:,:,i)')); % 数值 Hessian 常不对称，先对称化 ev = sort(abs(ev), 'descend'); ratio(i) = ev(1) / max(ev(2), 1e-12); end fprintf('各向异性比: 中位数 %.2f, 90分位 %.2f, 最大 %.2f\n', ... median(ratio), prctile_safe(ratio, 90), max(ratio));

判读标准很朴素：中位数在 1 到 3 之间、90 分位超过 20，说明场里存在明显的方向性结构，各向异性值得开；如果 90 分位只有 5 左右，老老实实用各向同性，省下的调试时间比省下的单元多。

3.2 度量张量的构造与等分布归一化

度量张量的定义是：一条边向量e在度量下的长度平方为e' * M * e，当这个值约等于 1 时，这条边就是"理想边长"。要让沿曲率大的方向边长变短，M在该方向的特征值就要大，因此M的特征值取 Hessian 特征值的倒数，方向取 Hessian 的特征向量。

% adapt/metric_from_hessian.m function M = metric_from_hessian(H, hmax, hmin, targetEdges, anisoCap) % H : 2x2xN，逐节点 Hessian % hmax/hmin : 允许的最大/最小边长 % targetEdges: 期望的总边数，用于行列式归一化 % anisoCap : 各向异性比上限，默认 1e4 if nargin < 5, anisoCap = 1e4; end n = size(H, 3); M = zeros(2, 2, n); for i = 1:n Hi = 0.5 * (H(:,:,i) + H(:,:,i)'); % 对称化 [V, D] = eig(Hi); d = abs(diag(D)); d = max(d, 1e-12); % 平坦区防零：否则边长趋于无穷 % 各向异性比截断：超过上限后单元退化，收益为负 d = max(d, max(d) / anisoCap); M(:,:,i) = V * diag(1 ./ d) * V'; end % 等分布归一化：常数场下总单元数由 det(M)^(1/2) 控制 for i = 1:n Mi = M(:,:,i); Mi = Mi / sqrt(det(Mi)) * targetEdges; % 2D 指数是 1/2，3D 是 1/3 [V, D] = eig(Mi); ev = diag(D); ev = min(ev, 1 / hmin^2); % 边长下限 ev = max(ev, 1 / hmax^2); % 边长上限 M(:,:,i) = V * diag(ev) * V'; end end

三处细节值得掰开。

为什么取倒数而不是直接用特征值。度量张量的物理含义是"把物理空间映射到参考空间"，映射后每条边都应该是单位长。Hessian 特征值大意味着该方向曲率大、误差大，需要更短的边，所以特征值取倒数是直接结论，不是调参技巧。

为什么必须做行列式归一化。不做归一化，M的整体尺度是自由的，你把所有特征值乘 2，边长全变短一半，单元数涨 4 倍，而相对分布一点没变。归一化的作用是把"单元数预算"这一维从分布问题里解耦出来：det(M)^(1/2)在 2D 下的积分近似等于总单元数，所以用targetEdges去除，就能让每轮的单元数落在你预期的量级。这一步要是不做，自适应跑两轮网格就爆炸，是新手最常见的翻车点。

为什么要有各向异性比上限。理论上比可以无穷大，实际上超过 1e4 后，单元雅可比的行列式在浮点精度下已经不可靠，刚度矩阵会出现接近奇异的行，稀疏求解器给的解带噪声，敏度算出来符号都可能错。我一般用 1e3 到 1e4，薄壁特别极端的时候才放到 1e5，且必须配合后面章节的网格质量检查。

3.3 网格操作三步：分裂、塌缩、翻转

有了度量，网格操作就是照着"边长应该接近 1"这个目标修正拓扑。标准三步走，顺序不能乱。

% adapt/adapt_loop.m —— 一轮自适应的骨架 for it = 1:opt.maxAdapt % 1) 在新网格上重解物理场 U = solve_dirichlet(K, F, fixed); % 2) 恢复 Hessian 并构造度量 H = recover_hessian(node, elem, U); M = metric_from_hessian(H, opt.hmax, opt.hmin, opt.targetEdges, opt.anisoCap); % 3) 边分裂：度量长度 > 1.3 判为长边 [node, elem, markSplit] = edge_split(node, elem, M, 1.3); % 4) 边塌缩：度量长度 < 0.5 判为短边，先做分裂再做塌缩可避免重复操作 [node, elem] = edge_collapse(node, elem, M, 0.5); % 5) 局部翻转：优化最小角，只翻不改变节点数 [node, elem] = edge_flip(node, elem, M); % 6) 收敛判据：被修改单元占比 changed = nnz(markSplit) / size(elem, 1); fprintf('adapt %d: nel = %d, changed = %.3f\n', it, size(elem,1), changed); if changed < opt.adaptTol, break; end % 7) 拓扑变了，几何量全部失效，重装配 [~, Ke] = fem_matrices(mat, elem); K = fem_assemble(node, elem, Ke); fixed = find(abs(node(:,1)) < 1e-9); % 坐标重算，不缓存旧索引 end

阈值 1.3 和 0.5 不是拍脑袋。等分布理论里，一条边的理想度量长度是 1，超过某个上界就该分裂、低于某个下界就该塌缩，1.3/0.5 这组值来自"分裂一次后两条子边都落在 0.65 到 1.0 之间"的几何约束——如果上界取 1.9，分裂后子边长度约 0.95，看似也行，但实际会让单元数增长过猛且新边立刻又接近上界，产生级联分裂。1.3 是经验上让每轮单元数增长控制在 30% 以内的保守值。

fixed那行注释值得强调：塌缩和分裂会改变节点编号，任何在自适应之前缓存的自由度索引都会指向错误的节点。约束集合必须在每轮之后按坐标重新计算，这是血泪经验，后面避坑章还会再提一次。

3.4 质量指标：什么时候该停下来

自适应不是跑得越久越好。度量构造有问题时，循环会持续产生退化单元而单元数不再增长，看起来"收敛了"，实际上网格已经烂了。每轮结束我固定看四个数。

四项里前两项是硬指标，后两项是观测量。我见过最典型的失败是：归一化那一步的指数写错（2D 用成 1/3），边长中位数一直停在 0.4 附近，循环永远不收敛，单元数每轮涨一倍，三轮之后内存就不够了。指数这个细节在不同维度的公式里不一样，抄代码的时候特别容易漏。

4. SIMP 主循环：密度、敏度过滤与自适应网格的耦合顺序

拓扑优化的主循环本身不复杂，密度法（SIMP）的公式在文献里翻来覆去就那几行，真正决定成败的是过滤放在哪一步、自适应什么时候插进来、移动限怎么调。这三件事的顺序错了，收敛曲线会呈现出"柔度先降后升、密度场变成雪花"的经典症状。

4.1 问题定义与参数取值

标准的最小柔度问题：

min c(rho) = F' U s.t. K(rho) U = F sum(Ve .* rho) / sum(Ve) <= volfrac rho_min <= rho <= 1

材料插值用 SIMP：E(rho) = Emin + rho^p * (E0 - Emin)。参数表如下，这张表我贴在显示器边上用了很久。

Emin = 1e-9这行看着无害，实际是数值稳定性的关键。取 1e-9 时 E0/Emin = 1e9，刚度矩阵条件数在 1e9 量级，双精度下还有 7 位有效数字可用；有些教程拿 1e-12，条件数到 1e12，稀疏求解器给的解末几位全是噪声，敏度跟着抖，优化过程表现为柔度曲线上蹿下跳。我一般从 1e-9 起步，只有在结构里出现大量"漂浮"的孤立单元时才考虑调小。

4.2 敏度分析与过滤的落点

% topopt/simp_loop.m function [rho, cHist] = simp_loop(node, elem, mat, bc, opt) nel = size(elem, 1); rho = opt.volfrac * ones(nel, 1); % 均匀初始化，比随机初值更稳 [~, Ke] = fem_matrices(mat, elem); % 满材料单元刚度，8x8xNel edof = elem_dof_table(elem); % 单元自由度索引表 H = density_filter(node, elem, opt.rmin); % 过滤矩阵只算一次，别放循环里 Ve = elem_volume(node, elem); % 单元体积，体积约束用 cHist = zeros(opt.maxIter, 1); for it = 1:opt.maxIter % ---- 1) 材料插值与装配 ---- E = mat.Emin + rho.^opt.penal * (mat.E0 - mat.Emin); K = fem_assemble_sparse(edof, Ke, E, 2*size(node,1)); U = solve_dirichlet(K, bc.F, bc.fixed); cHist(it) = full(bc.F' * U); % ---- 2) 单元应变能 -> 敏度 ---- ce = zeros(nel, 1); for e = 1:nel ue = U(edof(e, :)); ce(e) = ue' * Ke(:, :, e) * ue; % 单元应变能，永远 >= 0 end dc = -opt.penal * rho.^(opt.penal-1) .* (mat.E0 - mat.Emin) .* ce; % ---- 3) 敏度过滤（不是密度过滤）---- dc = H * (rho .* dc) ./ max(H * rho, 1e-3); % ---- 4) 体积敏度与 OC 更新 ---- dv = Ve; l1 = 0; l2 = 1e9; while (l2 - l1) > 1e-6 lmid = 0.5 * (l1 + l2); B = -dc ./ (lmid * dv); rhoNew = max(0, max(rho - opt.move, ... min(1, min(rho + opt.move, rho .* sqrt(B))))); if sum(Ve .* rhoNew) > opt.volfrac * sum(Ve) l1 = lmid; else, l2 = lmid; end end rho = rhoNew; % ---- 5) 收敛判据 ---- if it > 10 && max(abs(rho - rhoNew)) < 1e-3, break; end end end

这里有几个容易被忽略的判断。

敏度dc一定是负的。柔度对密度的导数是负值（加材料降低柔度），OC 更新里用-dc变成正值才能开方。如果符号搞反，第一轮迭代柔度就往上涨，后面再也降不下来。我一般会在循环里加一句assert(mean(dc) < 0)作为一次性检查，跑通几次后再删掉。

过滤放在敏度上，别放在密度上。密度过滤会把灰度过渡区整体拉宽，最终结构边界模糊；敏度过滤只平滑灵敏度场，黑白分离度明显更好。这是两类过滤最容易混淆的地方，很多教程默认给的是密度过滤，直接照搬会发现结构边界一圈灰。

过滤矩阵H只算一次。它只依赖节点坐标和 rmin，固定网格下是常量。我见过把它放在迭代循环里的实现，200 次迭代里重复计算 200 遍稀疏矩阵乘法，跑得慢到让人怀疑人生。

lmid是拉格朗日乘子。OC 的二分法在找一个满足体积约束的乘子，sqrt(B)里的开方来自rho * sqrt(-dc/(lambda*dv))这个解析式。二分上下界l1=0, l2=1e9覆盖了绝大多数情形，如果你的体积约束特别紧（volfrac < 0.1），把l2再放大一两个量级。

4.3 自适应网格什么时候插进来

这是整篇文章最需要强调的一句话：不要每一轮迭代都做网格自适应。原因是拓扑优化的目标函数在网格变化时会发生跳变——新网格上的离散柔度天然和旧网格不同，OC 里基于上一轮密度的移动限假设被打破，表现为柔度曲线在每次自适应后出现一个台阶，多次叠加就不收敛了。而过滤矩阵H、单元体积Ve、自由度索引edof全部依赖网格，一变就得重算，开销也不小。

我的做法是分阶段：

% 阶段一：固定网格上跑到大致收敛（比如 max|drho| < 5e-3） % 阶段二：进入外层循环，每 20~30 次 SIMP 迭代做一次网格自适应 for outer = 1:opt.maxAdapt [rho, cHist] = simp_loop(node, elem, mat, bc, optInner); % 内层跑 20~30 步 % 自适应前把密度场映射到新网格（按单元中心最近邻或面积加权） [node, elem, rhoNew] = adapt_and_project(node, elem, rho, opt); % 重置移动限：新网格上前几步要收着走 optInner.move = 0.1; optInner.rhoInit = rhoNew; end

adapt_and_project那一步是密度场投影：新单元的中心落在旧网格的某个单元里，直接取该旧单元密度，或者按面积加权平均。直接最近邻会引入阶梯状误差，面积加权更平滑但会额外产生一些中间密度——这两种取舍都在合理范围内，关键是投影后必须重算一遍sum(Ve.*rho)/sum(Ve)，把体积分数校准回volfrac，否则自适应会让体积约束漂移。

移动限的临时下调很有效。新网格上的敏度场整体偏移，用 0.2 的移动限容易一步跨过头；前两次迭代用 0.1，等密度场适应新网格后再恢复到 0.2，柔度曲线的台阶会小很多。

4.4 结果判读的三个信号

迭代结束后，我先看三样东西，而不是直接看结构图。

柔度收敛曲线应该是单调下降后趋于平台，末端有轻微振荡但幅度小于 1%。如果曲线在中段出现持续上升，几乎可以确定是过滤半径小于单元尺寸或者移动限过大。

体积分数闭合：sum(Ve.*rho)/sum(Ve)与volfrac的差应该在 1e-3 以内。差得大说明 OC 二分没有收敛，把l2上界放宽或者把二分容差从 1e-6 收紧到 1e-8。

灰度单元占比：统计0.1 < rho < 0.9的单元比例，健康值在 5% 以下。超过 15% 说明惩罚不够或者过滤过强，先调penal到 3.5 试试，别急着改网格。

5. 避坑排查：五个让算例当场翻车的现场记录

这套脚本我断断续续用了很久，下面五条是反复踩到、每次都浪费半天以上的。按"现象 → 原因 → 解决"写，方便对照。

5.1 函数明明在目录里，却报 Undefined function

现象：addpath之后执行rect_mesh报未定义，但文件确实存在。

原因：两种可能。一是 Linux 下文件名大小写敏感，文件叫Rect_Mesh.m，路径里写rect_mesh；二是包内存在同名函数被遮蔽，which只显示第一个命中，实际调用的是另一个目录下行为不同的版本。

解决：先which -all rect_mesh看命中数，多于一小时按路径顺序判断谁生效；大小写问题用ls逐级核对实际文件名。这个排查我一般会提前到第 2 章的自检脚本里做掉。

5.2 迭代 20 步柔度不降反升，密度场像雪花

现象：柔度曲线前 10 步下降，之后开始震荡上行，导出密度云图是黑白相间的碎斑。

原因：敏度过滤的权重矩阵H用了rmin小于 1.5 倍单元尺寸，棋盘格没有被抑制；或者敏度符号反了，dc算成了正值，OC 在往错误方向更新。

解决：先在第一次迭代打印mean(dc)，确认是负值；再把rmin提到 2 到 3 倍单元尺寸。如果rmin已经够大还抖，检查过滤公式是不是漏了./ max(H*rho, 1e-3)里的归一化分母，缺了它会让低密度区被过度平滑。

5.3 自适应跑两轮内存爆掉

现象：第一轮自适应正常，第二轮开始内存占用翻倍增长，最终Out of memory。

原因：塌缩操作后节点编号不连续，稀疏矩阵按max(node id)开维度，产生了大量空行；同时旧的K、edof、H在循环里没有释放，每轮都在累加。

解决：每次塌缩后做一次节点压缩，把废弃编号挤掉并重建elem的索引映射；循环内显式clear K edof H。另外装配函数返回后加一句issparse断言，我遇到过装配内部某条分支返回了full矩阵的版本，小算例看不出来，一万单元以上直接爆。

5.4 自适应之后约束加到空处，解出现刚体位移

现象：求解结果里位移量级到了 1e6，结构像飘在空中。

原因：塌缩操作把左端边界上的节点合并掉了，而fixed集合还是自适应之前缓存的旧索引，指向了已经不存在的自由度。刚度矩阵因此奇异，稀疏求解器给出的是最小二乘意义上的解，数值完全不可信。

解决：这条是硬规矩——任何网格操作之后，边界条件按坐标重新计算，不缓存索引。fixed = find(abs(node(:,1)) < tol)这行放进每一轮自适应之后，成本可以忽略。另外加个断言assert(numel(fixed) >= 2)，边界节点被塌缩光了能立刻发现。

5.5 Octave 下 patch 出图全黑，中文标注乱码

现象：MATLAB 下正常的密度云图，在 Octave 里整片黑色，坐标轴标签显示方块。

原因：默认图形后端在 headless 环境或缺少字体时渲染失败；facecolor用了'interp'在某些旧版本不被支持，退化成纯色填充。

解决：求解和出图解耦，先把结果存成.mat，再单独出图。出图时切后端graphics_toolkit('gnuplot')或者干脆用print -dpng导出后用外部工具查看。中文标注直接在脚本里换成英文，这条没什么优雅的解法，换来换去不如一开始就写英文标签。

6. 结果验收：四个数字判定这次优化到底能不能用

结构图好看不代表结果可信。我给自己定了一套验收流程，跑完优化不急着导图，先过四个数字。这套流程是从一次翻车里长出来的：那个算例的柔度曲线漂亮得不像话，导出结构给同事看，对方问了一句"网格加密一倍结构还长这样吗"，我答不上来，回去重跑，拓扑直接变了个样。

6.1 四个验收指标与判定阈值

% verify/acceptance.m —— 四个数字，全过才算可交付 function report = acceptance(node, elem, rho, bc, opt) report = struct(); % 1) 体积分数闭合 Ve = elem_volume(node, elem); report.volError = abs(sum(Ve .* rho) / sum(Ve) - opt.volfrac); % 2) 灰度单元占比 report.grayRatio = nnz(rho > 0.1 & rho < 0.9) / numel(rho); % 3) 网格无关性：整体加密一倍再跑 50 步，比较柔度 [node2, elem2] = uniform_refine(node, elem); rho2 = project_density(node, elem, rho, node2, elem2); c1 = eval_compliance(node, elem, rho, bc, opt); c2 = eval_compliance(node2, elem2, rho2, bc, opt); report.meshSensitivity = abs(c2 - c1) / c1; % 4) 边界离散误差：结构周长随自适应的变化率 report.perimeter = struct_boundary_length(node, elem, rho, 0.5); fprintf('volError = %.2e (阈值 1e-3)\n', report.volError); fprintf('grayRatio = %.4f (阈值 0.05)\n', report.grayRatio); fprintf('meshSensitivity = %.4f (阈值 0.02)\n', report.meshSensitivity); fprintf('perimeter = %.4f\n', report.perimeter); end

四个阈值的来历说清楚。volError的 1e-3 对应 OC 二分容差 1e-6 的量级，低于这个数说明二分收敛了；grayRatio的 5% 是经验值，超过 15% 说明惩罚不足或过滤过宽；meshSensitivity的 2% 是网格无关性的实际标准——注意这里比的是"加密后重新优化"的柔度，不是同一密度场在两套网格上评估，后者必然有差异，不能作为判据；perimeter没有固定阈值，看的是随自适应轮次的收敛趋势，如果周长在最后三轮还在大幅变化，说明边界还没被网格充分分辨。

6.2 网格无关性验证的具体操作

第三条最容易做错。常见错误是把优化后的密度场直接插值到细网格上，用细网格重算一次柔度，发现差 30%，然后归咎于网格依赖。这个做法有问题：插值会把灰度扩散，细网格上评估的其实是另一个结构。

正确做法是在细网格上从同一初始条件重新跑优化，得到rho2，再和rho1比较柔度。这多花一倍时间，但结论可靠。我一般只在最终交付前做一次，中间调试阶段用自适应循环的收敛性代替。

% 正确的网格无关性检查：重新优化，不是重新评估 opt2 = opt; opt2.volfrac = opt.volfrac; opt2.rmin = opt.rmin * 2; % 过滤半径随网格缩放，否则细网格上相对过滤变强 [node2, elem2] = uniform_refine(node, elem); [rho2, cHist2] = simp_loop(node2, elem2, mat, bc, opt2); c2 = cHist2(find(cHist2 > 0, 1, 'last')); fprintf('细网格柔度 %.6e vs 粗网格 %.6e, 偏差 %.2f%%\n', ... c2, c1, 100 * abs(c2 - c1) / c1);

rmin随网格缩放这一行是关键。过滤半径的物理意义是长度，网格加密一倍，如果rmin的单元数不变，物理长度就减半，相对过滤强度变了，结果当然对不上。这是网格无关性验证里最隐蔽的一处不一致。

6.3 一个我常用来省时间的习惯

出图之前先把密度场存成带坐标和单元连接的.mat，文件名带上参数组合（比如p3_rmin2_vol045_adapt3.mat）。理由很实际：同一批参数往往要跑好几组对比，中途改一个参数就得回滚重来，有存档就能直接加载对比而不必重跑。存档里我一定会存node、elem、rho、cHist、参数结构体这五项，缺任何一项，过两周自己都还原不出当时的条件。

还有个小技巧：把cHist和验收四个数字一起写进一个文本日志，每次跑完自动追加一行。攒到十几组之后，参数和结果的关系一眼就能看出来，比每次翻代码回忆参数值靠谱得多。

从那以后我每次交付拓扑优化结果，都会强制走一遍这四道验收，哪怕是给内部看的中间稿。柔度曲线漂亮、结构图好看，这两件事和结果可信之间没有必然联系，能替你兜底的只有那几个数字。希望帮到你。

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