填充结构拓扑优化的连续与离散变量法

Continuous and discrete variable methods for topology optimization of infill structures

摘要

针对填充结构拓扑优化问题，对比研究了连续与离散变量设计方法的性能差异。基于局部体积约束，构建了基于最优准则的固体各向同性材料惩罚模型（solid isotropic material with penalization，SIMP）连续变量优化模型和基于正则松弛算法的序列近似整数规划（sequential approximate integer programming，SAIP）离散变量优化模型。首先，通过理论推导介绍了两类方法应用于填充结构拓扑优化的基本列式。然后，利用数值算例从优化结果的性能、结构拓扑设计的清晰程度和计算效率三方面开展对比研究。最后，通过引入相同初始化策略、移动限制策略与收敛准则，消除算法外因干扰。结果表明，基于SIMP法的优化耗时约为基于SAIP法的50%，但其生成的拓扑设计结构边界存在灰度区域；SAIP法能够获得边界清晰、结构形态更优的拓扑设计，其优化结果目标函数值降低了4.2%。研究结果为填充结构拓扑设计优化提供了理论参考和依据。

Abstract

The performance differences between continuous and discrete variable design methods for topology optimization of infill structures were investigated. Based on local volume constraints， two optimization models were established： a continuous variable model using the SIMP method with the optimality criterion and a discrete variable model using the SAIP method with the canonical relaxation algorithm. Firstly， the fundamental formulations of these methods for topology optimization of infill structures were derived theoretically. Subsequently， numerical examples were utilized to conduct comparative studies from three aspects： optimization performance， clarity of structural topology design， and computational efficiency. Finally，by implementing identical initialization strategies， move limit strategies， and convergence criteria， the external interference in the algorithm was eliminated. The results demonstrate that the optimization time using the SIMP method is approximately 50% of that using the SAIP method； however， the topological designs generated by SIMP exhibit grayscale regions at structural boundaries. In contrast， the SAIP method yields topological designs with clearer boundaries and superior structural configurations， achieving about a 4% reduction in the objective function value. A theoretical reference and basis for the topology optimization design of infill structures is provided.

Graphical abstract

关键词

填充结构 / 拓扑优化 / 固体各向同性材料惩罚模型 / 序列近似整数规划 / 离散变量法 / 正则松弛算法

Key words

infill structures / topology optimization / SIMP / SAIP / discrete variable method / canonical relaxation algorithm

随着增材制造技术的飞速发展，复杂多孔结构的制造已成为可能，推动了轻量化与多功能集成设计的快速发展^［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所示。

为形成多孔结构，将全局体积约束替换为一组局部体积约束f_U，这些约束调节小邻域U内固体材料的存在，如式（1）所示。

f_U=\frac{V_U}{V_{U,\mathrm{T}\mathrm{o}\mathrm{t}}}-\phi \le \mathrm{0}（1）

式中：V_U 、V_U，Tot和\phi分别为邻域U内的材料体积、总体积和体积分数上限。

基于有限元分析，矩形设计域Ω可被分割成N_s 个矩形子域，每个子域Ω_s （s=1， 2， …， N_s ）进一步离散为R₁×R₂个正方形单元，总计N=N_s ×R₁×R₂个正方形单元，设计域的子域划分如图2所示。

全局体积约束表示为

