随着增材制造技术的飞速发展，复杂多孔结构的制造已成为可能，推动了轻量化与多功能集成设计的快速发展［1-2］。作为实现材料高效分布的先进设计方法，拓扑优化通过与增材制造技术融合，为高性能填充结构设计提供了全新途径。目前较为流行的连续变量方法如变密度SIMP法，通过密度惩罚策略生成拓扑结构，其产生的中间密度会导致结构边界存在模糊，难以满足增材制造对清晰几何轮廓的需求［3-4］；而离散变量方法，如SAIP、双向渐进结构优化（bi-directional evolutionary structural optimization，BESO）等，通过正则松弛或渐进优化策略，可直接获得离散化多孔构型，显著提高结构的可制造性［5-7］。近年来，Wu等［8］通过引入邻域体积分数上限控制材料分布，生成了仿骨小梁的多孔结构设计。针对多孔复合结构的跨尺度协同设计，Liu等［9］提出了基于动态聚类的多尺度拓扑优化框架。此外，针对如自支撑条件、最小尺寸限制等增材制造约束，动态体积约束、双场模型投影等技术进一步优化了填充结构的力学性能与工艺兼容性。高彤等［10］针对惯性载荷作用下的结构优化问题，提出了可变参数的材料特性的有理近似模型（rational approximation of material properties，RAMP）。高云凯等［11］针对低周疲劳优化需求，开发了基于疲劳寿命灵敏度的双向渐进结构优化方法。
在连续体拓扑优化领域，变密度SIMP法凭借其实施简洁、收敛速度快的特点得到广泛应用。焦洪宇等［12］讨论了基于变密度法的周期性拓扑优化。吴一帆等［13］提出一种灰度单元等效转换方法，用于解决SIMP法中的中间灰度问题。昌俊康等［14］提出了一种基于惩罚函数的改进密度法。廉睿超等［15］提出了分层双重惩罚方法，结合过滤和惩罚策略来抑制灰度单元。这些基于SIMP法的研究本质上仍依赖于连续变量松弛策略，难以彻底消除中间密度单元。
基于离散变量的SAIP法通过采用非0即1变量，在理论上可避免中间密度问题，但其计算效率与收敛稳定性仍有待验证［16］。现有研究缺乏对连续与离散变量方法的系统性对比分析，特别是对于填充结构的多尺度设计，两类方法的适用性差异尚不明确。填充拓扑优化生成的多孔结构能够在保证力学性能的前提下显著降低材料用量，通过引入可控孔隙，赋予材料多功能特性。
本文考虑填充结构拓扑问题，系统对比SIMP法与SAIP法在多孔结构设计中的差异。通过数值算例对比两种方法在目标函数收敛性、优化效率及可制造性等方面的表现，为多孔结构的高效设计及拓扑优化算法的选择与改进提供理论参考和依据。
1 局部体积约束
多孔结构有两个主要特征：一是其由相互交织和连接的子结构组成；二是其材料分布较非多孔结构更加均匀。多孔结构示意图如图1所示。
为形成多孔结构，将全局体积约束替换为一组局部体积约束fU，这些约束调节小邻域U内固体材料的存在，如式（1） 所示。
式中：VU 、VU ，Tot和 分别为邻域U内的材料体积、总体积和体积分数上限。
基于有限元分析，矩形设计域Ω可被分割成Ns 个矩形子域，每个子域Ωs （s=1， 2， …， Ns ）进一步离散为R 1×R 2个正方形单元，总计N=Ns ×R 1×R 2个正方形单元，设计域的子域划分如图2所示。
全局体积约束表示为
式中：V（ x ）为 x 所对应设计的体积；φ为规定的体积分数；V 0为设计域总体积。
子域局部体积约束表示为
式中：vs 、Vs 和V 0， s 分别为第s个子域的材料体积分数、材料体积和总体积； 为第s个子域内的单元相对密度。该约束明确了每个子域的体积分数应该不大于φ，因此全局体积约束也隐含其中。
2 基于SIMP法的填充结构拓扑优化
3 基于SAIP法的填充结构拓扑优化
3.1 问题列式
采用SAIP方法，在引入局部体积约束后考虑结构刚度优化的拓扑优化问题可表示为
式中： x 只可以取值为1（实心）或0（空洞）。其余符号与式（4） 中的含义相同。
3.2 离散变量灵敏度
对于离散变量方法，基于微分运算的灵敏度可以通过差分运算进行扩展［18］，即
优化问题中的目标函数可以表示为
这里将δc/δxi 记作bi，它将在3.3节关于正则松弛算法的说明中使用。
在其他单元密度不变的情况下，某个单元改变其密度时结构刚度矩阵的变化表示为Δ K，位移增量表示为δ U，因此有
联立两式得到
当xe =1时，Δ K 为-k 0和其他单元刚度矩阵取零时的组装，式（6） 和（9）分别变为
在式（10） 中，变化后的刚度矩阵 K +Δ K 是严格正定的，因此有
即 K -1（-Δ K ）的谱半径小于1，所以式（10） 中的第2式有以下展开式
在一般情况下，单元应变能应远小于总体结构应变能，则 K -1（-Δ K ）的谱半径远小于1。因此，根据式（12） 取δ U 的一阶近似
将式（13） 代入式（10） 得到
当xe =0时，Δ K 为k 0和其他单元刚度矩阵取零时的组装，对应向结构添加一个单元的情况。设想一种原结构所有位移自由度为零、新增位移自由度为1的特殊情况，有
即 K -1Δ K 的谱半径大于1，关于xe =1时的推导在xe =0时无法实施。本文采用类似SIMP法中的灵敏度计算公式［7］。
式中：x min=1×10-3，为一个极小值。式（16） 中第2式的灵敏度远小于第1式，这可以防止孤立单元在空隙区域内的积聚，更符合物理直觉。在计算过程中，对灵敏度的过滤进一步缓解了这个问题，即过滤后的灵敏度可以视为离散变量灵敏度的合理近似。
3.3 正则松弛算法
该算法的第一步是将离散变量约束等价为等式约束，并放松设计变量，即
结合式（7） ，引入拉格朗日乘子μ用于体积约束、拉格朗日乘子 σ 用于式（17） 的等式约束，构造问题（5）的拉格朗日函数。
根据密度设计变量的一阶条件
得到原始密度设计变量与对偶变量的关系
将式（20） 代入式（18） ，消去 x 、L（ x，μ， σ ）可以表示为如下的正则对偶函数。
然后，计算-P（μ， σ ）对σi 的二阶偏导。
令σi （i=1，2，…，N）非负，以保证这个对偶问题的凸性。在式（21） 中增加扰动以克服非唯一全局最优的问题［19］。
因此，有变量μ的卡鲁什-库恩-塔克（Karush-Kuhn-Tucker，KKT）条件为
式中：约束μ≥0的乘数记为κ。需要指出，体积约束是必须生效的约束，即μ≠0。
另外，变量 σ 的KKT条件可以表示为
式中：约束σi ≥0（i=1，2，…，N）的乘数记为γ i。式（24） 和式（25） 分别简化为
正则松弛算法首先初始化对偶变量μ，然后计算三次方程（27）的最大根，得到对偶变量σi （i=1，2，…，N）。然后通过关系式（20） 计算相应的设计变量。当设计变量满足收敛条件时得到问题（5）的近似解，否则使用式（26） 得到新的对偶变量μ继续迭代。本文使用目标函数值的相对变化作为收敛准则。
3.4 移动限制策略
对于SAIP方法，初始设计中所有单元被赋予的密度值为1，对应的体积分数值也为1，通过预定义体积减小因子χ逐渐减少材料的用量，即设定下一步迭代的材料用量为当前步的χ倍。设φ=0.5及体积减小200次，则体积减小因子为χ=0.51/200≈0.996 5。
通常情况下SIMP法初始设计中的所有单元会被赋予等于目标体积分数值的初始密度。本文统一了两种优化方法的初始化方式，即SIMP法也采用和SAIP法相同的移动限制策略，此时SIMP法的迭代计算时间会增加，因为该策略要求迭代逐步收敛至由初始体积分数至目标体积分数的一系列体积分数值。值得指出，在中间体积分数也收敛时，可以一次性得到一组不同体积分数对应的优化设计，如果只想获得目标体积分数的优化结果，则对中间体积分数的迭代收敛条件可以适当放松以减少计算时间。
4 数值算例和讨论
本节通过数值算例对基于SIMP法和SAIP法的填充结构拓扑优化设计进行对比研究。如图3所示，两种方法的差异主要体现在更新设计变量阶段：SIMP法基于最优准则实现材料分布的惩罚控制，而SAIP法通过正则松弛算法更新设计变量。该流程图直观呈现了两种方法在同一框架下的流程及在关键步骤的分支策略。
从问题定义与初始化开始，通过有限元分析和灵敏度计算循环迭代，将灵敏度分离至不同子区域进行优化。迭代后判断收敛性，满足条件则结束优化。在初始化中，采用移动限制策略后SIMP法的初始设计和SAIP法相同，所有单元密度值为1。
4.1 MBB梁
本节考虑一个简支MBB梁设计，MBB梁拓扑优化设计域和载荷条件如图4所示。梁在顶部边缘中间承受向下的点载荷为10 kN，考虑对称性，取MBB梁左半部分进行建模，仅施加总载荷的一半。弹性模量E=3×104 MPa，泊松比ν=0.3。设计域的尺寸L=480 mm，H=240 mm。对称处理后梁模型采用480×240的正方形网格离散，R 1=R 2=60。体积分数φ=0.5，过滤半径r min=2，即单元尺寸的2倍。
在未达到目标体积分数时，收敛条件为目标函数值的相对变化小于1×10-4，每次未收敛时，将该判定界限乘以1.5，收敛后重置。在达到目标体积分数后，收敛条件为目标函数值的相对变化小于恒定值1×10-4。
MBB梁填充结构拓扑优化结果如图5所示。从结构性能上比较，SAIP法优化结果目标函数值比SIMP法低3.1%，表明其优化结果的柔顺度更小、刚度更大；从结构形态上比较，两者都获得了填充结构拓扑设计，SIMP法的优化结果出现了灰度区域，需后处理以明确结构边界的材料分布，SAIP法能获得更加清晰的结构拓扑设计，其结构形态更优；从计算时间上比较，SIMP法用时590 s，SAIP法用时1 066 s，为SIMP法的1.8倍。MBB梁填充拓扑优化的目标函数值迭代历史如图6所示。两者的柔顺度都逐渐升高，这是因为按照移动限制策略，随着迭代次数增加，结构的材料用量逐渐减少。整个迭代过程中，同一体积分数下，SIMP法计算结果的柔顺度均高于SAIP法。
4.2 悬臂梁
对比SIMP法过滤半径为2、3.5、5的优化结果可见，较大的滤波半径可获得较低的孔隙率，尽管结构的整体外观相似，但细节有所差异。然而，SIMP法优化结果孔隙率降低后，得到的并不是如SAIP法优化结果那样低孔隙率、低柔顺度的设计结果，反而随过滤半径增大和孔隙率降低，柔顺度会增加。
5 结论
本文基于局部体积约束，将最优准则的SIMP连续变量方法和采用正则松弛算法的SAIP离散变量方法应用于填充结构拓扑优化设计。通过理论推导给出设计优化的基本列式和计算框架，并结合数值算例的分析得到以下结论：
1）SIMP法展现出计算效率优势，尤其适用于需要快速处理的大规模工程优化问题。
2）基于离散变量的SAIP法能够生成边界清晰的拓扑设计，而且目标函数值更优，在需要高精度几何表达的增材制造中更具实用性。
3）尽管两种方法均能满足多孔结构设计需求，但在实际应用中需综合考量计算资源限制及对结构设计几何精度的要求。