SIMP、BESO、Level Set：三大拓扑优化方法原理与避坑指南

简介：一套基于水平集方法的BESO拓扑优化Matlab开源代码，面向结构优化、工程设计与科研人员，用于在复杂载荷与约束条件下自动寻找最优结构布局。水平集方法擅长模拟形状演化与拓扑变化，BESO则通过迭代移动边界实现优化，与SIMP法相比可避免中间密度的灰色区域，得到边界清晰、易于制造的设计结果。压缩包内含10个.m源文件，整体仅5KB，但覆盖了完整优化流程：主程序调度迭代，刚度矩阵模块完成有限元分析，符号距离函数重初始化维持水平集稳定，演化方程更新驱动边界移动，另有材料信息与拓扑描述等辅助函数，模块划分明确，便于局部修改与算法对比。该实现将水平集灵活性与BESO高效性结合，适用于结构轻量化、刚度最大化等典型问题，是理解先进拓扑优化算法的优秀范例。目前已有306人学习浏览，适合具备Matlab与数值优化基础的学习者，基于此代码快速上手并针对实际需求进行二次开发。

1. 一个 .rar 打包三种拓扑优化方法：先分清 SIMP、BESO、Level Set 再动手

手头这个 Level Set Method.rar 的命名方式很典型：把三种主流拓扑优化实现一起打包，文件名里同时挂着 BESO、SIMP 和 Level Set Method。这不是整理者有收集癖，而是这个领域最常见的开源资源组织方式——三条路线分别代表拓扑优化的三种底层思路：SIMP 用连续密度加惩罚项逼近 0/1 分布，BESO 直接做离散单元的增删进化，Level Set 则把结构边界放进隐式函数的零等值面，让边界自由分裂、融合。你拿到这个压缩包时面对的其实是一个选型问题：同一批载荷和约束，方法选对事半功倍，选错会在棋盘格和振荡曲线里折腾两周。这篇笔记就把三条路线的原理、代码骨架、参数边界和踩坑记录拆开讲，按自己的网格规模和制造约束对号入座。

2. SIMP 拓扑优化：连续密度加惩罚，先跑通最小柔度问题

2.1 密度变量与惩罚机制：为什么 penal 默认取 3

SIMP（Solid Isotropic Material with Penalization，固体各向同性材料惩罚模型）是最容易上手的拓扑优化入口。它的设计变量是每个单元的伪密度 x_e，取值 0 到 1 之间，单元弹性模量按 E_e = E_min + x_e^penal * (E_0 - E_min) 插值。这里的 penal 是惩罚因子，典型取 3。目的是让中间密度（0.4 到 0.6 的模糊单元）在优化过程中被惩罚——因为中间密度的刚度收益低于材料成本，优化器会自动把它们推向 0 或 1，最终得到一个近似 0/1 的清晰拓扑。

penal 为什么是 3 而不是 2 或 4？取 2 时柔度函数接近线性，优化问题相对光滑，收敛快但中间密度大量残留；取 4 时边界干净，但非凸性变强，迭代容易卡在局部最优。3 是两者之间的平衡点，这也是绝大多数开源 SIMP 代码默认 penal=3 的原因。另一层的道理在于：penal 其实在改变设计空间的形状，太小的惩罚等于允许材料"半存在"，工程上没法解释；太大则是把连续优化硬掰成离散搜索，数值稳定性明显下降。

2.2 最小柔度问题核心循环：从 99 行代码移植的骨架

最经典的学习入口是 Sigmund 的 99 行 MATLAB 代码，后续大量 Python 开源实现基本保留同一结构：有限元求解位移场、算单元应变能灵敏度、过滤、OC 准则更新、检查收敛。这个核心循环我摘出来，改成更接近工程可读的 Python 版本：

