固体力学学报

方法名称	迭代次数	运行时间	柔顺度值	离散度值DR	离散度值M_nd
SIMPH	515	83.3	53.2	0.9385	0.0282
FPTO	169	49.8	50.1	0.5520	0.0861
SEMDOT	124	222.3	50.9	0.9072	0.0425
O-ICM	54	16.4	57.4	0.7893	0.1010
B-ICM	63	15.4	51.2	0.9392	0.0182

方法名称	迭代次数	运行时间	柔顺度值	离散度值DR	离散度值M_nd
SIMPH	515	83.3	53.2	0.9385	0.0282
FPTO	169	49.8	50.1	0.5520	0.0861
SEMDOT	124	222.3	50.9	0.9072	0.0425
O-ICM	54	16.4	57.4	0.7893	0.1010
B-ICM	63	15.4	51.2	0.9392	0.0182

方法名称	迭代次数	运行时间	柔顺度值	离散度值DR	离散度值M_nd
SIMPH	1262	311.6	30.8；30.8	0.9898	0.0027
FPTO	112	47.2	30.3；30.3	0.5379	0.0991
SEMDOT	122	314.5	30.4；30.4	0.9245	0.0299
O-ICM	78	33.4	32.1；32.1	0.8371	0.0774
B-ICM	46	17.8	30.6；30.6	0.9550	0.0135

方法名称	迭代次数	运行时间	柔顺度值	离散度值DR	离散度值M_nd
SIMPH	1262	311.6	30.8；30.8	0.9898	0.0027
FPTO	112	47.2	30.3；30.3	0.5379	0.0991
SEMDOT	122	314.5	30.4；30.4	0.9245	0.0299
O-ICM	78	33.4	32.1；32.1	0.8371	0.0774
B-ICM	46	17.8	30.6；30.6	0.9550	0.0135

结构拓扑优化ICM方法的B映射求解途径

PDF下载

彭细荣 ¹ , 隋允康 ²^,^** , 叶红玲 ² , 铁军 ³

固体力学学报 | 研究论文 2025,46(2): 162-176

收起

固体力学学报 | 研究论文 2025, 46(2): 162-176

结构拓扑优化ICM方法的B映射求解途径

全屏

彭细荣¹, 隋允康²^,^**, 叶红玲², 铁军³

作者信息

¹湖南城市学院土木工程学院，益阳，413000

²北京工业大学数学统计学与力学学院，北京，100022

³天津财经大学理工学院，天津，300222

通讯作者:

**E-mail：ysui@bjut.edu.cn.

Bi-mapping Solving Approach for the ICM Method of Structural Topology Optimization

Xirong Peng¹, Yunkang Sui²^,^**, Hongling Ye², Jun Tie³

Affiliations

¹School of Civil Engineering, Hunan City University, Yiyang, 413000

²College of Mathematics, Statistics and Mechanics, Beijing University of Technology, Beijing, 100022

³School of Polytechnic, Tianjin University of Finance and Economics, Tianjin, 300222

出版时间: 2025-04-23 doi: 10.19636/j.cnki.cjsm42-1250/o3.2024.052

文章导航

摘要

收起

本文基于结构拓扑优化的ICM（独立、连续和映射的简称）方法，发展出了双映射接力运作的求解途径，因“双”的英文简称为B，故称之为B-ICM求解途径. 该途径由两步组成：第一步是L（线性）映射作用于结构拓扑优化问题，使之成为离散模型，然后构造了约束函数；第二步是NL（非线性）映射作用于离散模型，使之成为连续模型，同时实现了单元拓扑变量由离散到连续的转换过程. 以往ICM方法的求解途径中，上述第一步只是起到理论推导的作用，构造约束函数同其它建模与求解等算法均包含在第二步里，因此属于“一步求解”途径. B-ICM虽然属于“两步求解”途径，但是优化模型寻优仍然沿用ICM方法惯用的序列对偶二次规划算法. 本文以位移约束的体积极小化结构拓扑优化问题为例，示例了上述建模及求解过程. 单载荷工况和多载荷工况的算例，均印证了本文的研究实现了预期构想，与目前旨在得到清晰拓扑的3种方法（1、考虑Heaviside投影的SIMP方法；2、浮动投影拓扑优化<简称为FPTO>方法；3、非惩罚的光滑边界材料分布拓扑优化<简称为SEMDOT>方法），以及ICM方法以往的求解途径，进行了迭代次数、清晰程度、寻优能力等方面的对比，结果表明B-ICM求解途径表现最好. 本文研究不仅丰富了ICM方法建模策略，推动了ICM方法求解途径的完善，也为解决模糊边界问题提供了一种优越的做法. 过往的连续体结构拓扑优化求解，因消除棋盘格和网格依赖性问题而采取的过滤操作导致最优拓扑构型产生模糊边界，而且过滤半径越大，则边界越模糊. 本文克服了这些令人堪忧的问题，可以成功地得到最优拓扑构型的清晰边界. 值得提及的是，本文研究的关键技术，可以移植到包括变密度方法等所有连续变量优化的方法中.

关键词

结构拓扑优化 / ICM方法 / 两步求解途径 / 线性映射 / 非线性映射 / 位移约束的体积极小化问题

Abstract

收起

