企业官网建设流程全解析

搞过结构优化的人应该都有体会：想快速验证一个思路、想看懂一篇拓扑优化论文的算法细节、想在课程作业里拿出一版能跑的结果，最容易卡住的往往不是理论，而是“理论到程序之间那层窗户纸”。今天这篇要拆的题目很典型——四边形元最小化应变能的二维拓扑优化，关键词就是四边形元、二维拓扑优化、应变能、Matlab源码。它本质上是把“材料怎么分布最合理”这个直觉问题，翻译成一道可以交给计算机迭代求解的数学问题：在体积约束下，用四节点矩形单元离散设计域，最小化结构的总应变能。这套流程在论文里叫SIMP拓扑优化，在工程里叫“给定材料、刚度最大化”的经典布局设计，在Matlab里则是从单元刚度矩阵到灵敏度更新的一整套可运行代码。

我见过太多人一上来就研究优化器细节，结果卡在刚度矩阵组装。这篇我想按自己的理解，把题目拆成“物理含义、算法原理、源码实现、调试排错”四个层次，最后再聊几个我踩过坑以后才明白的扩展方向。无论你是力学背景的本科生、做轻量化设计的工程师，还是刚接触拓扑优化的研究生，按照这个顺序往下读，应该能把程序跑通，并且知道每个变量、每次迭代到底在干什么。

1. 项目拆解与核心思路

1.1 “最小化应变能”到底在优化什么

先说人话版。应变能就是结构在外力作用下储存的能量，数值上等于外力做功。对一个线弹性结构，应变能越小，说明结构受力后变形越小、整体刚度越大。所以“最小化应变能”和“最大化刚度”是同一件事的两个说法，写出来就是：

$$\min_{\rho} ; c(\rho)=\mathbf{U}^{\top}\mathbf{K}(\rho)\mathbf{U}=\sum_{e=1}^{N} \rho_e^{p},\mathbf{u}_e^{\top}\mathbf{k}_0,\mathbf{u}_e$$

约束条件是平衡方程 $\mathbf{KU}=\mathbf{F}$、体积分数 $\sum \rho_e v_e \le fV_0$，以及每个单元密度的上下限 $0<\rho_{\min}\le\rho_e\le1$。

这里的 $\rho_e$ 是单元密度，也就是每个格子究竟放不放材料。理想情况下 $\rho_e$ 只能是0或1，但直接做离散0/1优化是NP难问题，没法解大网格，所以SIMP方法放了一个口子：允许中间密度存在，但用惩罚系数 $p$ 让中间密度在刚度贡献上“不划算”，最终迭代结果自然向0/1收敛。这个思路很聪明，相当于把“选哪个格子去掉”这种组合爆炸问题，变成了一个“每个格子密度填多少”的连续优化问题。而应变能这个目标函数是“自伴随”的，灵敏度计算不需要额外求解伴随方程，这让整个算法在Matlab里实现起来出奇地紧凑，也是为什么现在90%的入门拓扑优化代码都拿它开刀。

注意：很多教材里应变能写成 $c=\frac{1}{2}\mathbf{U}^{\top}\mathbf{KU}$，因为线弹性结构的应变能确实有个1/2系数。但简化程序时常把1/2省略，反正它只是个常数，不影响密度分布的最优结果。但一旦你确定不用1/2，灵敏度的公式里也不要加1/2，前后要自洽。

1.2 四边形元选型：比三角形强在哪儿，比八节点简单在哪儿

这个题目点名要用四边形元，非常合理。二维拓扑优化的网格无非三种选择：三节点三角形（常应变三角形，CST）、四节点四边形（Q4）、八节点四边形（Q8）。

三节点三角形最容易被方案“诱惑”上，因为三角形网格剖分灵活，能适应复杂几何边界。但它的最大问题是偏刚，一个单元内应变恒定，精度差，棋盘格现象非常严重，在拓扑优化里常常出现次优解。如果网格不均匀，结果还会受网格剖分方式影响——同一个设计域，换一种剖分方式能算出不一样的结构，这是很让人崩溃的。

Q4元则不同，它是双线性单元，应变在单元内线性变化，精度比CST高一个档次，而实现复杂度又没有Q8那么高。Q8带边中节点，更不容易剪切锁定，但每个单元的节点数翻倍，自由度编号、边界条件处理、后处理都更啰嗦。对一个二维拓扑优化的SIMP流程来说，Q4的组合性价比最高：既有足够的精度反映受力路径，代码量又能控制在核心程序几十行的规模。很多经典的99行拓扑优化程序用的就是Q4单元，这不是巧合，是无数人实践出来的选择。

