简介：水平集方法与BESO结合的拓扑优化开源实现，面向结构优化领域的研究人员，以及具备一定MATLAB和数值计算基础的工程师、学生，用于解决最小化负载传递路径、满足设计约束的结构布局寻优问题。压缩包内共10个文件，全部为.m脚本，涵盖主程序、刚度矩阵组装、水平集函数初始化与重初始化、边界演化、有限元分析等核心模块，覆盖完整优化迭代链路。整个包仅5KB，代码紧凑、便于逐行研读和修改调试。已有306人学习下载。通过该资源可直接运行基础算例，观察水平集方法如何清晰表达结构边界，规避SIMP法中中间状态“灰色区域”问题，并可根据设计需求调整载荷、体积约束等参数，实现定制化拓扑布局探索。对希望深入算法内部原理或进行二次开发的研究者，是一份轻量而实用的参考实现。

1. 拓扑优化选型从哪下手：一个标题串起 Level Set、BESO 与 SIMP 的全局

前阵子做支架轻量化，甲方要求刚度不降、重量砍掉三成。我第一反应是用 SIMP 跑一版密度法拓扑优化，但检索资料时看到不少免费分享把 Level Set Method、BESO、SIMP 三个关键词捆在一起，才发现这是一类问题的三条主线：拓扑优化里最常用的三种开源实现路径。这篇文章不从理论课本讲起，直接回答工程师最关心的四个问题——这三种方法差别在哪、怎么用开源代码快速跑通、参数怎么设、遇到棋盘格和局部极值时怎么收场。适合刚接手结构优化任务的机械工程师、CAE 工程师，也适合想做轻量化对比研究的学生。读完你应该能独立跑通一个 2D 算例，并且知道换方法时该改哪里。

2. 三种方法的关系：SIMP、BESO 与 Level Set 的数学模型和适用边界

用手算一个三杆桁架，你就能明白拓扑优化在干什么：给定载荷和支撑，材料该放哪里、该掏掉哪里，才能让结构最刚。SIMP、BESO、Level Set 是三种不同的“填空”策略，它们回答同一个问题——哪些材料对刚度贡献最大，哪些是纯负担。理解这三条线的差别，是后续所有调参的基础。

2.1 SIMP 的罚函数模型与灵敏度计算基础

SIMP（Solid Isotropic Material with Penalization）把每个单元当成一个可连续变化的密度变量 ρ，允许 0 到 1 的中间值，但优化目标是让中间值消失、解最终落到 0/1。为了实现这个目的，材料弹性模量被写成 E(ρ) = E_min + ρ^p (E_0 - E_min)。

p 就是惩罚因子，通常取 3。p 越大，中间密度单元的“性价比”越低，优化器会主动把它们推向 0 或 1。这里的 E_min 不能设成 0，否则刚度矩阵奇异，求解器直接报错。

灵敏度的计算是 SIMP 的关键环节。目标函数一般为柔度 C = F^T U，也就是结构在载荷下的应变能。柔度对某个单元密度的导数 ∂C/∂ρ_e 可以用伴随法求，最终表达式只包含该单元的位移、应力与惩罚因子 p。开源实现里最常见的 99 行 SIMP 代码就是一个完整闭环：有限元求解、灵敏度计算、OC 更新准则、灵敏度滤波，全部浓缩在不到一百行里，用 MATLAB 或 Octave 就能跑。

注意：SIMP 的中间密度只是数学工具，不是材料属性。后处理必须设定阈值（通常 0.5）把连续密度图转成黑白构型。

2.2 BESO：从连续密度到 0/1 离散解

BESO（Bi-directional Evolutionary Structural Optimization）走的是另一条路：它不是让密度连续变化，而是每步直接判定单元去留。单元要么存在（密度 1），要么删除（密度 0），没有中间态。“双向”指的是每一步既删单元也加单元：删除的依据是单元灵敏度，灵敏度低说明这块材料对刚度贡献小；同时允许被删的单元在后续迭代中重新长回来，前提是它现在的灵敏度足够高。这个“后悔药”机制让 BESO 比早期 ESO 稳定很多。

BESO 的灵敏度还有一个特色处理：历史平均。因为 0/1 离散化会让灵敏度震荡，常见做法是把当前灵敏度与上一步的灵敏度做加权平均，这个平滑操作直接影响收敛稳定性，是实现中最容易漏掉的细节。

BESO 不需要惩罚因子 p，因为不存在中间密度。核心参数是删除率（进化率 ER）和体积约束的逼近策略。对工程人员来说，BESO 输出天然是黑白结果，不需要做密度阈值化，这是它在工程落地中的一大优势。

