4. 拓扑优化基础与SIMP方法

拓扑优化要解决的核心问题是：在给定设计域、载荷与边界条件下，如何分配材料的有无，使结构在满足约束的同时让目标性能最优。连续体拓扑优化最常用的做法是把设计域离散成有限元网格，为每个单元分配一个密度变量 ρe，ρe=1 表示实体材料，ρe=0 表示孔洞。但直接使用 0/1 离散变量会导致组合爆炸，因此引入连续密度变量 0≤ρe≤1，再通过插值模型把中间密度“惩罚”到两端，最终得到接近 0/1 的清晰构型。

密度变量与SIMP插值模型

SIMP（Solid Isotropic Material with Penalization）方法的关键是把单元材料属性写成密度的幂函数。以刚度问题为例，单元弹性模量插值为：

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

其中 E0 是实体材料模量，Emin 是防止刚度矩阵奇异的小值（常取 E0×10-9），p 是惩罚因子。当 p>1 时，中间密度对应的刚度被压低，优化器会倾向于把密度推向 0 或 1。工程中 p 常取 3，也可采用连续化策略：前期 p=1 让问题更接近凸问题，随后逐步增大到 3，以改善收敛性。

热传导拓扑优化与刚度问题结构相似，只是把弹性模量替换为导热系数：

k_e(ρ_e) = k_min + ρ_e^p · (k_0 - k_min)

目标通常是最小化热柔度（即最大化散热能力），约束为材料用量。热柔度定义为 C = TTK T，其中 K 为热传导矩阵，T 为节点温度向量。热载荷可来自内部热源或边界热流，此时右端项与设计变量相关，灵敏度推导需额外考虑载荷对密度的导数。

灵敏度分析与优化求解

以最小化热柔度、体积约束为例，优化列式为：

min  C(ρ) = T^T K(ρ) T
s.t. K(ρ) T = Q
     V(ρ) = Σ ρ_e v_e ≤ V*
     0 ≤ ρ_e ≤ 1

利用伴随法可得灵敏度：

∂C/∂ρ_e = -T_e^T (∂k_e/∂ρ_e) T_e
        = -p ρ_e^(p-1) (k_0 - k_min) T_e^T k_e0 T_e

其中 Te 为单元节点温度向量，ke0 为单位导热系数对应的单元矩阵。灵敏度为负值，说明增加该单元密度会降低热柔度、提升散热性能。得到灵敏度后，常用优化准则法（OC）或移动渐近线法（MMA）更新密度。OC 更新格式可写为：

ρ_e^new = clamp( ρ_e · (-∂C/∂ρ_e / (λ v_e))^η , ρ_min, 1 )

其中 λ 为体积约束的拉格朗日乘子，η 为阻尼系数（常取 0.5），ρmin 为密度下限（如 0.001）。每步更新后需用二分法调整 λ，使体积约束满足。

滤波与投影技术

直接按单元密度更新会出现两个典型问题：网格依赖性和棋盘格。网格依赖性指细化网格后得到不同拓扑，说明问题不适定；棋盘格指相邻单元密度交替 0/1，形成类似棋盘的伪结构，其根源是有限元离散下的数值不稳定。解决方法是在灵敏度或密度上施加空间滤波。

密度滤波将单元密度替换为邻域加权平均：

ρ̃_e = Σ w(x_i,x_e) ρ_i / Σ w(x_i,x_e)
w(x_i,x_e) = max(0, r_min - ||x_i - x_e||)

rmin 为滤波半径，通常取 1.5～3 倍单元尺寸。滤波后密度过渡带变宽，可抑制棋盘格并缓解网格依赖。为进一步获得清晰边界，可叠加投影：

ρ̄_e = tanh(βη) + tanh(β(ρ̃_e - η)) / (tanh(βη) + tanh(β(1-η)))

其中 η 为阈值（常取 0.5），β 为锐度参数。优化初期 β 取小值（如 1），随后逐步增大到 16～32，使密度向 0/1 逼近，同时保持可导性。