A novel solving approach featuring bi-mapping, referred to as B-ICM, has been developed based on the ICM (independent, continuous, and mapping) method of structural topology optimization. B-ICM consists of two distinct steps: the first involves applying linear (L) mapping to the structural topology optimization problem, transforming it into a discrete model, and subsequently constructing the constraint function; the second entails implementing nonlinear (NL) mapping on this discrete model to create a continuous model, while converting continuous elemental topology variables into discrete ones. In contrast to the original ICM method, wherein the first step serves solely as a theoretical derivation, and the construction of constraint functions along with modeling and solving algorithms are all encompassed within the second step, which is categorized as a “one-step” approach, the B-ICM is classified as a “two-step” approach. Despite this distinction, it still employs the sequential dual quadratic programming algorithm commonly utilized in ICM methods for solving optimization models. We demonstrated this modeling and solving process using the structural topology optimization problem of volume minimization with displacement constraints. Results from both single-load and multi-load cases validated the effectiveness of our approach. We compared iteration count, clarity, and optimization capability across three methods for achieving distinct topologies: (1) the SIMP method considering Heaviside projection, (2) the floating projection topology optimization (FPTO) method, (3) the non-penalized method of smooth-edged material distribution for optimizing topology (SEMDOT), as well as the original ICM method. Results indicated that B-ICM outperforms these alternatives. This study not only enhances the modeling strategy and refines the solution approach of the ICM method, but also offers a superior technique for addressing blurry boundary problems. In continuum topology optimization, optimal topologies with blurry boundaries are typically generated through filtering operations designed to mitigate checkerboard patterns and mesh dependency issues. Notably, an increase in the filtering radius results in a more blurred boundary. Our study successfully addressed this challenge by achieving clear boundaries for optimal topologies. Importantly, the key techniques developed here are applicable to all continuous variable optimization methods, including the variable density method.

Key words

structural topology optimization / ICM (independent, continuous, and mapping) method / two-step solving approach / linear mapping / nonlinear mapping / volume minimization problem with displacement constraints

引用本文

彭细荣, 隋允康, 叶红玲, 铁军. 结构拓扑优化ICM方法的B映射求解途径. 固体力学学报, 2025 , 46 (2) : 162 -176 . DOI: 10.19636/j.cnki.cjsm42-1250/o3.2024.052

Xirong Peng, Yunkang Sui, Hongling Ye, Jun Tie. Bi-mapping Solving Approach for the ICM Method of Structural Topology Optimization[J]. Chinese Journal of Solid Mechanics, 2025 , 46 (2) : 162 -176 . DOI: 10.19636/j.cnki.cjsm42-1250/o3.2024.052

正文

收起

0　引言

收起

连续体结构拓扑优化是在指定的设计域内，设计材料的拓扑分布，以满足一定结构性能如位移、应力等约束条件下，使目标最优，如结构体积（或重量）极小化. 在采用有限元网格对结构进行离散、分析和优化的情况下，连续体结构拓扑优化问题本质上属于0-1离散优化问题，0代表相应子区域不分布材料，1代表相应子区域分布材料. 由于结构的有限元网格数通常很多，导致这类大规划的离散优化问题求解非常困难^[1]. 为解决这一困难，早期连续体结构拓扑优化的一类方法是将0-1之间的拓扑变量放松到[0，1]之间，从而可以利用成熟的基于连续变量的大规模优化求解算法. 如均匀化方法定义了微结构尺寸作为设计变量^[2]，变厚度法定义单元厚度为设计变量^[3]，变密度法定义单元人工密度为设计变量，采取惩罚函数使单元人工密度变量向0或1逼近，典型的如固体各向同性材料惩罚模型（SIMP：Solid Isotropic Material with Penalization）^[4]. 不同于上述各方法将拓扑变量挂靠在不同的物理量上，隋允康提出的ICM（独立、连续和映射）方法定义独立拓扑优化层次的设计变量，采用映射和逼近概念从数学层面建立基于连续变量的结构拓扑优化模型^[1,5-7]. 进化结构优化（ESO：Evolutionary Structural Optimization）^[8]以及后面发展的双向进化结构优化方法（BESO）^[9]则不建立优化模型，直接依据敏度数，采用启发式算法，对单元进行删除或保留. 上述各类基于单元设计变量的方法均存在棋盘格现象及网格依赖性问题^[10]，解决上述问题一种常用的方式是对单元敏度或单元设计变量进行过滤处理，但得到的最优拓扑在边界上存在大量灰度单元（即单元设计变量取值不为0或1，而取0到1之间的中间值，意味着此类单元需要圆整为0或1）^[11]. ESO或BESO得到的结果不存在灰度单元，但需要指定体积比约束条件，不同的体积比以及不同的删除策略对最优拓扑构型影响很大.

解决边界灰度单元过多的方法之一是采用Heaviside投影^[12-14]，但Heaviside投影方法中一次性太逼近于Heaviside函数将引起数值问题，表征对Heaviside函数逼近程度的参数需要渐次变化，从而导致优化迭代次数大幅增加. 另一种途径是采用基于离散变量的优化方法. Svanberg等将位移和应力约束等连续体拓扑优化问题采用基于离散变量的线性混合0-1规划算法求解，但只能求解小规模的离散变量拓扑优化问题^[15]. Svanberg和Werme等结合离散变量灵敏度和构造邻域搜索技术，采用分支定界法求解了中等规模的离散变量拓扑优化问题^[16]. Picelli借鉴SIMP和BESO方法中的高效灵敏度计算，构造一系列线性整数规划，使用商用的分支定界求解器求解了最小柔顺度等问题^[17]. Beckers采用序列近似规划构造非线性整数规划子问题，使用拉格朗日松弛算法求解了约束个数较少的大规模离散变量拓扑优化问题^[18]. 基于离散变量优化模型的方法还处于探索阶段，其对大规模问题的求解效率还有待提升.