\left\{\begin{array}{l}\mathit{x}=\left(x_{\mathrm{1}},x_{\mathrm{2}},\dots ,x_N\right)\\ V=V\left(\mathit{x}\right)-\phi V_{\mathrm{0}}\le \mathrm{0}\end{array}\right.（2）

式中：V（ x ）为 x 所对应设计的体积；φ为规定的体积分数；V₀为设计域总体积。

子域局部体积约束表示为

\left\{\begin{array}{cc}\begin{array}{l}{f}_s=v_s-\phi \le \mathrm{0}\\ v_s=\frac{V_s}{V_{\mathrm{0},s}}=\frac{\sum _{e\in s}{\rho }_e}{R_{\mathrm{1}}{R}_{\mathrm{2}}}\end{array}& ,s=\mathrm{1,2},\dots ,N_s\end{array}\right.（3）

式中：v_s 、V_s 和V_0，s 分别为第s个子域的材料体积分数、材料体积和总体积；{\rho }_e为第s个子域内的单元相对密度。该约束明确了每个子域的体积分数应该不大于φ，因此全局体积约束也隐含其中。

2 基于SIMP法的填充结构拓扑优化

变密度SIMP法作为一种普遍流行的连续体拓扑优化方法，其核心思想是通过引入惩罚因子p，将密度ρ相关的单元刚度定义为ρ^pk₀，其中，k₀是单元相对密度为1时的单元刚度，将中间密度推向0或1，获得优化的材料分布^［17］。采用最优准则法求解^［4］，以结构柔顺度最小化（刚度优化）为例，引入局部体积约束的拓扑优化问题，如式（4）所示^［4］。

\begin{array}{l}\begin{array}{cc}\mathrm{f}\mathrm{i}\mathrm{n}\mathrm{d}& \mathit{x}=\left(x_{\mathrm{1}},x_{\mathrm{2}},\dots ,x_N\right)\end{array}\\ \begin{array}{ll}\underset{\mathit{x}}{\mathrm{m}\mathrm{i}\mathrm{n}}& c\left(\mathit{x}\right)={\mathit{F}}^{\mathrm{T}}\mathit{U}\\ \mathrm{s}.\mathrm{t}.& \begin{array}{l}\mathit{F}=\mathit{K}\mathit{U}\\ \begin{array}{cc}{f}_s\le \mathrm{0},& s=\mathrm{1,2},\dots ,N_s\end{array}\\ \begin{array}{cc}\mathrm{0}\le x_i\le \mathrm{1},& i=\mathrm{1,2},\dots ,N\end{array}\end{array}\end{array}\end{array}（4）

式中：c为结构柔顺度； F 、 K 和 U 分别为载荷向量、刚度矩阵和位移向量。设计域被离散为N个单元。设计变量 x 为单元相对密度，其取值范围在0和1之间；局部体积约束f_s 在式（3）中给出。

3 基于SAIP法的填充结构拓扑优化

3.1 问题列式

采用SAIP方法，在引入局部体积约束后考虑结构刚度优化的拓扑优化问题可表示为

\begin{array}{l}\begin{array}{cc}\mathrm{f}\mathrm{i}\mathrm{n}\mathrm{d}& \mathit{x}=\left(x_{\mathrm{1}},x_{\mathrm{2}},\dots ,x_N\right)\end{array}\\ \begin{array}{ll}\underset{\mathit{x}}{\mathrm{m}\mathrm{i}\mathrm{n}}& c\left(\mathit{x}\right)={\mathit{F}}^{\mathrm{T}}\mathit{U}\\ \mathrm{s}.\mathrm{t}.& \begin{array}{l}\mathit{F}=\mathit{K}\mathit{U}\\ \begin{array}{cc}{f}_s\le \mathrm{0},& s=\mathrm{1,2},\dots ,N_s\end{array}\\ \begin{array}{cc}{x}_i\in \left\{\mathrm{0,1}\right\},& i=\mathrm{1,2},\dots ,N\end{array}\end{array}\end{array}\end{array}（5）

式中： x 只可以取值为1（实心）或0（空洞）。其余符号与式（4）中的含义相同。

3.2 离散变量灵敏度

对于离散变量方法，基于微分运算的灵敏度可以通过差分运算进行扩展^［18］，即

\begin{array}{l}\frac{\delta c}{\delta x_e}=\frac{c\left(x_e=\mathrm{1}\right)-c\left(x_e=\mathrm{0}\right)}{\mathrm{1}-\mathrm{0}}\\ =c\left(x_e=\mathrm{1}\right)-c\left(x_e=\mathrm{0}\right)\end{array}（6）

优化问题中的目标函数可以表示为

c={\displaystyle \sum _{i=\mathrm{1}}^N}\frac{\delta c}{\delta x_i}{x}_i={\displaystyle \sum _{i=\mathrm{1}}^N}{b}_i{x}_i（7）

这里将δc/δx_i 记作b_i，它将在3.3节关于正则松弛算法的说明中使用。

在其他单元密度不变的情况下，某个单元改变其密度时结构刚度矩阵的变化表示为Δ K，位移增量表示为δU，因此有

\left\{\begin{array}{l}\mathit{F}=\mathit{K}\mathit{U}\\ \mathit{F}=\left(\mathit{K}+\mathrm{\Delta }\mathit{K}\right)\left(\mathit{U}+\delta \mathit{U}\right)\end{array}\right.（8）

联立两式得到

\left(\mathit{K}+\mathrm{\Delta }\mathit{K}\right)\delta \mathit{U}=-\mathrm{\Delta }\mathit{K}\mathit{U}（9）

当x_e =1时，Δ K 为-k₀和其他单元刚度矩阵取零时的组装，式（6）和（9）分别变为

\left\{\begin{array}{l}\frac{\delta c}{\delta x_e}=\frac{\mathrm{1}}{\mathrm{2}}{\mathit{U}}^{\mathrm{T}}\mathit{K}\mathit{U}-\frac{\mathrm{1}}{\mathrm{2}}{\left(\mathit{U}+\delta \mathit{U}\right)}^{\mathrm{T}}\left(\mathit{K}+\mathrm{\Delta }\mathit{K}\right)\left(\mathit{U}+\delta \mathit{U}\right)\\ \delta \mathit{U}=-{\left(\mathrm{I}+{\mathit{K}}^{-\mathrm{1}}\mathrm{\Delta }\mathit{K}\right)}^{-\mathrm{1}}{\mathit{K}}^{-\mathrm{1}}\mathrm{\Delta }\mathit{K}\mathit{U}\end{array}\right.（10）

在式（10）中，变化后的刚度矩阵 K +Δ K 是严格正定的，因此有

（11）

即 K^-1（-Δ K ）的谱半径小于1，所以式（10）中的第2式有以下展开式

\delta \mathit{U}=\left(\mathrm{I}-{\mathit{K}}^{-\mathrm{1}}\mathrm{\Delta }\mathit{K}+{\left({\mathit{K}}^{-\mathrm{1}}\mathrm{\Delta }\mathit{K}\right)}^{\mathrm{2}}-\cdots \right){\mathit{K}}^{-\mathrm{1}}\left(-\mathrm{\Delta }\mathit{K}\right)\mathit{U}（12）

在一般情况下，单元应变能应远小于总体结构应变能，则 K^-1（-Δ K ）的谱半径远小于1。因此，根据式（12）取δU 的一阶近似

\delta \mathit{U}\approx {\mathit{K}}^{-\mathrm{1}}\left(-\mathrm{\Delta }\mathit{K}\right)\mathit{U}（13）

将式（13）代入式（10）得到

\frac{\delta c}{\delta x_e}\approx \frac{\mathrm{1}}{\mathrm{2}}{\mathit{U}}^{\mathrm{T}}\mathrm{\Delta }\mathit{K}\mathit{U}（14）

当x_e =0时，Δ K 为k₀和其他单元刚度矩阵取零时的组装，对应向结构添加一个单元的情况。设想一种原结构所有位移自由度为零、新增位移自由度为1的特殊情况，有

（15）

即 K^-1Δ K 的谱半径大于1，关于x_e =1时的推导在x_e =0时无法实施。本文采用类似SIMP法中的灵敏度计算公式^［7］。

\frac{\delta c}{\delta x_e}\approx \left\{\begin{array}{lll}-\frac{\mathrm{1}}{\mathrm{2}}{\mathit{U}}^{\mathrm{T}}{k}_{\mathrm{0}}\mathit{U},& & x_e=\mathrm{1}\\ -x_{\mathrm{m}\mathrm{i}\mathrm{n}}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{U}}^{\mathrm{T}}{k}_{\mathrm{0}}\mathit{U},& & x_e=x_{\mathrm{m}\mathrm{i}\mathrm{n}}\end{array}\right.（16）

式中：x_min=1×10^-3，为一个极小值。式（16）中第2式的灵敏度远小于第1式，这可以防止孤立单元在空隙区域内的积聚，更符合物理直觉。在计算过程中，对灵敏度的过滤进一步缓解了这个问题，即过滤后的灵敏度可以视为离散变量灵敏度的合理近似。

3.3 正则松弛算法

该算法的第一步是将离散变量约束等价为等式约束，并放松设计变量，即

\begin{array}{ccc}{x}_i\left(x_i-\mathrm{1}\right)=\mathrm{0},& \mathrm{0}\le x_i\le \mathrm{1},& i=\mathrm{1,2},\dots ,N\end{array}（17）

结合式（7），引入拉格朗日乘子μ用于体积约束、拉格朗日乘子 σ 用于式（17）的等式约束，构造问题（5）的拉格朗日函数。

\begin{array}{l}L\left(\mathit{x},\mu ,\mathit{\sigma }\right)={\displaystyle \sum _{i=\mathrm{1}}^N}{b}_i{x}_i+\mu \left({\displaystyle \sum _{i=\mathrm{1}}^N}{v}_i{x}_i-\phi V_{\mathrm{0}}\right)\\ +{\displaystyle \sum _{i=\mathrm{1}}^N}{\sigma }_i{x}_i\left(x_i-\mathrm{1}\right)\end{array}（18）

根据密度设计变量的一阶条件

\begin{array}{cc}\frac{\partial L\left(\mathit{x},\mu ,\mathit{\sigma }\right)}{\partial x_i}=\mathrm{0},& i=\mathrm{1,2},\dots ,N\end{array}（19）

得到原始密度设计变量与对偶变量的关系

\begin{array}{cc}{x}_i=\frac{{\sigma }_i-\mu v_i-b_i}{\mathrm{2}{\sigma }_i},& i=\mathrm{1,2},\dots ,N\end{array}（20）

将式（20）代入式（18），消去 x 、L（ x，μ， σ ）可以表示为如下的正则对偶函数。

P\left(\mu ,\mathit{\sigma }\right)=-\mu \phi V_{\mathrm{0}}-{\displaystyle \sum _{i=\mathrm{1}}^N}\frac{{\left(b_i+\mu v_i-{\sigma }_i\right)}^{\mathrm{2}}}{\mathrm{4}{\sigma }_i}（21）

然后，计算-P（μ， σ ）对σ_i 的二阶偏导。

\frac{{\partial }^{\mathrm{2}}\left(-P\left(\mu ,\mathit{\sigma }\right)\right)}{\partial {\sigma }_i^{\mathrm{2}}}=\frac{{\left(b_i+\mu v_i\right)}^{\mathrm{2}}}{\mathrm{2}{\sigma }_i^{\mathrm{3}}}（22）

令σ_i （i=1，2，…，N）非负，以保证这个对偶问题的凸性。在式（21）中增加扰动以克服非唯一全局最优的问题^［19］。

P_{\beta }\left(\mu ,\mathit{\sigma }\right)=-\mu \phi V_{\mathrm{0}}-{\displaystyle \sum _{i=\mathrm{1}}^N}\left(\frac{{\left(b_i+\mu v_i-{\sigma }_i\right)}^{\mathrm{2}}}{\mathrm{4}{\sigma }_i}+\frac{{\sigma }_i^{\mathrm{2}}}{\mathrm{4}\beta }\right)（23）

因此，有变量μ的卡鲁什-库恩-塔克（Karush-Kuhn-Tucker，KKT）条件为

\left\{\begin{array}{l}\frac{\partial P_{\beta }\left(\mu ,\mathit{\sigma }\right)}{\partial \mu }-\mathrm{\kappa }=\mathrm{0}\\ \mu \ge \mathrm{0}\\ \mu \mathrm{\kappa }=\mathrm{0}\\ \mathrm{\kappa }\ge \mathrm{0}\end{array}\right.（24）

式中：约束μ≥0的乘数记为κ。需要指出，体积约束是必须生效的约束，即μ≠0。

另外，变量 σ 的KKT条件可以表示为

\left\{\begin{array}{l}\frac{\partial P_{\beta }\left(\mu ,\mathit{\sigma }\right)}{\partial {\sigma }_i}-{\mathrm{\gamma }}_i=\mathrm{0}\\ {\sigma }_i\ge \mathrm{0}\\ {\sigma }_i{\mathrm{\gamma }}_i=\mathrm{0}\\ {\mathrm{\gamma }}_i\ge \mathrm{0}\end{array}\right.（25）

式中：约束σ_i ≥0（i=1，2，…，N）的乘数记为γ _i。式（24）和式（25）分别简化为

\mu =\frac{{\displaystyle \sum _{i=\mathrm{1}}^N}\frac{\left({\sigma }_i-b_i\right)v_i}{\mathrm{2}{\sigma }_i}-\phi V_{\mathrm{0}}}{{\displaystyle \sum _{i=\mathrm{1}}^N}\frac{v_i^{\mathrm{2}}}{\mathrm{2}{\sigma }_i}}（26）

{\sigma }_i^{\mathrm{3}}+\frac{\beta }{\mathrm{2}}{\sigma }_i^{\mathrm{2}}-\frac{\beta }{\mathrm{2}}{\left(b_i+\mu v_i\right)}^{\mathrm{2}}=\mathrm{0}（27）

正则松弛算法首先初始化对偶变量μ，然后计算三次方程（27）的最大根，得到对偶变量σ_i （i=1，2，…，N）。然后通过关系式（20）计算相应的设计变量。当设计变量满足收敛条件时得到问题（5）的近似解，否则使用式（26）得到新的对偶变量μ继续迭代。本文使用目标函数值的相对变化作为收敛准则。

3.4 移动限制策略

对于SAIP方法，初始设计中所有单元被赋予的密度值为1，对应的体积分数值也为1，通过预定义体积减小因子χ逐渐减少材料的用量，即设定下一步迭代的材料用量为当前步的χ倍。设φ=0.5及体积减小200次，则体积减小因子为χ=0.5^1/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×10⁴ MPa，泊松比ν=0.3。设计域的尺寸L=480 mm，H=240 mm。对称处理后梁模型采用480×240的正方形网格离散，R₁=R₂=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 悬臂梁

图5中优化后的结构左上角有一处凸起结构，它对结构整体的刚度贡献很小，但占用了较多的材料用量，这是由于局部体积约束给各个子域规定了体积分数。为消减这种结构，本小节考虑初始设计包含椭圆边界的悬臂梁结构，其拓扑优化设计域和载荷条件如图7所示。设计域尺寸L=480 mm、H=240 mm；初始设计包括椭圆短轴左侧的矩形部分和右半椭圆，短轴长度为H、半长轴长度W=180 mm。梁在右边缘中点承受向下的集中载荷 F =10 kN。弹性模量E=3×10⁴ MPa，泊松比ν=0.3。该梁模型的离散采用480×240的正方形网格，R₁=R₂=60。体积分数设置为φ=0.5，过滤半径r_min=2。

悬臂梁填充结构拓扑优化结果如图8所示。SAIP法优化结果目标函数值比SIMP法低4.2%。SAIP法用时1 266 s，为SIMP法用时629 s的2.01倍。悬臂梁填充拓扑优化的目标函数值迭代历史如图9所示。

另外，图8中的结果显示，SAIP法优化结果的孔隙率低于SIMP法，这一结论同样适用于图5中的结果。为探讨孔隙率对优化结果的影响，利用密度法调整过滤半径来获取不同孔隙率^［3］，基于SIMP法的悬臂梁填充拓扑优化结果如图10所示。

对比SIMP法过滤半径为2、3.5、5的优化结果可见，较大的滤波半径可获得较低的孔隙率，尽管结构的整体外观相似，但细节有所差异。然而，SIMP法优化结果孔隙率降低后，得到的并不是如SAIP法优化结果那样低孔隙率、低柔顺度的设计结果，反而随过滤半径增大和孔隙率降低，柔顺度会增加。

5 结论

本文基于局部体积约束，将最优准则的SIMP连续变量方法和采用正则松弛算法的SAIP离散变量方法应用于填充结构拓扑优化设计。通过理论推导给出设计优化的基本列式和计算框架，并结合数值算例的分析得到以下结论：

1）SIMP法展现出计算效率优势，尤其适用于需要快速处理的大规模工程优化问题。

2）基于离散变量的SAIP法能够生成边界清晰的拓扑设计，而且目标函数值更优，在需要高精度几何表达的增材制造中更具实用性。

3）尽管两种方法均能满足多孔结构设计需求，但在实际应用中需综合考量计算资源限制及对结构设计几何精度的要求。

参考文献

[1] 王辰， 刘义畅， 陆宇帆， 等. 考虑增材制造填充结构强度的拓扑优化方法［J］.上海交通大学学报， 2024， 58（3）： 333-341.

[2] Garaigordobil A， Ansola R， Querin O M， et al. Infill topology optimization of porous structures with discrete variables by the sequential element rejection and admission method［J］.Engineering Optimization， 2023， 55（3）： 457-475.

[3] Sigmund O. A 99 line topology optimization code written in Matlab［J］.Structural & Multidiscip-linary Optimization， 2001， 21（2）： 120-127.

[4] Andreassen E， Clausen A， Schevenels M， et al. Efficient topology optimization in MATLAB using 88 lines of code［J］.Structural & Multidisciplinary Optimization， 2011， 43（1）： 1-16.

[5] Liang Y， Cheng G. Topology optimization via sequential integer programming and Canonical relaxation algorithm［J］.Computer Methods in Applied Mechanics and Engineering， 2019， 348（10）： 64-96.

[6] Liang Y， Cheng G. Further elaborations on topo-logy optimization via sequential integer programming and Canonical relaxation algorithm and 128-line MATLAB code［J］.Structural and Multidisciplinary Optimization，2020， 61（1）： 411-431.

[7] Liang Y， Cheng G. Review of discrete variable topology optimization by sequential approximate integer programming［J］.Engineering Optimization， 2025， 57（1）： 130-160.

[8] Wu J， Aage N， Westermann R， et al. Infill Optimization for additive manufacturing—approaching bone-like porous structures［J］.IEEE Transactions on Visualization and Computer Graphics， 2018， 24（2）： 1127-1140.

[9] Liu J， Zou Z， Li Z， et al. A clustering-based multiscale topology optimization framework for efficient design of porous composite structures［J］.Computer Methods in Applied Mechanics and Engineering， 2025， 439： 117881.

[10] 高彤， 张卫红， 朱继宏. 惯性载荷作用下结构拓扑优化［J］. 力学学报，2009，41（4）：530-541.

[11] 高云凯， 张锁， 袁泽.考虑疲劳性能的驾驶室拓扑优化设计［J］. 汽车工程，2023，45（3）：468-476.

[12] 焦洪宇， 周奇才， 李文军， 等. 基于变密度法的周期性拓扑优化［J］.机械工程学报，2013，49（13）： 132-138.

[13] 吴一帆， 郑百林， 何旅洋， 等. 结构拓扑优化变密度法的灰度单元等效转换方法［J］.计算机辅助设计与图形学学报，2017，29（4）： 759-767.

[14] 昌俊康， 段宝岩. 连续体结构拓扑优化的一种改进变密度法及其应用［J］. 计算力学学报， 2009， 26（2）： 188-192.

[15] 廉睿超， 敬石开， 何志军， 等. 拓扑优化变密度法的灰度单元分层双重惩罚方法［J］. 计算机辅助设计与图形学学报， 2020，32（8）：1349-1356.

[16] 梁缘. 基于序列近似整数规划的通用高性能离散变量拓扑优化新方法 ［D］. 大连：大连理工大学， 2021.

[17] Rozvany G I N， Zhou M， Birker T. Generalized shape optimization without homogenization［J］. Structural Optimization， 1992， 4（3）： 250-252.

[18] Svanberg K， Werme M. Topology optimization by a neighbourhood search method based on efficient sensitivity calculations［J］. International Journal for Numerical Methods in Engineering， 2010， 67（12）： 1670-1699.

[19] Gao D Y. Solutions and optimality criteria to box constrained nonconvex minimization problems［J］.Journal of Industrial and Management Optimization， 2007， 3（2）： 293-304.

基金资助

国家自然科学基金(12002218)

辽宁省高校基本科研业务费项目(LJ212410143080)

工业装备结构分析优化与CAE软件全国重点实验室项目(GZ22108)