1. 悬臂梁柔度最小化拓扑优化
问题定义与建模对象
本例以二维悬臂梁为对象：左端固定，右端下边缘施加集中载荷 F = 1 N（向下）。设计域尺寸取 120 mm × 40 mm，离散为 120 × 40 个四边形单元，单元边长 1 mm。材料弹性模量 E₀ = 1 MPa（归一化单位），泊松比 ν = 0.3，体积分数上限 f = 0.4。
柔度最小化等价于刚度最大化，其数学形式为：
min c(ρ) = UᵀKU = Σₑ (ρₑ)ᵖ · uₑᵀk₀uₑ
s.t. KU = F
V(ρ)/V₀ = Σₑ ρₑ·vₑ / V₀ ≤ f
0 < ρ_min ≤ ρₑ ≤ 1其中 ρₑ 为单元相对密度，p = 3 为 SIMP 惩罚因子，ρ_min = 0.001 防止刚度矩阵奇异。体积约束通过拉格朗日乘子 λ 引入，最优性条件为各单元灵敏度与体积灵敏度之比相等。
SIMP 材料插值与灵敏度
SIMP 方法将单元刚度插值为 Eₑ = ρₑᵖ · E₀，使中间密度在物理上"不经济"，从而驱动设计趋向 0/1 分布。柔度对设计变量的灵敏度为：
∂c/∂ρₑ = -p · ρₑ^(p-1) · uₑᵀk₀uₑ体积灵敏度为 ∂V/∂ρₑ = vₑ。采用优化准则（OC）更新：
ρₑ_new = ρₑ · ( -∂c/∂ρₑ / (λ · ∂V/∂ρₑ) )^η其中 η = 0.5 为阻尼系数，λ 由二分法在每步迭代中调整以满足体积约束。密度变化限制在 ±0.2 以内，避免振荡。
关键实现与参数选择
有限元求解采用四节点四边形单元，每个单元 8 个自由度。全局刚度矩阵组装后施加固定端约束（左端所有节点 uₓ = uᵧ = 0），右端下角节点施加集中力。为避免棋盘格，可选用灵敏度过滤或密度过滤，过滤半径取 r_min = 2.5 倍单元边长。
OC 更新中的拉格朗日乘子 λ 通过二分法在 [10⁻⁹, 10⁹] 区间内搜索，使体积约束恰好满足。每步迭代后检查柔度相对变化，当 |cₖ - cₖ₋₁| / cₖ < 10⁻³ 且连续 5 步满足时终止。
收敛判据建议：除柔度变化外，可同时监控体积分数偏差
|V/V₀ - f| < 10⁻⁴ 和设计变量最大变化量 max|Δρₑ| < 0.01，三者同时满足时结果更可靠。
结果判读与验证
- 构型特征：最优材料分布应呈现从固定端向右端载荷点延伸的桁架状结构，上下缘形成连续传力路径，中间区域材料被删除。
- 柔度值：最终柔度应显著低于均匀密度初始设计（约下降 60%–80%），具体数值取决于网格密度和过滤半径。
- 体积分数：最终体积分数应精确等于
f = 0.4，偏差不超过10⁻⁴。 - 网格无关性：将网格加密至
240 × 80后重新求解，构型拓扑应基本一致，柔度值变化在合理范围内。
常见问题：若出现棋盘格或网格依赖，检查过滤半径是否过小；若迭代振荡不收敛，减小密度变化上限或增大阻尼系数
η；若体积约束始终无法满足，检查二分法区间和灵敏度符号是否正确。