与前述基于单元定义设计变量方法不同，另一类结构拓扑优化方法是定义边界相关参数为设计变量，通过边界演化来实现结构拓扑寻优. 如水平集方法以水平集函数的边界演化实现材料的最优分布^[19]. 可移动变形组件法^[20]、可移动变形孔洞法^[21]、移动可变形杆件法^[22]和几何投影法^[23,24]等，最近的如Zhang等提出的将可移动孔洞（MMV）与边界元法结合的方法^[25]，均是以几何组件的移动及边界形状实现拓扑变化. 这类方法可得到显式定义的几何边界，但通常其结构拓扑仍映射到有限元的网格单元上，其边界单元仍是以灰度单元的形式进行模拟计算，而且最优拓扑受初始的几何参数集影响.

相比于基于边界演化的方法，基于单元设计变量的方法可自由生成拓扑，如Zhang等将SIMP法与神经网络技术结合实现在结构拓扑中考虑人文的个性化设计^[26]. 但最优拓扑边界模糊，存在大量灰度单元. Huang结合基于单元设计变量法和基于边界设计法两类方法的优点，研究了一种浮动投影拓扑优化（Floating Projection Topology Optimization，简写为FPTO）方法，可以得到清晰拓扑，只有与边界线相交的单元处于灰度状态，而且其材料插值不需要如变密度法一样采取惩罚材料插值模型如SIMP插值，而仅需要采用自然的线性材料插值方案^[27,28]. Fu等基于SIMP方法的框架提出一种光滑边界材料分布拓扑优化（Smooth-Edged Material Distribution for Optimizing Topology简写为SEMDOT）方法，使用单元体分比作为设计变量，早期仍用SIMP的材料惩罚模型，之后可以不使用材料惩罚模型，而是直接用线性材料模型. 该方法同样能得到只有边界相交的单元处于灰度状态，边界光滑的最优拓扑^[29,30].

本文工作基于结构拓扑优化ICM方法^[5-7]的概念及建模思想，在优化的求解途径解法上进行了调整，第一步通过“L（线性）映射”表达离散优化模型时，用离散拓扑变量构造了约束函数；第二步通过“NL（非线性）映射”转换离散模型为连续模型时，在将单元拓扑变量连续化的同时，使目标和约束条件皆成为连续拓扑变量的显函数. 以B-ICM命名这个B（双）映射求解途径.

尽管B-ICM的第二步同以往ICM方法的运作完全相同，但是约束函数近似显式化的表达算法却留在第一步里，尤其在迭代寻优中，这种保留显示了数值求解的合理性. 在达到同样最优拓扑离散效果的要求下，B-ICM收敛的更快.

1　结构拓扑优化离散模型的L映射和建模

收起

ICM方法的理论体系其实就是两步走出来的：第一步采用阶跃函数的反函数跨栏函数，把物理量表示为离散拓扑变量的严格函数关系，形成了准确的离散优化模型；第二步采用过滤函数逼近跨栏函数，离散优化模型转换为近似的连续优化模型，此时离散的0/1拓扑变量扩展为[0，1]之间的连续拓扑变量. ICM方法的第一步只是理论推导过程，求解集中在第二步上，包括构造约束函数与其它建模与求解等算法，因此遵循“一步求解”策略.

B-ICM属于“两步求解”途径，第一步视为“L映射”，置离散优化模型于求解的相对独立阶段. 对结构设计域进行有限元离散，对应于每一个单元，定义一个独立的拓扑层次的设计变量

，以

表示对应单元删除，

表示对应单元保留. 由于仅取0或1两个离散值，称之为离散拓扑设计变量.

单元各物理量与离散拓扑变量

的关系可表示为乘积形式：

(1)

其中

(2)

称完成离散线性变换运作的式（1）为“L映射”，L是Linear的缩写.

由此，物理量只能取到两个状态：0或

，称

为结构单元对应的固有量. 物理量可以为单元的体积（或重量）、单元刚度矩阵、单元质量矩阵、单元材料对应的弹性模量等等.

以单元体积和结构总体积为例：

(3)

其中，

为单元固有体积，N为有限元网格中单元总数.

单元刚度矩阵和结构总刚度矩阵分别为：

(4)

其中，

为单元固有刚度矩阵. 对矩阵而言，这里的求和符号意指单元矩阵组装成整体矩阵的过程.

将约束函数表示为函数向量形式为G=（G₁，…，G_M）^T，由此可表示基于离散拓扑变量的优化模型：

(5)

其中，

为离散拓扑变量向量，目标函数为结构总体积，M为约束总数，N为设计变量总数.

以位移约束为例，说明约束函数的显式化的过程. 式（5）中的位移约束

，按惯用表达，记为

. 根据莫尔定理，j自由度的位移为^[5,6]：

(6)

其中，

表示位移实载荷工况下单元i的节点位移向量，

与

分别表示实、虚载荷工况下单元i的节点力向量. 将式（4）代入式（6）得到：

(7)

其中

(8)

很可惜，实、虚载荷下式（7）中两个节点力向量皆为离散拓扑变量的隐函数，因而尚未得到位移以离散拓扑变量表达的显函数. 不过，此式提供启发：若式中分子在迭代中为常数，则位移成为非零的拓扑变量倒函数. 为此，下面在一次寻优迭代中进行讨论，第V次迭代时，式（7）的分子为：

(9)

尽管超静定结构的内力随单元面积的变化而改变，但是在每一次迭代寻优中，按照暂时“静定化假设”，得知各单元内力不随单元改变，这种处理无论在截面层次还是拓扑层次的优化，皆屡试不爽. 换句话说，在第v次迭代寻优中，式（7）中的约束系数若用式（9），则可视为常数. 于是，由式（7）得到离散拓扑变量的如下显函数：

(10)

为减少计算量，根据式（6），式（9）可转换为：

(11)

其中下标j对应于虚工况V，其为只在j自由度方向作用单位荷载的虚工况^[5,6]. 具体编程时，可以自行选上述的一种. 总结一下，位移约束推导的可用结果是式（8）、式（10）和式（11）三式.

2　化离散模型为连续模型的NL映射处理

