做结构拓扑优化的人对“应力”这两个字往往又爱又恨。柔顺性目标是个全局量灵敏度可以用自伴随公式一行写完可应力是局部量一个应力集中点就能决定整个结构是否失效。我前阵子在做3D应力敏感度分析与拓扑优化时选择用p-范数把“全局应力”光滑化再用伴随方法在Matlab里做有限元分析和灵敏度推导。论文里的推导看着干净真落到代码上才发现3D模型带来的内存问题、B矩阵组装、伴随方程右侧向量、还有低密度单元应力的处理每一步都可能让数值结果面目全非。这篇文章把能跑通的一套实现和我踩过的坑都整理出来希望能给同样用Matlab做应力相关拓扑优化的朋友一个参考。1. 为什么应力不能像柔顺性那样省心地进入优化循环1.1 柔顺性为什么是“自伴随”的幸运儿拓扑里最常优化的目标函数是柔顺性 C F^T u约束条件是平衡方程 K u F。对柔顺性求密度灵敏度时大多数人背过这个结果dC/dρ_e -u_e^T (dK_e/dρ_e) u_e。这里不需要额外解伴随方程原因是目标函数对位移的导数恰好等于载荷向量F而K^{-1}F就是位移u本身。这个性质叫自伴随我第一次学的时候觉得是天上掉馅饼因为整个灵敏度求解几乎零成本。但应力目标没有这个好运气。应力是局部量必须先通过某种方式聚合成一个全局测度再对密度求导目标函数对位移的导数也不再是F而是一个由单元应力导出的向量。也就是说每做一次优化迭代至少要解两次线性方程组一次用于位移一次用于伴随向量。这个成本在二维网格上还不明显一旦进入3D网格自由度成倍增长求解次数直接决定程序能不能在合理时间内跑完。1.2 应力目标的三个“老熟人”麻烦第一个麻烦是局部性与全局性的冲突。结构某一点的应力超标理论上就可能导致失效。但连续体拓扑优化不可能逐点控制应力只能用单元或高斯点应力代替再把这些局部的值用一个光滑标量聚合起来。聚得不好优化器就会“钻空子”比如把高应力区域的材料密度改成中间值从而骗过全局测度。第二个麻烦是最大应力函数不可导。max函数是一个非光滑函数在最大值对应的单元切换时导数不连续。梯度优化算法在这种目标下会出现强烈的锯齿震荡甚至直接发散。所以工程上几乎不会直接做min max问题而是用p-范数、KS函数这类光滑聚合函数去逼近最大应力。第三个麻烦更隐蔽叫奇异解。拓扑优化里的密度为0或接近0的孔洞区域理论上可以承受任意大的真实应力但因为这些区域几乎没有材料有限元应力已经被密度惩罚压得很低所以优化器会倾向于把高应力区域的“中间密度”保留成灰色而不是形成清晰的0/1拓扑。要处理它通常需要引入应力松弛项。网上很多纯代码实现不谈这个问题实际跑起来就会看到灰色单元一片一片地冒出来。1.3 p-范数如何“光滑地逼近最大应力”p-范数聚合应力最简单也最常用。假设每个单元的von Mises应力是σ_vm,e那么全局应力测度可以写成σ_PN ( Σ_e (σ_vm,e)^p )^(1/p)当p趋向无穷大时这个测度收敛到最大的单元应力当p2时它更接近均方根水平。换句话说p值是“选代表”的尺度p越小所有单元应力都在平均拉票p越大最大应力的单元几乎独裁。拓扑优化里常用p6到p12既保证光滑性又能体现危险部位的影响。这里还要注意一个问题如果不除单元数σ_PN会随网格加密而增大因为求和项数变多。实际实现时我会建议把目标函数写成“每个单元平均”的形式即g ( (1/N_ele) Σ_e (σ_vm,e)^p )^(1/p)这样网格粗细对目标量级影响小一些也好设定p值和收敛判断。2. 伴随法应力灵敏度推导三项都要算漏一项整个梯度就错2.1 离散模型与符号先把离散问题写清楚。设计变量是单元密度ρ_eSIMP插值下单元弹性模量为 E_e ρ_e^penal E0penal一般取3。有限元平衡方程为 K(ρ) u F其中单元刚度矩阵是k_e ρ_e^penal k0_ek0_e是密度为1时的单元刚度矩阵。单元应力按本构关系和几何关系得到σ_e D_e B_e u_e ρ_e^penal D0 B_e u_e等效应力用von Mises公式计算可以写成二次型s_e σ_vm,e sqrt( σ_e^T V σ_e )这里V是一个6×6的常系数矩阵它把应力分量组合成Mises等效应力。目标函数定义为 p-范数全局应力g ( Σ_e (s_e)^p )^(1/p)对密度求导最关键的一点是s_e既通过u间接依赖密度也通过E_e直接依赖密度。因此总导数有三部分来源少算任何一个有限差分验证都过不了。2.2 链式法则与伴随方程把g对ρ_e求导展开dg/dρ_e ∂g/∂ρ_e (∂g/∂u)^T (du/dρ_e)其中∂g/∂u是一个全局向量由所有单元贡献我把它记为a。再看平衡方程对密度求导K (du/dρ_e) (dK/dρ_e) u 0所以 du/dρ_e -K^{-1} (dK/dρ_e) u。代回去得到dg/dρ_e ∂g/∂ρ_e - a^T K^{-1} (dK/dρ_e) u如果直接算K^{-1}向量成本太高。于是定义伴随向量λ满足K^T λ a如果K对称就写成K λ a。把λ代入最终灵敏度为dg/dρ_e ∂g/∂ρ_e - λ_e^T (dK_e/dρ_e) u_e这就是伴随法的核心。和柔顺性灵敏度相比形式上只多了一个∂g/∂ρ_e但实质性变化是a不再等于载荷向量F必须单独组装。2.3 三项灵敏度的具体形式第一项是直接依赖∂g/∂ρ_e。它来自应力σ_e对密度ρ_e的直接依赖因为D_e ρ_e^penal D0。逐单元计算时∂s_e/∂ρ_e (1/s_e) σ_e^T V · (penal ρ_e^{penal-1} D0 B_e u_e)再乘以全局测度导数∂g/∂s_e g^{1-p} s_e^{p-1}。第二项是a的组装。a ∂g/∂u对单元e贡献是g^{1-p} s_e^{p-1} · ∂s_e/∂u_e而∂s_e/∂u_e (1/s_e) σ_e^T V D_e B_e。这些单元贡献要叠加到全局向量a的对应自由度上不能只放单元内部导数值。第三项是刚度矩阵变化项。dK_e/dρ_e penal ρ_e^{penal-1} k0_e把它乘上λ_e和u_e就是伴随项对灵敏度的贡献。这三项中第二项组装最容易写错。因为∂s_e/∂u_e是一个1×24行向量对应Hex8单元的24个自由度而在Matlab里组装a时要把它按edof索引加到对应的行位置而不是用矩阵赋值覆盖。2.4 用有限差分验证梯度这一步省不了不管推导多自信我都会先在2D小网格上做一次有限差分验证。做法很简单对某个单元密度加一个小扰动δ重新求解K u F重新计算目标函数g得到数值梯度再和解析梯度比较。% 验证灵敏度的小模板 delta 1e-6; dofs size(K,1); U K \ F; [g, dg] computeObjSens(U, rho); % 解析灵敏度 rho_pert rho; rho_pert(e) rho(e) delta; [U_pert] solveFE(rho_pert); g_pert computeObj(rho_pert, U_pert); dg_num (g_pert - g) / delta; rel_error abs(dg(e) - dg_num) / abs(dg_num);3D模型里自由度大扰动一个单元也意味着要重新组装K和求解一次所以验证时可以选一个小网格比如10×10×5把所有单元都测一遍画一个误差分布图。我通常要求相对误差在1e-3以下。如果伴随项漏掉一个人误差会成几十倍地跳。3. 3D有限元与伴随求解Matlab里真正让人头疼的地方3.1 自由度规模与稀疏矩阵组装3D六面体网格的规模比2D大一个量级。以40×20×10的网格为例节点数是(41×21×11)自由度大约28万。如果拿满矩阵存KMatlab直接内存溢出即使没溢出求解也慢得让人怀疑人生。所以第一步就是必须用稀疏矩阵并且在组装时用i、j、s三向量一次性构造K不要在一个循环里不断做K(edof,edof)K(edof,edof)KE。稀疏组装的套路在三维里依然成立。先根据每个单元的edof向量生成64×64的索引矩阵然后把所有单元的i、j、s拼成三个长向量最后调用一次sparse(i,j,s,ndof,ndof)。这样组装速度比边循环边赋值快很多内存也可控。求解fx时我建议直接用K\F而不建议用inv(K)*F。inv会破坏稀疏结构内存和耗时都不可接受。如果网格规模更大可以用预处理共轭梯度pcg配一个不完全Cholesky预处理因子。对于几万自由度的小算例K\F足够稳定。3.2 B矩阵、D矩阵与单元应力计算3D应力的计算核心是Hex8单元的B矩阵。下面给出一个中心高斯点的B矩阵实现节点坐标按标准八节点顺序排列function [B, detJ] hex8_B_center(node_coord) % node_coord: 8x3 矩阵按 Hex8 标准节点顺序排列 % 顺序下表面逆时针1-4上表面逆时针5-8 dN 1/8 * ... [-1 1 1 -1 -1 1 1 -1; -1 -1 1 1 -1 -1 1 1; -1 -1 -1 -1 1 1 1 1]; J dN * node_coord; % 3x3 Jacobian detJ det(J); invJ J \ eye(3); B zeros(6, 24); for i 1:8 g invJ * dN(:, i); idx (i-1)*3 1 : i*3; B(1:3, idx) diag(g); B(4, idx) [g(2), g(1), 0]; % gamma_xy B(5, idx) [0, g(3), g(2)]; % gamma_yz B(6, idx) [g(3), 0, g(1)]; % gamma_xz end物理应力是σ D * B * u_e其中D是3D各向同性弹性矩阵。注意如果用中心高斯点单元刚度积分的体积权重是detJ乘以8因为三维高斯点权重是2×2×28。很多人在小网格上忘记这个8柔顺性灵敏度会差一个体积因子。von Mises的组合矩阵V可以这样定义Vmat [1.0 -0.5 -0.5 0 0 0; -0.5 1.0 -0.5 0 0 0; -0.5 -0.5 1.0 0 0 0; 0 0 0 3 0 0; 0 0 0 0 3 0; 0 0 0 0 0 3];然后把每个单元的σ_vm存起来用于后续p-范数目标计算和灵敏度组装。3.3 伴随求解就是“再解一次线性方程组”伴随方程K λ a的求解本质上是把K当成系数矩阵再解一次线性方程。由于K对称可以复用同一个因式分解结果。如果用的是K\FMatlab会自动做分解如果网格规模大到需要手动控制可以使用分解一次[L, U, P] lu(K); u U \ (L \ (P * F)); lambda U \ (L \ (P * a));这样做的好处是位移求解和伴随求解可以共用同一个LU分解第二次求解量只剩回代。3D优化迭代中主求解和伴随求解是每次迭代的两大成本省一次分解很关键。3.4 过滤与棋盘格要在3D里重新写一遍2D的灵敏度过滤很容易写3D则需要考虑z方向的邻居单元。过滤半径rmin通常取1.5倍单元尺寸如果网格是均匀立方体可以直接以索引差作为距离。一个简单的实现是把所有单元的坐标中心预先算出来再对每个单元搜索距离小于rmin的邻居。这个操作在网格较小时直接用三重循环也行但网格大了会很慢建议用k近邻搜索或者分桶预处理。过滤是拓扑优化里抑制棋盘格的最基本手段。应力目标下棋盘格比柔顺性目标更严重因为应力对局部材料分布极其敏感。不加过滤的程序跑出来的结构看起来像被小尺度噪声覆盖应力云图也一团乱麻。4. 代码骨架3D拓扑优化主循环、目标灵敏度与OC更新4.1 主程序结构下面是我常用的主循环框架。目标函数是最小化p-范数全局应力体积约束固定优化更新用OC准则。为了让代码清晰我把有限元求解、目标灵敏度计算、过滤器、OC更新都拆成独立函数。function topopt3D_stress(nelx, nely, nelz, volfrac, penal, p, rmin, move) % 3D拓扑优化以p-范数应力和最小为目标的简化实现 Vol ones(nelx*nely*nelz, 1) * volfrac; rho Vol; change 1; hist []; while change 0.01 % 有限元求解 [U, K, edofMat] solve_FE(nelx, nely, nelz, rho, penal); % 目标函数与解析灵敏度 [g, dg] compute_Stress_Obj_Sens(nelx, nely, nelz, ... U, rho, penal, p, edofMat); % 灵敏度过滤 dg filter3D(nelx, nely, nelz, rmin, rho, dg); % OC更新 rho_new OC_update(rho, dg, volfrac, move); change max(abs(rho_new - rho)); rho rho_new; hist(end1, :) [g, change]; fprintf(iter %d, g %.4e, change %.4f\n, ... length(hist), g, change); end end这里我故意没有把solve_FE写全因为标准的3D有限元单元刚度矩阵、边界条件设置等和常规3D拓扑优化完全一样。核心差异集中在目标函数和灵敏度函数里也就是下面一小节。4.2 目标函数与灵敏度函数计算目标和灵敏度的函数是整个程序的关键。我按前面推导的三个部分来写先算每个单元应力并聚合目标再组装伴随右侧向量a接着求解伴随方程最后叠加上直接项和刚度项。function [g, dg] compute_Stress_Obj_Sens(nelx, nely, nelz, ... U, rho, penal, p, edofMat) ndof 3 * (nelx1) * (nely1) * (nelz1); ne nelx * nely * nelz; svm zeros(ne, 1); sig_store zeros(ne, 6); Bstore cell(ne, 1); detJstore zeros(ne, 1); E0 1.0; nu 0.3; D0 isotropicD(E0, nu); % 第一遍计算单元应力、von Mises应力、目标函数值 for e 1:ne node_coord getNodeCoord(e, nelx, nely, nelz); [B, detJ] hex8_B_center(node_coord); D rho(e)^penal * D0; Ue U(edofMat(e, :)); sigma D * B * Ue; svm(e) sqrt(sigma * Vmat * sigma); sig_store(e, :) sigma; Bstore{e} B; detJstore(e) detJ; end g (sum(svm.^p))^(1/p); % 第二遍组装伴随右侧向量a以及直接项 a zeros(ndof, 1); dg_direct zeros(ne, 1); for e 1:ne s svm(e); coef g^(1-p) * s^(p-1); B Bstore{e}; Ue U(edofMat(e, :)); sigma sig_store(e, :); D0B D0 * B; if s 1e-12 dsdU (sigma * Vmat) / s * (rho(e)^penal * D0B); dsdrho (sigma * Vmat) / s * ... (penal * rho(e)^(penal-1) * D0B * Ue); else dsdU zeros(1, 24); dsdrho 0; end a(edofMat(e, :)) a(edofMat(e, :)) coef * dsdU; dg_direct(e) coef * dsdrho; end % 伴随求解 lambda K \ a; % 注意这里的K需要从solve_FE传进来或重复组装 % 第三遍叠加伴随项 dg dg_direct; for e 1:ne B Bstore{e}; k0e B * D0 * B * detJstore(e) * 8; Ue U(edofMat(e, :)); lambda_e lambda(edofMat(e, :)); dg(e) dg(e) - lambda_e * (penal * rho(e)^(penal-1) * k0e * Ue); end end这段代码里K没有直接传入实际使用时要把它从solve_FE里带出来。也可以改成在函数内部重新组装一次K再解伴随方程但那样会浪费一次组装时间。更好的写法是在主循环里把K传到compute_Stress_Obj_Sens里或者像我前面说的把LU分解结果也传进去。4.3 OC更新与过滤函数OC更新在柔顺性最小化里已经很成熟应力目标下也可以用但要注意两个地方一是灵敏度可能出现正值因为增加材料未必在所有位置都降低全局应力二是目标函数非凸OC更新很容易震荡。我一般把move限制在0.1到0.2之间并且加一个对灵敏度的裁剪如果-dg/dv小于等于0就把它设成一个很小的正值避免sqrt里出现负数。function rho_new OC_update(rho, dg, volfrac, move) l1 0; l2 1e9; dv ones(size(rho)); while abs(l2 - l1) 1e-5 lmid 0.5 * (l1 l2); s -dg ./ dv / lmid; s(s 0) 1e-6; rho_new max(0, max(rho - move, min(1, min(rho move, rho .* sqrt(s))))); if sum(rho_new) - volfrac * length(rho) 0 l1 lmid; else l2 lmid; end end end过滤函数本身和柔顺性优化里的标准实现几乎一样只是三维邻居索引需要额外处理z方向。需要提醒的是过滤可以施加在灵敏度上也可以施加在密度场上。如果做密度过滤需要在灵敏度里再乘一个滤波矩阵的转置我在教程里先用灵敏度过滤代码最简单也最容易复现。5. 算例与典型结果从3D悬臂梁看应力优化行为5.1 模型设置与载荷我用一个3D L形支架来做算例。这类结构的内凹角会产生强烈的应力集中是检验应力相关拓扑优化的经典模型。网格取30×20×10共6000个六面体单元左端面固定右上方的局部位置施加向下的载荷体积约束0.3。材料参数用无量纲的E01泊松比0.3。p值取8惩罚系数3过滤半径1.5OC移动步长0.15。在这个规模下每轮迭代需要解两次稀疏线性方程Matlab在我的机器上大概需要两到三分钟完成全过程。如果内存紧张可以把网格缩小到20×15×8结果趋势是一样的。5.2 优化结果解读优化后的拓扑会比较直观内凹角附近会保留下更多材料形成一个近似圆角的过渡区域远离传力路径的角落被挖空结构整体看起来比柔顺性优化更“粗壮”因为应力目标倾向于把材料铺设在受力路径上而不是单纯追求最小柔顺性。p-范数应力值在优化过程中通常先快速下降随后趋于平稳。但因为应力目标不是凸问题收敛曲线会有一定波动。我建议不要只看最后一步要看整个历史的趋势。如果曲线在某一代后反复横跳就减小move值或者改用MMA更新。5.3 收敛曲线与运行时间下面是典型的收敛行为前20代目标值下降幅度最大20到60代缓慢下降60代之后基本稳定。虽然OC更新的数学保证比MMA弱但在这个算例里配合过滤和小步长依然能得到清晰的0/1结构。需要注意的是p越大目标函数越接近“最大应力”曲线抖动越明显我的经验是p8在稳定性和信息量之间比较平衡。这种3D算例最耗时的环节依次是稀疏组装、位移求解、伴随求解。如果用LU分解一次性处理两个方程组总时间能比重复求解少百分之三十左右。代码里一定要避免在循环内反复调用inv或者full否则6000单元的网格就能让Matlab卡到无法忍受。6. 调参与避坑p值、惩罚、棋盘格、奇异解6.1 p值的选择p值直接影响优化结果的性格。p2时所有应力被“平均”看待优化器倾向于把应力分布磨平但最大应力可能并不可控p20时几乎等于在追着最大应力跑更新量很容易剧烈震荡。不同文献的取值范围差别很大但我自己实测下来p在8到12之间最顺手。可以先跑一个p8的版本看趋势再逐步增大。如果目标函数里用了平均形式的p-范数除以单元数那么p值带来的量级变化会小很多收敛判断也更稳定。否则网格加密后同样的p值目标值会大不少容易被误判成没有收敛。6.2 SIMP惩罚与应力松弛应力优化里最微妙的是低密度单元的应力。SIMP模型下单元刚度是ρ^penal应力也是ρ^penal所以低密度单元的内力反馈很弱应力值很低。这看起来是合理的但它让优化器可以把高应力单元“稀释”成低密度从而绕过应力惩罚。真实材料不存在这种取巧行为孔洞区域只要有一点材料局部应力可能非常大。常用的处理办法是应力松弛比如对单元应力再除一个密度相关项σ_physical,e ρ_e^{penal-1} σ_e这样密度下降时应力不会同比例下降低密度区域的应力贡献会被放大从而破坏灰色单元绕过约束的路径。这篇文章的代码里没有加这项是为了先讲清楚伴随方法和p-范数基础实际工程中建议加上。加了之后目标函数和灵敏度的推导需要同步改直接项里会多出(d(ρ^{penal-1})/dρ)的因子。6.3 过滤半径与move限制过滤半径rmin小于1.5个单元尺寸时棋盘格压制不彻底rmin太大结构又会失去细节。3D网格中我习惯用1.5到2.0倍单元尺寸。密度场滤波比灵敏度滤波更平滑但代码实现稍重需要构造线性滤波矩阵建议先用灵敏度滤波入门遇到严重棋盘格再升级。OC更新对move值非常敏感。move0.2在柔顺性优化里很常见但在应力目标下容易震荡我通常从0.15起步如果目标历史平滑就慢慢加到0.2。每轮更新后还要检查体积约束是否满足否则二分会失效。6.4 一个小型的检查清单我每次跑新的3D应力算例前都会按下面这个清单快速过一遍第一确认B矩阵的行列式detJ没有出现负数说明节点顺序错了第二用一个均匀密度场算一次目标函数看是否和手算量级一致第三用有限差分验证随机选出的几个单元灵敏度第四观察第一次迭代的密度更新方向如果高应力区域的材料反而减少说明某个符号搞反了第五检查p-范数目标值是否随网格加密而异常增大。这套检查听起来琐碎但能省下好几天的debug时间。尤其是第一次从2D转到3D时B矩阵的维度、edof索引范围、稀疏组装行数任何一个地方错一点程序都会在某个角落给出离谱的灵敏度而优化器会把错误当成正确方向一路跑偏。最后再分享两个小技巧一是灵敏度验证时扰动密度不要取太小1e-6比较合适因为3D应力目标对舍入误差更敏感二是如果3D网格大到求解吃力先把z方向厚度压到两层退化成一个“准2D”模型来调通所有代码再恢复完整厚度。把这两步做好3D应力敏感度分析和拓扑优化就不会再像一开始那样劝退人了。希望这份记录对你有用。
3D应力灵敏度分析与拓扑优化：p-范数与伴随方法的Matlab实现
做结构拓扑优化的人对“应力”这两个字往往又爱又恨。柔顺性目标是个全局量灵敏度可以用自伴随公式一行写完可应力是局部量一个应力集中点就能决定整个结构是否失效。我前阵子在做3D应力敏感度分析与拓扑优化时选择用p-范数把“全局应力”光滑化再用伴随方法在Matlab里做有限元分析和灵敏度推导。论文里的推导看着干净真落到代码上才发现3D模型带来的内存问题、B矩阵组装、伴随方程右侧向量、还有低密度单元应力的处理每一步都可能让数值结果面目全非。这篇文章把能跑通的一套实现和我踩过的坑都整理出来希望能给同样用Matlab做应力相关拓扑优化的朋友一个参考。1. 为什么应力不能像柔顺性那样省心地进入优化循环1.1 柔顺性为什么是“自伴随”的幸运儿拓扑里最常优化的目标函数是柔顺性 C F^T u约束条件是平衡方程 K u F。对柔顺性求密度灵敏度时大多数人背过这个结果dC/dρ_e -u_e^T (dK_e/dρ_e) u_e。这里不需要额外解伴随方程原因是目标函数对位移的导数恰好等于载荷向量F而K^{-1}F就是位移u本身。这个性质叫自伴随我第一次学的时候觉得是天上掉馅饼因为整个灵敏度求解几乎零成本。但应力目标没有这个好运气。应力是局部量必须先通过某种方式聚合成一个全局测度再对密度求导目标函数对位移的导数也不再是F而是一个由单元应力导出的向量。也就是说每做一次优化迭代至少要解两次线性方程组一次用于位移一次用于伴随向量。这个成本在二维网格上还不明显一旦进入3D网格自由度成倍增长求解次数直接决定程序能不能在合理时间内跑完。1.2 应力目标的三个“老熟人”麻烦第一个麻烦是局部性与全局性的冲突。结构某一点的应力超标理论上就可能导致失效。但连续体拓扑优化不可能逐点控制应力只能用单元或高斯点应力代替再把这些局部的值用一个光滑标量聚合起来。聚得不好优化器就会“钻空子”比如把高应力区域的材料密度改成中间值从而骗过全局测度。第二个麻烦是最大应力函数不可导。max函数是一个非光滑函数在最大值对应的单元切换时导数不连续。梯度优化算法在这种目标下会出现强烈的锯齿震荡甚至直接发散。所以工程上几乎不会直接做min max问题而是用p-范数、KS函数这类光滑聚合函数去逼近最大应力。第三个麻烦更隐蔽叫奇异解。拓扑优化里的密度为0或接近0的孔洞区域理论上可以承受任意大的真实应力但因为这些区域几乎没有材料有限元应力已经被密度惩罚压得很低所以优化器会倾向于把高应力区域的“中间密度”保留成灰色而不是形成清晰的0/1拓扑。要处理它通常需要引入应力松弛项。网上很多纯代码实现不谈这个问题实际跑起来就会看到灰色单元一片一片地冒出来。1.3 p-范数如何“光滑地逼近最大应力”p-范数聚合应力最简单也最常用。假设每个单元的von Mises应力是σ_vm,e那么全局应力测度可以写成σ_PN ( Σ_e (σ_vm,e)^p )^(1/p)当p趋向无穷大时这个测度收敛到最大的单元应力当p2时它更接近均方根水平。换句话说p值是“选代表”的尺度p越小所有单元应力都在平均拉票p越大最大应力的单元几乎独裁。拓扑优化里常用p6到p12既保证光滑性又能体现危险部位的影响。这里还要注意一个问题如果不除单元数σ_PN会随网格加密而增大因为求和项数变多。实际实现时我会建议把目标函数写成“每个单元平均”的形式即g ( (1/N_ele) Σ_e (σ_vm,e)^p )^(1/p)这样网格粗细对目标量级影响小一些也好设定p值和收敛判断。2. 伴随法应力灵敏度推导三项都要算漏一项整个梯度就错2.1 离散模型与符号先把离散问题写清楚。设计变量是单元密度ρ_eSIMP插值下单元弹性模量为 E_e ρ_e^penal E0penal一般取3。有限元平衡方程为 K(ρ) u F其中单元刚度矩阵是k_e ρ_e^penal k0_ek0_e是密度为1时的单元刚度矩阵。单元应力按本构关系和几何关系得到σ_e D_e B_e u_e ρ_e^penal D0 B_e u_e等效应力用von Mises公式计算可以写成二次型s_e σ_vm,e sqrt( σ_e^T V σ_e )这里V是一个6×6的常系数矩阵它把应力分量组合成Mises等效应力。目标函数定义为 p-范数全局应力g ( Σ_e (s_e)^p )^(1/p)对密度求导最关键的一点是s_e既通过u间接依赖密度也通过E_e直接依赖密度。因此总导数有三部分来源少算任何一个有限差分验证都过不了。2.2 链式法则与伴随方程把g对ρ_e求导展开dg/dρ_e ∂g/∂ρ_e (∂g/∂u)^T (du/dρ_e)其中∂g/∂u是一个全局向量由所有单元贡献我把它记为a。再看平衡方程对密度求导K (du/dρ_e) (dK/dρ_e) u 0所以 du/dρ_e -K^{-1} (dK/dρ_e) u。代回去得到dg/dρ_e ∂g/∂ρ_e - a^T K^{-1} (dK/dρ_e) u如果直接算K^{-1}向量成本太高。于是定义伴随向量λ满足K^T λ a如果K对称就写成K λ a。把λ代入最终灵敏度为dg/dρ_e ∂g/∂ρ_e - λ_e^T (dK_e/dρ_e) u_e这就是伴随法的核心。和柔顺性灵敏度相比形式上只多了一个∂g/∂ρ_e但实质性变化是a不再等于载荷向量F必须单独组装。2.3 三项灵敏度的具体形式第一项是直接依赖∂g/∂ρ_e。它来自应力σ_e对密度ρ_e的直接依赖因为D_e ρ_e^penal D0。逐单元计算时∂s_e/∂ρ_e (1/s_e) σ_e^T V · (penal ρ_e^{penal-1} D0 B_e u_e)再乘以全局测度导数∂g/∂s_e g^{1-p} s_e^{p-1}。第二项是a的组装。a ∂g/∂u对单元e贡献是g^{1-p} s_e^{p-1} · ∂s_e/∂u_e而∂s_e/∂u_e (1/s_e) σ_e^T V D_e B_e。这些单元贡献要叠加到全局向量a的对应自由度上不能只放单元内部导数值。第三项是刚度矩阵变化项。dK_e/dρ_e penal ρ_e^{penal-1} k0_e把它乘上λ_e和u_e就是伴随项对灵敏度的贡献。这三项中第二项组装最容易写错。因为∂s_e/∂u_e是一个1×24行向量对应Hex8单元的24个自由度而在Matlab里组装a时要把它按edof索引加到对应的行位置而不是用矩阵赋值覆盖。2.4 用有限差分验证梯度这一步省不了不管推导多自信我都会先在2D小网格上做一次有限差分验证。做法很简单对某个单元密度加一个小扰动δ重新求解K u F重新计算目标函数g得到数值梯度再和解析梯度比较。% 验证灵敏度的小模板 delta 1e-6; dofs size(K,1); U K \ F; [g, dg] computeObjSens(U, rho); % 解析灵敏度 rho_pert rho; rho_pert(e) rho(e) delta; [U_pert] solveFE(rho_pert); g_pert computeObj(rho_pert, U_pert); dg_num (g_pert - g) / delta; rel_error abs(dg(e) - dg_num) / abs(dg_num);3D模型里自由度大扰动一个单元也意味着要重新组装K和求解一次所以验证时可以选一个小网格比如10×10×5把所有单元都测一遍画一个误差分布图。我通常要求相对误差在1e-3以下。如果伴随项漏掉一个人误差会成几十倍地跳。3. 3D有限元与伴随求解Matlab里真正让人头疼的地方3.1 自由度规模与稀疏矩阵组装3D六面体网格的规模比2D大一个量级。以40×20×10的网格为例节点数是(41×21×11)自由度大约28万。如果拿满矩阵存KMatlab直接内存溢出即使没溢出求解也慢得让人怀疑人生。所以第一步就是必须用稀疏矩阵并且在组装时用i、j、s三向量一次性构造K不要在一个循环里不断做K(edof,edof)K(edof,edof)KE。稀疏组装的套路在三维里依然成立。先根据每个单元的edof向量生成64×64的索引矩阵然后把所有单元的i、j、s拼成三个长向量最后调用一次sparse(i,j,s,ndof,ndof)。这样组装速度比边循环边赋值快很多内存也可控。求解fx时我建议直接用K\F而不建议用inv(K)*F。inv会破坏稀疏结构内存和耗时都不可接受。如果网格规模更大可以用预处理共轭梯度pcg配一个不完全Cholesky预处理因子。对于几万自由度的小算例K\F足够稳定。3.2 B矩阵、D矩阵与单元应力计算3D应力的计算核心是Hex8单元的B矩阵。下面给出一个中心高斯点的B矩阵实现节点坐标按标准八节点顺序排列function [B, detJ] hex8_B_center(node_coord) % node_coord: 8x3 矩阵按 Hex8 标准节点顺序排列 % 顺序下表面逆时针1-4上表面逆时针5-8 dN 1/8 * ... [-1 1 1 -1 -1 1 1 -1; -1 -1 1 1 -1 -1 1 1; -1 -1 -1 -1 1 1 1 1]; J dN * node_coord; % 3x3 Jacobian detJ det(J); invJ J \ eye(3); B zeros(6, 24); for i 1:8 g invJ * dN(:, i); idx (i-1)*3 1 : i*3; B(1:3, idx) diag(g); B(4, idx) [g(2), g(1), 0]; % gamma_xy B(5, idx) [0, g(3), g(2)]; % gamma_yz B(6, idx) [g(3), 0, g(1)]; % gamma_xz end物理应力是σ D * B * u_e其中D是3D各向同性弹性矩阵。注意如果用中心高斯点单元刚度积分的体积权重是detJ乘以8因为三维高斯点权重是2×2×28。很多人在小网格上忘记这个8柔顺性灵敏度会差一个体积因子。von Mises的组合矩阵V可以这样定义Vmat [1.0 -0.5 -0.5 0 0 0; -0.5 1.0 -0.5 0 0 0; -0.5 -0.5 1.0 0 0 0; 0 0 0 3 0 0; 0 0 0 0 3 0; 0 0 0 0 0 3];然后把每个单元的σ_vm存起来用于后续p-范数目标计算和灵敏度组装。3.3 伴随求解就是“再解一次线性方程组”伴随方程K λ a的求解本质上是把K当成系数矩阵再解一次线性方程。由于K对称可以复用同一个因式分解结果。如果用的是K\FMatlab会自动做分解如果网格规模大到需要手动控制可以使用分解一次[L, U, P] lu(K); u U \ (L \ (P * F)); lambda U \ (L \ (P * a));这样做的好处是位移求解和伴随求解可以共用同一个LU分解第二次求解量只剩回代。3D优化迭代中主求解和伴随求解是每次迭代的两大成本省一次分解很关键。3.4 过滤与棋盘格要在3D里重新写一遍2D的灵敏度过滤很容易写3D则需要考虑z方向的邻居单元。过滤半径rmin通常取1.5倍单元尺寸如果网格是均匀立方体可以直接以索引差作为距离。一个简单的实现是把所有单元的坐标中心预先算出来再对每个单元搜索距离小于rmin的邻居。这个操作在网格较小时直接用三重循环也行但网格大了会很慢建议用k近邻搜索或者分桶预处理。过滤是拓扑优化里抑制棋盘格的最基本手段。应力目标下棋盘格比柔顺性目标更严重因为应力对局部材料分布极其敏感。不加过滤的程序跑出来的结构看起来像被小尺度噪声覆盖应力云图也一团乱麻。4. 代码骨架3D拓扑优化主循环、目标灵敏度与OC更新4.1 主程序结构下面是我常用的主循环框架。目标函数是最小化p-范数全局应力体积约束固定优化更新用OC准则。为了让代码清晰我把有限元求解、目标灵敏度计算、过滤器、OC更新都拆成独立函数。function topopt3D_stress(nelx, nely, nelz, volfrac, penal, p, rmin, move) % 3D拓扑优化以p-范数应力和最小为目标的简化实现 Vol ones(nelx*nely*nelz, 1) * volfrac; rho Vol; change 1; hist []; while change 0.01 % 有限元求解 [U, K, edofMat] solve_FE(nelx, nely, nelz, rho, penal); % 目标函数与解析灵敏度 [g, dg] compute_Stress_Obj_Sens(nelx, nely, nelz, ... U, rho, penal, p, edofMat); % 灵敏度过滤 dg filter3D(nelx, nely, nelz, rmin, rho, dg); % OC更新 rho_new OC_update(rho, dg, volfrac, move); change max(abs(rho_new - rho)); rho rho_new; hist(end1, :) [g, change]; fprintf(iter %d, g %.4e, change %.4f\n, ... length(hist), g, change); end end这里我故意没有把solve_FE写全因为标准的3D有限元单元刚度矩阵、边界条件设置等和常规3D拓扑优化完全一样。核心差异集中在目标函数和灵敏度函数里也就是下面一小节。4.2 目标函数与灵敏度函数计算目标和灵敏度的函数是整个程序的关键。我按前面推导的三个部分来写先算每个单元应力并聚合目标再组装伴随右侧向量a接着求解伴随方程最后叠加上直接项和刚度项。function [g, dg] compute_Stress_Obj_Sens(nelx, nely, nelz, ... U, rho, penal, p, edofMat) ndof 3 * (nelx1) * (nely1) * (nelz1); ne nelx * nely * nelz; svm zeros(ne, 1); sig_store zeros(ne, 6); Bstore cell(ne, 1); detJstore zeros(ne, 1); E0 1.0; nu 0.3; D0 isotropicD(E0, nu); % 第一遍计算单元应力、von Mises应力、目标函数值 for e 1:ne node_coord getNodeCoord(e, nelx, nely, nelz); [B, detJ] hex8_B_center(node_coord); D rho(e)^penal * D0; Ue U(edofMat(e, :)); sigma D * B * Ue; svm(e) sqrt(sigma * Vmat * sigma); sig_store(e, :) sigma; Bstore{e} B; detJstore(e) detJ; end g (sum(svm.^p))^(1/p); % 第二遍组装伴随右侧向量a以及直接项 a zeros(ndof, 1); dg_direct zeros(ne, 1); for e 1:ne s svm(e); coef g^(1-p) * s^(p-1); B Bstore{e}; Ue U(edofMat(e, :)); sigma sig_store(e, :); D0B D0 * B; if s 1e-12 dsdU (sigma * Vmat) / s * (rho(e)^penal * D0B); dsdrho (sigma * Vmat) / s * ... (penal * rho(e)^(penal-1) * D0B * Ue); else dsdU zeros(1, 24); dsdrho 0; end a(edofMat(e, :)) a(edofMat(e, :)) coef * dsdU; dg_direct(e) coef * dsdrho; end % 伴随求解 lambda K \ a; % 注意这里的K需要从solve_FE传进来或重复组装 % 第三遍叠加伴随项 dg dg_direct; for e 1:ne B Bstore{e}; k0e B * D0 * B * detJstore(e) * 8; Ue U(edofMat(e, :)); lambda_e lambda(edofMat(e, :)); dg(e) dg(e) - lambda_e * (penal * rho(e)^(penal-1) * k0e * Ue); end end这段代码里K没有直接传入实际使用时要把它从solve_FE里带出来。也可以改成在函数内部重新组装一次K再解伴随方程但那样会浪费一次组装时间。更好的写法是在主循环里把K传到compute_Stress_Obj_Sens里或者像我前面说的把LU分解结果也传进去。4.3 OC更新与过滤函数OC更新在柔顺性最小化里已经很成熟应力目标下也可以用但要注意两个地方一是灵敏度可能出现正值因为增加材料未必在所有位置都降低全局应力二是目标函数非凸OC更新很容易震荡。我一般把move限制在0.1到0.2之间并且加一个对灵敏度的裁剪如果-dg/dv小于等于0就把它设成一个很小的正值避免sqrt里出现负数。function rho_new OC_update(rho, dg, volfrac, move) l1 0; l2 1e9; dv ones(size(rho)); while abs(l2 - l1) 1e-5 lmid 0.5 * (l1 l2); s -dg ./ dv / lmid; s(s 0) 1e-6; rho_new max(0, max(rho - move, min(1, min(rho move, rho .* sqrt(s))))); if sum(rho_new) - volfrac * length(rho) 0 l1 lmid; else l2 lmid; end end end过滤函数本身和柔顺性优化里的标准实现几乎一样只是三维邻居索引需要额外处理z方向。需要提醒的是过滤可以施加在灵敏度上也可以施加在密度场上。如果做密度过滤需要在灵敏度里再乘一个滤波矩阵的转置我在教程里先用灵敏度过滤代码最简单也最容易复现。5. 算例与典型结果从3D悬臂梁看应力优化行为5.1 模型设置与载荷我用一个3D L形支架来做算例。这类结构的内凹角会产生强烈的应力集中是检验应力相关拓扑优化的经典模型。网格取30×20×10共6000个六面体单元左端面固定右上方的局部位置施加向下的载荷体积约束0.3。材料参数用无量纲的E01泊松比0.3。p值取8惩罚系数3过滤半径1.5OC移动步长0.15。在这个规模下每轮迭代需要解两次稀疏线性方程Matlab在我的机器上大概需要两到三分钟完成全过程。如果内存紧张可以把网格缩小到20×15×8结果趋势是一样的。5.2 优化结果解读优化后的拓扑会比较直观内凹角附近会保留下更多材料形成一个近似圆角的过渡区域远离传力路径的角落被挖空结构整体看起来比柔顺性优化更“粗壮”因为应力目标倾向于把材料铺设在受力路径上而不是单纯追求最小柔顺性。p-范数应力值在优化过程中通常先快速下降随后趋于平稳。但因为应力目标不是凸问题收敛曲线会有一定波动。我建议不要只看最后一步要看整个历史的趋势。如果曲线在某一代后反复横跳就减小move值或者改用MMA更新。5.3 收敛曲线与运行时间下面是典型的收敛行为前20代目标值下降幅度最大20到60代缓慢下降60代之后基本稳定。虽然OC更新的数学保证比MMA弱但在这个算例里配合过滤和小步长依然能得到清晰的0/1结构。需要注意的是p越大目标函数越接近“最大应力”曲线抖动越明显我的经验是p8在稳定性和信息量之间比较平衡。这种3D算例最耗时的环节依次是稀疏组装、位移求解、伴随求解。如果用LU分解一次性处理两个方程组总时间能比重复求解少百分之三十左右。代码里一定要避免在循环内反复调用inv或者full否则6000单元的网格就能让Matlab卡到无法忍受。6. 调参与避坑p值、惩罚、棋盘格、奇异解6.1 p值的选择p值直接影响优化结果的性格。p2时所有应力被“平均”看待优化器倾向于把应力分布磨平但最大应力可能并不可控p20时几乎等于在追着最大应力跑更新量很容易剧烈震荡。不同文献的取值范围差别很大但我自己实测下来p在8到12之间最顺手。可以先跑一个p8的版本看趋势再逐步增大。如果目标函数里用了平均形式的p-范数除以单元数那么p值带来的量级变化会小很多收敛判断也更稳定。否则网格加密后同样的p值目标值会大不少容易被误判成没有收敛。6.2 SIMP惩罚与应力松弛应力优化里最微妙的是低密度单元的应力。SIMP模型下单元刚度是ρ^penal应力也是ρ^penal所以低密度单元的内力反馈很弱应力值很低。这看起来是合理的但它让优化器可以把高应力单元“稀释”成低密度从而绕过应力惩罚。真实材料不存在这种取巧行为孔洞区域只要有一点材料局部应力可能非常大。常用的处理办法是应力松弛比如对单元应力再除一个密度相关项σ_physical,e ρ_e^{penal-1} σ_e这样密度下降时应力不会同比例下降低密度区域的应力贡献会被放大从而破坏灰色单元绕过约束的路径。这篇文章的代码里没有加这项是为了先讲清楚伴随方法和p-范数基础实际工程中建议加上。加了之后目标函数和灵敏度的推导需要同步改直接项里会多出(d(ρ^{penal-1})/dρ)的因子。6.3 过滤半径与move限制过滤半径rmin小于1.5个单元尺寸时棋盘格压制不彻底rmin太大结构又会失去细节。3D网格中我习惯用1.5到2.0倍单元尺寸。密度场滤波比灵敏度滤波更平滑但代码实现稍重需要构造线性滤波矩阵建议先用灵敏度滤波入门遇到严重棋盘格再升级。OC更新对move值非常敏感。move0.2在柔顺性优化里很常见但在应力目标下容易震荡我通常从0.15起步如果目标历史平滑就慢慢加到0.2。每轮更新后还要检查体积约束是否满足否则二分会失效。6.4 一个小型的检查清单我每次跑新的3D应力算例前都会按下面这个清单快速过一遍第一确认B矩阵的行列式detJ没有出现负数说明节点顺序错了第二用一个均匀密度场算一次目标函数看是否和手算量级一致第三用有限差分验证随机选出的几个单元灵敏度第四观察第一次迭代的密度更新方向如果高应力区域的材料反而减少说明某个符号搞反了第五检查p-范数目标值是否随网格加密而异常增大。这套检查听起来琐碎但能省下好几天的debug时间。尤其是第一次从2D转到3D时B矩阵的维度、edof索引范围、稀疏组装行数任何一个地方错一点程序都会在某个角落给出离谱的灵敏度而优化器会把错误当成正确方向一路跑偏。最后再分享两个小技巧一是灵敏度验证时扰动密度不要取太小1e-6比较合适因为3D应力目标对舍入误差更敏感二是如果3D网格大到求解吃力先把z方向厚度压到两层退化成一个“准2D”模型来调通所有代码再恢复完整厚度。把这两步做好3D应力敏感度分析和拓扑优化就不会再像一开始那样劝退人了。希望这份记录对你有用。