选用四边形元还有一个容易被忽略的好处：规则矩形网格下，单元形状整齐，雅可比行列式均为常数，高斯积分计算稳定，滤波半径的处理也方便。网格拓扑关系用一个二维索引(ely, elx)就能描述清楚，Matlab里用矩阵而不是稀疏索引来组织单元会顺畅得多。

1.3 先定边界再定优化：标准算例的选择

拓扑优化有个特点：同样一套算法，边界条件和载荷位置不同，结果完全两样。所以复现代码时，不要一上来就设计自己想象中的工况，先用几个标准测试算例验证程序对不对，再改边界。

最常见的三个算例：

• MBB梁：设计域长宽比2:1，底部左下角约束竖直位移，右下角约束水平和竖直位移，顶部中点受竖直向下的集中力。得到的结果是经典的两根斜撑桥式结构，几乎每一篇拓扑优化论文都会拿它做对比图。
• 左端固支悬臂梁：设计域固定左边界，右下角节点受竖直向下的力。结果通常是类似树枝分叉的悬臂结构。
• 短悬臂梁：设计域右侧两个角点局部固定，左侧中部施加载荷，结果会看到类似螺栓连接抗拉板的结构。

标准算例的最大价值，是结果有公开对比图。如果你的程序跑出来和论文里那张经典图差不多，那基本说明单元刚度矩阵、灵敏度、OC更新都没错。我第一次复现时就是因为直接换了个怪异载荷，结果看起来怎么都“不对劲”，后来才发现是灵敏度符号搞反了。所以我的建议是：源码到手先跑MBB梁，跑对了再改工况。

2. SIMP插值、灵敏度分析与OC更新的完整推导

2.1 SIMP材料模型：让“中间密度”无处遁形

SIMP全称Solid Isotropic Material with Penalization，直译是“带惩罚的各向同性实体材料模型”。它的核心是给每个单元赋予一个“虚拟杨氏模量”：

$$E(\rho_e)=E_{\min}+\rho_e^{p},(E_0-E_{\min})$$

$E_0$ 是实体材料模量，$E_{\min}$ 是一个很小的数（通常取 $10^{-9}E_0$），目的是避免空单元导致刚度矩阵奇异。$p$ 是惩罚系数，一般取3，这是实践出来的经验值：$p=1$ 时问题退化成线性材料插值，结果充满灰度中间密度；$p$ 太大（比如7以上）则优化收敛困难，容易陷入局部最优。$p=3$ 的妙处在于，密度0.5的单元刚度贡献只有 $0.5^3=12.5%$，也就是说放着50%的材料，只换来12.5%的“战斗力”，优化器很快就会觉得这是浪费，把密度推向两侧。

为什么泰勒展开？因为目标是让连续解尽可能接近0/1离散解。可以这样理解：你给每个单元开一扇门，门开多大是连续的，但SIMP让“半开门”的性价比极差，于是最终门要么几乎全开、要么几乎全关，这就模拟了一个“去材料”的过程。

2.2 灵敏度推导：如何把“改哪个单元”算出来

灵敏度就是“目标函数对每个单元密度的导数”。没有它，优化器不知道往哪个方向改材料分布。推导的核心是链式法则。令 $c(\rho)=\mathbf{U}^{\top}\mathbf{KU}$，对 $\rho_e$ 求导：

$$\frac{\partial c}{\partial \rho_e} = \frac{\partial \mathbf{U}^{\top}}{\partial \rho_e}\mathbf{KU} + \mathbf{U}^{\top}\frac{\partial \mathbf{K}}{\partial \rho_e}\mathbf{U} + \mathbf{U}^{\top}\mathbf{K}\frac{\partial \mathbf{U}}{\partial \rho_e}$$

由于 $\mathbf{KU}=\mathbf{F}$，且载荷 $\mathbf{F}$ 与密度无关，对等式两边求导得 $\frac{\partial \mathbf{K}}{\partial \rho_e}\mathbf{U} + \mathbf{K}\frac{\partial \mathbf{U}}{\partial \rho_e}=0$，于是前两项和第三项正好抵消（这就是自伴随特性的由来），剩下：

$$\frac{\partial c}{\partial \rho_e} = \mathbf{U}^{\top}\frac{\partial \mathbf{K}}{\partial \rho_e}\mathbf{U} = -p\rho_e^{p-1}\mathbf{u}_e^{\top}\mathbf{k}_0\mathbf{u}_e$$