收起

上一节虽然叙述了L映射后得到离散拓扑优化模型的过程，但是离散问题的求解存在组合爆炸困难，须将其转换能够近似逼近它的连续化问题.

式（1）与式（2）其实揭示了ICM方法中跨栏函数^[5-7]的作用，用过滤函数可以逼近它，见图1.

采用过滤函数f（t_i）逼近跨栏函数

，取得一石二鸟的效果：取0/1的离散变量拓广为闭区间[0，1]上的连续变量，离散模型式（5）扩展为如下连续模型：

(12)

其中f_v（t_i）为目标函数的过滤函数，f_g（t_i）为约束函数的过滤函数.

以上做法是ICM方法里熟知的由离散到连续的变换，为了行文方便，本文命名为NL映射，NL是非线性的英文缩写.

此时结构拓扑优化问题的位移约束，由式（10）和式（11）连续化为如下两式：

(13)

(14)

其中，f_k（t_i）是单元刚度矩阵的过滤函数.

这里需要做一个重要说明：进行NL映射后，L映射阶段的式（8）就不需要单独予以考虑了，原因是连续变量在0与1的闭区间上，接近0的数值可以用一个很小的正数例如取

替代，这不仅可以防止有限元分析时出现奇异，而且优化模型中待删单元对应的项也不会出现数值发散. 因此，原本是难以预判的主动集计算，相应的式（8）就可以免去考虑了，这是十分有利的处理.

本文篇幅所限，B映射之后的求解，将结合位移约束的结构拓扑优化问题，在下一节予以介绍.

3　位移约束结构体积极小化问题求解

收起

累计前面的阐述，顺序进行了包括L映射和NL映射的B映射之后，位移约束下结构总体积最小的结构拓扑优化问题，可以表示为如下显式模型：

(15)

其中f_v（t_i）为单元体积过滤函数.

令

(16)

则

(17)

于是式（15）转化为：

(18)

其中，

将目标函数进行二阶Taylor展开，其第i项为：

(19)

其中，目标函数一阶和二阶导数为：

(20)

(21)

由式（16）得：

(22)

从式（22）求二阶导数：

(23)

将式（22）代入式（20），同时将式（22）和式（23）代入式（21），分别得：

(24)

(25)

略去式（19）中对最优解没有影响的常数项，并记：

(26)

(27)

即得到式（18）的二阶近似优化模型：

(28)

前述过滤函数可取专著[5,6]中介绍的任意一种形式，例如当取为幂函数形式时，单元及刚度矩阵的过滤函数为：

(29)

此时，对应的目标函数中的系数可求得为：

(30)

(31)

对式（28）可以采用ICM方法惯用的序列对偶二次规划算法求解^[5,6].

收敛准则取为：

(32)

其中ε为设定的收敛限值，通常可取ε=1×10^-4.

从本质看，L映射是离散-离散的过程，而NL映射是离散-连续的过程. 尽管对B映射模型的求解得到了很接近离散的最优拓扑构型，但是还缺一个从连续-离散逆映射的一锤定音过程. 为此，这里简要地介绍一下本文方法中的逆映射做法，包括边界圆整化技术.

第一步采用加权平均的方式，用单元拓扑设计变量得出节点拓扑变量，即：

(33)

式中

表示节点拓扑变量，

表示单元拓扑变量，N_ie是与i节点相连的单元数，α_j是与i节点相连的单元j的夹角.

第二步是按节点拓扑分布场的等值线（或面）将区域划分成删除区域和保留区域. 等值线（或面）阀值简单处理可按最优体积比不变的原则确定，也可按使主动性能约束取等式的条件来确定，但需要做有限元分析，前者按最优体积比不变处理则不需要做有限元分析，但可能会使性能约束与连续模型时的结果有略微差别.

以下算例中，过滤函数取为幂函数形式，指数分别取为α_v=1，α_k=10. 最优拓扑边界按最优体积比不变的原则确定等值线（或面）阀值.

4　算例与讨论

收起

以下算例将下列5种方法得到的结果进行对比：

（1）采用Heaviside投影的SIMP方法^[12-14]，惩罚函数采用幂函数形式，指数取3，简记为SIMPH；

（2）FPTO方法^[27,28]，物理量与设计变量间采用线性关系式，简记为FPTO；

（3）采用非惩罚的SEMDOT方法^[30]，物理量与设计变量间采用线性关系式，简记为SEMDOT；

（4）ICM一步解法^[5,6]，过滤函数采用幂函数形式，单元刚度矩阵的过滤函数幂指数取3，简记为OICM；

（5）B-ICM方法，过滤函数采用幂函数形式，单元刚度矩阵的过滤函数幂指数取10.

鉴于前三种方法一般求解体积约束下结构柔顺度极小化问题，而本文遵循ICM方法一直秉承的理念，以结构性能为约束结构体积极小为目标建模. 两种模型在某些情况下是具有互逆关系^[31]，为方便比较，在设置ICM和B-ICM解法的位移约束条件时，通过确定适当的位移约束值，使之最优体积目标值正好与前三种方法设置的体积约束值相等.

运行时间的比较是基于所有方法均运行在同一硬、软件环境之下的前提之下. 柔顺度值的比较基于所有最优拓扑结构的结构体积比均相同的前提之下. 离散度值用于衡量最优拓扑结构时拓扑变量分布值的离散程度. 文献[32]定义

(34)

描述离散程度，但其值越小，代表离散程度超高，与离散程度成反比关系. 本文采取如下更为直观并合乎常理的离散度定义，其值越大，代表离散程度超高，与离散程度成正比关系.

定义

(35)

为设计变量值逼近0（或