2.3 Level Set：边界隐式表达与拓扑变化的代价

Level Set Method 的思路换了个坐标系：不再问每个单元用不用，而是用一个函数 φ(x) 在空间里画一条零水平集作为结构的边界。φ > 0 的区域是材料，φ < 0 的区域是空洞，边界就是 φ = 0 的等值面。

优化过程中，结构边界在速度场作用下移动，速度由灵敏度推导而来。数学上要解 Hamilton-Jacobi 方程：∂φ/∂t + V·∇φ = 0。这个方程的本质是让边界沿法线方向运动，V 是法向速度，正负决定边界向外扩张还是向内收缩。

Level Set 最大的优势是边界天然光滑，不会有棋盘格这类密度场伪影；拓扑变化（孔洞的生成与消失）由水平集的拓扑演化自然处理。缺点也很现实：方程是时间依赖的，需要 CFL 条件控制时间步，实现复杂度比 SIMP 高一个量级；φ 场需要周期性重新初始化为符号距离函数，否则边界运动失真。

2.4 三种方法的选型逻辑：什么时候用哪个

选型上没有绝对的“最优方法”，取决于你要交付什么。如果只出概念构型，SIMP 半小时就能出图；结果要直接丢进五轴加工路径规划，Level Set 的光滑边界能省大量模型修复时间；BESO 则处在两者之间，适合想快速看到 0/1 方案又不想被棋盘格纠缠的团队。我的建议是先把 SIMP 跑通，因为它包含了灵敏度、滤波、优化准则更新这套核心思维。换到 BESO 和 Level Set 时，你至少能看懂每一步在干什么。

3. 用开源代码跑通第一个拓扑优化算例：以 99 行 SIMP 为例

拓扑优化开源实现很多，但绝大多数教学代码都能追溯到 Ole Sigmund 的 99 行 SIMP 程序。它短、无依赖、结果直观，适合作为第一个跑通的算例。

3.1 环境准备与最小代码骨架

你需要一个支持矩阵运算的环境。GNU Octave 是开源且与 MATLAB 语法高度兼容的选择；有 MATLAB 授权就用 MATLAB，跑法一致。环境就绪后，把下面的最小骨架存成 top99.m：

% top99.m - 99行SIMP拓扑优化最小闭环骨架（Octave兼容） function top99(nelx, nely, volfrac, penal, rmin) % 输入: nelx/nely为单元网格数, volfrac为体积约束, % penal为惩罚因子, rmin为滤波半径 x(1:nely, 1:nelx) = volfrac; % 均匀初始密度 for loop = 1:200 % 1) 组装刚度矩阵并求解位移场 [U, KE] = FE_solve(nelx, nely, x, penal); % 2) 计算每个单元的柔度灵敏度 dc = sensitivity(nelx, nely, KE, U, x, penal); % 3) 灵敏度滤波，抑制棋盘格 dc = filter_sens(nelx, nely, rmin, x, dc); % 4) OC准则更新密度场 x = update_oc(nelx, nely, volfrac, x, dc, loop); end end

这个骨架把问题拆成四个子函数：有限元求解、灵敏度计算、滤波、OC 更新。拆开的好处是，后续换成 BESO 或 Level Set 时，你只需要替换“更新准则”这一个子函数，前面三步可以复用。常见做法是把四个子函数写在同一个 .m 文件里方便调试，Octave 下直接把它拖进命令行运行即可，不需要额外装工具包。

3.2 悬臂梁算例的实现：主循环与灵敏度

真正 99 行程序里最核心的是 OC 更新准则。它本质上是寻找满足体积约束的拉格朗日乘子，让密度按灵敏度比例变化：

% OC更新准则：二分法求乘子，并更新密度 function xnew = update_oc(nelx, nely, volfrac, x, dc, loop) l1 = 0; l2 = 1e9; move = 0.2; while (l2 - l1) / (l1 + l2) > 1e-3 lmid = 0.5 * (l1 + l2); xnew = max(0, max(x - move, ... min(1, min(x + move, x .* sqrt(-dc / lmid))))); if sum(xnew(:)) > volfrac * nelx * nely l1 = lmid; else l2 = lmid; end end end

逻辑说明：这里的 dc 是柔度灵敏度，lmid 是虚构的“材料价格”。密度更新公式中的 sqrt 项来自 OC 准则推导，保证每一步都朝更优方向走；move 限制每步密度最多变化 0.2，防止迭代震荡。二分法循环结束后得到满足体积约束的新密度分布。

