简介这是一份基于水平集方法的BESO边界元形状优化拓扑优化开源MATLAB代码包面向机械、航空航天等领域的结构优化研究者和工程师适合需要实现结构拓扑优化算法、探究水平集与边界演化方法的用户。压缩包内共包含10个文件全部为.m脚本大小仅5KB核心模块包括主程序main.m、刚度矩阵计算stiffnessMatrix.m、水平集重初始化reinit.m、演化更新evolve.m以及有限元分析FE.m等完整覆盖了从材料信息定义到拓扑迭代的主要流程。已有306人学习下载。通过阅读和运行这些源码可以直观理解水平集函数如何描述结构边界并驱动拓扑演化也可结合SIMP等经典方法对比分析并针对具体工况修改惩罚因子、体积约束或边界移动准则便于科研教学和二次开发。1. 打开 Level Set Method.rar 之前先分清三件事一个以 “Level Set Method.rar” 命名的压缩包里面同时躺着 BESO 和三套 SIMP 相关脚本这在拓扑优化的开源分享里太常见了打包人按自己的目录习惯命名没打算按教科书分类。拓扑优化领域绕不开三条技术路线——SIMP 靠密度惩罚逼近实体分布BESO 靠渐进增删直接生成 0/1 设计水平集方法则把结构边界当作隐式曲面去追踪。把这三件事捋清楚一个包里的代码才能按顺序跑起来而不是换来一个 “BESO 结果怎么全是灰度” 的困惑。这篇文章就拆这三条线概念边界、最小可跑实现、必调参数、复现时五个高频翻车点以及最后怎么把结果拿去做工程输出。2. SIMP 与 BESO 的关系为什么同一个包里同时出现两套密度逻辑2.1 最小柔度问题所有拓扑优化共用的“黑匣子入口”不管用哪种方法连续体拓扑优化的目标几乎都是同一个在给定的材料体积约束下最小化结构柔度。柔度 C 的定义是外力做的功离散到有限元网格上写作C 0.5 * f^T · u 0.5 * u^T · K · u其中 K 是全局刚度矩阵u 是位移解。体积约束写作 V / V0 volfracvolfrac 是允许的材料体积占比。设计变量在 SIMP 里是每个单元的密度 rho0 到 1在 BESO 里是每个单元的存在状态 x0 或 1在水平集方法里则是一个隐式函数的取值。目标函数本身不是难点难点全在灵敏度每个单元的材料增减会让目标函数变化多少。这个梯度信息是所有方法的“黑匣子入口”因为后续的删材料、填材料、推进边界全靠它决定方向。灵敏度算错了后面所有结果都是错的而且错得很隐蔽。2.2 SIMP密度惩罚如何把中间密度推向 0/1SIMPSolid Isotropic Material with Penalization的基本假设是单元弹性模量与密度之间满足幂函数关系E(rho) E_min rho^p * (E0 - E_min)当惩罚因子 p 3 时rho 0.5 的单元实际模量只有实体单元的 12.5%中间密度的“性价比”被压得很低优化器会倾向于把单元推向 0 或 1。这就是 SIMP 名字里 Penalization 的含义。但 SIMP 的惩罚再强最终结果里仍然会残留大量 0.3 到 0.7 的灰色单元尤其是载荷边界附近。这类灰色区域没法直接进 CAM 或切片软件必须做后处理二值化。我一般会把 SIMP 看作“先找拓扑骨架”的工具它的优点是数学形式干净、梯度信息完整、收敛行为相对平稳缺点是灰度多、边界毛糙、制造前必须再做一次阈值切分。2.3 BESO双向渐进增删如何在 SIMP 之后收尸BESOBi-directional Evolutionary Structural Optimization的思路完全不同设计变量就是 0/1每个迭代步直接删掉一批灵敏度最低的实体单元同时激活一批灵敏度最高的空单元。这个“双向”是核心老版 ESO 只删不加常常把一些本应该重新出现的材料路径永久删掉最后得到次优甚至断裂的结构。BESO 的典型迭代流程是先做一次有限元分析计算每个单元的应变能作为灵敏度再对灵敏度做空间滤波防止棋盘格然后根据当前体积约束把低灵敏度的实体单元置 0把高灵敏度的空单元置 1。灵敏度滤波是整个 BESO 的“后悔药”没有它单像素的棋盘格会让结果完全不可用。在开源实现里BESO 的有限元求解骨架几乎都源自 SIMP 那套经典代码区别只在于设计变量更新那一段SIMP 用梯度投影BESO 用排序硬阈值。所以一个压缩包里同时出现 SIMP 和 BESO不是重复而是前者负责灵敏度分析后者负责 0/1 化。2.4 三条路线的选型对比方法设计变量边界表示灰度问题后处理成本计算开销SIMPrho 连续密度灰度斜面严重需要阈值二值化低成熟求解器可复用BESO0/1 离散状态阶梯状边界无需要平滑边界低但每代要重新排序Level Set隐式函数 phi零等值面无边界可直接提取高需处理重新初始化2.5 选型建议什么场景用哪种做概念设计、想快速验证载荷路径用 SIMP因为实现简单、参数少做制造约束研究或者增材制造准备用 BESO因为 0/1 结果离 STL 更近做边界质量要求高的设计比如气动外形、流道用 Level Set 或 BESO Level Set 组合。实际开源项目里最常见的组合是BESO 跑出 0/1 骨架再把骨架转成符号距离函数交给 Level Set 做几十轮边界光顺。这个组合我在后面会给出具体做法。3. 从 0 写一个 BESO 拓扑优化二维骨架代码与四个必调参数3.1 先定问题悬臂梁的边界条件与设计域复现 BESO我建议不要一上来就选 MBB 梁。MBB 梁的滚动支座和对称边界条件容易搞错悬臂梁的边界更直观左侧全固定右下角节点受一个向下的集中力。设计域取 60 x 20 个正方形平面应力单元目标体积占比 0.5进化率 ER 0.02灵敏度滤波半径 rmin 2.4单位单元尺寸。3.2 最小可跑实现Python NumPy SciPy完整骨架代码如下这个版本我把单元刚度矩阵用 2x2 高斯积分现算不依赖任何记忆中的 8x8 常量矩阵避免抄错。import numpy as np import scipy.sparse as sp import scipy.sparse.linalg as spla # ---------- 参数 ---------- nelx, nely 60, 20 # 网格尺寸 volfrac 0.5 # 目标体积占比 er 0.02 # 进化率每次最多改变的材料比例 rmin 2.4 # 灵敏度滤波半径单元尺寸 E0, Emin, nu 1.0, 1e-6, 0.3 n_el nelx * nely ndof 2 * (nelx 1) * (nely 1) # ---------- 单元与自由度映射 ---------- def node_id(ex, ey): return ex * (nely 1) ey edof np.zeros((n_el, 8), dtypeint) for ey in range(nely): for ex in range(nelx): n1 node_id(ex, ey) # 左下 n2 node_id(ex 1, ey) # 右下 n3 node_id(ex 1, ey 1) # 右上 n4 node_id(ex, ey 1) # 左上 e ey * nelx ex edof[e] [2*n1, 2*n11, 2*n2, 2*n21, 2*n3, 2*n31, 2*n4, 2*n41] # ---------- 平面应力 Q4 单元刚度矩阵 ---------- def unit_K(nu): D E0 / (1 - nu**2) * np.array([ [1, nu, 0], [nu, 1, 0], [0, 0, (1 - nu) / 2]]) k np.zeros((8, 8)) gp np.sqrt(1 / 3) for gx in [-gp, gp]: for gy in [-gp, gp]: # 形函数对局部坐标的导数单位正方形 [-1,1]^2 dN_local 0.25 * np.array([ [-(1 - gy), (1 - gy), (1 gy), -(1 gy)], [-(1 - gx), -(1 gx), (1 gx), (1 - gx)]]) # 物理空间导数单元边长 1故乘 2 dN_phys 2.0 * dN_local B np.zeros((3, 8)) B[0, 0::2] dN_phys[0] B[1, 1::2] dN_phys[1] B[2, 0::2] dN_phys[1] B[2, 1::2] dN_phys[0] k B.T D B * 0.25 # detJ 0.25 return k KE0 unit_K(nu) # ---------- 全局刚度矩阵的自由度索引 ---------- I np.repeat(edof, 8, axis0).flatten() # 行索引 J np.tile(edof, (1, 8)).flatten() # 列索引 # ---------- 载荷与约束 ---------- # 右下角节点受向下单位力 load_node node_id(nelx, 0) F np.zeros(ndof) F[2 * load_node 1] -1.0 # 左侧全部固定 fixed [] for ey in range(nely 1): n node_id(0, ey) fixed [2 * n, 2 * n 1] fixed np.array(fixed) # ---------- 灵敏度滤波 ---------- def sens_filter(sens, rmin): sf np.zeros_like(sens) r int(np.ceil(rmin)) for ey in range(nely): for ex in range(nelx): e ey * nelx ex sw, sv 0.0, 0.0 for dx in range(-r, r 1): for dy in range(-r, r 1): nx, ny ex dx, ey dy if 0 nx nelx and 0 ny nely: w max(0.0, rmin - np.sqrt(dx**2 dy**2)) sw w sv w * sens[ny * nelx nx] sf[e] sv / max(sw, 1e-12) return sf # ---------- 主循环 ---------- x np.ones(n_el) # 设计变量1 实体0 空 vol_cur 1.0 c_hist [] for it in range(200): # 组装当前刚度矩阵BESO 中 x 为 0/1乘 x 即 0/1 模量 E_full Emin x**3 * (E0 - Emin) K sp.coo_matrix( (np.repeat(E_full, 64) * np.tile(KE0.flatten(), n_el), (I, J)), shape(ndof, ndof) ).tocsc() # 施加边界条件 K K.tolil() K[fixed, :] 0 K[:, fixed] 0 K[fixed, fixed] 1 K K.tocsr() F_bc F.copy() F_bc[fixed] 0 U spla.spsolve(K, F_bc) # 单元应变能作为灵敏度 u_e U[edof] # (n_el, 8) sens 0.5 * np.sum(u_e * (u_e KE0), axis1) sens sens_filter(sens, rmin) # 目标函数 c_hist.append(0.5 * F_bc U) # ---------- 设计变量更新 ---------- vol_tar max(vol_cur * (1 - er), volfrac) if vol_cur vol_tar: # 阶段一只删按最低灵敏度删到本代目标体积 n_del int(vol_cur * n_el) - int(vol_tar * n_el) n_del max(n_del, 1) solid np.where(x 1)[0] order np.argsort(sens[solid]) x[solid[order[:n_del]]] 0 else: # 阶段二双向增删保持体积不变 n_move max(int(er * n_el), 1) solid np.where(x 1)[0] void np.where(x 0)[0] del_th np.sort(sens[solid])[min(n_move, len(solid) - 1)] add_th np.sort(sens[void])[max(-n_move, -len(void))] x[solid[sens[solid] del_th]] 0 x[void[sens[void] add_th]] 1 vol_cur x.mean() # ---------- 收敛判断 ---------- if len(c_hist) 10: d abs(c_hist[-1] - c_hist[-10]) / max(abs(c_hist[-10]), 1e-12) if abs(vol_cur - volfrac) 0.01 and d 0.005: print(fconverged at iteration {it}) break3.3 代码逻辑说明单元刚度矩阵用 2x2 高斯积分生成原因是四边形等参元的刚度矩阵没有简单一致的常数形式手工抄一个 8x8 矩阵容易错写成公式反而可检查。全局刚度矩阵用 COO 格式按自由度索引组装I和J两个数组把每个单元的 8x8 刚度块映射到全局这是稀疏矩阵组装的标准做法。固定自由度的处理是直接置零、对角线置 1对二维小规模问题足够快。灵敏度取单元应变能物理含义是“这个单元对结构刚度贡献了多少”。删除低应变能的实体、添加高应变能的空单元就是在保留承载路径的同时填掉冗余区域。滤波半径 rmin 直接影响棋盘格。半径小于 1.5 倍单元尺寸时单像素的交替结构无法被抑制大于 3 时拓扑会丢失细长杆件。2.4 是一个在 60x20 网格上表现可靠的经验值。进化率 er 控制每代最多改变多少材料占比。er 太大比如 0.1体积下降太快灵敏度的排序信息还没反映到新拓扑上结果会是断裂的碎杆er 太小比如 0.005收敛速度慢对三维网格尤其不划算。3.4 从二维到三维要注意的差异二维代码推到三维三个地方会直接卡住一是稀疏矩阵组装不能再依赖柯西稠密索引建议直接用scipy.sparse的csr_array构建二是滤波函数如果还用双层循环几十万单元会让每次迭代慢到不可接受要改用scipy.ndimage.convolve或skimage做平滑三是求解器从spsolve换成cg或minres配合 ILU 预条件。三维 BESO 的内存瓶颈通常在滤波而不是有限元求解。4. Level Set Method 的工程化从偏微分方程到一段能跑的更新骨架4.1 符号距离函数把“边界”变成等值面水平集方法的核心是把结构表示成一个标量函数 phi(x)约定 phi 0 表示实体phi 0 表示孔洞phi 0 的等值面就是结构边界。phi 通常被初始化为符号距离函数也就是点到边界的最短带符号距离。这个表示方式最大的好处是拓扑变化不需要人为处理边界交叉、合并、断裂都只是等值面的正常演变。而 BESO 和 SIMP 每步都要显式判断哪些单元被删、哪些单元被填。边界还能直接提供法向量和曲率法向量是 nabla(phi) / |nabla(phi)|曲率是散度。这两个几何量在做应力约束或尺寸控制时非常有用这是密度方法很难直接给的。4.2 Hamilton-Jacobi 演化与速度场来源水平集的演化方程是经典的 Hamilton-Jacobi 方程d(phi)/dt Vn * |nabla(phi)| 0Vn 是边界沿法向的运动速度。在拓扑优化里Vn 从目标函数对 phi 的变分导数推导出来。直觉上边界向内收缩会减少材料体积向外扩张会增加材料体积所以 Vn 和柔度灵敏度直接相关。这里有一个工程上必须解决的问题phi 经过多轮演化后会变得非常陡峭或非常平缓不再保持符号距离函数性质导致 |nabla(phi)| 失真边界运动速度被算错。常见的补救措施是定期做重新初始化reinitialization最常见的是逐点把 phi 拉回“距离边界多远”的语义。这个操作很吃数值功底也是很多开源代码里最难看懂的部分。4.3 参数化水平集把偏微分方程变成有限维优化纯水平集方法要解 Hamilton-Jacobi 方程需要 CFL 条件控制时间步长还要周期性重新初始化对结构工程师并不友好。一种更工程化的替代是参数化水平集方法把 phi 展开成一组径向基函数的线性组合phi(x) sum(alpha_i * g_i(x))其中 g_i 取高斯函数alpha_i 是需要优化的系数。这样原来的偏微分方程演化变成了对有限个系数 alpha 的梯度下降可以直接用和 SIMP 相同的灵敏度分析框架。一个可运行的更新骨架如下# 控制点与网格点设计域 60x20 单元物理尺寸 1x1 gx np.linspace(0, 1, 30) gy np.linspace(0, 1, 10) ctrl np.array(np.meshgrid(gx, gy)).T.reshape(-1, 2) grid np.array( np.meshgrid(np.linspace(0, 1, nelx), np.linspace(0, 1, nely)) ).T.reshape(-1, 2) eps 0.12 # 高斯基函数的支撑半径 def rbf_basis(ctrl, grid, eps): diff grid[:, None, :] - ctrl[None, :, :] return np.exp(-np.sum(diff**2, axis-1) / eps**2) G rbf_basis(ctrl, grid, eps) # (n_grid, n_ctrl) alpha np.ones(ctrl.shape[0]) * 0.3 # 初始正数全部实体 def pls_step(alpha, grad_c, G, dt, reg1e-4): # grad_c: 目标函数对每个网格点 phi 的梯度 dalpha -dt * (G.T grad_c / G.shape[0] reg * G.T (G alpha)) return alpha dalpha # 每次迭代后phi G alpha # phi 0 的网格为实体phi 等值面为边界这个骨架里最关键的是支撑半径 eps。eps 太小每个基函数只能影响周围几个单元优化结果退化成逐点噪声eps 太大边界表达能力下降细节特征丢失。对 60x20 的网格eps 在 0.1 到 0.15 之间比较合理。参数化水平集的代价是边界表达能力受限于控制点密度。控制点 30x10边界最多能表达 30x10 分辨率下的拓扑细节这比 BESO 的逐单元 0/1 要粗。所以更常见的做法不是让水平集从头优化而是先跑 BESO 拿到骨架再转成 phi 做几十轮边界光顺。4.4 BESO 与 Level Set 的组合先骨架后光顺我在开源项目里见过的比较可靠的落地流水线是两阶段式。第一阶段跑 BESO得到 0/1 单元分布第二阶段把实体单元转成 phi 的初始符号距离函数然后让水平集在固定拓扑或弱拓扑变化下做边界光顺把阶梯状的单元边界变成平滑的等值面。这样组合的好处是BESO 负责全局搜索水平集负责边界质量。BESO 单跑的结果放进流体仿真或应力分析里阶梯边界会造成局部应力集中平滑后数值结果更接近真实。5. 拓扑优化复现避坑五个让结果翻车的参数陷阱5.1 棋盘格花纹灵敏度滤波半径没到位现象结果里实体单元和空单元呈棋盘状交错黑白相间像象棋盘一样。目标函数还很漂亮地下降但实际上完全没有工程意义。原因单元灵敏度和单元本身一一对应没有空间连续性。相邻单元的灵敏度差异很大时优化器会利用这种不连续形成单像素交替的高刚度路径。解决把灵敏度滤波半径 rmin 至少调到单元尺寸的 1.5 倍以上。在我这个骨架代码里 rmin 2.4 是在 60x20 网格上验证过的经验值。注意滤波时权重必须除以权重和否则灵敏度的量纲会漂移收敛后的拓扑偏刚或偏柔。5.2 结果里灰度一堆BESO 跑成了 SIMP现象x 数组里 0.3 到 0.7 的中间值占了一大半设计变量根本没有收敛到 0/1。原因最常见是把 SIMP 的密度投影逻辑直接搬进 BESO而 BESO 的设计变量更新逻辑是“排序后硬阈值”不是“对 rho 做梯度下降”。另一类原因是惩罚因子 p 设置过小比如 p 1中间密度没有惩罚SIMP 本身就收不到 0 和 1。解决BESO 的单元模量直接用 E Emin x * (E0 - Emin)不要用 x 的幂次每次更新后可以对 x 做一个硬阈值投影x 0.5 置 0x 0.5 置 1。这个二值化会让目标函数曲线出现小的跳变属于正常现象不用担心。5.3 加载点周围出现蘑菇头状实心块现象集中力作用的节点附近糊了一大块实体像个蘑菇头而真正有意义的承载杆件在别处被删得七零八落。原因集中力在节点处产生应力奇异性局部应变能极高灵敏度排序永远把这块区域排在“不能删”的最前面。优化器宁可牺牲远端结构也要保住载荷点的局部材料。解决把载荷作用点周边 2 到 3 层单元设为非设计域这部分单元固定为实体不参与增删。另一个做法是把单个节点力分散到相邻 3 个节点上降低奇异性。这个坑在三维模型里更严重因为网格更密局部奇异的影响范围更集中。5.4 目标函数一直下降但拓扑每代都在变现象柔度曲线平滑下降体积也已到达目标值可是每次迭代输出的拓扑图案完全不一样有时候右边多一条杆下一轮又消失了。原因收敛判据只看了目标函数变化没看拓扑变化。柔度在最优解附近本来就平缓可能只差 0.1%拓扑却还在大范围调整。还有一个隐藏原因双向增删阶段的阈值取单代灵敏度排序代与代之间的噪声会让边界反复横跳。解决收敛条件改成两条同时满足——体积误差小于 1%且连续 10 代目标函数相对变化小于 0.5%。如果还震荡对灵敏度做历史加权平均把上一代的灵敏度按 0.5 权重混合进来。这个技巧在论文里叫 averaging ratio工程上直接好用。5.5 从二维扩到三维内存先炸求解器再卡住现象把 60x20 换成 60x40x20 的三维网格程序在组装刚度矩阵时内存直接耗尽或者用稠密求解器几分钟算不出来。原因三维网格的自由度数是二维的几十倍以上如果全局刚度矩阵还按稠密数组存内存直接爆炸滤波函数还是双层循环每步都要做几百万次距离计算慢是必然的。解决三维代码必须满足三个条件刚度矩阵用 scipy.sparse 的 CSR/CSC 格式滤波改用 ndimage 的均匀滤波或高斯滤波替代逐单元循环求解器从 spsolve 换成迭代法比如共轭梯度配 ILU 预条件。我踩过最深的一次翻车就是把二维滤波直接平移结果一次迭代跑了二十分钟。6. 把 BESO 结果变成能出图的边界等值线提取与一个小习惯BESO 输出的 x 是 0/1 矩阵画成像素图能用但拿去做 CAD 和有限元复核就得先提边界。最常见做法是用 marching squares 提取 0.5 等值线一个函数就能完成from skimage import measure import numpy as np x2d x.reshape(nely, nelx).astype(float) contours measure.find_contours(x2d, 0.5) # contours 是多段折线直接导出 DXF 前先用高斯平滑去掉单像素毛刺 from scipy.ndimage import gaussian_filter x_smooth gaussian_filter(x2d, sigma0.8) contours_smooth measure.find_contours(x_smooth, 0.5)一个有价值的验证技巧是把 BESO 的拓扑结果做一次镜像再跑一遍有限元对比镜像前后的柔度。偏差小于 1% 说明边界条件和网格映射处理得没有方向性错误偏差大先检查载荷和约束的自由度索引。对称约束是个更强的手段把设计域设为全宽、载荷对称优化结果若不关于中心线对称说明灵敏度滤波的窗口扫描方向或网格编号有 bug。这个检查比盯着彩色云图可靠得多。我自己的一个习惯是任何拓扑优化结果先提取 0.5 等值线看一眼再回去看彩色灵敏度图。只看灵敏度云图很容易被灰度骗到以为结构快要收敛了。边界折线一拉出来单像素断裂、悬浮材料、锐角尖点全都藏不住。这个习惯帮我省了不少后悔药。希望帮到你。本文还有配套的精品资源点击获取 拓扑优化三大主流方法：SIMP、BESO与水平集工程实践解析
简介这是一份基于水平集方法的BESO边界元形状优化拓扑优化开源MATLAB代码包面向机械、航空航天等领域的结构优化研究者和工程师适合需要实现结构拓扑优化算法、探究水平集与边界演化方法的用户。压缩包内共包含10个文件全部为.m脚本大小仅5KB核心模块包括主程序main.m、刚度矩阵计算stiffnessMatrix.m、水平集重初始化reinit.m、演化更新evolve.m以及有限元分析FE.m等完整覆盖了从材料信息定义到拓扑迭代的主要流程。已有306人学习下载。通过阅读和运行这些源码可以直观理解水平集函数如何描述结构边界并驱动拓扑演化也可结合SIMP等经典方法对比分析并针对具体工况修改惩罚因子、体积约束或边界移动准则便于科研教学和二次开发。1. 打开 Level Set Method.rar 之前先分清三件事一个以 “Level Set Method.rar” 命名的压缩包里面同时躺着 BESO 和三套 SIMP 相关脚本这在拓扑优化的开源分享里太常见了打包人按自己的目录习惯命名没打算按教科书分类。拓扑优化领域绕不开三条技术路线——SIMP 靠密度惩罚逼近实体分布BESO 靠渐进增删直接生成 0/1 设计水平集方法则把结构边界当作隐式曲面去追踪。把这三件事捋清楚一个包里的代码才能按顺序跑起来而不是换来一个 “BESO 结果怎么全是灰度” 的困惑。这篇文章就拆这三条线概念边界、最小可跑实现、必调参数、复现时五个高频翻车点以及最后怎么把结果拿去做工程输出。2. SIMP 与 BESO 的关系为什么同一个包里同时出现两套密度逻辑2.1 最小柔度问题所有拓扑优化共用的“黑匣子入口”不管用哪种方法连续体拓扑优化的目标几乎都是同一个在给定的材料体积约束下最小化结构柔度。柔度 C 的定义是外力做的功离散到有限元网格上写作C 0.5 * f^T · u 0.5 * u^T · K · u其中 K 是全局刚度矩阵u 是位移解。体积约束写作 V / V0 volfracvolfrac 是允许的材料体积占比。设计变量在 SIMP 里是每个单元的密度 rho0 到 1在 BESO 里是每个单元的存在状态 x0 或 1在水平集方法里则是一个隐式函数的取值。目标函数本身不是难点难点全在灵敏度每个单元的材料增减会让目标函数变化多少。这个梯度信息是所有方法的“黑匣子入口”因为后续的删材料、填材料、推进边界全靠它决定方向。灵敏度算错了后面所有结果都是错的而且错得很隐蔽。2.2 SIMP密度惩罚如何把中间密度推向 0/1SIMPSolid Isotropic Material with Penalization的基本假设是单元弹性模量与密度之间满足幂函数关系E(rho) E_min rho^p * (E0 - E_min)当惩罚因子 p 3 时rho 0.5 的单元实际模量只有实体单元的 12.5%中间密度的“性价比”被压得很低优化器会倾向于把单元推向 0 或 1。这就是 SIMP 名字里 Penalization 的含义。但 SIMP 的惩罚再强最终结果里仍然会残留大量 0.3 到 0.7 的灰色单元尤其是载荷边界附近。这类灰色区域没法直接进 CAM 或切片软件必须做后处理二值化。我一般会把 SIMP 看作“先找拓扑骨架”的工具它的优点是数学形式干净、梯度信息完整、收敛行为相对平稳缺点是灰度多、边界毛糙、制造前必须再做一次阈值切分。2.3 BESO双向渐进增删如何在 SIMP 之后收尸BESOBi-directional Evolutionary Structural Optimization的思路完全不同设计变量就是 0/1每个迭代步直接删掉一批灵敏度最低的实体单元同时激活一批灵敏度最高的空单元。这个“双向”是核心老版 ESO 只删不加常常把一些本应该重新出现的材料路径永久删掉最后得到次优甚至断裂的结构。BESO 的典型迭代流程是先做一次有限元分析计算每个单元的应变能作为灵敏度再对灵敏度做空间滤波防止棋盘格然后根据当前体积约束把低灵敏度的实体单元置 0把高灵敏度的空单元置 1。灵敏度滤波是整个 BESO 的“后悔药”没有它单像素的棋盘格会让结果完全不可用。在开源实现里BESO 的有限元求解骨架几乎都源自 SIMP 那套经典代码区别只在于设计变量更新那一段SIMP 用梯度投影BESO 用排序硬阈值。所以一个压缩包里同时出现 SIMP 和 BESO不是重复而是前者负责灵敏度分析后者负责 0/1 化。2.4 三条路线的选型对比方法设计变量边界表示灰度问题后处理成本计算开销SIMPrho 连续密度灰度斜面严重需要阈值二值化低成熟求解器可复用BESO0/1 离散状态阶梯状边界无需要平滑边界低但每代要重新排序Level Set隐式函数 phi零等值面无边界可直接提取高需处理重新初始化2.5 选型建议什么场景用哪种做概念设计、想快速验证载荷路径用 SIMP因为实现简单、参数少做制造约束研究或者增材制造准备用 BESO因为 0/1 结果离 STL 更近做边界质量要求高的设计比如气动外形、流道用 Level Set 或 BESO Level Set 组合。实际开源项目里最常见的组合是BESO 跑出 0/1 骨架再把骨架转成符号距离函数交给 Level Set 做几十轮边界光顺。这个组合我在后面会给出具体做法。3. 从 0 写一个 BESO 拓扑优化二维骨架代码与四个必调参数3.1 先定问题悬臂梁的边界条件与设计域复现 BESO我建议不要一上来就选 MBB 梁。MBB 梁的滚动支座和对称边界条件容易搞错悬臂梁的边界更直观左侧全固定右下角节点受一个向下的集中力。设计域取 60 x 20 个正方形平面应力单元目标体积占比 0.5进化率 ER 0.02灵敏度滤波半径 rmin 2.4单位单元尺寸。3.2 最小可跑实现Python NumPy SciPy完整骨架代码如下这个版本我把单元刚度矩阵用 2x2 高斯积分现算不依赖任何记忆中的 8x8 常量矩阵避免抄错。import numpy as np import scipy.sparse as sp import scipy.sparse.linalg as spla # ---------- 参数 ---------- nelx, nely 60, 20 # 网格尺寸 volfrac 0.5 # 目标体积占比 er 0.02 # 进化率每次最多改变的材料比例 rmin 2.4 # 灵敏度滤波半径单元尺寸 E0, Emin, nu 1.0, 1e-6, 0.3 n_el nelx * nely ndof 2 * (nelx 1) * (nely 1) # ---------- 单元与自由度映射 ---------- def node_id(ex, ey): return ex * (nely 1) ey edof np.zeros((n_el, 8), dtypeint) for ey in range(nely): for ex in range(nelx): n1 node_id(ex, ey) # 左下 n2 node_id(ex 1, ey) # 右下 n3 node_id(ex 1, ey 1) # 右上 n4 node_id(ex, ey 1) # 左上 e ey * nelx ex edof[e] [2*n1, 2*n11, 2*n2, 2*n21, 2*n3, 2*n31, 2*n4, 2*n41] # ---------- 平面应力 Q4 单元刚度矩阵 ---------- def unit_K(nu): D E0 / (1 - nu**2) * np.array([ [1, nu, 0], [nu, 1, 0], [0, 0, (1 - nu) / 2]]) k np.zeros((8, 8)) gp np.sqrt(1 / 3) for gx in [-gp, gp]: for gy in [-gp, gp]: # 形函数对局部坐标的导数单位正方形 [-1,1]^2 dN_local 0.25 * np.array([ [-(1 - gy), (1 - gy), (1 gy), -(1 gy)], [-(1 - gx), -(1 gx), (1 gx), (1 - gx)]]) # 物理空间导数单元边长 1故乘 2 dN_phys 2.0 * dN_local B np.zeros((3, 8)) B[0, 0::2] dN_phys[0] B[1, 1::2] dN_phys[1] B[2, 0::2] dN_phys[1] B[2, 1::2] dN_phys[0] k B.T D B * 0.25 # detJ 0.25 return k KE0 unit_K(nu) # ---------- 全局刚度矩阵的自由度索引 ---------- I np.repeat(edof, 8, axis0).flatten() # 行索引 J np.tile(edof, (1, 8)).flatten() # 列索引 # ---------- 载荷与约束 ---------- # 右下角节点受向下单位力 load_node node_id(nelx, 0) F np.zeros(ndof) F[2 * load_node 1] -1.0 # 左侧全部固定 fixed [] for ey in range(nely 1): n node_id(0, ey) fixed [2 * n, 2 * n 1] fixed np.array(fixed) # ---------- 灵敏度滤波 ---------- def sens_filter(sens, rmin): sf np.zeros_like(sens) r int(np.ceil(rmin)) for ey in range(nely): for ex in range(nelx): e ey * nelx ex sw, sv 0.0, 0.0 for dx in range(-r, r 1): for dy in range(-r, r 1): nx, ny ex dx, ey dy if 0 nx nelx and 0 ny nely: w max(0.0, rmin - np.sqrt(dx**2 dy**2)) sw w sv w * sens[ny * nelx nx] sf[e] sv / max(sw, 1e-12) return sf # ---------- 主循环 ---------- x np.ones(n_el) # 设计变量1 实体0 空 vol_cur 1.0 c_hist [] for it in range(200): # 组装当前刚度矩阵BESO 中 x 为 0/1乘 x 即 0/1 模量 E_full Emin x**3 * (E0 - Emin) K sp.coo_matrix( (np.repeat(E_full, 64) * np.tile(KE0.flatten(), n_el), (I, J)), shape(ndof, ndof) ).tocsc() # 施加边界条件 K K.tolil() K[fixed, :] 0 K[:, fixed] 0 K[fixed, fixed] 1 K K.tocsr() F_bc F.copy() F_bc[fixed] 0 U spla.spsolve(K, F_bc) # 单元应变能作为灵敏度 u_e U[edof] # (n_el, 8) sens 0.5 * np.sum(u_e * (u_e KE0), axis1) sens sens_filter(sens, rmin) # 目标函数 c_hist.append(0.5 * F_bc U) # ---------- 设计变量更新 ---------- vol_tar max(vol_cur * (1 - er), volfrac) if vol_cur vol_tar: # 阶段一只删按最低灵敏度删到本代目标体积 n_del int(vol_cur * n_el) - int(vol_tar * n_el) n_del max(n_del, 1) solid np.where(x 1)[0] order np.argsort(sens[solid]) x[solid[order[:n_del]]] 0 else: # 阶段二双向增删保持体积不变 n_move max(int(er * n_el), 1) solid np.where(x 1)[0] void np.where(x 0)[0] del_th np.sort(sens[solid])[min(n_move, len(solid) - 1)] add_th np.sort(sens[void])[max(-n_move, -len(void))] x[solid[sens[solid] del_th]] 0 x[void[sens[void] add_th]] 1 vol_cur x.mean() # ---------- 收敛判断 ---------- if len(c_hist) 10: d abs(c_hist[-1] - c_hist[-10]) / max(abs(c_hist[-10]), 1e-12) if abs(vol_cur - volfrac) 0.01 and d 0.005: print(fconverged at iteration {it}) break3.3 代码逻辑说明单元刚度矩阵用 2x2 高斯积分生成原因是四边形等参元的刚度矩阵没有简单一致的常数形式手工抄一个 8x8 矩阵容易错写成公式反而可检查。全局刚度矩阵用 COO 格式按自由度索引组装I和J两个数组把每个单元的 8x8 刚度块映射到全局这是稀疏矩阵组装的标准做法。固定自由度的处理是直接置零、对角线置 1对二维小规模问题足够快。灵敏度取单元应变能物理含义是“这个单元对结构刚度贡献了多少”。删除低应变能的实体、添加高应变能的空单元就是在保留承载路径的同时填掉冗余区域。滤波半径 rmin 直接影响棋盘格。半径小于 1.5 倍单元尺寸时单像素的交替结构无法被抑制大于 3 时拓扑会丢失细长杆件。2.4 是一个在 60x20 网格上表现可靠的经验值。进化率 er 控制每代最多改变多少材料占比。er 太大比如 0.1体积下降太快灵敏度的排序信息还没反映到新拓扑上结果会是断裂的碎杆er 太小比如 0.005收敛速度慢对三维网格尤其不划算。3.4 从二维到三维要注意的差异二维代码推到三维三个地方会直接卡住一是稀疏矩阵组装不能再依赖柯西稠密索引建议直接用scipy.sparse的csr_array构建二是滤波函数如果还用双层循环几十万单元会让每次迭代慢到不可接受要改用scipy.ndimage.convolve或skimage做平滑三是求解器从spsolve换成cg或minres配合 ILU 预条件。三维 BESO 的内存瓶颈通常在滤波而不是有限元求解。4. Level Set Method 的工程化从偏微分方程到一段能跑的更新骨架4.1 符号距离函数把“边界”变成等值面水平集方法的核心是把结构表示成一个标量函数 phi(x)约定 phi 0 表示实体phi 0 表示孔洞phi 0 的等值面就是结构边界。phi 通常被初始化为符号距离函数也就是点到边界的最短带符号距离。这个表示方式最大的好处是拓扑变化不需要人为处理边界交叉、合并、断裂都只是等值面的正常演变。而 BESO 和 SIMP 每步都要显式判断哪些单元被删、哪些单元被填。边界还能直接提供法向量和曲率法向量是 nabla(phi) / |nabla(phi)|曲率是散度。这两个几何量在做应力约束或尺寸控制时非常有用这是密度方法很难直接给的。4.2 Hamilton-Jacobi 演化与速度场来源水平集的演化方程是经典的 Hamilton-Jacobi 方程d(phi)/dt Vn * |nabla(phi)| 0Vn 是边界沿法向的运动速度。在拓扑优化里Vn 从目标函数对 phi 的变分导数推导出来。直觉上边界向内收缩会减少材料体积向外扩张会增加材料体积所以 Vn 和柔度灵敏度直接相关。这里有一个工程上必须解决的问题phi 经过多轮演化后会变得非常陡峭或非常平缓不再保持符号距离函数性质导致 |nabla(phi)| 失真边界运动速度被算错。常见的补救措施是定期做重新初始化reinitialization最常见的是逐点把 phi 拉回“距离边界多远”的语义。这个操作很吃数值功底也是很多开源代码里最难看懂的部分。4.3 参数化水平集把偏微分方程变成有限维优化纯水平集方法要解 Hamilton-Jacobi 方程需要 CFL 条件控制时间步长还要周期性重新初始化对结构工程师并不友好。一种更工程化的替代是参数化水平集方法把 phi 展开成一组径向基函数的线性组合phi(x) sum(alpha_i * g_i(x))其中 g_i 取高斯函数alpha_i 是需要优化的系数。这样原来的偏微分方程演化变成了对有限个系数 alpha 的梯度下降可以直接用和 SIMP 相同的灵敏度分析框架。一个可运行的更新骨架如下# 控制点与网格点设计域 60x20 单元物理尺寸 1x1 gx np.linspace(0, 1, 30) gy np.linspace(0, 1, 10) ctrl np.array(np.meshgrid(gx, gy)).T.reshape(-1, 2) grid np.array( np.meshgrid(np.linspace(0, 1, nelx), np.linspace(0, 1, nely)) ).T.reshape(-1, 2) eps 0.12 # 高斯基函数的支撑半径 def rbf_basis(ctrl, grid, eps): diff grid[:, None, :] - ctrl[None, :, :] return np.exp(-np.sum(diff**2, axis-1) / eps**2) G rbf_basis(ctrl, grid, eps) # (n_grid, n_ctrl) alpha np.ones(ctrl.shape[0]) * 0.3 # 初始正数全部实体 def pls_step(alpha, grad_c, G, dt, reg1e-4): # grad_c: 目标函数对每个网格点 phi 的梯度 dalpha -dt * (G.T grad_c / G.shape[0] reg * G.T (G alpha)) return alpha dalpha # 每次迭代后phi G alpha # phi 0 的网格为实体phi 等值面为边界这个骨架里最关键的是支撑半径 eps。eps 太小每个基函数只能影响周围几个单元优化结果退化成逐点噪声eps 太大边界表达能力下降细节特征丢失。对 60x20 的网格eps 在 0.1 到 0.15 之间比较合理。参数化水平集的代价是边界表达能力受限于控制点密度。控制点 30x10边界最多能表达 30x10 分辨率下的拓扑细节这比 BESO 的逐单元 0/1 要粗。所以更常见的做法不是让水平集从头优化而是先跑 BESO 拿到骨架再转成 phi 做几十轮边界光顺。4.4 BESO 与 Level Set 的组合先骨架后光顺我在开源项目里见过的比较可靠的落地流水线是两阶段式。第一阶段跑 BESO得到 0/1 单元分布第二阶段把实体单元转成 phi 的初始符号距离函数然后让水平集在固定拓扑或弱拓扑变化下做边界光顺把阶梯状的单元边界变成平滑的等值面。这样组合的好处是BESO 负责全局搜索水平集负责边界质量。BESO 单跑的结果放进流体仿真或应力分析里阶梯边界会造成局部应力集中平滑后数值结果更接近真实。5. 拓扑优化复现避坑五个让结果翻车的参数陷阱5.1 棋盘格花纹灵敏度滤波半径没到位现象结果里实体单元和空单元呈棋盘状交错黑白相间像象棋盘一样。目标函数还很漂亮地下降但实际上完全没有工程意义。原因单元灵敏度和单元本身一一对应没有空间连续性。相邻单元的灵敏度差异很大时优化器会利用这种不连续形成单像素交替的高刚度路径。解决把灵敏度滤波半径 rmin 至少调到单元尺寸的 1.5 倍以上。在我这个骨架代码里 rmin 2.4 是在 60x20 网格上验证过的经验值。注意滤波时权重必须除以权重和否则灵敏度的量纲会漂移收敛后的拓扑偏刚或偏柔。5.2 结果里灰度一堆BESO 跑成了 SIMP现象x 数组里 0.3 到 0.7 的中间值占了一大半设计变量根本没有收敛到 0/1。原因最常见是把 SIMP 的密度投影逻辑直接搬进 BESO而 BESO 的设计变量更新逻辑是“排序后硬阈值”不是“对 rho 做梯度下降”。另一类原因是惩罚因子 p 设置过小比如 p 1中间密度没有惩罚SIMP 本身就收不到 0 和 1。解决BESO 的单元模量直接用 E Emin x * (E0 - Emin)不要用 x 的幂次每次更新后可以对 x 做一个硬阈值投影x 0.5 置 0x 0.5 置 1。这个二值化会让目标函数曲线出现小的跳变属于正常现象不用担心。5.3 加载点周围出现蘑菇头状实心块现象集中力作用的节点附近糊了一大块实体像个蘑菇头而真正有意义的承载杆件在别处被删得七零八落。原因集中力在节点处产生应力奇异性局部应变能极高灵敏度排序永远把这块区域排在“不能删”的最前面。优化器宁可牺牲远端结构也要保住载荷点的局部材料。解决把载荷作用点周边 2 到 3 层单元设为非设计域这部分单元固定为实体不参与增删。另一个做法是把单个节点力分散到相邻 3 个节点上降低奇异性。这个坑在三维模型里更严重因为网格更密局部奇异的影响范围更集中。5.4 目标函数一直下降但拓扑每代都在变现象柔度曲线平滑下降体积也已到达目标值可是每次迭代输出的拓扑图案完全不一样有时候右边多一条杆下一轮又消失了。原因收敛判据只看了目标函数变化没看拓扑变化。柔度在最优解附近本来就平缓可能只差 0.1%拓扑却还在大范围调整。还有一个隐藏原因双向增删阶段的阈值取单代灵敏度排序代与代之间的噪声会让边界反复横跳。解决收敛条件改成两条同时满足——体积误差小于 1%且连续 10 代目标函数相对变化小于 0.5%。如果还震荡对灵敏度做历史加权平均把上一代的灵敏度按 0.5 权重混合进来。这个技巧在论文里叫 averaging ratio工程上直接好用。5.5 从二维扩到三维内存先炸求解器再卡住现象把 60x20 换成 60x40x20 的三维网格程序在组装刚度矩阵时内存直接耗尽或者用稠密求解器几分钟算不出来。原因三维网格的自由度数是二维的几十倍以上如果全局刚度矩阵还按稠密数组存内存直接爆炸滤波函数还是双层循环每步都要做几百万次距离计算慢是必然的。解决三维代码必须满足三个条件刚度矩阵用 scipy.sparse 的 CSR/CSC 格式滤波改用 ndimage 的均匀滤波或高斯滤波替代逐单元循环求解器从 spsolve 换成迭代法比如共轭梯度配 ILU 预条件。我踩过最深的一次翻车就是把二维滤波直接平移结果一次迭代跑了二十分钟。6. 把 BESO 结果变成能出图的边界等值线提取与一个小习惯BESO 输出的 x 是 0/1 矩阵画成像素图能用但拿去做 CAD 和有限元复核就得先提边界。最常见做法是用 marching squares 提取 0.5 等值线一个函数就能完成from skimage import measure import numpy as np x2d x.reshape(nely, nelx).astype(float) contours measure.find_contours(x2d, 0.5) # contours 是多段折线直接导出 DXF 前先用高斯平滑去掉单像素毛刺 from scipy.ndimage import gaussian_filter x_smooth gaussian_filter(x2d, sigma0.8) contours_smooth measure.find_contours(x_smooth, 0.5)一个有价值的验证技巧是把 BESO 的拓扑结果做一次镜像再跑一遍有限元对比镜像前后的柔度。偏差小于 1% 说明边界条件和网格映射处理得没有方向性错误偏差大先检查载荷和约束的自由度索引。对称约束是个更强的手段把设计域设为全宽、载荷对称优化结果若不关于中心线对称说明灵敏度滤波的窗口扫描方向或网格编号有 bug。这个检查比盯着彩色云图可靠得多。我自己的一个习惯是任何拓扑优化结果先提取 0.5 等值线看一眼再回去看彩色灵敏度图。只看灵敏度云图很容易被灰度骗到以为结构快要收敛了。边界折线一拉出来单像素断裂、悬浮材料、锐角尖点全都藏不住。这个习惯帮我省了不少后悔药。希望帮到你。本文还有配套的精品资源点击获取