📚 课程目录第2章拓扑优化数学模型与SIMP方法

2. 拓扑优化数学模型与SIMP方法

拓扑优化之所以能从“凭经验挖孔”变成可重复、可验证的工程方法，核心在于它被写成了一个数学规划问题。OptiStruct 的拓扑优化求解器并不直接“猜”哪里该留材料，而是把连续体离散成有限元网格，给每个单元分配一个设计变量，然后在目标与约束之间反复迭代。理解这套数学模型，是判断结果是否可信、参数是否合理的前提。

设计变量：每个单元的相对密度

在 OptiStruct 的密度法（Density Method）拓扑优化中，设计变量是每个单元的相对密度 ρ_e，取值范围通常为 0 到 1。ρ_e=1 表示该单元保留实体材料，ρ_e=0 表示该单元被挖空，中间值则是数值上的“灰度”单元。设计变量向量可写为：

设计变量：  ρ = [ρ₁, ρ₂, …, ρₙ]ᵀ ,  0 < ρ_min ≤ ρₑ ≤ 1
其中 n 为设计域单元总数，ρ_min 为下限（常取 0.001，防止刚度矩阵奇异）

这里有一个工程上极易被忽略的细节：设计变量绑定的是设计域内的单元，而非整个模型。在 OptiStruct 中通过 DESVAR 卡片或 topology 面板定义设计域，非设计域（如安装点、载荷引入区、配合面）的单元密度被固定为 1，不参与优化。因此建模阶段就要想清楚：哪些区域允许被挖空，哪些区域必须保留。若把螺栓连接区误划入设计域，优化器很可能把它掏空，导致后续无法装配。

目标函数：柔度最小化与频率最大化

最常见的拓扑优化目标是柔度最小化（Compliance Minimization），等价于刚度最大化。柔度定义为外力在位移上做的功：

目标函数：  min  C(ρ) = Uᵀ K U = Σₑ ρₑ^p · uₑᵀ k₀ uₑ
其中 U 为全局位移向量，K 为全局刚度矩阵，
uₑ 为单元位移向量，k₀ 为实体单元刚度矩阵

柔度越小，结构在给定载荷下变形越小，刚度越高。OptiStruct 默认的拓扑优化目标就是最小化加权柔度（MINIMIZE COMPLIANCE）。除柔度外，还可选择最大化第一阶固有频率、最小化质量、最小化最大应力等目标。目标函数的选择直接决定优化器“往哪个方向努力”，例如以频率为目标时，优化器会倾向于把材料布置在能提高整体刚度和模态刚度的路径上，结果形态与柔度最小化往往明显不同。

约束条件：体积分数与制造约束

拓扑优化若只给目标不给约束，最优解必然是“全部填满”，因此必须引入约束。最核心的是体积分数约束：

约束条件：  V(ρ) / V₀ = Σₑ ρₑ · vₑ / V₀ ≤ f
其中 V₀ 为设计域初始体积，vₑ 为单元体积，f 为目标体积分数（如 0.3）

体积分数 f 是拓扑优化中最重要的工程旋钮。f=0.3 意味着最终只保留设计域 30% 的材料。f 越小，结构越轻，但柔度越大、应力越高，且更容易出现细长杆件和灰度单元。实际项目中通常从 0.3 左右起步，再根据刚度、应力和工艺可行性逐步调整。除体积约束外，OptiStruct 还支持位移约束、应力约束、频率约束，以及拔模、挤压、最小/最大成员尺寸等制造约束。制造约束不是数学上的“锦上添花”，而是让优化结果能被加工出来的必要条件。

SIMP 材料插值模型与惩罚因子

如果直接把单元刚度按 ρ_e 线性缩放，即 K_e=ρ_ek₀，优化结果会充满大量中间密度单元，无法判断该留还是该删。SIMP（Solid Isotropic Material with Penalization）方法通过引入惩罚因子 p 来抑制中间密度：