参数说明：move 取值 0.1~0.3 常见，越大收敛越快但越容易震荡。二分法容差 1e-3 对工程精度足够，再小只会增加循环次数。悬臂梁算例中网格用 60×30，左端固定、右下角施加集中力，这是拓扑优化最经典的验证场景之一。

3.3 参数说明：volfrac、penal、rmin 三个核心参数

这三个参数互相牵连：penal 太高会让灵敏度对中间密度过于敏感，rmin 太小则棋盘格清理不干净。我常用的调参顺序是固定 volfrac，把 rmin 从 1.2 逐步加到 3，找到没有棋盘格且传力路径连续的最小值，再微调 penal。注意 rmin 的单位是单元尺寸而不是绝对长度，网格一变，同样的 rmin 对应的物理尺度就变了，网格加密时要同步调整。

3.4 结果判读：目标函数曲线与拓扑构型

跑完迭代先看目标函数（柔度）的迭代曲线。健康的曲线是前二十步快速下降，之后缓慢收敛，最后几乎水平。如果曲线像锯齿一样上下跳动，说明 move 或滤波半径不合适。

再看最终密度云图：黑白分明、传力路径清晰是正常状态；大片灰色说明 penal 不够或没收敛；细碎斑点、交替的实心与空心方块就是棋盘格，需要加大 rmin。做这一步时不要用肉眼看色阶，直接把密度矩阵按 0.5 阈值二值化，再统计黑白单元是否成片分布，这样判断更可靠。

注意：2D 算例任何时候都不要直接拿最终构型去加工。拓扑优化结果是概念设计，圆角、装配接口、尺寸细节都要结构工程师重新绘制，这一点在第五、六章还会展开。

4. BESO 与 Level Set 的开源实现路径：从迭代逻辑到工程文件

跑通 SIMP 之后，你会觉得拓扑优化不过如此。直到换 BESO 或 Level Set，才发现选方法比写代码更需要经验。这一章给出两条路径的迭代骨架和工程导出方法。

4.1 BESO 的迭代逻辑与 Python 实现骨架

BESO 的核心循环比 SIMP 更直观，因为它每一步都在做删和加的决策。用 Python 写一个骨架如下：

import numpy as np def beso_topology(nelx, nely, volfrac, er, rmin, max_iter=100): # er: 进化率，控制每步删除/重生的单元比例 x = np.ones((nely, nelx)) # 初始满材料 alpha_prev = np.zeros((nely, nelx)) # 上一步灵敏度，用于历史平均 target_volume = volfrac * nelx * nely vol_current = 1.0 for loop in range(max_iter): alpha = sensitivity_analysis(nelx, nely, x) # FEA + 灵敏度 alpha = filter_sensitivity(nelx, nely, rmin, x, alpha) alpha_hist = 0.5 * (alpha + alpha_prev) # 历史平均，抑制振荡 vol_target = max(target_volume, vol_current - er * nelx * nely) alpha_flat = np.sort(alpha_hist.flatten()) ncut = int((1.0 - vol_target / (nelx * nely)) * nelx * nely) ncut = min(max(ncut, 0), nelx * nely - 1) threshold_low = alpha_flat[ncut] threshold_high = alpha_flat[nelx * nely - 1 - ncut] delete_mask = (x == 1) & (alpha_hist < threshold_low) add_mask = (x == 0) & (alpha_hist > threshold_high) x[delete_mask] = 0 x[add_mask] = 1 vol_current = np.sum(x) / (nelx * nely) alpha_prev = alpha_hist.copy() return x

逻辑说明：threshold_low决定哪些单元删除，threshold_high决定哪些空单元重生，两个阈值从排序后的灵敏度序列里取，本质是控制每步修改的单元比例。历史平均那行是 BESO 收敛的关键，很多复现翻车都是漏了这行，导致目标函数震荡到底。这里的sensitivity_analysis和filter_sensitivity可以复用你在 SIMP 中写好的版本，只是 BESO 的灵敏度对 0/1 密度做滤波时，权重需要稍作调整。

参数说明：er 常见取 0.02~0.05，太小收敛慢，太大每步改动过猛，容易丢失重要传力路径。滤波半径 rmin 的作用与 SIMP 相同，但因为密度是离散的，它对噪声更敏感，建议从 SIMP 的推荐值开始再往上加 0.5。工程上我不会用纯 Python 从零写有限元求解，矩阵组装与求解放进去会非常慢；常见做法是复用开源 FEA 内核解位移场，BESO 骨架只负责灵敏度计算与单元删除/重生逻辑。

4.2 Level Set 的 Hamilton-Jacobi 更新与边界处理