）或1的单元，称为离散单元，离散度DR定义为满足式（35）的离散单元总数除以总单元数. 后面算例中所比较的及图形中所描述的均是DR值. 但在表1和表3中同时也列举了M_nd值.

4.1　算例1：单工况二维悬臂梁结构拓扑优化设计

此算例为经典算例，算例参数参考文献[29]. 如图2所示：基结构为150×100的平面体，厚度为1，材料弹性模量为1，泊松比为0.3. 一集中载荷F=1作用于右边界中心位置. 左边界采用固定支承. 有限元网格为150×100的边长为1的正方形单元，过滤半径取为2.5.

SIMPH、FPTO和SEMDOT三种方法建立的模型为体积比约束下柔顺度极小的优化模型，体积比约束值设为0.3. 为了比较，后面两种O-ICM和B-ICM解法建立的是位移约束下结构总体积极小化的优化模型，位移约束力作用方向上位移大于或等于一定值. 为方便比较，使所有最优结构体积比相同（即O-ICM和B-ICM解法的目标函数值为0.3），通过试凑确定后两种方法的位移约束值分别为OICM的-57.4及B-ICM的-51.2. 依据本算例条件及文献31的结论可知，本算例中结构柔顺度等于力作用方向的位移绝对值，基于此原因，后面直接将O-ICM和B-ICM解法最优点的位移值转化为柔顺度，再与其它计算方法进行比较.

参与的5种方法得到计算结果比较见表1所示，最优拓扑见表2所示. 值得指出的是：FPTO和SEMDOT两种方法声称得到清晰且光滑边界的最优拓扑，其有限元分析模型仍是基于单元拓扑变量（或单元人工密度）连续分布的模型，故此，对这两种方法补充输出单元拓扑变量分布图，对其它3种方法则依据文献[29，30]的方式，补充输出清晰且光滑边界的最优拓扑.

从表1的结果对比中可以看到，声称可得到清晰拓扑的4种方法SIMPH、FPTO、SEMDOT和BICM中，本文提出的B-ICM解法离散度值最大，迭代次数最少，运行时间最短，而柔顺度值相差并不大，柔顺度最小的是FPTO，本文方法比其结果高2.2%，但运行时间省220.8%，离散度值高0.3872. 从表2所列的最优拓扑图比较来看，除SIMPH方法图形不相同外，其它4种方法图形基本相同. BICM解法与原来的O-ICM解法相比，边界更清晰，图形几乎一样，前两个结点位置更向结构中点靠近一些，FPTO得到的最优拓扑的柔顺度最小，其结构的前两个结点位置也是最靠近结构的水平中轴线.

4.2　算例2：多工况二维悬臂梁结构拓扑优化设计

此算例亦为经典算例，算例参数参考文献[14]. 如图3所示150×150矩形平板，材料弹性模量E=1，泊松比0.3. 平板左边界固定约束，受两种荷载工况作用，工况1：集中荷载F₁=1竖直向上作用于右上角；工况2：集中载荷F₂=1竖直向下作用于右下角. 有限元网格为150×100的边长为1的正方形单元，过滤半径取为2.5.

与算例1采用相同的方式进行比较. SIMPH、FPTO和SEMDOT三种方法设定体积比约束0.4. 后面两种O-ICM和B-ICM解法通过试算出位移约束值，使最优目标体积比也正好等于0.4. 经过计算，O-ICM方法的位移约束值分别F₁作用力方向位移小于或等于32.1和F₂作用力方向位移大于或等于-32.1；B-ICM解法的位移约束值分别F₁作用力方向位移小于或等于30.6和F₂作用力方向位移大于或等于-30.6.

用5种方法计算得到的结果分别如表3和表4所示. 从表3和表4中的结果比较可以看到，在多荷载工况下，也是B-ICM方法求解效率最高，离散度也非常高，仅SIMPH方法的离散度值比它更高，这是因为SIMPH方法求解这个问题时原程序设定的0.01收敛值下无法正常收敛，放松收敛限值为0.03时，也经1262次迭代才收敛，致使其逼近Heaviside函数的参数非常大，使得结果的离散度值很高. 结构的柔顺度值相差都不大，所得到的最优拓扑是几乎相同的. 值得指出的是，其它方法在处理多荷载工况问题时，是将各工况的柔顺度加权处理成单目标问题，其加权系数是会对最优拓扑有影响的，在复杂的多工况问题时，难以合理地确定各工况间的加权系数. 而O-ICM和B-ICM解法都是以位移为约束建立模型，不需要进行加权处理，所建立的优化模型更为合理.

B-ICM解法的目标值体积比和离散度值的迭代历史曲线，以及典型步骤的拓扑变量分布图形如图4所示，从图中可以看到，迭代寻优过程非常平稳、高效，目标值体积比快速趋于稳定，离散度值快速增大，最优拓扑在第15迭代时已经非常清晰，离散度值为0.8971，体积比为41.43%. 表明本文所提出的B-ICM解法非常高效且寻优能力很强，得到的最优拓扑边界清晰.

为更清晰地看到处理模糊边界的效果，采取较大的过滤半径6.0，用O-ICM和B-ICM得到的最优拓扑如图5所示. 在较大的过滤半径下，O-ICM经119次迭代收敛得到的最优拓扑边界模糊，离散度值为0.6714，而B-ICM经59次迭代得到的最优拓扑边界非常清晰，离散度值为0.8965.

4.3　算例3：三维悬臂梁结构拓扑优化设计

