搞过结构优化的人应该都有体会想快速验证一个思路、想看懂一篇拓扑优化论文的算法细节、想在课程作业里拿出一版能跑的结果最容易卡住的往往不是理论而是“理论到程序之间那层窗户纸”。今天这篇要拆的题目很典型——四边形元最小化应变能的二维拓扑优化关键词就是四边形元、二维拓扑优化、应变能、Matlab源码。它本质上是把“材料怎么分布最合理”这个直觉问题翻译成一道可以交给计算机迭代求解的数学问题在体积约束下用四节点矩形单元离散设计域最小化结构的总应变能。这套流程在论文里叫SIMP拓扑优化在工程里叫“给定材料、刚度最大化”的经典布局设计在Matlab里则是从单元刚度矩阵到灵敏度更新的一整套可运行代码。我见过太多人一上来就研究优化器细节结果卡在刚度矩阵组装。这篇我想按自己的理解把题目拆成“物理含义、算法原理、源码实现、调试排错”四个层次最后再聊几个我踩过坑以后才明白的扩展方向。无论你是力学背景的本科生、做轻量化设计的工程师还是刚接触拓扑优化的研究生按照这个顺序往下读应该能把程序跑通并且知道每个变量、每次迭代到底在干什么。1. 项目拆解与核心思路1.1 “最小化应变能”到底在优化什么先说人话版。应变能就是结构在外力作用下储存的能量数值上等于外力做功。对一个线弹性结构应变能越小说明结构受力后变形越小、整体刚度越大。所以“最小化应变能”和“最大化刚度”是同一件事的两个说法写出来就是$$\min_{\rho} ; c(\rho)\mathbf{U}^{\top}\mathbf{K}(\rho)\mathbf{U}\sum_{e1}^{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这是实践出来的经验值$p1$ 时问题退化成线性材料插值结果充满灰度中间密度$p$ 太大比如7以上则优化收敛困难容易陷入局部最优。$p3$ 的妙处在于密度0.5的单元刚度贡献只有 $0.5^312.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_em) \ \min(1, \rho_em), \min(1, \rho_em) \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.21.5想更稳健可以取2。这个参数不是物理参数纯粹是数值调节旋钮但这恰恰是拓扑优化“玄学”的一部分。3. Matlab源码复现从参数表到主循环的实现细节3.1 环境准备与参数含义这套程序不需要额外的Matlab工具箱纯基础矩阵运算就能跑。版本方面2016b以上皆可更早的版本也基本没问题只要支持函数句柄和稀疏矩阵运算就行。网络上那些复杂的“OOP架构多算法融合图像处理系统”“2026b下载”之类的关键词和本程序无关别被带偏这里的核心就是一个topopt_main.m主脚本加若干函数文件。拿到源码后先看参数表。无论哪个版本下面这几个参数是标配参数含义典型取值nelx水平方向单元数60120nely垂直方向单元数2040volfrac允许的材料体积分数0.30.5penalSIMP惩罚系数3rmin滤波半径1.5倍的单元边长ft滤波类型标记1或2跑第一个算例时建议先用小网格比如60×20验证流程迭代一两百步也只要几秒钟速度够快你可以随时打印中间结果观察变化。等把逻辑摸透了再根据机器配置加大网格和迭代步数。3.2 网格、自由度映射与边界条件的装配Q4元的自由度映射是所有装配的基础这里值得花几分钟彻底理解。对规则矩形网格用(elx, ely)定位单元其中elx1..nelx是水平编号ely1..nely是垂直编号。每个节点有两个自由度x方向和y方向的平动位移。单元四个角节点的全局节点编号按“从左下角顺时针”的顺序可以这样计算左下角节点n1 (nely1)*(elx-1) ely右下角节点n2 (nely1)*elx ely右上角节点n3 (nely1)*elx ely 1左上角节点n4 (nely1)*(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*(nely1)*nelx 1, 1) -1; % 顶部中点施加向下的单位力 fixeddofs [1:2*(nely1), 2*(nely1)*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梁本身就是对称结构优化域取一半镜像回去就得到完整结果计算量直接减半。这个方法很多论文里都在用是合法的工程技巧不是偷懒。下面把最常见的问题整理成一张速查表方便你在调试时对照现象可能原因快速排查解决建议收敛曲线震荡不降灵敏度符号错误打印首次迭代的dc是否全为负检查负号和1/2系数一致性结果全是棋盘格滤波半径太小 / 滤波未生效查看第5次迭代密度图提高rmin到1.5刚度矩阵奇异Emin0或约束不足检查K行列式或求解警告设Emin1e-9E0检查固定自由度结果灰度太重penal太低 / 迭代不够看密度直方图分布penal3迭代300步以上迭代不收敛移动极限太大观察每步密度最大变化量move0.2必要时降到0.1速度慢循环组装刚度矩阵用profile定位耗时函数向量化sparse组装5. 从教程代码到工程应用的扩展路径5.1 从2D到3D三个主要改动二维程序跑通之后很多人会想往三维扩展。如果只是把网格从矩形扩展成立方体有三个地方必须同步调整。自由度映射变成每个节点3个自由度单元刚度矩阵从8×8变成24×24高斯积分从2×2变成2×2×2多了一重循环载荷和边界的定义复杂度也上升一个量级一个面节点的自由度数很大固定面时不能犯“只固定一个点”的老毛病。三维程序的体积约束灵敏度仍然很简单因为单元体积因子依然可以直接代入。但三维问题的规模膨胀很快100×30×20网格就是60000个单元稀疏矩阵求解成为真正的性能瓶颈这时就不是入门阶段需要考虑的问题了。5.2 加约束、加制造限制、换优化器下一步怎么走SIMPOC只适合解决最简单的问题。一旦你开始加入应力约束、位移约束、多工况载荷、自支撑约束OC的简单形式就不够用了需要换用MMAMethod of Moving Asymptotes或SQP等通用优化器。不过即使换优化器整个有限元分析框架和灵敏度推导思路完全不变OC程序里的那一套“组装→求解→求灵敏度”的流程依然可以作为基础框架复用。增材制造约束是近几年很热门的方向要求设计结果没有悬垂结构或者所有材料都能沿某个方向打印。这本质上是给密度场加了一个方向性的几何约束SIMP框架仍然适用只是目标函数和约束集合更复杂。还有多材料拓扑优化密度从标量变成“每种材料的体积分数向量”目标函数不变约束变成了多个体积分数约束的和OC的拉格朗日乘子也要变成多变量的版本。5.3 一点点个人经验程序跑通只是起点真正有用的是理解每一步背后的物理。我实际使用中最大的感受是不要一上来就追求大网格和高精度先用小网格把物理趋势看清楚再慢慢加细。拓扑优化对网格密度非常敏感同一个体积分数在60×20网格下会得到清晰的单层结构在200×60下可能就分层了这不是程序bug而是“优化结果包含更细的特征”的表现。另外滤波半径和网格尺寸是一对相关参数网格加细了滤波半径如果不跟着按比例调整结果特征尺寸会变因此“网格无关性”并不是天然成立的需要刻意控制滤波半径所对应的物理尺寸。最后分享一个小技巧拿到任何拓扑优化源码第一件事不是看公式而是把网格缩到特别小比如30×10跑一遍把每步迭代的目标函数值打印出来。如果这个小网格能跑出经典算例的基本形态那大网格大概率没问题。反过来大网格直接跑出了问题你连在哪个环节都定位不清楚。这个习惯救了我很多次希望也能帮到你。