★ FEATURED ARTICLE 简介拓扑优化基础实现合集面向刚开始接触结构优化与 MATLAB 编程的初学者。既有 99 行、88 行等精简版本演示二维拓扑优化主流程也纳入优化准则法OC与移动渐近线法MMA两种典型算法并引入 3D 扩展程序展示从二维到三维的材料分布设计覆盖网格划分、荷载与约束设定、材料属性分配等基本环节。压缩包共含 9 个文件全部为 .m 源码整体仅 22KB便于逐行阅读和快速修改MMA 相关子程序与主程序分离适合拆解算法迭代与收敛细节。资源已被 3195 人学习下载入门阶段参考价值较高。通过运行和对比 OC 与 MMA 的优化结果读者能直观理解不同策略在刚度最大化问题上的收敛差异3D 版本又为进一步研究三维拓扑优化提供了可复用的代码底座适合用于课程作业、论文复现或算法改造。1. 拓扑优化程序到底在写什么99行、88行和MMA/OC算法的关系拓扑优化听起来像是力学论文里才有的东西但你手头这套99行、88行程序却能把你从公式堆里捞出来。它用几十行MATLAB代码实现了完整的结构拓扑优化流程给一个设计域、边界条件和体积约束程序自己决定哪里保留材料、哪里掏空。核心就两个算法——OC优化准则法和MMA移动渐近线法。OC在经典99行代码里作为主力更新策略MMA则在88行版本或扩展版中提供更稳定的收敛尤其适合带约束、多工况的问题。适合谁刚接触拓扑优化的研究生、想快速复现优化效果的结构工程师还有想拿现成代码做二次开发的硬件工程师。读懂它你就掌握了从“任意初值”到“清晰材料分布”的完整链路。2. 技术栈与程序结构从材料插值到灵敏度过滤2.1 SIMP材料插值为什么用惩罚因子让中间密度消失在讲代码之前先把技术栈锁死。这套拓扑优化程序建立在“变密度法SIMP材料插值”上也就是把每个有限单元的材料密度x当成设计变量x1表示有材料x0表示空洞。如果直接让x取0到1的连续值优化完后会得到一大片灰色区域比如0.4、0.6的中间密度那没法加工。SIMP的思路是把单元弹性模量写成密度的幂函数E( x ) x^penal * E0其中penal是惩罚因子一般取3E0是实体材料的杨氏模量。在99行代码中单元刚度矩阵在组装时都要乘以x^penal比如下面的核心语句E E0 * x(e)^penal; ke E / (1 - nu^2) * [ ... ]; % 四节点平面应力单元刚度矩阵逻辑说明当x0.5时若penal3E只有实体材料的12.5%这样结构会把半密度区域视为“很软”的浪费材料优化器自然倾向于把x推向0或1。参数说明penal不能太小太小了灰色区域多也不能太大太大会导致收敛硬、容易出现局部最优。教学代码里penal3是经验值二开时可以在2.5~4之间试。这里还涉及材料参数弹性模量E0、泊松比nu四节点平面应力单元默认厚度为1。你把E0随意改一个量级对拓扑形态几乎没影响但会影响位移数值所以如果后面要做应力约束需要保持单位制一致。2.2 OC与MMA的迭代逻辑差异OC和MMA是两种层级的优化器。OC的全称是优化准则法它的出发点是拉格朗日乘子下的KKT条件。对于“最小化柔度单个体积约束”这种经典组合可以推导出一个非常紧凑的启发式更新公式再用二分法找个乘子满足体积约束。OC的优点是代码短、每轮迭代计算量小缺点是一旦你加入应力约束、位移约束或多工况它的推导公式就不再成立基本得重新推导。MMA则走的是另一条路移动渐近线法。它不依赖具体物理问题的闭式解而是在每个迭代点用“移动渐近线”构造一组凸近似子问题然后求解该子问题得到下一步设计变量。它天然支持多约束因为每个约束函数都可以各自做近似。缺点是子程序复杂需要维护很多内部状态比如上渐近线、下渐近线、外推因子等。一顿操作下来代码量从10行变成几百行但稳定性和泛化性都更好。维度OCMMA适用场景单约束、单工况多约束、多工况、含应力约束代码量约10行数百行每轮计算量小中等对非线性约束适应差好收敛稳定性快但可能振荡慢但稳参数对比看下来99行版本选OC88行版本默认也是OC但程序结构容易换成MMA。后面第4章会单独说怎么切。2.3 目录和文件清单一组.m文件如何组织这套程序拿到手目录一般很干净。顶层脚本以top99.m或top88.m为主里面包含了设计域网格划分、边界条件、载荷定义、有限元循环、灵敏度过滤和优化器更新。旁边大概率放着一个mma.m子函数这是MMA算法的核心黑匣子里面是几百行近似子问题的求解逻辑。有些版本会额外带KE.m或compound.m等辅助文件用来生成单元刚度或画出迭代曲线。. ├── top99.m # OC主程序经典99行 ├── top88.m # OC主程序改进过滤和边界处理 └── mma.m # MMA子程序供88行版本扩展调用运行入口就是主脚本。你不需要额外安装任何工具箱只要MATLAB基础环境就能跑。执行完脚本工作区里会留下x密度矩阵、c柔度历史曲线和U位移向量。你可以在命令行敲plot(1:loop, c_history); xlabel(迭代步); ylabel(柔度);就能看到收敛曲线。整个项目就是这样一个“有限元求解灵敏度过滤优化器更新”的三层结构理解这个三层结构后面改代码不会迷路。2.4 99行和88行主程序差异与扩展方向99行版本的主要目的是教学代码刻意压缩在一屏之内把SIMP、灵敏度过滤和OC更新都堆在主循环里。88行版本是它的升级重点优化了边界条件的索引方式和过滤逻辑。具体说99行版本把节点自由度编号写成了硬编码的数组88行版本则用edofMat和freeDofs单独计算改载荷工况时不容易错。另一个区别是88行版本把灵敏度过滤封装成了一个函数主循环结构更清晰适合继续加模块。很多人拿着一套代码就开始改结果发现99行版本改不动因为所有变量都耦合在一起。我建议二开从88行版本入手。它对SIMP惩罚项、过滤、体积约束都留了很清晰的入口你只需要替换某个子函数就能切换优化算法。这个细节在后面几章会反复用到。3. 核心模块拆解99行代码里最关键的三段逻辑3.1 主循环有限元求解与目标函数计算99行代码的主循环只有几段但每一段都有讲究。首先要初始化网格参数nelx和nely是设计域横向和纵向的单元数volfrac是总体体积占比penal是SIMP惩罚因子rmin是过滤半径。这些是全局变量几乎每个函数都要用到。随后进入迭代循环每次循环的第一步是组装全局有限元方程。这里有个效率技巧先用稀疏矩阵的索引向量iK和jK把所有单元贡献收集到一个一维数组里再用sparse一次性构造出全局刚度阵避免在循环里反复扩写矩阵。sK sparse(iK, jK, sKCoef, ndof, ndof); sK(freedofs, freedofs) sK(freedofs, freedofs); U zeros(ndof,1); U(freedofs) sK(freedofs, freedofs) \ F(freedofs); c sum(F .* U);逻辑说明先通过预分配的索引把单元刚度组装成大K然后只保留自由度为freedofs的区域做一次线性求解。这里F是外力向量和节点自由度一一对应U是节点位移向量。目标函数柔度c等于外力点乘位移也就是结构应变能的两倍最小化它就是最大化刚度。参数说明ndof是总自由度数等于2*(nelx1)*(nely1)freedofs是剔除了固定边界的自由度编号集合。如果这里边界索引算错得到的位移是错的后续所有灵敏度都跟着错。在99行代码里边界条件通常这样定义fixeddofs [2*floor(nelx/2)*(nely1)1 : 2*(nely1)*(nelx1)]; freedofs setdiff(1:2*(nelx1)*(nely1), fixeddofs); F zeros(2*(nelx1)*(nely1),1); F(2*(nely1)*floor(nelx/2)1,1) -1;这段的作用是把MBB梁的右下角约束住再在左下角施加竖向载荷。floor(nelx/2)是为了定位中间节点的索引如果你把nelx改成奇数这个索引会偏进而导致载荷位置不在设计域中央。这也是为什么我建议nelx用偶数的原因之一。3.2 灵敏度计算与过滤解决棋盘格的玄学柔度对设计变量的灵敏度是优化的核心。对单工况线性结构柔度导数可以写成dc -penal * x(e)^(penal-1) * U_e * KE * U_e;这里U_e是当前单元的节点位移向量KE是单元刚度矩阵。负号表示增加材料会降低柔度也就是让结构变硬。理论上拿这个灵敏度直接做梯度下降也行但拓扑优化有个经典的玄学现象如果你不过滤灵敏度优化结果会布满棋盘格。这不是数值误差而是问题本身存在棋盘模式的最优解实际没法加工。过滤的方式是把每个单元的灵敏度替换成以它为中心、半径rmin范围内所有单元灵敏度的加权平均。权重为rmin - |距离|离得越近影响越大。下面是99行和88行版本里常见的过滤实现function dcn filter_sensitivities(nelx,nely,rmin,x,dc) dcn zeros(nely,nelx); for i 1:nelx for j 1:nely sumw 0.0; dcn_sum 0.0; for k max(i-rmin,1):min(irmin,nelx) for l max(j-rmin,1):min(jrmin,nely) w rmin - sqrt((i-k)^2 (j-l)^2); sumw sumw w; dcn_sum dcn_sum w * dc(l,k); end end dcn(j,i) dcn_sum / sumw; end end end逻辑说明对每个单元先找出过滤半径覆盖的邻居单元按距离计算权重然后加权求和并归一化。dcn是过滤后的灵敏度优化器要用的是dcn而不是原dc。参数说明rmin的单位是单元尺寸的倍数通常取1.5到3如果你把rmin设成0过滤就是白做结果会像电视雪花一样稀碎。除了灵敏度过滤有些扩展版本还会做密度过滤即对设计变量本身做类似卷积。密度过滤比灵敏度过滤更物理因为最终结果会得到更清晰的尺寸控制。但88行版本默认只做灵敏度过滤因为它代码量小且对MBB这类经典算例已经够用。3.3 OC更新准则拉格朗日乘子二分法OC更新是99行代码的压轴戏。它的推导过程可以简化为把“目标函数体积约束”写成拉格朗日函数求每个设计变量的偏导并令其为零然后变形得到一个启发式更新公式。最终的形式是x_new max(0.001, max(x - move, min(1, min(x move, x * sqrt(-dc / lmid)))))这里lmid是拉格朗日乘子move是设计变量每步允许改变的最大幅度。用sqrt而不是线性项是因为柔度灵敏度与密度之间存在幂律关系这个启发式在SIMP框架下很好用。0.001和1是密度下限和上限避免刚度阵奇异。为了满足体积约束需要对lmid做二分搜索。因为体积是lmid的单调函数lmid增大更新后的密度整体减小体积也随之减小。所以可以放心二分function xnew OC(nelx, nely, x, volfrac, dc) l1 0; l2 100000; move 0.2; while (l2 - l1 1e-4) lmid 0.5 * (l2 l1); xnew max(0.001, max(x - move, min(1, min(x move, x .* sqrt(-dc ./ lmid))))); if sum(sum(xnew)) - volfrac * nelx * nely 0 l1 lmid; else l2 lmid; end end end逻辑说明循环里先按当前lmid生成新密度然后判断这个新密度的总体积是否超过目标体积。超过说明lmid太小把下界抬高否则把上界压低。迭代几百次后区间长度小于1e-4就停此时体积约束误差在可接受范围内。参数说明move0.2是经典值限制了每步最多变化0.2防止迭代震荡。l2100000看起来很大但对柔度灵敏度来说足够。当你的单元数和边界条件变化后如果二分找不到最优解先从l2入手把它调到1e6找不到再说。4. MMA算法接入从OC到移动渐近线的改造4.1 MMA子问题的近似模型MMA之所以叫“移动渐近线”是因为它在每一轮迭代里都要构造一组渐近线把非线性原问题近似成一系列凸子问题。这个概念不复杂你有一根能移动的“渐近线”靠近设计变量时近似更激进离得远时更保守。子问题的目标函数通常是设计变量的线性项倒数列项的和约束函数也类似然后解这个显式子问题比直接解原问题容易得多。代码里mma.m一般会接收以下输入xval 当前设计变量 xold1, xold2 前两轮设计变量 fval, fvalold 当前和前一轮目标函数值 gval, gvalold 当前和前一轮约束值 dfdx, gfdx 目标函数灵敏度和约束灵敏度输出一个更新后的设计变量xnew和更新后的约束值。这些历史输入是MMA的灵魂因为它要用xold来估计曲率用fvalold来估计变化趋势所以如果你每次调用都重新传空数组MMA就退化成一种很慢的梯度法。内部结构大致是这样function [xnew, fnew, gnew] mma_update(x, fval, gval, dfdx, dgdx, pars) % 1. 根据 x 和 xold 更新上下渐近线 % 2. 构造显式子问题 % 3. 求解子问题得到 xnew end逻辑说明设计变量在上下渐近线之间被限制渐近线位置随迭代移动因此叫移动渐近线。参数说明这个子函数虽然叫mma.m但它不是MATLAB内置函数必须从资源包里带过来否则运行时会报“未定义函数或变量 mma”。4.2 如何把OC换成MMA而不破坏主流程在88行版本里把OC替换成MMA其实不用大改主循环。你只要在迭代步里维护好设计变量、目标函数和约束数组。比如设两个约束体积约束和位移上限约束g1 sum(x(:)) / (nelx*nely) - volfrac; g2 max(U(freedofs)) / disp_limit - 1; g [g1; g2]; [x_new, fnew, gnew] mma_update(x, fval, g, dfdx, dgdx, pars);调用示意。注意MMA要求你提供所有的约束值并且对每个约束都要给一组灵敏度dgdx所以你要对限制位移的节点位移做一次伴随求解。如果不做伴随只给体积约束的灵敏度那么g2就没法被正确优化程序会报尺寸不匹配。我见过很多人在这里翻车以为把dc当成dfdx传进去就算接好MMA结果忘了dgdx第二维长度要与约束数量一致。mma.m内部会严格检查维度报错通常很明确但问题在于你很难定位是哪一列写错。4.3 参数设置与收敛性对比MMA的内部参数按经验来设置。初始渐近线偏移asymp_init可以取0.1到0.5偏移越小近似模型离当前点越近收敛越慢但更稳外推因子albefa在1.05到1.2之间控制渐近线移动幅度太大容易振荡移动限制move取0.1到0.3。如果你发现MMA目标函数下降太快但振铃严重就把albefa调小、把move调小。参数含义经验范围备注asymp_init初始渐近线偏移0.1 ~ 0.5越小越稳albefa渐近线外推因子1.05 ~ 1.2太大易振荡move设计变量移动限制0.1 ~ 0.3越小越保守asymp_dec第2轮渐近线缩放0.7 ~ 0.9控制历史权重对比OC和MMA跑同一个MBB梁算例OC大概30~40步收敛到稳定值MMA大概40~80步。OC每步只调用一个二分MMA每步要解一个子问题所以总时间OC更快。但当你增加了一个位移约束后OC的更新公式需要重新推导而MMA只需要在约束数组里多塞一行。这就是为什么你看到的多数扩展拓扑优化程序都选择MMA。5. 运行与避坑从MATLAB报错到结果异常排查5.1 运行前准备版本、工具箱、主程序入口运行这套程序不需要额外工具箱我建议MATLAB R2020以后。原因倒不是语法不兼容而是早期版本对sparse和绘图函数的显示有些差异容易让你误判优化效果。把主脚本和mma.m放在同一目录在命令行敲top99就能跑。默认算例是MBB梁设计域60x20个单元体积比例0.5惩罚因子3过滤半径1.5整个迭代大概几十秒。如果你看到命令行疯狂刷数字那是循环里打印了目标函数和体积比不用慌。提示如果提示缺少mma函数先执行which mma确保它在你当前路径上。不要同时打开top99和top88它们共享一些变量名混着跑容易覆盖工作区。5.2 常见问题一迭代后拓扑左右不对称现象MBB梁算出来左右不对称密度云图明显歪了。原因设计域列数为奇数中间单元不能成对出现或者载荷加载在单一节点上产生不平衡的力流。解决把nelx改成偶数比如60、80如果结构要求严格对称可以只算半个设计域另一边镜像复制。注意镜像时要同时考虑载荷和约束的镜像。5.3 常见问题二密度图全是棋盘格现象优化结果出现密密麻麻的黑白相间小格子像棋盘。原因灵敏度过滤没有生效rmin设成0或过滤半径远小于单元尺寸。解决把rmin调到至少1.5倍单元边长并确认过滤循环里的坐标距离单位正确。如果你把rmin改成10结果会变成大块实心连细节都没有这也说明过滤半径不宜过大。5.4 常见问题三MMA迭代不收敛目标函数来回跳现象换成MMA后柔度曲线不下降甚至反弹。原因历史迭代参数没存对或者渐近线初始值太靠近设计变量导致近似模型过于激进。解决放宽MMA的初始渐近线偏移把移动限制调小到0.1附近并检查每次调用时输入的设计变量必须是上次返回的值不能重新初始化。还有个血泪经验MMA子程序内部状态是全局变量如果你在不同算例之间没清除它们上一次的收敛信息会污染下一次。可以在每次运行前执行clear mma或重启脚本。5.5 常见问题四体积约束一直不满足最终体积偏高或偏低现象迭代结束时体积与设定体积差了几个百分点。原因OC二分法的上界l2太小或者volfrac乘以nelx*nely时用了浮点误差累积。解决把l2从100000改成1e6并在体积检查时加一个小容差比如1e-3。如果用的是MMA检查体积约束的灵敏度是否被归一化g和dgdx单位不统一也会让体积无法收敛。5.6 常见问题五优化结果中出现了孤立的细小杆件现象结果中有细长杆或漂浮材料块像是断掉的桁架。原因过滤半径太小或设计变量下限0.001被截断过多次。解决把rmin稍微调大或者对最终密度做后处理小于阈值的单元强制置0。对于实际加工还要检查最小杆件厚度这个在第6章验证部分会提到。6. 二次开发加一个约束条件并验证拓扑结果6.1 在88行框架里加入体积约束的缩放如果你只是想调调体积占比把主脚本顶部的volfrac0.5改成0.3就行。但改成0.3后结构会更纤细可能产生细小杆件改成0.7后优化会倾向于保留大量材料拓扑纹理不明显。我通常会把惩罚因子从3提高到4避免中间密度被摸鱼。这一步几乎零风险跑一次就能看出效果。6.2 改载荷工况从单工况到多工况加权实际项目中结构可能要承受多个方向载荷。你可以把多个载荷向量F1、F2分别求解位移然后目标函数改为加权柔度c alpha*c1 (1-alpha)*c2灵敏度也要加权。这样能生成类似“同时抗竖直压力和横向剪切”的拓扑构型。c alpha * (F1*U1) (1-alpha) * (F2*U2); dc alpha * dc1 (1-alpha) * dc2;逻辑说明alpha是第一个工况的权重取值范围0到1。每个工况的位移必须独立求解灵敏度也必须用各自工况的位移来计算。如果你把两个载荷向量直接相加再求解那只是单工况得不到同时承受多方向载荷的拓扑。参数说明如果两个载荷大小量级不同建议先归一化载荷再调alpha否则权重会被量纲吃掉。6.3 验证方法目标函数曲线与对称性检查改完之后不要只看最后一张密度图。我会把每一步c打印出来画成曲线确认在稳定下降。然后检查材料分布是否满足最小尺寸约束比如最小杆件厚度至少有3个单元否则实际加工无法制造。另一个常用方法是用对称几何验证把设计域设为对称跑完看结果是否对称。从那以后我每次改完这套拓扑优化代码都会强制走一遍“体积约束→灵敏度→收敛曲线→对称性”四步检查才敢把结果拿去做详细设计。这四步看着简单但能滤掉大部分参数设置翻车。这套代码的价值不在“99行”这个数字而在于你把它当成试验台用OC理解优化迭代再用MMA拓宽问题类型。希望帮到你。本文还有配套的精品资源点击获取
阅读完成 · 觉得有帮助？