注意，这里的“共振折衷”很重要：我们完全不需要额外求解一组伴随方程，只用在收敛后按单元提取位移向量，再乘上单元刚度矩阵，就能算出灵敏度。好多初学者第一次推到这里都会愣一下：为什么这么简单？对，拓扑优化的入门代码之所以短，很大程度就是因为选了应变能这个自伴随目标。如果你换成应力约束或频率约束，就不能这么偷懒了，这也是为什么后续扩展时问题会复杂得多。

体积约束的灵敏度更简单：$\partial V/\partial \rho_e = v_e$，因为每个单元的体积只和自身密度线性相关。

我想特别提醒一个细节：如果目标函数带1/2系数，那么灵敏度同样要带1/2。我在调试时见过有人程序里目标函数写c = U'*K*U，灵敏度却套用了教材里带1/2的公式，结果灵敏度整体放大两倍，OC更新里的拉格朗日乘子也得跟着瞎调，收敛曲线“上蹿下跳”，就是不知道问题在哪。

2.3 OC准则与拉格朗日乘子的二分搜索

有了灵敏度之后，怎么把密度从旧值更新到新值？入门程序最常用 Optimality Criteria（优化准则法，简称OC）。OC不是万能的，但在单约束、最小柔度问题里效果极好，更新非常稳定。

先构造拉格朗日函数 $\mathcal{L}=c + \lambda\left(\sum\rho_e v_e - fV_0\right)$，对 $\rho_e$ 求导并令其为0，得到：

$$B_e = \frac{-\partial c/\partial \rho_e}{\lambda,\partial V/\partial \rho_e} = \frac{-\partial c/\partial \rho_e}{\lambda v_e}$$

这个 $B_e$ 可以理解为“材料放在这个单元上的边际收益与边际体积成本的比值”。OC的更新策略是让密度向 $B_e^{\eta}$ 方向移动，$\eta$ 通常取0.5，阻尼系数 $m$ 通常取0.2，即每次迭代每个单元密度最大只能变化0.2。写成三段式：

$$\rho_e^{\text{new}} = \begin{cases} \max(\rho_{\min}, \rho_e - m), & \rho_e B_e^{\eta} \le \max(\rho_{\min}, \rho_e-m) \ \rho_e B_e^{\eta}, & \max(\rho_{\min}, \rho_e-m) < \rho_e B_e^{\eta} < \min(1, \rho_e+m) \ \min(1, \rho_e+m), & \min(1, \rho_e+m) \le \rho_e B_e^{\eta} \end{cases}$$

这里 $\lambda$ 是拉格朗日乘子，它不是一个由公式直接算出来的量，而是需要靠二分法搜索：给定一个 $\lambda$，用上面的公式更新所有密度，统计总体积，如果超过目标体积 $fV_0$，说明 $\lambda$ 太小，放大；如果低于目标体积，说明 $\lambda$ 太大，缩小。如此往复几十次，就能找到一个让体积约束刚好满足的 $\lambda$。

实操心得：二分搜索的迭代次数设50次足够，多了没意义，少了体积约束会飘。上界设 $10^{9}$ 是一个安全值，但如果你看到程序里 $\lambda$ 一直顶着上界不动，第一反应不应该是加大上界，而应该回头检查目标函数正负号——我就在这个问题上浪费过一整天。

2.4 棋盘格问题与密度滤波

所有初学拓扑优化的人都会遇到同一个画面：迭代结果像国际象棋棋盘，高密度和低密度单元交错出现。这其实是数学上合理的解，但工程上不可制造，而且应变能会虚低，结构看起来很“密集”，本质上却是一些细条相互铰接的脆弱骨架。

Q4单元能缓解棋盘格，但不能根除。最有效的工程手段是滤波（filter）。核心思想很简单：一个单元的“真实密度”（或灵敏度）不应该只看它自己，而要看它周围一个半径 $r_{\min}$ 内所有单元的加权平均。比如密度滤波的公式：

$$\tilde{\rho}e = \frac{\sum{j \in N_e} w_{ej} v_j \rho_j}{\sum_{j \in N_e} w_{ej} v_j}, \quad w_{ej}=\max(0, r_{\min}-\text{dist}(e,j))$$