算例参数参考文献[28]. 如图6所示60×40×20三维长方体，材料弹性模量E=1，泊松比0.3. 长方体左端面固定约束，右下角边线中点受一集中力F=1作用. 有限元网格为60×40×20的边长为1的正方体单元，过滤半径取为3.0. 位移约束集中力作用方面向下位移（位移竖直向上为正）大于或等于-10. 经54次迭代收敛，最优拓扑值分布及光滑边界的最优拓扑结构如图7所示. 最优体积比为12.43%，位移约束值为-9.99，满足约束条件. 离散度值为0.9244. 目标体积比和离散度值迭代历史曲线如图8所示，约束迭代曲线如图9所示. 从图形可以看到得到的最优拓扑值离散度很高，结构边界清晰. 迭代过程平稳收敛. 相同问题，文献[28]用FPTO方法，以体积比15%为约束条件，极小化结构柔顺度，经164次迭代收敛，本文方法迭代次数少110次.

5　结语

收起

本文基于ICM方法，强调了两次映射，展示了建模的两步途径：

在第一步的L映射之后，得到了迭代中每次建立的离散优化模型，显示了离散拓扑变量的线性作用——体现在离散线性函数的准确乘积上.

在第二步的NL映射之后，得到了迭代中每次寻优的连续优化模型，显示了连续拓扑变量的非线性作用——体现在过滤函数对跨栏函数的逼近上.

求解算法仍然沿用ICM方法惯用的序列对偶二次规划算法，可以得到边界清晰的拓扑结构，克服了ICM原有的“一步解法”采用过滤操作消除棋盘格和网格依赖性问题时所产生的模糊边界的不足.

以位移约束结构体积极小化的连续体结构拓扑优化问题为例，示例了上述建模及求解过程. 并以单载荷工况和多载荷工况的算例，与3种旨在得到清晰拓扑的方法如SIMPH、FPTO和SEMDOT进行了计算效率、清晰程度和寻优能力等综合对比，也与O-ICM进行了对比，结果表明本文所提出的B-ICM解法综合比较表现最优.

本文进行了大量的数值计算及比较，限于篇幅不能全面展开，简要提示如下几点：

（1）SIMPH、FPTO和SEMDOT方法均需要指定体积比约束，而体积比约束事先是无法确知的，以上各方法在处理多荷载工况问题时，各工况间的加权系数也难以合理确定. 而本文的方法沿用ICM方法一直沿用的建模传统，以体积极小为目标，结构性能为约束，该建模的优势在文献[31，33-35]等中已经进行了详细的介绍，本文不予赘述.

（2）B-ICM解法收敛条件简单，只需要指定目标函数的相对误差限；过滤函数的参数取值范围较宽泛，数值实验表明，采用幂函数形式的过滤函数时，其参数α_k可以从3到20范围内取值均不会出现数值问题；寻优算法也不需要设置每迭代步设计变量的运动极限. 简单的参数设置为程序的稳健性提供了保障. 而SIMPH、FPTO和SEMDOT等各方法均需要设置复杂的组合收敛条件或渐近的参数变化，其采用的MMA算法或准则算法也都是需要设置设计变量的运动极限，程序参数的设置比本文所提出方法更为复杂. 此外，不同于SIMP等方法中常用的一阶近似的MMA寻优算法，ICM方法一直使用高效的序列对偶二阶近似的算法，如序列对偶二次规划算法（SDQP）^[1,5-7]或对偶显式模型DPEM方法^[36,37]，此亦是本文方法效率更高的原因之一.

（3）逆映射采取两步做法，第一步用各个汇交单元的拓扑变量值加权插值得到该汇交节点的拓扑变量值；第二步按拓扑变量场的等值线（或面）确定最优拓扑边界，可以得到清晰光滑的拓扑构型.

（4）本文研究的关键技术，可以移植到包括变密度方法等所有连续变量优化模型对应的方法中，即使相关模型并不是独立的拓扑优化层次，而是挂靠在较低的优化层次上，也必定会产生两个效果：其一是相关的映射函数可以比原来的非线性程度大大提高，其二是迭代收敛将会加快. 这些优点本文在对于ICM方法的改进上都做了详细阐述，并且已见成效，其缘由皆归属于双映射的结果，并不因为依赖于ICM方法的本身. 也就是说，双映射是一个普适的处理技术，本文只是以ICM方法做了成功的实验而已. 相信读者能够理解，限于篇幅，本文不予赘述了.

基金

收起

湖南省自然科学基金项目(2022JJ30113)
湖南省教育厅重点项目(21A0507)
大连理工大学工业装备结构分析国家重点实验室开放基金项目(GZ23104)

参考文献

收起

文献

收起

参考文献引证文献

排序方式：

[1]

隋允康. 建模·变换·优化：结构综合方法新进展[M]. 大连: 大连理工大学出版社, 1996.

(Sui

Y K

. Modelling, Transformation, and Optimization-New Developments of Structural Synthesis Method[M]. Dalian University of Technology Press, Dalian, China. 1996. (in Chinese))

[2]

Bendsoe

, Kikuchi

. Generating optimal topologies in structural design using a homogenization method[J]. Computer Methods in Applied Mechanics and Engineering, 1988, 71(2): 197-224.

[3]

程耿东, 张东旭. 受应力约束的平面弹性体的拓扑优化[J]. 大连理工大学学报, 1995, 35(1): 1-9.

(Cheng

G D

, Zhang

D X

. Topological optimization of plane elastic continuum with stress constraints[J]. Journal of Dalian University of Technology, 1995, 35(1): 1-9. (in Chinese))

[4]

Mlejnek

H P

. Some aspects of the genesis of structures[J]. Structural optimization, 1992, 5: 64-69.

[5]

隋允康, 叶红玲. 连续体结构拓扑优化的ICM方法[M]. 北京: 科学出版社, 2013.

(Sui

Y K

, Ye

H L

. Continuum Topology Optimization Methods ICM[M]. Beijing: Science Press, 2013. (in Chinese))

[6]

Sui

Y K

, Peng

X R

