3. 有限元分析与拓扑优化的耦合

拓扑优化并不是独立于有限元分析之外的一套算法，它本质上是把有限元求解器嵌入到一个迭代优化框架中，反复调用“分析—求导—更新”三个环节。理解这条耦合链路，是判断优化结果是否可信、迭代为何不收敛、灵敏度为何异常的基础。本节按数据在链路中的流动顺序展开：结构离散化如何产生设计变量与单元刚度，刚度矩阵如何组装并求解位移，灵敏度如何从位移场反推，伴随法如何把成千上万个设计变量的求导压缩为一次附加求解，最后串成完整的优化循环。

结构离散化：设计变量与有限元网格的绑定

拓扑优化的第一步是把连续设计域离散为有限元网格，同时把每个单元的相对密度 ρ_e（0 到 1 之间）作为设计变量。以二维悬臂梁为例，设计域取 120 mm × 40 mm，用 120×40 的四节点平面应力单元划分，共 4800 个单元，即 4800 个设计变量。每个单元采用 SIMP 插值模型：

E_e = E_min + ρ_e^p · (E_0 - E_min)
p = 3，E_0 = 210 GPa（钢），E_min = 1e-3 · E_0

这里 ρ_e=1 表示实体材料，ρ_e=0 表示空洞，惩罚因子 p=3 使中间密度在刚度上“不划算”，从而把结果推向 0/1 分布。离散化阶段必须同时确定：单元类型（平面应力/壳/实体）、网格密度、设计域与非设计域（如加载点、约束区通常锁定为实体）、以及过滤半径 r_min（常取 1.5～3 倍单元尺寸，本例取 3 mm）。网格越细，可分辨的最小结构特征越小，但设计变量与计算量同步增长。

刚度矩阵组装与位移求解

每个单元的刚度矩阵由形函数与材料矩阵积分得到，单元刚度写为 k_e = ∫ B^T D(ρ_e) B dΩ，其中 D 随 ρ_e 变化。组装时按单元自由度编号把 k_e 累加到全局刚度矩阵 K 的对应位置：

K = Σ_e  A_e^T · k_e(ρ_e) · A_e
K · U = F

A_e 是单元到全局自由度的映射矩阵。由于 ρ_e 在迭代中不断变化，K 每次迭代都要重新组装并重新分解。对 4800 单元、约 5000 个自由度的模型，直接求解尚可；当自由度达到百万级时，必须采用稀疏存储与迭代求解器（如 PCG），否则单次分析就会成为瓶颈。求解得到位移场 U 后，目标函数柔度 C = F^TU = U^TKU 即可计算，它等价于结构总应变能，最小化柔度就是最大化刚度。

灵敏度分析：目标对设计变量的导数

优化器需要知道每个 ρ_e 变化时目标 C 如何变化，即灵敏度 ∂C/∂ρ_e。对柔度问题，利用平衡方程 K U = F 对 ρ_e 求导，可得解析表达式：

∂C/∂ρ_e = -U_e^T · (∂k_e/∂ρ_e) · U_e
∂k_e/∂ρ_e = p · ρ_e^(p-1) · (E_0 - E_min) · k_e^0

其中 U_e 是该单元的自由度位移向量，k_e⁰ 是单位弹性模量下的单元刚度。这个式子说明：灵敏度只依赖单元自身的位移与刚度，不需要对整个 K 求逆，计算量很小。但直接使用会带来两个问题——棋盘格（相邻单元密度交替 0/1）和网格依赖性（网格越细结构越细碎）。因此要对灵敏度做过滤：

∂C~/∂ρ_e = Σ_f w_f · ρ_f · (∂C/∂ρ_f) / (ρ_e · Σ_f w_f)
w_f = max(0, r_min - dist(e,f))

过滤半径 r_min 决定了结构的最小特征尺寸，是控制可制造性的关键参数。过滤后灵敏度再交给优化器（如 OC 准则法或 MMA）更新密度场。

伴随法：把大规模求导压缩为一次附加求解

当目标或约束不是简单的柔度，而是应力、频率、位移等一般响应 g(ρ, U) 时，直接对 g 求全导数会遇到 ∂U/∂ρ_e 这一项，而 U 的维度等于自由度数，逐个求导不可行。伴随法引入拉格朗日乘子 λ，构造：

L = g(ρ, U) + λ^T (K U - F)
dL/dρ_e = ∂g/∂ρ_e + λ^T (∂K/∂ρ_e) U
K^T λ = -∂g/∂U   （伴随方程）

由于 K 对称，伴随方程与原平衡方程共用同一个刚度矩阵，只需一次额外求解即可得到 λ，随后所有设计变量的导数都能用 λ^T(∂K/∂ρ_e)U 批量算出。这就是伴随法的核心价值：把“对每个设计变量各求一次导”变成“解一次伴随方程，再统一回代”。在应力约束、频率约束、多工况问题中，伴随法几乎是唯一可行的灵敏度计算方式。

优化循环流程与收敛判断

把上述环节串起来，就得到标准的拓扑优化循环。以最小柔度、体积分数约束为例，目标体积分数 f=0.4，即最终保留 40% 的材料：

1. 初始化：ρ_e=f，设置过滤半径、惩罚因子、迭代上限（如 100 步）。
2. 有限元分析：组装 K(ρ)，求解 KU=F，得到 U 与柔度 C。
3. 灵敏度分析：计算 ∂C/∂ρ_e，做灵敏度过滤。
4. 优化更新：用 OC 或 MMA 更新 ρ_e，并满足体积约束 Σρ_eV_e=f·V₀。
5. 收敛判断：若密度最大变化量 max|ρ_e^k+1-ρ_e^k| < 0.01 且连续若干步目标变化小于 0.1%，则停止；否则返回第 2 步。

实际迭代中常见的现象是：前 10～20 步柔度快速下降，密度场从均匀灰色逐渐分化出传力路径；30 步后结构拓扑基本稳定，仅边界微调；若出现目标震荡或体积约束反复越界，通常要检查过滤半径是否过小、惩罚因子是否过大、或优化器步长是否激进。收敛后还需做后处理：以 ρ_e≥0.5 提取实体区域，重新划分光滑网格，再做一次验证性有限元分析，确认柔度、应力与制造约束满足要求。