import numpy as np def simp_core(nelx, nely, volfrac, penal=3.0, rmin=1.5, max_iter=200): """SIMP 最小柔度拓扑优化核心循环。 nelx/nely: 水平/垂直单元数；volfrac: 体积分数上限。 """ n = nelx * nely x = volfrac * np.ones(n) # 均匀初值，避免冷启动成型困难 change = 1.0 loop = 0 while change > 0.01 and loop < max_iter: loop += 1 # 1) 由 x 插值单元刚度：E_e = E_min + x**penal * (E_0 - E_min) # 2) 组装全局 K，施加边界条件，求解 K u = f（4 节点四边形单元） # 3) 单元应变能计入灵敏度：dc_e = -penal * x_e**(penal-1) * u_e^T k0 u_e # 4) 灵敏度过滤（rmin 卷积，抑制棋盘格） dc = filter_sens(dc, x, nelx, nely, rmin) # 5) OC 准则更新设计变量 x_new = oc_update(x, dc, volfrac, move=0.2) change = float(np.max(np.abs(x_new - x))) x = x_new if loop % 20 == 1: print(f"iter {loop:3d} change={change:.4f} volfrac={x.mean():.4f}") return x.reshape(nely, nelx)

这里有限元组装被省略，因为不同问题（悬臂梁、MBB 梁、L 形梁）的载荷和边界位置各不相同，单元刚度矩阵组装本身是最容易抄错的部分，不值得在循环里反复展开。灵敏度的符号是负的：柔度越小越好，所以灵敏度为负值时密度倾向于增大，理解这一点后续排查 BESO 增删方向时会很有用。

OC 准则是 SIMP 最常见的更新器，代码单独拆出来：

def oc_update(x, dc, volfrac, move=0.2): """OC 准则：拉格朗日乘子二分法更新密度。""" l1, l2 = 0.0, 1e5 # 乘子二分区间 while (l2 - l1) / (l1 + l2) > 1e-4: lmid = 0.5 * (l1 + l2) x_new = np.maximum(0.0, np.maximum( x - move, np.minimum(1.0, np.minimum( x + move, x * np.sqrt(-dc / lmid))))) if np.sum(x_new) - volfrac * x.size > 0: l1 = lmid else: l2 = lmid return x_new

参数含义：move 限制单步密度最大变化量，默认 0.2，防止迭代跳变；l1/l2 是拉格朗日乘子二分区间，l2 取 1e5 保证上限覆盖多数问题。如果换 MMA（移动渐近线法）求解器，更新逻辑不同但灵敏度输入完全一样——SIMP 的灵敏度公式和求解器解耦，这是它生态丰富的重要原因。

2.3 四个必调参数：penal、rmin、volfrac、move

任何一份 SIMP 开源代码，核心参数就这么几个，我整理成一张常用参数表：

rmin 和网格尺寸的关系必须注意：rmin 通常取单元尺寸的 1.5~3 倍，它直接决定结果的最小特征尺寸，也决定了结构能不能被 3D 打印。换句话说，rmin 是 SIMP 里唯一同时影响数值稳定性和制造约束的参数。很多开源代码里 rmin 的默认单位是"单元个数"，改为绝对长度时换算错一位，结果就会从清晰桁架变成一团糊。

volfrac 设太小时常见的问题是结构出现悬空的孤岛——材料块和支撑完全断开。这是因为过滤让高密度区域连成片，但体积约束太紧时优化器会保留局部高刚度区而牺牲连通性。工程上我一般先跑 0.4 看整体构型，再往下压到目标值，而不是直接用最终体积分数冷启动。

2.4 从教学代码到工程：多工况加权与对称约束

教学代码默认只有一个工况和一个载荷点，真实产品往往同时承受静载、振动和温度载荷。多工况的标准做法是把每个工况的柔度按权重线性叠加，灵敏度也同步加权。我一般把目标写成 sum(w_i * C_i)，w_i 按工况重要性给，而不是平均分配。这个改动在代码里只影响目标函数和灵敏度累加两处，收益却非常明显——零件不会在某一个工况下被优化得特别好，却在另一个工况下直接断裂。

另一个实用技巧是善用对称约束。MBB 梁的标准算例就是只算半边模型，网格规模直接减半，迭代速度接近翻倍。不过对称约束会隐藏非对称的失稳模式，如果是压杆或薄壁结构，完全对称的载荷下也可能出现反对称屈曲，这时候强制对称反而会得到偏危险的设计。我自己的习惯是：先用对称模型快速迭代选型，锁定拓扑后放开对称做一次全模型校验。

