📚 课程目录第3章连续体拓扑优化的数学表达

3. 连续体拓扑优化的数学表达

连续体拓扑优化要解决的核心问题，可以概括为一句话：在给定的设计域内，决定材料“哪里留、哪里去”，使结构在满足材料用量限制的前提下刚度最大。这句话听起来直观，但要交给计算机求解，必须把它翻译成一套严格的数学表达——设计变量如何定义、材料属性如何随设计变量变化、用什么指标衡量刚度、体积约束如何写、最终优化问题如何组织。本节就围绕这五个要素，建立可计算的问题描述。

一、设计域与设计变量：把“留不留材料”变成连续可导的变量

连续体拓扑优化首先需要一个设计域（design domain）Ω，即允许布置材料的空间范围。它通常由结构的受力条件、装配边界和工艺限制共同确定。设计域内还要区分两类区域：一类是不可设计域，例如必须保留的连接孔、安装面、载荷作用区；另一类是可设计域，材料可以自由分布。

最常用的离散方式是有限元网格。把设计域划分为 N 个单元，每个单元赋予一个设计变量 ρ_e。在经典的SIMP（Solid Isotropic Material with Penalization）方法中，ρ_e 取连续值：

ρ_e ∈ [ρ_min, 1]，通常 ρ_min = 0.001

ρ_e ≈ 1  → 该单元为实体材料
ρ_e ≈ 0  → 该单元近似为空（孔洞）

这里的关键在于：如果直接用 0/1 离散变量，问题会变成组合爆炸的整数规划，几乎无法求解。引入连续变量后，优化问题变得可导，可以用梯度类算法高效求解。ρ_min 取一个很小的正数而不是 0，是为了避免刚度矩阵奇异。

二、材料插值：让中间密度“不划算”

连续变量带来一个新问题：ρ_e = 0.5 这种中间密度在物理上意味着什么？如果直接按线性插值，中间密度单元会获得与其体积不成比例的高刚度，优化结果就会充满灰色区域，无法得到清晰的黑白拓扑。SIMP 模型用幂函数惩罚来解决这个问题：

E_e(ρ_e) = ρ_e^p · E_0

E_0：实体材料的弹性模量
p ：惩罚因子，通常取 p = 3

当 p > 1 时，中间密度单元的刚度被压低：ρ_e = 0.5、p = 3 时，E_e 仅为实体材料的 12.5%，而它却占用了 50% 的体积。这样，优化器会主动把材料推向 0 或 1 两端，最终得到接近黑白的清晰拓扑。p 越大，惩罚越强、结果越“黑白”，但过大容易导致收敛困难或陷入局部最优，工程中常从 p = 1 逐步增加到 p = 3，即连续化策略。

三、目标函数与体积约束：刚度最大化的标准写法

“刚度最大”在数学上通常等价于柔度最小。柔度定义为外力做功：

c(ρ) = U^T K U = Σ_e u_e^T k_e(ρ_e) u_e

U ：整体位移向量
K ：整体刚度矩阵
u_e、k_e：单元位移与单元刚度矩阵

柔度越小，结构在同样载荷下变形越小，刚度越大。体积约束则限制材料用量：

V(ρ) = Σ_e ρ_e · v_e ≤ f · V_0

v_e ：单元体积
V_0 ：设计域总体积
f  ：体积分数上限，如 f = 0.4

把目标、约束和设计变量组织起来，就得到连续体拓扑优化的标准数学表达：

min_ρ   c(ρ) = U^T K U
s.t.    K U = F                    （平衡方程）
        Σ_e ρ_e v_e ≤ f · V_0      （体积约束）
        ρ_min ≤ ρ_e ≤ 1            （设计变量上下限）

其中 K U = F 是有限元平衡方程，它作为等式约束保证结构在物理上成立。这个问题的本质是：在体积预算内，寻找使柔度最小的材料分布。

四、灵敏度分析：告诉优化器往哪走

有了目标函数和约束，还需要知道“改变某个单元的密度，柔度会怎么变”，这就是灵敏度。对 SIMP 模型，柔度对设计变量的导数有简洁形式：

∂c/∂ρ_e = -p · ρ_e^(p-1) · u_e^T k_e^0 u_e

k_e^0 ：实体单元的刚度矩阵

负号表示增加材料会降低柔度（提高刚度）。体积约束的灵敏度则是 ∂V/∂ρ_e = v_e。有了这两组导数，就可以用OC（Optimality Criteria，最优准则）或MMA（Method of Moving Asymptotes，移动渐近线）等算法迭代更新设计变量。OC 更新格式通常写作：

ρ_e^new = ρ_e · ( -∂c/∂ρ_e / (λ · v_e) )^η

λ ：体积约束的拉格朗日乘子
η ：阻尼系数，常取 0.5，用于稳定迭代

λ 通过二分法调整，使更新后的体积恰好满足约束。迭代直到设计变量变化量或目标函数变化量小于收敛容差（如 0.01）为止。

五、一个可复现的示例设定

下面给出一个二维悬臂梁拓扑优化的示例参数，用于说明上述数学表达如何落地。这些数值属于教学示例设定，不是实测结果。

项目	取值	说明
设计域	80 × 40 mm	矩形区域，左端固定
网格	80 × 40 单元	单元尺寸 1 mm
材料	E₀ = 210 GPa，ν = 0.3	钢的典型值
载荷	F = 100 N，右端中点向下	集中力
体积分数	f = 0.4	最多使用 40% 材料
惩罚因子	p = 3	SIMP 惩罚
过滤半径	r_min = 2.5 mm	抑制棋盘格
收敛容差	0.01	设计变量最大变化量

按上述设定迭代，通常几十到上百步后，柔度会从初始均匀分布的较高值下降到稳定值，材料逐渐向主传力路径聚集，形成类似桁架的清晰构型。判断结果是否可用，主要看三点：拓扑是否黑白分明、体积是否恰好满足约束、柔度是否已收敛平稳。

SIMP 的惩罚因子 p 与过滤半径 r_min 是两个最需要关注的参数：p 决定结果的“黑白程度”，r_min 决定最小结构尺寸并抑制数值不稳定。二者需要配合调整，不能只改其一。

六、完整问题描述的逻辑链

把本节内容串起来看，连续体拓扑优化的数学表达遵循一条清晰的逻辑链：设计域定义“在哪里优化”，设计变量定义“优化什么”，材料插值建立变量与物理属性之间的桥梁，目标函数回答“优化成什么样”，体积约束限定“能用多少材料”，灵敏度与更新算法则给出“怎么一步步算出来”。这六个环节缺一不可，共同构成一个可计算、可迭代、可收敛的优化问题。

理解这套数学表达，是后续使用 OptiStruct、ANSYS、COMSOL 或自编 MATLAB/Python 程序进行拓扑优化的基础。无论工具如何变化，设计变量、插值模型、目标、约束和灵敏度这五个要素始终是问题描述的核心。