SIMP热传导拓扑优化：从密度变量到清晰构型 设计域离散 ρ_e ∈ [0,1]，单元密度变量 SIMP插值模型 k_e = k_min + ρ_e^p (k_0 - k_min) p = 3（惩罚因子） k_min = k_0 × 10⁻⁹ 中间密度被惩罚 → 趋向0/1 连续化：p从1逐步增至3 有限元热分析 K(ρ) T = Q C = Tᵀ K T 热柔度目标 体积约束 Σρ_e v_e ≤ V* 灵敏度分析 ∂C/∂ρ_e = -p ρ_e^(p-1)· (k_0-k_min)· T_eᵀ k_e0 T_e 伴随法求解 滤波与投影：抑制棋盘格与网格依赖 密度滤波：ρ̃_e = Σw·ρ_i / Σw，w = max(0, r_min - ||x_i-x_e||) 投影：ρ̄_e = [tanh(βη)+tanh(β(ρ̃_e-η))] / [tanh(βη)+tanh(β(1-η))] r_min ≈ 1.5~3倍单元尺寸；β从1逐步增至16~32；η=0.5 优化准则法（OC）更新密度 ρ_e^new = clamp( ρ_e · (-∂C/∂ρ_e / (λ v_e))^η , ρ_min, 1 ) 二分法调整λ满足体积约束；η=0.5；ρ_min=0.001 迭代直至收敛 收敛判据 max|ρ_e^new - ρ_e| < 0.01 或目标函数变化 < 0.1%

热传导拓扑优化的基本实现

下面给出一个可直接参考的实现骨架，以二维矩形域、左侧固定温度、右侧热流输入、内部均匀热源为例。设计域尺寸 100 mm × 60 mm，网格 100×60，单元尺寸 1 mm，材料导热系数 k0=200 W/(m·K)，体积分数上限 0.4，滤波半径 rmin=3 mm，惩罚因子 p=3，投影参数 β 从 1 逐步增至 32。

% 初始化
rho = volfrac * ones(nely, nelx);
beta = 1; p = 3; rmin = 3; eta = 0.5;
loop = 0; change = 1;

while change > 0.01 && loop < 200
    loop = loop + 1;
    % 1. 密度滤波
    rho_tilde = density_filter(rho, rmin);
    % 2. 投影
    rho_bar = projection(rho_tilde, beta, eta);
    % 3. 组装热传导矩阵并求解
    K = assemble_K(rho_bar, p, k0, kmin);
    T = K \ Q;
    % 4. 计算热柔度与灵敏度
    C = T' * K * T;
    dc = sensitivity(T, rho_bar, p, k0, kmin);
    % 5. 滤波灵敏度
    dc = density_filter(dc, rmin);
    % 6. OC更新
    [rho, lambda] = oc_update(rho, dc, volfrac, v_e, eta_oc);
    % 7. 收敛判断
    change = max(abs(rho(:) - rho_old(:)));
    rho_old = rho;
    % 8. 逐步增大beta
    if mod(loop, 20) == 0 && beta < 32
        beta = beta * 2;
    end
end

运行上述流程后，可观察到的典型现象是：前 20 步 β 较小时，密度场呈现较宽的灰色过渡带，结构轮廓模糊；随着 β 增大，过渡带收窄，最终形成由实体单元和孔洞单元组成的清晰散热通道。判断结果是否合理，可检查三项指标：体积分数是否稳定在 0.4 附近；热柔度是否单调下降并趋于平稳；密度场中是否出现大面积棋盘格或孤立单元。若棋盘格明显，应增大 rmin 或检查滤波实现是否正确；若网格依赖明显，可在同一体积分数下对比 50×30 与 100×60 两套网格的拓扑是否一致。

实践提示：热传导拓扑优化中，若热源或热流边界与设计变量相关，灵敏度需增加载荷导数项。对于固定热流边界，该项为零；对于内部热源，若热源强度不随密度变化，也可忽略。但若采用“热源密度”模型，则必须补充 ∂Q/∂ρe 项，否则灵敏度不完整，优化结果会偏离预期。
常见误区：惩罚因子 p 并非越大越好。p 过大（如 p=5）会使中间密度刚度过低，优化器过早陷入局部极小，得到不自然的细长结构；p 过小则中间密度残留多，边界模糊。建议从 p=1 开始，每 20～30 步增加 0.5，直至 p=3，兼顾收敛性与清晰度。

总结来说，SIMP 方法通过密度变量、幂函数插值、灵敏度分析、滤波与投影四个环节，把连续体拓扑优化转化为可求解的非线性规划问题。热传导问题与刚度问题在数学结构上高度一致，只需替换材料属性矩阵和目标函数定义。掌握这套流程后，可进一步扩展到热-力耦合、多材料以及考虑对流边界的热拓扑优化场景。

← 上一章 📖 返回目录 下一章 →