Level Set 的代码量比前两个方法大得多，核心在把 φ 场更新实现正确。下面是一步时间推进的简化片段：

import numpy as np def update_level_set(phi, velocity, dt, dx): # phi: 符号距离函数, 材料区域 phi > 0 # velocity: 法向速度场, 由灵敏度推导 # dt: 时间步, dx: 网格间距 grad_x = np.gradient(phi, axis=1) / dx grad_y = np.gradient(phi, axis=0) / dx norm = np.sqrt(grad_x**2 + grad_y**2 + 1e-10) phi_new = phi - dt * velocity * norm return phi_new

逻辑说明：这个式子对应 ∂φ/∂t + V·|∇φ| = 0，速度场正负决定边界收缩或扩张。实际做时 velocity 不是经验值，而是由形状灵敏度推导出的场量。这里用中心差分近似梯度，严格实现应该用上风格式保证数值稳定，否则边界会出现伪振荡。另一个容易忽略的点是，φ 更新多轮后会偏离符号距离函数（|∇φ| ≠ 1），边界运动会失真，所以每 10~20 步要做一次重新初始化。

参数说明：时间步 dt 必须满足 CFL 条件，即 dt ≤ dx / max(|velocity|)，多数实现取上限的一半。网格间距 dx 在 2D 单元尺寸为 1 时也取 1，3D 同样以单元尺寸为单位。Level Set 对初始设计极其敏感：实心块初始场很难生成新孔洞，容易掉进局部极值。常见做法是先用随机小孔作为初始拓扑，或者用 SIMP 跑一个粗解，插值成 φ 的初始场，再交给 Level Set 细化边界。

4.3 把拓扑结果导出成 CAD/CAE 可用的几何

SIMP 输出的密度场只是一堆 0~1 的数，要变成工程文件就必须提取等值面。MATLAB/Octave 里用 isosurface，Python 里用 scikit-image 的 marching_cubes：

from skimage import measure import numpy as np # density_map: SIMP 输出的三维密度场, 取 0.5 等值面为材料边界 verts, faces, normals, values = measure.marching_cubes( density_map, level=0.5, spacing=(1.0, 1.0, 1.0) ) # 写出 ASCII STL，通用性最好 with open("topology.stl", "w") as f: f.write("solid topology\n") for tri in faces: f.write(f"facet normal {normals[tri[0]][0]:.6f} " f"{normals[tri[0]][1]:.6f} {normals[tri[0]][2]:.6f}\n") f.write("outer loop\n") for vid in tri: f.write(f"vertex {verts[vid][0]:.6f} " f"{verts[vid][1]:.6f} {verts[vid][2]:.6f}\n") f.write("endloop\nendfacet\n") f.write("endsolid topology\n")

逻辑说明：marching_cubes 把密度场的 0.5 等值面提取成三角形网格，spacing 参数对应 3D 网格尺寸。不做缩放的话，导出的 STL 单位是网格单元尺寸，导入 CAD 前要乘以实际模型尺寸。法向量这里取了三角形第一个顶点的法向，属于简化写法，追求水密模型时建议用跨立算三角形法向。

参数说明：level 阈值取 0.5 是把中间密度归为材料，想保守就把阈值抬到 0.6，材料区域更实。STL 导出后一定要检查孤立碎片和水密性，拓扑优化结果经常在低灵敏度区域出现少量漂浮单元，建议顺手加一步拉普拉斯平滑，再用网格检查工具做一次封闭性检测。

5. 拓扑优化落地避坑：棋盘格、局部极值与制造约束的 5 个注意点

开源代码能跑通只是第一步，真正消耗时间的是边缘问题。以下五类问题按“现象、原因、解决”给出判断路径，每一条都在实际项目中反复出现过。

5.1 棋盘格：看似合理实则脆弱的交错结构

现象：最终构型里交替出现黑白方块，整体呈棋盘状。这种结构数值上刚度不错，但实际加工后强度很差，甚至没法按图加工。

原因：有限元离散导致的高频数值振荡。棋盘格单元的灵敏度交替高低，密度场出现高频波动，这类不稳定构型在某些情况下也是局部最优，程序不会自动消除。

解决：灵敏度滤波是最有效的手段，把 rmin 从 1 提高到 2~3，棋盘格基本消失；细分网格或改用高阶单元有帮助，但成本比滤波高。经验上 rmin 小于 1.5 时棋盘格风险显著上升，低于这个值时不要轻易相信结果。

5.2 网格依赖性与滤波半径的矛盾