距离越近权重越大，距离超过 $r_{\min}$ 就完全忽略。这样就强行抹去了那些尺度小于滤波半径的棋盘格特征。灵敏度滤波的做法类似，只是把滤波作用在灵敏度 $\partial c/\partial \rho_e$ 上，而不是密度场上。

$r_{\min}$ 的取值直接影响结果特征尺寸。太小，棋盘格压不住；太大，结构变成一大坨，细节全丢。参考经验是：规则网格下，$r_{\min}$ 取单元边长的1.5倍左右比较合适，想得到更清晰的结构可以取1.2～1.5，想更稳健可以取2。这个参数不是物理参数，纯粹是数值调节旋钮，但这恰恰是拓扑优化“玄学”的一部分。

3. Matlab源码复现：从参数表到主循环的实现细节

3.1 环境准备与参数含义

这套程序不需要额外的Matlab工具箱，纯基础矩阵运算就能跑。版本方面2016b以上皆可，更早的版本也基本没问题，只要支持函数句柄和稀疏矩阵运算就行。网络上那些复杂的“OOP架构多算法融合图像处理系统”“2026b下载”之类的关键词和本程序无关，别被带偏，这里的核心就是一个topopt_main.m主脚本加若干函数文件。

拿到源码后，先看参数表。无论哪个版本，下面这几个参数是标配：

跑第一个算例时，建议先用小网格（比如60×20）验证流程，迭代一两百步也只要几秒钟，速度够快，你可以随时打印中间结果观察变化。等把逻辑摸透了，再根据机器配置加大网格和迭代步数。

3.2 网格、自由度映射与边界条件的装配

Q4元的自由度映射是所有装配的基础，这里值得花几分钟彻底理解。对规则矩形网格，用(elx, ely)定位单元，其中elx=1..nelx是水平编号，ely=1..nely是垂直编号。每个节点有两个自由度（x方向和y方向的平动位移）。单元四个角节点的全局节点编号，按“从左下角顺时针”的顺序，可以这样计算：

• 左下角节点：n1 = (nely+1)*(elx-1) + ely
• 右下角节点：n2 = (nely+1)*elx + ely
• 右上角节点：n3 = (nely+1)*elx + ely + 1
• 左上角节点：n4 = (nely+1)*(elx-1) + ely + 1

然后单元自由度向量就是：

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

这个向量要和单元刚度矩阵中“局部自由度排列顺序”严格一致，否则组装的全局刚度矩阵就是乱的。很多改网格时改出bug，问题都出在这里。

边界条件的处理上，经典做法是定义fixeddofs（固定自由度索引数组）和F（力向量）。MBB梁的写法通常是：

F(2*(nely+1)*nelx + 1, 1) = -1; % 顶部中点施加向下的单位力 fixeddofs = [1:2*(nely+1), 2*(nely+1)*nelx + 2]; % 左边界全部固定 + 右下角竖向固定

施加载荷时注意，集中力必须落在节点上。如果想把力加在单元中间，要么细分网格，要么改用分布载荷近似，初学者最容易忽略这一点。

3.3 Q4单元刚度矩阵与全局组装

Q4单元的刚度矩阵按标准等参元流程计算：对每个高斯积分点，计算形函数导数、雅可比矩阵、B矩阵，再累加：

$$\mathbf{k}e = \int{-1}^{1}\int_{-1}^{1} \mathbf{B}^{\top}\mathbf{D}\mathbf{B},t,\det\mathbf{J},d\xi d\eta$$

实际代码可以写得很紧凑。平面应力问题的本构矩阵是：

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

四个高斯积分点取坐标 $\pm 1/\sqrt{3}$，每个权重为1，按2×2循环累加即可。

全局刚度矩阵的组装，关键性能技巧是不要用循环一个个叠加。先建好三个大向量iK、jK、sK，分别记录每个元素贡献的行索引、列索引和数值，全部收集完以后一次性调用sparse(iK, jK, sK)生成稀疏矩阵。这样比三重循环快几十倍，网格规模从40×20变成120×40时，程序依然能在秒级完成一次迭代。

组装完成后，固定自由度的处理用“划去行列”的标准做法。有些程序为了省事把固定自由度对应的行列清零、对角线置1，这样会导致刚度矩阵变成病态；更推荐的做法是只解自由部分：

U(freedofs, 1) = K(freedofs, freedofs) \ F(freedofs, 1);

未约束的自由度位移保持0，之后的U'*K*U计算不受影响。

3.4 OC更新与体积约束的代码级实现

