★ FEATURED ARTICLE 做结构优化的人应该都清楚SIMP框架下的二维拓扑优化看起来只是几十行代码的事但真正跑出一个边界干净、无棋盘格、迭代稳定的结果里面有相当多容易被忽略的细节。今天我把一套基于四边形元、以最小化应变能也就是柔度为目标函数的二维拓扑优化Matlab实现完整梳理一遍这套代码来自我最近整理的源码合集里的15185期网格划分、有限元求解、灵敏度分析、OC更新、敏度过滤全链路都有。如果你正在写自己的拓扑优化程序或者想把一套现成源码改到自己的边界条件里这篇文章应该能帮你少走不少弯路。我默认你至少知道有限元的基本概念读过一点结构优化的入门资料但不需要特别精通因为下面每一步我都会讲清楚“为什么这样做”。标题里那句“四边形元最小化应变能”拆开其实就是三件事用什么单元去离散设计域、以什么物理量作为优化目标、以及怎么迭代更新单元的“有/无”状态。把这三条线理清了这套源码基本就吃透了。1. 从“最小化应变能”说起目标函数为什么选柔度而不是应力1.1 柔度最小化的物理意义结构在载荷下会发生变形内部储存的应变能越少说明结构抵抗变形的能力越强也就是刚度越大。所以“最小化应变能”和“最大化刚度”在数学上是同一个问题。给定材料用量和边界条件我们希望用尽可能少的材料获得尽可能大的刚度这就是拓扑优化最经典的“最小柔度”问题。目标函数写成c(ρ) U^T K U Σ_e E_e(ρ_e) · u_e^T k0 u_e其中U是全局位移向量K是全局刚度矩阵ρ是单元密度向量。从有限元角度看柔度值本质上就是外力在最终位移场上做的功单位是能量量纲所以也有人直接管它叫“应变能最小化”。为什么不是直接以应力为目标这是很多初学者会问的。应力是局部量某个单元应力大不代表整体结构不安全而且应力约束会引入成千上万个局部约束处理起来需要应力惩罚、凝聚函数这类额外手段对入门阶段并不友好。柔度是全局能量量积分意义明确、数学性质好灵敏度形式非常简洁因此成为拓扑优化最常用的目标函数。1.2 离散0/1到连续密度的松弛如果每个单元只用0或1表示“无材料”或“有材料”那就是一个组合优化问题单元数量一上千基本无法直接求解。SIMP方法的核心是把密度放宽到[0,1]连续区间再通过惩罚因子使中间密度在优化过程中变得“不划算”。SIMP插值公式E(ρ) E_min ρ^p (E_0 - E_min)当惩罚因子p3时密度0.5的单元实际弹性模量大约只有E0的12.5%远低于线性插值的50%。这样一来优化器把材料分配到中间密度上带来的刚度收益很低最终结果会自发地向0/1两端靠拢。密度松弛解决的是“能不能算”的问题惩罚因子解决的是“结果像不像结构”的问题两者缺一不可。注意E_min不能取0否则无材料单元的刚度矩阵为零整体K奇异无法求解。常规取E_min1e-9既不会显著影响最优解又保证数值稳定。这一条几乎每次调试都会遇到。2. 四边形网格与有限元求解的细节Q4单元刚度矩阵的组装与自由度约定2.1 为什么选四边形Q4单元二维拓扑优化里常见单元有两种常应变三角形单元CST和四边形四节点单元Q4。CST单元的位移场是线性的单元内应变为常数模拟弯曲时表现偏“刚硬”位移结果偏小精度相对差。Q4单元用双线性位移插值精度高一个层次对弯曲变形也更友好。更重要的是在相同网格尺寸下Q4单元能得到更平滑的灵敏度场。因为拓扑优化每一轮都要根据当前位移场计算每个单元的灵敏度位移场的精度直接决定了梯度方向的质量Q4比CST的梯度噪声小迭代过程中密度分布也更规整。另外同一个四边形区域只需要一个Q4单元而三角形网格通常要拆成两个单元Q4的整体自由度也更少计算代价更低。如果是三维问题我一般建议直接上六面体单元思路与Q4完全一致只是自由度编号更复杂。2.2 单元刚度矩阵与高斯积分Q4单元刚度矩阵可以在自然坐标系下用2×2高斯积分得到。经典Matlab教学实现里8×8的单元刚度矩阵被直接展开成符号表达式代码长这样function [KE] Q4Stiffness(E0, nu) k [ 1/2-nu/6, 1/8nu/8, -1/4-nu/12, -1/83*nu/8, ... -1/4nu/12, -1/8-nu/8, nu/12, 1/8-3*nu/8 ]; KE E0/(1-nu^2) * [ k(1) k(2) k(3) k(4) k(5) k(6) k(7) k(8) k(2) k(1) k(8) k(7) k(6) k(5) k(4) k(3) k(3) k(8) k(1) k(6) k(7) k(4) k(5) k(2) k(4) k(7) k(6) k(1) k(8) k(3) k(2) k(5) k(5) k(6) k(7) k(8) k(1) k(2) k(3) k(4) k(6) k(5) k(4) k(3) k(2) k(1) k(8) k(7) k(7) k(4) k(5) k(2) k(3) k(8) k(1) k(6) k(8) k(3) k(2) k(5) k(4) k(7) k(6) k(1) ]; end这段代码本质上是数值积分后的符号结果。拿到15185期源码时我发现它的单元刚度就是从这里来的验证后发现和数值积分结果完全一致直接用没有问题。要注意的是这段代码假设单元是单位正方形边长为1如果你的设计域单元不是单位尺寸需要做等参变换或直接改用数值积分。2.3 网格编号、自由度映射和总体刚度组装结构化网格的编号有固定套路。假设设计域有nelx×nely个单元节点数为(nelx1)×(nely1)。我习惯用“列优先”编号节点索引和坐标的对应关系是node (nely1)*(elx-1) ely其中elx、ely分别从1开始。每个节点有2个自由度所以全局自由度编号是 2node-1 和 2node。四个节点的自由度向量可以预先算成一个edofMat表避免每轮迭代重复查索引。总体刚度矩阵用稀疏矩阵组装K sparse(2*(nelx1)*(nely1), 2*(nelx1)*(nely1)); for ely 1:nely for elx 1:nelx Ke E(ely, elx) .* KE0; K(edofMat(ely,elx,:), edofMat(ely,elx,:)) ... K(edofMat(ely,elx,:), edofMat(ely,elx,:)) Ke; end end这里E(ely, elx)就是SIMP插值后的单元弹性模量。需要特别强调如果不用sparse而使用满矩阵160×80的网格自由度约为2.6万满矩阵存储需要5.4GB以上拓扑优化要迭代上百轮根本没法跑。稀疏矩阵在这里不是一个优化技巧而是可行性的前提。2.4 边界条件载荷和位移约束的处理拓扑优化中边界条件通常固定在某几个节点上。求解时只保留自由节点的位移分量free setdiff(1:ndof, fixed); U(free) K(free, free) \ F(free);这里用Matlab反斜杠直接做稀疏分解。二维问题规模不大自由度几万到几十万直接分解完全够用。如果网格特别大可以换PCG迭代求解但注意更新密度时要复用同一套预条件子否则每轮重新分解非常浪费。我踩过的一个坑是载荷向量F的长度要和总自由度一致并且载荷必须落在free自由度上。如果把力施加到固定节点相当于约束直接承受了载荷该节点位移恒为0附近单元的灵敏度会出现异常尖峰优化结果往往会在受力点周围堆一大团材料怎么看都不对。3. 灵敏度推导与OC更新密度场每一轮迭代的“方向感”来源3.1 自伴随问题柔度灵敏度的简洁形式设计变量是单元密度ρ_e目标函数柔度c U^T K U其中U是平衡方程K U F的解隐式依赖于ρ。如果直接把c对ρ_e求偏导需要处理位移项的变化但柔度问题恰好有一个特殊性质利用平衡方程和刚度矩阵的对称性可以得到dc/dρ_e -u_e^T (dK_e/dρ_e) u_e这个形式里完全没有位移导数项也就是所谓的“自伴随性质”。因为柔度对每个单元密度的灵敏度只需要该单元自身的位移结果不需要额外求解伴随方程。这也是为什么几乎所有拓扑优化教程都用柔度作为第一个目标函数——换成应力约束就没有这么便宜了。负号代表增加材料会降低柔度增大刚度符合物理直觉。3.2 灵敏度计算的Matlab实现在SIMP插值下dK_e/dρ_e p·ρ_e^(p-1)·(E0 - Emin)·KE0_unit于是每个单元的灵敏度可以写成ce reshape(sum((U(edofMat) * KE0) .* U(edofMat), 2), nely, nelx); dc -p * (E0 - Emin) * rho.^(p-1) .* ce;其中edofMat是单元自由度矩阵KE0是弹性模量为1时的单元刚度矩阵。这行代码的本质是计算每个单元节点位移向量与单元刚度矩阵的双线性型也就是单元应变能。向量化写法避免了逐单元循环在Matlab里提速非常明显。我早期实现用for循环60×20网格跑100轮迭代要二十多秒改成向量化后基本缩短到1秒以内。另外体积约束的灵敏度dv ones(nely, nelx)因为四边形单元归一化后体积相同。如果用了非均匀网格记得替换成实际单元面积。3.3 OC更新准则为什么用它而不是梯度下降拓扑优化常用两类更新算法移动渐近线法MMA和最优化准则法OC。OC法本质上基于拉格朗日乘子对体积约束做KKT推导得到单个设计变量的启发式更新规则B_e -dc_e / (λ · dv_e)ρ_new max(0, ρ - move) 如果 ρ·√B ≤ max(0, ρ - move) ρ_new min(1, ρ move) 如果 ρ·√B ≥ min(1, ρ move) 否则 ρ_new ρ·√B其中move是步长上限一般取0.2。OC的优势在于每个变量更新只涉及一个标量迭代极快非常适合体积约束这个单一约束场景。MMA能处理多约束但参数多入门阶段容易调崩。lambda的求解用二分法l1 0; l2 1e5; while (l2-l1)/(l2l1) 1e-8 lambda 0.5*(l1l2); rho_new OC_update(rho, dc, dv, volfrac, lambda, move); if sum(rho_new(:)) volfrac*nelx*nely l1 lambda; else l2 lambda; end end源码15185期在OC部分的实现基本就是这个路子但它把二分终止条件固定为1e-8在某些算例上会导致内部循环次数偏多。我后面会提一个改进方法。4. 棋盘格问题与灵敏度过滤最容易复现的坑及过滤半径选择4.1 棋盘格是怎么出现的把设计域离散成单元后如果没有额外处理优化结果经常出现黑白相间的棋盘状纹理。根因是有限元离散本身对高频模式的刚度有过高估计优化器会利用这种数值伪影在相邻单元间制造交替0/1分布来“骗”低柔度。这既不是真实结构也不是制造友好的方案必须抑制。最常用的抑制方法是灵敏度过滤把每个单元的灵敏度替换为其邻域内灵敏度的加权平均值从而滤掉高频梯度。权重函数取线性衰减w(i,j) max(0, rmin - dist(i,j))其中dist是两单元中心的欧氏距离rmin是过滤半径。代码里预计算H矩阵H sparse(nelx*nely, nelx*nely); Hs zeros(nelx*nely, 1); for i 1:nelx*nely for j 1:nelx*nely delta dist(i,j); w max(0, rmin - delta); H(i,j) w; Hs(i) Hs(i) w; end end dc(:) (H * dc(:)) ./ Hs(:);由于H只依赖网格和过滤半径可以在迭代开始前一次性算完。双重循环在超大网格上会比较慢实际工程中可以利用网格的规则性只搜索邻域单元避免全局配对。但对入门算例来说上面的稀疏矩阵预计算方式已经够用。4.2 过滤半径怎么选过滤半径的物理意义是结构最小特征尺寸的下限。rmin太小过滤效果不足棋盘格依然会出现rmin太大结构细节被“糊”掉边界变得臃肿。我的经验是结构化四边形网格下rmin取1.2到2.4倍单元尺寸常见算例用1.5倍。对60×20的MBB梁单元尺寸为1rmin1.5时结果比较干净。提示如果改用密度过滤也就是对密度场本身做加权平均而不是过滤灵敏度边界会更光滑但灵敏度公式要相应改成分式求导还要在迭代中额外存储滤波导数。15185期源码默认的是灵敏度过滤改动最小初学建议先用这个。4.3 为什么过滤后还会出现灰度带过滤会让边界单元的灵敏度被周围混合导致密度停留在0.2到0.8之间形成灰度带。迭代中期出现灰度是正常的惩罚因子会随着迭代把它们逐渐推向0或1。但如果200轮之后边界仍是明显的大片灰度通常有两个原因惩罚因子太小比如p2或者move limit太大导致密度反复横跳。想进一步压灰度可以在过滤后加Heaviside投影把低于阈值的密度压向0、高于阈值的压向1。但投影参数需要小心调节太激进会让边界出现锯齿。初学阶段先接受灰度边界把主流程跑通更重要。5. 用15185期源码跑通经典算例参数设定、迭代流程与结果验证5.1 初始参数与设计域我拿到的15185期源码主程序入口是runTopOpt.m默认算例是MBB梁设计域60×20左下角节点固定铰支右下角节点滚动支撑上边中点施加竖直向下的单位集中力。初始参数如下表参数值说明nelx, nely60, 20单元数volfrac0.5材料体积约束penal3.0SIMP惩罚因子rmin1.5灵敏度过滤半径E01材料弹性模量nu0.3泊松比Emin1e-9空单元模量move0.2OC更新步长上限maxiter100最大迭代次数初始设计变量直接取volfrac的均匀密度场。这样做的好处是初始解已经满足体积约束OC二分法第一轮就能收敛。如果取随机初始值前几轮体积约束可能出现抖动而且随机种子不同最终拓扑可能不一样复现性差。5.2 主迭代循环的完整步骤一次迭代的流程是根据当前密度rho计算每个单元的弹性模量E Emin rho.^penal * (E0 - Emin)组装总刚度K求解位移U只保留free自由度计算目标函数柔度c U * K * U输出当前迭代信息计算目标灵敏度dc和体积灵敏度dv对dc做灵敏度过滤用OC公式更新rho二分法求lambda满足体积约束绘制密度云图检查收敛情况。主循环的骨架代码for iter 1:maxiter E Emin rho.^penal * (E0 - Emin); K assembleGlobal(E, KE0, edofMat); U(free) K(free, free) \ F(free); c(iter) U * K * U; dc computeSensitivity(U, rho, KE0, edofMat, penal, E0, Emin); dc filterSensitivity(dc, H, Hs); rho OC_update(rho, dc, dv, volfrac, move); if iter 10 abs(c(iter)-c(iter-1))/c(iter-1) 1e-3 break; end end收敛判据建议用最近两轮目标函数相对变化量小于千分之一。绝对变化量的阈值在不同规模问题上不好统一相对变化量基本与网格规模无关。实际运行中这个MBB算例一般50到80轮就能收敛到稳定拓扑。5.3 经典结果长什么样MBB梁优化完成后结构会呈现对称的两根类桁架斜撑结合中间主梁的形状左右镜像材料集中在主要传力路径上其余区域密度接近0。把结果与文献中的参考形状对比主要传力路径基本一致差异只在边界灰度和细节尖角。我多次实验发现一个规律过滤半径调大到rmin3结果会更“胖”构件更粗壮rmin1.2时结构更纤细但容易出现尖角。复现算例时不要执着于某一组精确参数只要体积分数一致、传力路径清晰主链路基本就算跑通了。想验证程序的正确性还可以对比目标函数曲线——正常情况是单调下降、后期趋于平稳如果曲线上下跳动剧烈说明参数需要调整。6. 源码改造与常见异常排查我的调试记录6.1 改边界条件和载荷最容易踩的坑拿到源码第一件事通常是换自己的算例比如把MBB梁改成悬臂梁。很多人直接改受力节点坐标结果发现优化结果乱七八糟出现局部大片材料。我排查过的常见问题有三个。第一固定自由度集合和载荷自由度冲突。载荷必须落在free集合内一旦加在固定节点上该节点位移恒为0灵敏度在该单元附近会出现异常尖峰结构会在受力点堆料。第二载荷向量F的维度没同步。某些改法只在固定节点处增删位移没同步调FMatlab会报维度不匹配这个比较低级但很容易忽略。第三节点编号方向。列优先编号和行优先编号差别很大如果载荷节点位置写反优化结果会上下颠倒虽然柔度值一样但结构形状就不对了。建议在改算例之前先画一下设计域节点编号图确认载荷和支撑节点的全局编号再动代码。6.2 OC二分法不收敛lambda区间的坑OC更新里的体积约束方程是非线性的二分法需要保证lambda区间覆盖真实解。15185期源码默认l10、l21e5对大多数算例够用但如果你把move limit调大、又把体积分数设得很小比如0.2有时会出现二分区间不够大的情况体积约束始终无法满足。我的做法是每轮迭代前动态扩展上限先算一个粗略估计lambda_est max(10 * max(-dc ./ dv))然后以它为基准设置l2的初始值。这样无论体积分数高低都能快速定位实测可以减少不少二分次数整体迭代速度提升10%左右。另外OC更新公式中可以对灵敏度做历史平均来抑制振荡。如果发现目标函数曲线上下跳动给灵敏度做一个加权平均是稳妥的缓解手段。15185期源码没有做这个阻尼所以大步长下偶尔能看到目标函数跳动这是正常的可以适当把move降到0.1再观察。6.3 图像输出与结果后处理的小技巧源码里默认用imagesc显示密度场但imagesc会上下翻转y轴我习惯加一句axis equal和axis off让图片按真实比例显示。如果想导出清晰位图可以用imwrite配合colormap gray把密度矩阵归一化后保存方便写实验报告。还可以在结果中统计实际体积分数和0.5阈值以上的材料占比。两个数差异过大说明边界灰度比例高一般迭代充分时两者应接近到5%以内如果差异超过10%就得检查是不是收敛过早或者惩罚因子设置有问题。这些经验在我自己调试其他算例时反复用到尤其是从源码包起步做改造的场景。只要你把目标函数、灵敏度、过滤、更新这几块的关系理清换算例、换边界条件、甚至换成三维六面体单元都是按同样的逻辑去改不需要推倒重来。
阅读完成 · 觉得有帮助？