提示：过滤半径的单位在不同开源代码里约定不同，有的是单元尺寸倍数，有的是绝对长度。读代码前先确认，否则 rmin 参数在两个项目之间迁移时会得到完全不同的拓扑。

3. BESO 拓扑优化：离散 0/1 变量的双向进化与灵敏度过滤

3.1 双向进化的逻辑：ESO 的后悔药

BESO（Bi-directional Evolutionary Structural Optimization，双向渐进结构优化）的思路和 SIMP 完全不同：每个单元要么存在（1）要么删除（0），没有中间状态。每一轮先做有限元分析，算出每个单元的灵敏度，然后同时做两件事——删掉灵敏度最低的实单元，添加灵敏度最高的空单元（空单元指被删除但仍在设计域内的格子）。"双向"就体现在这里：传统 ESO 只删不加，删过头就再也回不来；BESO 的添加机制相当于给了后悔药，允许中间过程修正过失。

这个机制让 BESO 对初始设计不敏感。从满材料开始，或者从均匀随机分布开始，最终都能收敛到相近的拓扑，这是它和 Level Set 的最大差异之一。代价是离散变量的优化问题非凸，算法没有严格的收敛证明，工程上靠体积约束和柔度变化组合判断何时收手。开源代码里 BESO 的实现通常比 SIMP 复杂一些，核心在于灵敏度过滤和增删策略，而不是有限元部分——两者共用同一套求解器。

3.2 灵敏度计算与三次平均过滤：BESO 的稳定性命门

BESO 灵敏度公式在最小柔度问题里和 SIMP 同源，都是单元应变能密度的导数，但处理方式有本质区别。SIMP 的过滤是作用于连续密度场的低通滤波；BESO 的过滤则必须针对 0/1 离散状态做特殊处理，常见做法是"单元灵敏度映射到节点、节点邻域加权平均、再映射回单元"的流程，并且通常重复三次。

def beso_sensitivity_filter(dc, x, nelx, nely, rmin): """BESO 灵敏度过滤：节点映射 + rmin 加权平均 + 映射回单元。 与 SIMP 的密度过滤不同，操作对象是离散单元的灵敏度。 """ # 1) 单元灵敏度映射到节点，每个节点累加相邻单元灵敏度 node_sens = np.zeros((nely + 1, nelx + 1)) # 2) 以 rmin 为半径做节点加权平均（含镜像边界处理） node_sens = average_with_radius(node_sens, rmin) # 3) 节点灵敏度映射回单元，取单元四个角点平均值 dc_filtered = map_to_element(node_sens) # 重复三次，抑制单个单元噪声对删/增判断的影响 return dc_filtered

三次平均的目的不是让图更平滑，而是让删单元的判断不因为单个单元的数值波动翻车。BESO 对灵敏度噪声比 SIMP 敏感得多——每个单元只有两个状态，一个噪声单元就可能让一整块材料被误删。过滤半径 rmin 的取值逻辑和 SIMP 一样，但注意 BESO 里 rmin 通常需要比 SIMP 略大，因为离散状态下的滤波效果更弱。

3.3 进化率 ER 与最大添加率 ARmax：增删节奏的配合

BESO 的增删策略直接决定能否收敛。两个参数必须配合：ER（Evolutionary Rate）是每轮体积减少比例，典型值 0.02；ARmax 是单轮最多添加单元数的上限比例，典型值 0.05。

def beso_add_remove(x, dc, volfrac, ER=0.02, ARmax=0.05): """BESO 增删：按灵敏度阈值批量翻转 0/1。""" n = x.size v_cur = x.sum() v_target = max(v_cur * (1.0 - ER), volfrac * n) # 逐步逼近目标体积 dc_sorted_idx = np.argsort(dc) # 升序，前面是低灵敏度 n_remove = int(v_cur - v_target) # 本轮要删的单元数 n_add = int(min(v_target - v_cur + n_remove, ARmax * n)) for idx in dc_sorted_idx[:n_remove]: # 删：灵敏度最低的实单元 if x[idx] == 1: x[idx] = 0 for idx in dc_sorted_idx[::-1][:n_add]: # 加：灵敏度最高的空单元 if x[idx] == 0: x[idx] = 1 return x, v_target