OC更新是整个迭代过程的最核心循环。核心逻辑可以写成这样：

function xnew = OC_update(x, volfrac, dc, move) l1 = 0; l2 = 1e9; % 拉格朗日乘子二分区间 for i = 1:50 mid = 0.5 * (l1 + l2); xnew = x .* sqrt(-dc ./ (mid * ones(size(x)))); xnew = max(0, max(x - move, xnew)); xnew = min(1, min(x + move, xnew)); if sum(xnew(:)) > volfrac * numel(x) l1 = mid; else l2 = mid; end end end

这个函数里x .* sqrt(...)对应的就是 $\rho_e B_e^{0.5}$，其中 $B_e = -dc / (\lambda v_e)$，假设各单元体积相同可以约掉，只保留单元数numel(x)。max和min的嵌套是为了同时满足密度上下限和移动极限move。二分法的出口是体积约束尽可能接近目标值，所以主程序里还需要根据当前总体积反推实际体积分数，作为收敛判断的一部分。

主循环的结构也就自然清晰了：

for iter = 1:maxIter % 1. 按当前密度重新组装全局刚度矩阵K % 2. 求解位移 U % 3. 计算目标函数 c = U'*K*U % 4. 计算灵敏度 dc，并做滤波 % 5. OC更新得到新密度 xnew % 6. 判断收敛：|c_new - c_old| / c_old < tol，若满足则停止 end

这里很多人会犯一个错误：只检查目标函数变化，不检查密度变化。目标函数在小数值上波动很小但密度还在大范围移动的情况确实存在。稳妥做法是同时跟踪change = norm(xnew - x, Inf)，当这个最大值变化小于比如0.001时，再判断收敛。

3.5 后处理与收敛曲线

拓扑优化的后处理其实很简单，但做对了能帮你判断程序是否正常。密度分布图推荐用：

colormap(gray); imagesc(-reshape(x, nely, nelx)); axis equal; axis off; drawnow;

imagesc(-x)会让高密度区域显示为黑色，低密度区域显示为白色，这是拓扑优化社区里约定俗成的可视化方式。收敛曲线则用semilogy或者loglog画目标函数随迭代的变化，正常收敛曲线应该前期急剧下降、后期平缓趋稳，如果看到曲线来回震荡或者缓慢爬升，基本可以判定算法逻辑有问题。

这套15185期源码拿到手后，我建议你先别急着改算例，打开主脚本看清楚主循环的顺序：每一次迭代里，刚度矩阵是用更新前的密度还是更新后的密度？通常应该用上一轮的当前密度。顺序一旦颠倒，优化过程就会变得非常不稳定。

4. 调试实录：拓扑优化最常见的五个坑

4.1 收敛曲线永远不降

这是最常见的翻车现场。程序能跑，密度也在变，但目标函数先是暴跌一段，然后开始震荡甚至回升。问题十有八九出在灵敏度符号上。记住，我们的目标是让应变能变小，所以随着某个单元密度增加、刚度变大，结构的应变能应该下降，也就是 $\partial c/\partial \rho_e$ 应该是个负数。如果你的灵敏度算出来是正数，OC更新就会反向操作，把材料搬到“帮倒忙”的位置。排查方法很简单：在第一次迭代时，手动打印前几个单元的灵敏度数值，检查它们是否全部小于等于0。如果出现明显正数，检查是不是少了个负号，或者目标函数带1/2而灵敏度没带。

4.2 结果全是棋盘格

棋盘格在细网格和小滤波半径下特别容易出现。解决方法依次排查三件事：滤波半径rmin是否太小，建议至少1.2个单元边长；滤波后参与刚度组装的到底是滤波前的原始密度还是滤波后的密度，很多半吊子程序只有灵敏度滤波，密度本身不滤波，导致棋盘格依然存在；第三个是滤波权重公式是否写成了“距离越远权重越大”。这个错误很隐蔽，因为结果仍然有滤波效果，只是结构形态奇奇怪怪，需要你用单步调试追踪w矩阵的数值分布才能发现。

4.3 求解器报矩阵奇异

如果出现 “Matrix is singular” 或者求解出来的位移出现NaN、Inf，第一反应查Emin的设置。Emin设成0是绝对不行的，空单元的刚度矩阵为零，全局矩阵缺秩。正确做法是设一个很小的比例系数，比如 $E_{\min}=10^{-9}E_0$。第二反应查固定自由度：一个结构必须至少约束住刚体位移才能求解。悬臂梁只固定左端一根梁的端点是不行的，网格边界至少固定一整列节点。第三反应：检查载荷是否被放在被固定自由度上，如果有冲突节点同时被施加力和约束，求解器也会给出匪夷所思的结果。