, Modeling, Solving and Application for Topology Optimization of Continuum Structures: ICM Method Based on Step Function[M]. Cambridge: Elsevier Inc., 2018.

[7]

隋允康, 彭细荣. 结构拓扑优化ICM方法新进展——概念深化及理论拓展[M]. 北京: 科学出版社, 2023.

(Sui

Y K

, Peng

X R

. New Progress in ICM Method for Structural Topology Optimization-Concepts Deepening and Theoretical Expansion[M]. Beijing: Science Press, 2023. (in Chinese))

[8]

Xie

Y M

, Steven

G P

. A simple evolutionary procedure for structural optimization[J]. Computers & Structures, 1993, 49: 885-896.

[9]

Yang

X Y

, Xie

Y M

, Steven

G P

. Bidirectional evolutionary method for stiffness optimization[J]. AIAA Journal. 1999, 37(11): 1483-1488.

[10]

Sigmund

, Petersson

. Numerical instabilities in topology optimization: a survey on procedures dealing with checkerboards, mesh-dependencies and local minima[J]. Structural Optimization, 1998, 16(1): 68-75.

[11]

Bourdin

. Filters in topology optimization[J]. International Journal for Numerical Methods in Engineering, 2001, 50(9): 2143-2158.

[12]

Sigmund

. On the design of compliant mechanicms using topology optimization[J]. Mechanics of Structures and Machines, 1997, 25(3): 485-526.

[13]

Guest

, Prevost

, Belytschko

. Achieving minimum length scale in topology optimization using nodal design variables and projection functions[J]. International Journal for Numerical Methods in Engineering, 2004, 61(2): 238-254.

[14]

Andreassen

, Clausen

, Lazarov

B S

, Sigmund

. Efficient topology optimization in MATLAB using 88 lines of code[J]. Structural Multidisciplinary Optimization, 2011, 43(1): 1-16.

[15]

Stolpe

, Svanberg

. Modelling topology optimization problems as linear mixed 0-1 programs[J]. International Journal for Numerical Methods in Engineering, 2003, 57(5): 723-739.

[16]

Svanberg

, Werme

. Sequential integer programming methods for stress constrained topology optimization[J]. Structural and Multidisciplinary Optimization, 2007, 34(4): 277-299.

[17]

Sivapuram

, Picelli

. Topology optimization of binary structures using Integer Linear Programming[J]. Finite Elements in Analysis and Design, 2018, 139: 49-61.

[18]

Beckers

. Topology optimization using a dual method with discrete variables[J]. Structural Optimization, 1999, 17(1): 14-24.

[19]

Osher

, Sethian

. Fronts propagating with curvature dependent speed-algorithms based on hamilton-jacobi formulations[J]. J Comput Phys, 1988, 79(1): 12-49.

[20]

Guo

, Zhang

W S

, Zhong

W L

. Doing topology optimization explicitly and geometrically-A new moving morphable components based framework[J]. Journal of Applied Mechanics, 2014, 81(8): 081009.

[21]

Zhang

W S

, Chen

J S

, Zhu

X F

, et al. Explicit three dimensional topology optimization via moving morphable void (MMV) approach[J]. Computer Methods in Applied Mechanics and Engineering, 2017, 322: 590-614.

[22]

Hoang

V N

, Jang

G W

. Topology optimization using moving morphable bars for versatile thickness control[J]. Computer Methods in Applied Mechanics and Engineering, 2017, 317: 153-173.

[23]

Norato

J A

. A geometry projection method for the optimal distribution of short fiber reinforcements[C]. In: ASME 2015 International Design Engineering Technical Conferences and Computers and Information in Engineering Conference, American Society of Mechanical Engineers.

[24]

Zhang

, Norato

J A

, Gain

A L

, Lyu

. A geometry projection method for the topology optimization of plate structures[J]. Structural and Multidisciplinary Optimization, 2016, 54(5): 1173-1190.

[25]

Zhang

W S

, Wang

, Du

Z L

, Liu

, Youn

S K

, Guo

. Machine-learning assisted topology optimization for architectural design with artistic flavor[J]. Computer Methods in Applied Mechanics and Engineering, 2023, 413: 116041.

[26]

Zhang

W S

, Tian

H H

, Sun

, Feng

W Z

. Topology optimization for pressure loading using the boundary element-based moving morphable void approach[J]. Advances in Engineering Software, 2024, 195: 103689.

[27]

Huang

X D

. Smooth topological design of structures using the floating projection[J]. Engineering Structures, 2020, 208: 110330.

[28]

Huang

X D

, On smooth or 0/1 designs of the fixed-mesh element-based topology optimization[J]. Advances in Engineering Software, 2020, 151: 102942.

[29]

Y F

, Rolfe

, Chiu

L N

, Wang

, Huang

X D

, Ghabraie

. SEMDOT: smooth-edged material distribution for optimizing topology algorithm[J]. Adv. Eng. Softw., 2020, 150: 102921.

[30]

Y F

, Long

, Rolfe

. On non-penalization SEMDOT using discrete variable sensitivities[J]. Journal of Optimization Theory and Applications, 2023, 198: 644-677.

[31]

彭细荣, 隋允康. 对连续体结构拓扑优化合理模型的再探讨[J]. 固体力学学报, 2016, 37(02): 181-191.

(Peng

X R

, Sui

Y K

. A further discussion on rational topology optimization models for continuum structures[J]. Chinese Journal of Solid Mechanics., 2016, 37(02): 181-191. (in Chinese))

[32]

Sigmund

. Morphology-based black and white filters for topology optimization[J]. Structural Multidisciplinary Optimization, 2007, 33(4/5): 401-424.

[33]

彭细荣, 隋允康. 应该为结构拓扑优化的设计变量正名和处理方法顺言[J]. 北京工业大学学报, 2016, 42(12): 27-37.