ER 决定每轮删掉多少体积，ARmax 控制最多加回多少。两个参数不匹配的典型情况：ER 大、ARmax 小，结构只会缩不会长，失去"双向"意义,退化成带添加限制的 ESO；ER 小、ARmax 大，结构在目标体积附近来回抖动，目标函数振荡。我一般固定 ER=0.02、ARmax=0.05 起调，最后的收敛判据设为：体积达到目标，且连续 10 轮柔度变化小于 0.1%。

这里有一个经常被忽略的细节：删除单元后，空单元必须保留在设计域里参与下一轮灵敏度计算。如果代码把空单元从模型中彻底移除，BESO 就退化成 ESO，只能删不能加，柔度必然持续上升且无法回头。检查方法是在迭代中打印"存在单元数"和"参与计算单元数"两个量，正常情况下前者小于后者，如果二者相等，说明空单元被错误排除。

3.4 与 SIMP 的取舍：什么时候选 BESO

SIMP 的结果天然带有 0.01 到 0.99 的过渡单元，进入制造前需要阈值化处理；BESO 本身就是 0/1，省掉了后处理。但 BESO 的迭代轮数通常是 SIMP 的 2 到 3 倍，每轮都要完整求解有限元，网格大时成本明显更高。我的经验是：网格规模很小、需要直接对接制造的 0/1 单元模型，用 BESO；大规模快速迭代，先用 SIMP 粗筛拓扑，再局部精修。

另一个现实问题是 BESO 的灵敏度符号约定。目标是最小柔度时，灵敏度取负值，删除"灵敏度最低"即对刚度贡献最小的单元。如果你拿到的开源代码里出现删除灵敏度最高单元的逻辑，目标函数一定在上升。检查方法很简单：完整跑一轮，对比首末柔度，上升就说明增删方向反了，把排序反向即可。

4. Level Set Method：把边界演化和拓扑成核分开看

4.1 零等值面、Hamilton-Jacobi 方程与速度场

Level Set Method 和 SIMP/BESO 最大的区别：结构的边界不显式存储为节点坐标或多边形网格，而是藏在一个标量函数 φ(x) 的零等值面里。设计域内 φ < 0 是材料域，φ > 0 是空域，φ = 0 就是结构边界。边界移动通过求解 Hamilton-Jacobi 方程完成：

∂φ/∂t + Vn |∇φ| = 0

Vn 是边界法向速度，由形状导数推导。对于最小柔度问题，最常见做法是取 Vn = ξ_e - λ，ξ_e 是单元应变能密度，λ 是拉格朗日乘子，通过二分法求得以保证体积守恒。这个形式和 SIMP 的灵敏度只差一个符号约定，但意义完全不同——SIMP 在修改单元密度，Level Set 在移动一条真实几何边界。这也是为什么 Level Set 的结果边界最光滑：边界本身就是一条连续曲线，不存在 SIMP 那种锯齿状过渡带。

4.2 演化实现与 CFL 条件：dt 和重新初始化怎么配

实现层级上，常见的开源代码流程是：初始化 φ 为符号距离函数 → 有限元求解 → 计算应变能密度 → 构造速度场 Vn → 用迎风格式离散 ∇φ → 更新 φ → 周期性重新初始化 φ。核心演化代码：

def level_set_evolve(phi, strain_energy, lam, dt, dx): """水平集一步演化。 phi: 符号距离函数（负值=材料域，正值=空域） strain_energy: 单元应变能密度（有限元后处理得到） lam: 拉格朗日乘子，保证体积守恒 """ # 1) 速度场：Vn = strain_energy - lam vn = strain_energy - lam # 2) 迎风差分计算 |grad(phi)| grad_norm = upwind_gradient_norm(phi, dx) # 3) 演化：phi^{n+1} = phi^n - dt * vn * |grad(phi)| phi_new = phi - dt * vn * grad_norm # 4) 周期性重新初始化 phi 为符号距离函数 # （每 5~10 步执行一次，不是每步都做） return phi_new

