3. 拓扑优化数学基础

拓扑优化的本质，是把“材料该放在哪里”这一工程直觉，翻译成一个可以被计算机反复求解的数学问题。要真正读懂优化结果、判断结果是否可信，必须理解五个核心要素：设计变量定义“能改什么”，目标函数定义“想变多好”，约束条件定义“不能越过的红线”，灵敏度分析回答“往哪个方向改最有效”，优化求解框架则把前四者组织成可迭代收敛的闭环。本节以增材制造中常见的悬臂梁刚度最大化问题为主线，把这五个要素串起来。

一、设计变量：把连续体离散成可决策的密度场

连续体拓扑优化最常用的做法是密度法（SIMP）。将设计域用有限元网格离散为 N 个单元，每个单元赋予一个伪密度变量 ρ_e，取值在 [ρ_min, 1] 之间：ρ_e=1 表示实体材料，ρ_e=ρ_min（常取 0.001）表示近似空洞。设计变量向量即 x = [ρ₁, ρ₂, …, ρ_N]^T。

为抑制中间密度、逼近 0-1 设计，单元弹性模量按幂律插值：

E_e(ρ_e) = E_min + ρ_e^p · (E_0 − E_min)

E_0   = 210 GPa   （实体材料弹性模量，示例取钢）
E_min = 1e-9 · E_0（空洞单元，防止刚度矩阵奇异）
p     = 3         （惩罚因子，常用 3）

惩罚因子 p 越大，中间密度单元的刚度被压得越低，优化越倾向产生黑白分明的结果；但 p 过大易导致收敛震荡。实践中常采用连续化策略：前若干次迭代用 p=1，随后逐步升至 p=3。

二、目标函数与约束条件：把工程需求写成数学式

悬臂梁刚度最大化的常用写法是最小化柔度（即最大化刚度）：

min  c(x) = U^T K U = Σ_e E_e(ρ_e) · u_e^T k_0 u_e
s.t. K U = F                     （平衡方程）
     V(x)/V_0 = Σ_e ρ_e v_e / V_0 ≤ f   （体积分数约束）
     0 < ρ_min ≤ ρ_e ≤ 1

其中 U 为全局位移向量，K 为全局刚度矩阵，F 为载荷向量，v_e 为单元体积，f 为允许的体积分数（如 0.4 表示最多用 40% 材料）。平衡方程是等式约束，体积分数是不等式约束，密度上下限是边界约束。增材制造场景下还需追加工艺约束，例如最小特征尺寸（控制悬垂与可打印性）、悬垂角约束（限制无支撑悬垂面）、连通性约束（避免封闭空腔残留粉末）。这些约束通常以附加不等式或过滤算子的形式并入同一优化模型。

三、灵敏度分析：告诉求解器往哪走

灵敏度即目标函数与约束对设计变量的偏导。对柔度最小化问题，利用伴随法可得到简洁的解析式：

∂c/∂ρ_e = −p · ρ_e^(p−1) · (E_0 − E_min) · u_e^T k_0 u_e
∂V/∂ρ_e = v_e

式中 u_e^Tk₀u_e 是单元在单位弹性模量下的应变能，物理含义清晰：应变能越大的单元对整体刚度贡献越大，越应保留材料。灵敏度为负值，说明增大该单元密度会降低柔度（提升刚度）。伴随法只需一次额外有限元求解即可获得全部 N 个偏导，计算量与设计变量数目无关，这是密度法能处理百万级单元的关键。

由于网格离散会带来棋盘格与网格依赖，需对灵敏度进行过滤（如密度过滤或灵敏度过滤），把单元灵敏度与其邻域加权平均，从而获得光滑、网格无关的设计。

四、优化求解框架：从初值到收敛的闭环

密度法通常采用移动渐近线法（MMA）或最优准则法（OC）更新设计变量。以 OC 为例，更新格式为：

ρ_e^new = clip( ρ_e · ( −∂c/∂ρ_e / (λ · v_e) )^η , ρ_min, 1 )

η  = 0.5   阻尼系数，控制步长
λ          拉格朗日乘子，由体积约束二分法确定
clip       将结果限制在 [ρ_min, 1]

每轮迭代的流程是：有限元求解位移 → 计算目标与约束 → 伴随法求灵敏度 → 灵敏度过滤 → OC/MMA 更新密度 → 检查收敛。收敛判据常用“连续两次迭代设计变量最大变化量小于 0.01”或“目标函数相对变化小于 1e-4”。

五、结果判读与常见问题

收敛后应检查三件事：一是体积分数是否恰好落在约束边界（通常应贴紧 f，说明约束起作用）；二是柔度是否单调下降并趋于平稳；三是拓扑是否出现细碎杆件、孤立材料或大面积中间密度。若出现棋盘格，说明过滤半径过小；若出现网格依赖，说明未做过滤或过滤半径随网格变化；若结果无法打印，说明缺少最小尺寸或悬垂约束。