SIMP 插值：  Eₑ(ρₑ) = ρₑ^p · E₀
             Kₑ(ρₑ) = ρₑ^p · k₀
其中 E₀ 为实体材料弹性模量，p 为惩罚因子，通常取 p = 3

当 p>1 时，中间密度单元（如 ρ=0.5）的刚度被大幅压低：0.5³=0.125，仅为实体刚度的 12.5%，而它却占用了 50% 的体积。这种“高成本、低刚度”的惩罚使优化器主动把中间密度推向 0 或 1，从而得到接近黑白的清晰拓扑。p 的取值需要权衡：p 太小（如 1.5）惩罚不足，灰度单元多；p 太大（如 5）虽能加速黑白化，但目标函数非凸性增强，容易陷入局部最优，且迭代震荡加剧。OptiStruct 默认 p=3，这是工程上经过大量验证的折中值。

SIMP 的惩罚机制本质上是一种“材料成本惩罚”：中间密度单元在体积上按 ρ 计费，在刚度上却只按 ρ^p 贡献。p 越大，这种“性价比”越差，优化器就越不愿意保留中间密度。

完整数学模型与求解流程

把上述要素合在一起，OptiStruct 拓扑优化的标准数学模型为：

find      ρ = [ρ₁, ρ₂, …, ρₙ]ᵀ
minimize  C(ρ) = Uᵀ K U = Σₑ ρₑ^p · uₑᵀ k₀ uₑ
subject to  K U = F
            Σₑ ρₑ · vₑ / V₀ ≤ f
            0 < ρ_min ≤ ρₑ ≤ 1 ,  e = 1, 2, …, n

求解时，OptiStruct 采用基于灵敏度的优化算法（如 MMA 或 SQP）。每次迭代先做有限元分析得到位移场，再计算目标函数和约束对设计变量的灵敏度 ∂C/∂ρₑ，然后沿可行方向更新密度场，直到收敛。灵敏度分析是这套方法高效的关键：它告诉优化器“把哪个单元的材料减少一点，对柔度影响最小”，从而避免盲目搜索。

参数设置与结果判读要点

在 OptiStruct 中，上述数学模型通过 OPTIMIZATION 面板和 DESVAR、DRESP1、DCONSTR 等卡片落地。一个典型的拓扑优化卡片设置如下：

DESVAR, 1, TOPO, 0.001, 1.0, 1.0
DVPREL1, 1, PSHELL, 1, 1, , 1.0, 0.0
DRESP1, 1, COMP, COMP, , , , , , , 1
DCONSTR, 1, 1, , 0.3
OBJECTIVE, MINIMIZE, 1

其中 DESVAR 定义设计变量上下限，DVPREL1 把设计变量与单元属性关联，DRESP1 定义柔度响应，DCONSTR 施加体积分数 0.3 的约束，OBJECTIVE 指定最小化该响应。求解完成后，OptiStruct 输出密度云图，通常以 0.5 为阈值判断保留与挖除：ρ_e≥0.5 视为实体，ρ_e≤0.5 视为空洞。但阈值并非绝对，工程上更应关注 0.3～0.7 之间的灰度单元比例——若灰度单元过多，说明惩罚不足或网格过粗，需要调整 p 值或细化网格。

惩罚因子 p 并非越大越好。p 过大时，目标函数非凸性增强，优化器可能收敛到较差的局部最优，表现为结果对初始网格或载荷步敏感、迭代曲线震荡。建议先用 p=3 跑通，再根据灰度单元比例微调。

最后需要强调：拓扑优化给出的是概念级材料分布，不是最终零件。密度云图中的细长杆件、尖角、悬空区域在制造和后续重构中往往需要修正。因此，理解设计变量、目标、约束、SIMP 插值与惩罚因子之间的数学关系，不仅是为了设置参数，更是为了在结果不合理时知道该动哪个旋钮——是放宽体积分数、调整惩罚因子，还是补充制造约束。这套数学模型是后续结构重构与工程化处理的起点。