时间步长 dt 必须满足 CFL 条件：dt ≤ dx / max(|Vn|)，dx 是单元尺寸。很多 Level Set 开源实现翻车都翻在 dt 上——取太大，φ 的梯度在几次迭代后严重变形，零等值面开始抖动，最终的拓扑像被啃过一样；取太小，迭代轮数成倍增长。稳妥做法是每步动态计算 max(|Vn|)，把 dt 压到 0.5 倍 CFL 上限。

提示：固定步长的 Level Set 代码在网格加密后几乎必翻车，因为 CFL 条件随 dx 变化。改用动态 dt 后，同一套代码可以从粗网格直接换到细网格而不用重新调参。

重新初始化是另一个常见坑。符号距离函数在演化中会逐渐失去"距离"性质，导致边界移动速度被扭曲，必须每 5 到 10 步重做一次，常见做法是重新求解一个 Eikonal 方程把 φ 拉回符号距离形式。频率太高会引入额外的数值耗散，频率太低则边界变形明显，这个平衡得靠观察拓扑演化来调。

4.3 初始洞的数量决定拓扑自由度

Level Set 方法有一个物理特性：拓扑变化（洞的生成和消失）只能通过边界的合并与分裂实现，边界本身不能凭空成核。也就是说，初始 φ 里画了多少个洞，最终拓扑的"孔洞自由度"就被限制了多少。一个实心方块初始设计，最终只能得到外轮廓变形，内部空腔基本不可能出现——这是 Level Set 和 SIMP/BESO 最大的行为差异。

解决方法很朴素：初始设计域里人为均匀布置小洞。常见做法是用 3x3 或 5x5 的规则孔阵，洞的半径取单元尺寸的 1 到 2 倍。洞太多会增加非凸性，迭代容易卡在奇怪的局部最优；洞太少会丢失有效的拓扑解。我一般从 3x3 起步，看目标函数和拓扑演化方向再决定加密布局。另有一个现代变体叫拓扑导数，在速度场里加入成核项允许新洞生成，实现复杂度明显更高，初学不建议直接碰。

4.4 三种方法选型对比：网格、制造成本与实现代价

选型逻辑我总结成一句话：快速验证用 SIMP，制造对接要 0/1 网格用 BESO，边界形状精度和应力集中敏感问题用 Level Set。三个方法共享同一套有限元求解器，几何和载荷可以直接切换方法做对比，这也是这类开源压缩包存在的意义——不用重复造基础轮子。

5. 避坑清单：SIMP、BESO、Level Set 各自的翻车现场

5.1 棋盘格：SIMP 的经典数值黑匣子

现象：优化结果像国际象棋棋盘，每个单元和邻居颜色相反，目标函数看着很低，但结构完全没法制造。

原因：单元密度和柔度之间存在非物理的交替模式，有限元离散在这种模式下会虚低应变能，优化器被骗了。这是网格依赖性的直接表现，网格加密后结果可能完全不同，同一套参数换网格就出问题。

解决：灵敏度过滤或密度过滤，rmin 从单元尺寸的 1.5 倍起步。过滤不够就加大 rmin；如果加大到 3 倍仍有棋盘格，检查过滤代码是不是写成了只对当前迭代的 x 过滤而不是对灵敏度过滤——这两个效果差别很大，很多移植版翻在这里。

5.2 柔度锯齿振荡：BESO 增删比例失调

现象：柔度曲线不下降反而上下跳动，甚至每轮都在上涨，迭代 50 轮后拓扑还在剧烈变化。

原因：增删比例失调。常见是 ER 设 0.02 但 ARmax 设 0.1，添加过猛；或者三次平均过滤没做，单元噪声让删单元的位置反复横跳。另一个可能是空单元被物理删除，BESO 退化为 ESO，结构只缩不长。

解决：把 ARmax 压到 0.05 以内，确认三次平均过滤生效。空单元必须保留在设计域内参与灵敏度计算，检查"参与计算单元数"是否大于"存在单元数"。按这个顺序排查，九成以上的振荡问题都能定位。