现象：同一问题用 60×30 和 120×60 网格算，得到的传力路径完全不同。网格越细，细节越碎，甚至出现许多须状细杆。

原因：拓扑优化天然具有网格依赖性，数学上的最优拓扑随离散精度变化。没有正则化手段时，网格越细越容易生成细小特征。

解决：滤波半径相当于给设计加上了最小特征尺寸约束。把 rmin 与绝对尺寸挂钩而不是与网格数挂钩，比如 rmin = 实际最小壁厚 / 单元尺寸，网格加密时保持绝对半径不变，结果才会趋同。工程实践里，先用粗网格定拓扑骨架，再加密网格修细节，能节省大量无效计算。

5.3 局部极值：不同初始设计得到不同结果

现象：同样的载荷与边界条件，换一个初始密度分布，最终拓扑差异很大。SIMP 相对温和，Level Set 对初始拓扑尤其敏感。

原因：拓扑优化本质是非凸优化，存在大量局部极值。优化器一旦落入局部低谷就很难爬出来，结果因此依赖初始猜测。

解决：用多种初始设计各跑一遍，比较目标函数值再做选择；或者采用先粗后细的策略，先用粗网格得到宏观拓扑，再插值成细网格的初始设计，能明显抑制随机局部极值。注意，目标函数最低的拓扑在工程上不一定最好，数值优化倾向于生成细长构件，实际加工很难保证强度，方案对比时要把可制造性放进评价标准。

5.4 制造约束：最小尺寸、拔模与对称

现象：拓扑结果里有大量小于工艺极限的细小杆件，或单一方向无法脱模的倒扣结构，直接加工就是废品。

原因：纯柔度最小化目标函数里不含制造信息，优化器不会主动避开加工禁忌。

解决：常见做法是在目标函数里加制造约束惩罚项，对小于最小尺寸的特征加惩罚，对拔模方向上的悬垂结构加投影约束。开源实现里这类约束很少内置，多数团队靠后处理阶段做几何修复兜底。如果允许，把结构对称性体现为设计变量的对称映射，能同时收敛优化空间并保证安装一致性。

5.5 从 2D 到 3D：内存与计算时间暴涨

现象：2D 30×60 算几十秒，换成 3D 100×50×10 要跑几小时甚至几十小时，笔记本内存直接被吃满。

原因：3D 的未知量相对 2D 是数量级放大，刚度矩阵求解与灵敏度计算随之膨胀，迭代次数也要多很多才能收敛。

解决：先用对称面和周期边界条件把模型减半；用粗网格扫出大拓扑再局部细化；迭代求解器比直接求解器在 3D 模型上省大量内存。另外别在工位上跟 3D 模型硬扛，临时租个高性能实例算，拓扑优化并行效率通常不错，算完把结果拉回本地做后处理。

6. 收敛性验证与参数调优：一个可复用的实验流程

6.1 用目标函数与体积分数曲线判断收敛

跑完算例，第一件事是看两条曲线：柔度曲线和体积分数曲线。柔度曲线应连续下降，最后二十步变化小于 0.5% 才算收敛。体积分数曲线是平缓下降的阶梯，因为 OC 准则每步都在逼近目标体积。如果体积分数已经达标但柔度还在大幅下降，说明仍有大量中间密度在变，程序没到稳定状态。让程序多跑一百步，并打印每步的密度变化率，变化率低于 0.01 就是收敛信号。

6.2 参数扫描：把 rmin 和 penal 的组合做成小实验

我每次换新结构时，不直接跑精算，而是先做一个 3×3 参数网格实验：rmin 取 1.5、2.0、2.5，penal 取 2.5、3.0、3.5，其余固定。每个组合只跑几十步，比较目标函数和拓扑形态。这个小实验花不了几分钟，却能告诉你这个结构对参数敏感不敏感。如果三种 rmin 下拓扑差异很大，说明算例对滤波半径敏感，后处理要格外小心；如果差异很小，设计比较稳健，可以放心定参数。

6.3 与制造端对接的最终检查

拓扑优化的终点不是一张漂亮的密度图，而是能进入工艺链条的几何模型。导出 STL 之前做三道检查：模型水密性，看有没有开口边界；最小特征尺寸，看有没有小于工艺极限的细杆；拔模方向，按加工方向做一次投影检测。现在我的习惯是所有算例都保留原始参数脚本和导出模型文件，方便追溯。哪怕当时觉得“就一个小支架，以后不会再碰”，三个月后要解释结果来源时，能直接翻出参数和生成文件，省下无数沟通成本。这个习惯帮我避过很多次责任不清的扯皮，希望帮到你。

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