📚 课程目录第2章拓扑优化数学基础

2. 拓扑优化数学基础

拓扑优化之所以能从“凭经验挖孔”走向“由算法自动生成结构”，根本原因在于它被严格地写成了一个数学规划问题。理解这一节，关键不是记住几个公式，而是建立一种映射思维：把工程需求翻译成设计变量、目标函数和约束条件，再交给求解器迭代。下面按“变量—目标—约束—插值—标准形式”的顺序展开，并在最后给出一个可手算的示例。

一、设计变量：优化的“可调旋钮”

在连续体拓扑优化中，最常用的做法是密度法：把设计域离散成有限元网格，每个单元赋予一个伪密度变量 ρ_e，取值范围通常为 [0,1]。ρ_e=1 表示该单元保留材料，ρ_e=0 表示该单元为空，中间值则是数值求解过程中的过渡状态。设计变量向量可写为：

x = [ρ₁, ρ₂, …, ρ_N]ᵀ,  0 ≤ ρ_e ≤ 1,  e = 1…N

其中 N 为设计域单元总数。需要注意，ρ 是“伪密度”，并不直接等于真实材料密度，它只是控制单元刚度贡献的数学量。设计变量的数量等于单元数，这也是拓扑优化计算量大的直接原因：一个十万单元的模型就有十万个设计变量。

二、目标函数：优化的“方向盘”

目标函数决定算法朝哪个方向改进。轻量化设计中最常见的是柔度最小化（等价于刚度最大化）：

min  C(x) = Uᵀ K U = Σ_e ρ_e^p · u_eᵀ k₀ u_e

式中 U 为整体位移向量，K 为整体刚度矩阵，u_e 为单元位移，k₀ 为实体单元刚度矩阵。柔度越小，结构在载荷下变形越小，刚度越高。若以减重为目标，则写成质量最小化 min M(x)=Σρ_ev_eρ₀，但纯质量最小会趋向“什么都不留”，因此必须配合刚度或位移约束。工程中常把柔度作为目标、体积作为约束，这样物理意义清晰、数值也稳定。

三、约束条件：不可逾越的边界

约束把“数学上可行”收窄到“工程上可用”。典型约束包括：

体积约束：Σρ_ev_e ≤ V*，即材料用量不超过给定上限，V* 常取设计域体积的 30%～50%。
位移约束：关键节点位移 u_i ≤ u_allow，用于控制刚度或精度。
应力约束：σ_e ≤ σ_allow，防止局部屈服，但应力约束非线性强、求解代价高。
制造约束：最小/最大成员尺寸、对称性、拔模方向等，保证结果可加工。

约束数量与类型直接决定求解难度。体积约束是单标量约束，处理最方便；应力约束是逐单元约束，数量可达数万，通常需要凝聚或伴随法处理。

四、材料插值模型：SIMP 与 RAMP

如果直接让 ρ_e 线性控制刚度，中间密度单元会大量出现，结果无法制造。为此引入材料插值模型，对中间密度进行“惩罚”，使其在经济上不划算，从而把解推向 0 或 1。

SIMP（Solid Isotropic Material with Penalization）是最常用的模型，其单元刚度插值为：

E_e(ρ_e) = E_min + ρ_e^p (E₀ − E_min)

其中 E₀ 为实体材料弹性模量，E_min 为防止刚度矩阵奇异而设的极小值（如 E₀×10⁻⁹），p 为惩罚因子，通常取 p=3。p>1 时，ρ=0.5 的单元刚度只有 0.5³=0.125 倍实体刚度，即“用了 50% 材料却只贡献 12.5% 刚度”，优化器自然不愿保留它。p 越大惩罚越强，但过大易导致收敛震荡或局部最优。

RAMP（Rational Approximation of Material Properties）采用有理函数形式：

E_e(ρ_e) = E_min + ρ_e/(1 + q(1 − ρ_e)) · (E₀ − E_min)

q 为惩罚参数，常取 q=5～10。RAMP 在 ρ 接近 0 时惩罚较温和，接近 1 时逐渐增强，因此对中间密度的“驱赶”更平滑，在应力约束或几何非线性问题中有时比 SIMP 更稳定。两者本质都是构造一个凸性可控的插值函数，让连续松弛问题的解尽量逼近 0/1 离散解。

五、优化问题标准形式与求解闭环

把上述要素组合起来，连续体拓扑优化的标准形式为：

find    x = [ρ₁, ρ₂, …, ρ_N]ᵀ
min     C(x) = UᵀKU
s.t.    V(x)/V₀ ≤ f
        KU = F
        0 ≤ ρₑ ≤ 1,  e = 1…N

求解时通常采用“有限元分析—灵敏度计算—变量更新”的迭代闭环。以柔度最小化为例，目标对设计变量的灵敏度为：

∂C/∂ρₑ = −p ρₑ^(p−1) (E₀ − E_min) uₑᵀ k₀ uₑ

该值为负，说明增加材料会降低柔度（提高刚度），符合直觉。优化准则法（OC）据此构造更新格式，例如：

ρₑ_new = ρₑ · ( (∂C/∂ρₑ) / (λ ∂V/∂ρₑ) )^η

其中 λ 为体积约束的拉格朗日乘子，η 为阻尼系数（常取 0.5），用于抑制震荡。每轮迭代后检查柔度变化量，当 |C_k+1−C_k|/C_k < 0.1% 且体积约束满足时，认为收敛。

六、一个可手算的示例

设设计域为 100 mm × 100 mm 的悬臂梁，左端固定，右端中点受集中力 F=1000 N，材料 E₀=210 GPa，ν=0.3，体积分数上限 f=0.4，惩罚因子 p=3。若某单元当前 ρ=0.5，则其刚度贡献系数为 0.5³=0.125；若优化器把它提升到 ρ=0.8，贡献系数变为 0.512，提升约 4 倍，而材料只增加 60%。这种“非线性收益”正是 SIMP 能把材料推向两端的数学动力。实际运行中，前 20～30 次迭代柔度快速下降，之后进入微调，最终在体积分数 40% 附近得到类似桁架的传力路径。

判断结果是否可信，可检查三点：体积约束是否恰好贴上限（说明材料被充分利用）；是否出现大量灰色单元（说明惩罚不足或网格过粗）；传力路径是否与载荷方向、边界条件一致。若灰色单元过多，可适当增大 p 或细化网格后重算。

本节的核心可以概括为一句话：设计变量定义“能改什么”，目标函数定义“往哪改”，约束定义“不能越界”，插值模型定义“怎么把连续解逼成可制造解”。四者齐备，拓扑优化才成为一个可求解、可复现、可判断的工程计算问题。