5.3 Level Set 出不了内部空腔

现象：优化结果只有外轮廓变形，内部该有的筋、孔全没有，看起来像"缩水版"结构。

原因：初始 φ 是实心方块，没有种洞，水平集边界无法成核。这是方法本身的几何限制，不是求解器 bug。

解决：初始化时在域内均匀布置小洞，半径取 1 到 2 倍单元尺寸，起始布局用 3x3。同时检查重新初始化频率，如果每两轮就强制重初始化，小洞容易被抹平，适当放宽到 5 到 10 轮一次。如果想彻底解决成核问题，再调研拓扑导数方案。

5.4 rmin 和网格尺寸的连锁反应

现象：SIMP 结果要么糊成一片、要么碎成渣，怎么调都不像经典梁结构；BESO 那边则是删单元的位置飘忽不定。

原因：rmin 和网格尺寸矛盾。rmin 取 0.5 倍单元尺寸时过滤几乎失效，棋盘格和碎条状结构出现；rmin 取 5 倍以上时所有特征被抹平，结构变成一大坨。

解决：先固定网格，再扫 rmin：分别试 1.5、2.0、3.0 倍单元尺寸，看拓扑变化趋势。工程上我更愿意先用小 rmin 拿到细节丰富的拓扑，再在制造约束阶段检查最小杆径。rmin 数值不追求精确，它在结果里直接对应最小特征尺寸，要和打印工艺的壁厚要求对齐。

5.5 体积约束失效：二分法爆炸与增删数量错位

现象：代码不报错，但体积分数一路飘，不到 5 轮就超出 volfrac 两倍；BESO 那边则是目标体积到了但结构还在持续删除。

原因：SIMP 的 OC 更新里拉格朗日乘子二分区间没初始化好，或者 move 限制和二分条件冲突；BESO 是增删数量计算失误，用了当前体积减目标体积的绝对值，导致目标和实际永远对不上。

解决：SIMP 检查 l1/l2 初始区间，l2 取 1e5 是常规起点；BESO 检查 v_target 的计算逻辑，每轮输出当前体积、目标体积、增删数量三个值，人工核对第一轮就应该建立单调递减关系。这个排查就是打印中间量，拓扑优化没有太多黑匣子，大多数问题在数据流里原形毕露。

6. 从跑通到可信：基准算例、调参次序与 3D 化的现实

6.1 三个基准算例建立参照系

接手任何一份新的开源拓扑优化代码，先别急着跑自己的模型。MBB 梁（两端支撑、中点载荷）、悬臂梁（一端固定、另一端加载）、L 形梁（内角应力集中）是三个公认基准。在这三个算例上跑到文献经典图对得上，再谈工程问题。我习惯用对称约束只算半边 MBB 梁，既加速又能对照教科书结果，同时验证代码的对称性处理是否正确——错误代码算出来的半边模型往往不对称。

6.2 调参次序：先网格后算法，先粗后细

调参次序固定为：先锁定网格和载荷，再调体积分数，最后才碰算法参数。penal、ER、rmin 这些一次只动一个，每次只改 20% 到 50% 比例，记录拓扑外形和柔度曲线。算法参数同时动两个，出了问题不知道是谁造成的，这是最常见的调参翻车方式。用这个次序，三方法从同一个全局刚度矩阵出发，对比才有意义。

3D 化时的额外注意点是内存。三维网格单元数随尺寸立方增长，稀疏矩阵求解成为瓶颈，99 行代码的二维逻辑到三维可以直接扩，但内存可能先爆炸。常见做法是先 2D 设计、拉伸成 3D，再在关键截面做局部 3D 优化，避免一上来就全模型 3D。

这些年处理过不少这类开源拓扑优化压缩包，最大的体会是：不要一开始就追求最好看的拓扑，先接受一个能复现基准算例的结果，再逐步逼近。参数玄学大多出在过滤和增删逻辑上，而不是有限元求解器；把三个方法的边界条件理清楚，这个 .rar 里的代码就能变成一套真正顺手的工具。希望帮到你。

本文还有配套的精品资源，点击获取