4.4 结果灰度太重，看不出清晰构型

如果你看到最终密度分布整个图都是灰色的云状图，没有鲜明的黑白色块，原因不外乎三个。惩罚系数penal太低，比如用了1，那几乎必然全是灰色；迭代次数太少，密度还没来得及分离成0和1；体积分数定得过高，比如0.8以上，因为大部分区域都要放材料，灰色区域自然多。正常参数下，迭代到后期应该出现明显的黑白分界，少数灰色单元聚集在材料边界过渡带，这才算正常。如果你希望边界更锐利，可以把penal从3提到4，或者加一轮“去灰度”后处理——但注意，更强的惩罚会让优化更容易陷入次优解，3这个值不是随便定的。

4.5 计算太慢怎么办

对二维问题，最能拉低速度的通常是刚度矩阵组装的循环写法。用sparse向量化组装后，另一个瓶颈就是求解大型稀疏线性方程组。如果网格规模到了200×100，别指望普通台式机的\运算能秒级完成。我一般会先把网格调到能看趋势的大小（比如60×20），确认拓扑形态对路，再适当放大。还可以利用对称性只优化一半，比如MBB梁本身就是对称结构，优化域取一半，镜像回去就得到完整结果，计算量直接减半。这个方法很多论文里都在用，是合法的工程技巧，不是偷懒。

下面把最常见的问题整理成一张速查表，方便你在调试时对照：

5. 从教程代码到工程应用的扩展路径

5.1 从2D到3D：三个主要改动

二维程序跑通之后，很多人会想往三维扩展。如果只是把网格从矩形扩展成立方体，有三个地方必须同步调整。自由度映射变成每个节点3个自由度，单元刚度矩阵从8×8变成24×24；高斯积分从2×2变成2×2×2，多了一重循环；载荷和边界的定义复杂度也上升一个量级，一个面节点的自由度数很大，固定面时不能犯“只固定一个点”的老毛病。三维程序的体积约束灵敏度仍然很简单，因为单元体积因子依然可以直接代入。但三维问题的规模膨胀很快：100×30×20网格就是60000个单元，稀疏矩阵求解成为真正的性能瓶颈，这时就不是入门阶段需要考虑的问题了。

5.2 加约束、加制造限制、换优化器：下一步怎么走

SIMP+OC只适合解决最简单的问题。一旦你开始加入应力约束、位移约束、多工况载荷、自支撑约束，OC的简单形式就不够用了，需要换用MMA（Method of Moving Asymptotes）或SQP等通用优化器。不过，即使换优化器，整个有限元分析框架和灵敏度推导思路完全不变，OC程序里的那一套“组装→求解→求灵敏度”的流程依然可以作为基础框架复用。

增材制造约束是近几年很热门的方向：要求设计结果没有悬垂结构，或者所有材料都能沿某个方向打印。这本质上是给密度场加了一个方向性的几何约束，SIMP框架仍然适用，只是目标函数和约束集合更复杂。还有多材料拓扑优化，密度从标量变成“每种材料的体积分数向量”，目标函数不变，约束变成了多个体积分数约束的和，OC的拉格朗日乘子也要变成多变量的版本。

5.3 一点点个人经验

程序跑通只是起点，真正有用的是理解每一步背后的物理。我实际使用中最大的感受是：不要一上来就追求大网格和高精度，先用小网格把物理趋势看清楚，再慢慢加细。拓扑优化对网格密度非常敏感，同一个体积分数，在60×20网格下会得到清晰的单层结构，在200×60下可能就分层了，这不是程序bug，而是“优化结果包含更细的特征”的表现。另外，滤波半径和网格尺寸是一对相关参数：网格加细了，滤波半径如果不跟着按比例调整，结果特征尺寸会变，因此“网格无关性”并不是天然成立的，需要刻意控制滤波半径所对应的物理尺寸。

最后分享一个小技巧：拿到任何拓扑优化源码，第一件事不是看公式，而是把网格缩到特别小（比如30×10）跑一遍，把每步迭代的目标函数值打印出来。如果这个小网格能跑出经典算例的基本形态，那大网格大概率没问题。反过来，大网格直接跑，出了问题你连在哪个环节都定位不清楚。这个习惯救了我很多次，希望也能帮到你。