(Peng

X R

, Sui

Y K

. Name correction for design variable of structural topology optimization and presentation of its corresponding treatment method[J]. Journal of Beijing University of Technology, 2016, 42(12): 27-37. (in Chinese))

[34]

隋允康, 彭细荣, 叶红玲, 铁军. 互逆规划理论及其用于建立结构拓扑优化的合理模型[J]. 力学学报, 2019, 51(06): 1940-1948.

(Sui

Y K

, Peng

X R

, Ye

H L

, Tie

. Reciprocal programming theory and its application to establish a reasonable model of structural topology optimization[J]. Chinese Journal of Theoretical and Applied Mechanics, 2019, 51(06): 1940-1948. (in Chinese))

[35]

铁军, 隋允康, 彭细荣. 互逆规划的拓宽和深化及其在结构拓扑优化的应用[J]. 力学学报, 2020, 52(06): 1822-1837.

(Tie

, Sui

Y K

, Peng

X R

. Widening and deepening of reciprocal programming and its application to structural topology optimiation[J]. Chinese Journal of Theoretical and Applied Mechanics, 2020, 52(06): 1822-1837. (in Chinese))

[36]

Sui

Y K

, Peng

X R

. Explicit model of dual programming and solving method for a class of separable convex programming problems[J]. Engineering Optimization, 2019, 51(9): 1604-1625.

[37]

彭细荣, 隋允康. 化可分离凸规划为对偶规划显式模型的普适解法研究比较[J]. 力学学报, 2023, 55(10): 2417-2426.

(Peng

X R

, Sui

Y K

. Research on universal solution of transforming separable convex programing to dual programming with explict model[J]. Chinese Journal of Theoretical and Applied Mechanics, 2023, 55(10): 2417-2426. (in Chinese))

2025年第46卷第2期

PDF下载

引用本文

BibTeX

文章信息

doi: 10.19636/j.cnki.cjsm42-1250/o3.2024.052

接收时间：2024-11-03
首发时间：2026-03-20
出版时间：2025-04-23

补充材料

相关文章

文章信息

作者

出版历史

收稿日期：2024-11-03

基金

湖南省自然科学基金项目(2022JJ30113)

湖南省教育厅重点项目(21A0507)

大连理工大学工业装备结构分析国家重点实验室开放基金项目(GZ23104)

作者信息

¹湖南城市学院土木工程学院，益阳，413000

²北京工业大学数学统计学与力学学院，北京，100022

³天津财经大学理工学院，天津，300222

通讯作者:

**E-mail：ysui@bjut.edu.cn.

参考文献

分享链接

https://castjournals.cast.org.cn/joweb/gtlxxb/CN/10.19636/j.cnki.cjsm42-1250/o3.2024.052

分享至

全文二维码

扫描看全文

引用本文

BibTeX

本文的引用情况

2种不同金属材料的力学参数

科 Family	属数 Number of genus	种数 Number of species	占总种数比例 Percentage of total species (%)	属 Genus	种数 Number of species	占总种数比例 Percentage of total species (%)
鹅膏菌科Amanitaceae	2	11	5.26	鹅膏菌属 Amanita	10	4.78
小菇科 Mycenaceae	2	12	5.74	丝盖伞属 Inocybe	5	2.39
多孔菌科 Polyporaceae	8	14	6.70	蜡蘑属 Laccaria	5	2.39
红菇科 Russulaceae	3	23	11.00	小皮伞属 Marasmius	6	2.87
				小菇属 Mycena	11	5.26
				光柄菇属 Pluteus	5	2.39
				红菇属 Russula	17	8.13
				栓菌属 Trametes	5	2.39

关闭全屏

BibTeX
EndNote
RefWorks
TxT

方法名称	迭代次数	运行时间	柔顺度值	离散度值DR	离散度值M_nd
SIMPH	515	83.3	53.2	0.9385	0.0282
FPTO	169	49.8	50.1	0.5520	0.0861
SEMDOT	124	222.3	50.9	0.9072	0.0425
O-ICM	54	16.4	57.4	0.7893	0.1010
B-ICM	63	15.4	51.2	0.9392	0.0182

方法名称

迭代次数

运行时间

柔顺度值

离散度值DR

离散度值M_nd

SIMPH

515

83.3

53.2

0.9385

0.0282

FPTO

169

49.8

50.1

0.5520

0.0861

SEMDOT

124

222.3

50.9

0.9072

0.0425

O-ICM

16.4

57.4

0.7893

0.1010

B-ICM

15.4

51.2

0.9392

0.0182

方法名称

单元拓扑值分布

光滑边界最优结构

SIMPH

FPTO

SEMDOT

O-ICM

B-ICM

方法名称	迭代次数	运行时间	柔顺度值	离散度值DR	离散度值M_nd
SIMPH	1262	311.6	30.8；30.8	0.9898	0.0027
FPTO	112	47.2	30.3；30.3	0.5379	0.0991
SEMDOT	122	314.5	30.4；30.4	0.9245	0.0299
O-ICM	78	33.4	32.1；32.1	0.8371	0.0774
B-ICM	46	17.8	30.6；30.6	0.9550	0.0135

方法名称

迭代次数

运行时间

柔顺度值

离散度值DR

离散度值M_nd

SIMPH

1262

311.6

30.8；30.8

0.9898

0.0027

FPTO

112

47.2

30.3；30.3

0.5379

0.0991

SEMDOT

122

314.5

30.4；30.4

0.9245

0.0299

O-ICM

33.4

32.1；32.1

0.8371

0.0774

B-ICM

17.8

30.6；30.6

0.9550

0.0135

方法名称

单元拓扑值分布

光滑边界最优结构

SIMPH

FPTO

SEMDOT

O-ICM

B-ICM