基于改进四叉树和比例边界有限元法的自适应设计域拓扑优化方法

王浩 王江北 罗浩东 王立文

摘要：

针对大型结构拓扑优化计算成本高和固体各向同性材料惩罚模型（SIMP）在优化后结构边界处求解精度低的问题，提出了一种基于SIMP法的自适应设计域（ADD）拓扑优化方法。将改进四叉树法应用在拓扑优化的过程中，通过自动划分不同等级的网格单元来减少网格数量、减轻计算负担并提高边界处求解精度；采用比例边界有限元法（SBFEM）实时计算划分后结构的有限元信息，解决了不同等级网格间悬挂节点的问题。所提方法可在初始网格相对较少的情况下得到更加精确的结果，大幅度地降低了计算成本。数值算例结果表明，所提方法在最终结构边界处精度相同的情况下，计算时间最快可缩短为原来的1/100，可以为后续降低大型结构拓扑优化的计算成本提供参考。

关键词：拓扑优化；改进四叉树法；比例边界有限元法；网格自适应

中图分类号：TH122

DOI：10.3969/j.issn.1004132X.2024.05.016

开放科学（资源服务）标识码（OSID）：

An Adaptive Design Domain Topology Optimization Method Based on

Improved Quadtree and SBFEM

WANG Hao1，2  WANG Jiangbei3  LUO Haodong3  WANG Liwen4

1.College of Safety Science and Engineering，Civil Aviation University of China，Tianjin，300300

2.Engineering Techniques Training Center，Civil Aviation University of China，Tianjin，300300

3.College of Aeronautical Engineering，Civil Aviation University of China，Tianjin，300300

4.Institute of Science and Technology lnnovation，Civil Aviation University of China，Tianjin，300300

Abstract： Aiming at the problems of high computational cost of large-scale structure topology optimization and low solving accuracy of solid isotropic material with penalization（SIMP） at the optimized structure boundary， an adaptive design domain（ADD） topology optimization method was proposed based on SIMP. The improved quadtree method was applied in the processes of topology optimization to reduce the computational burden and improve the solution accuracy at the boundaries through different levels of grid cell. The SBFEM was used to calculate the finite element information of the partitioned structures in real-time， solving the problem of hanging nodes between different level cells. The proposed method might obtain more accurate results with fewer initial grids and significantly reduce the computational cost. Numerical example results demonstrate that， with the same final structural boundary accuracy， the proposed method may reduce the computation time to 1/100 of the original at most， offering a reference for reducing the computational cost of large-scale structural topology optimization in subsequent processes.

Key words： topology optimization; improved quadtree method; scaled boundary finite element method（SBFEM）; grid adaptation

收稿日期：20230914

基金項目：国家自然科学基金民航联合研究基金重点项目（U2133202）；四川省重大科技专项（2021ZDZX0001）；中国民航大学研究生科研创新项目（2022YJS058）

0  引言

拓扑优化方法是一种给定约束和目标函数，在规定设计域内对结构材料进行寻优的算法，可以对结构零件进行轻量化设计。目前，该方法已被广泛地应用于航空航天、汽车工程、机械工程、海洋工程［1-6］等领域。早期由于计算条件的限制，主要应用于小型结构件，随着有限元的发展，越来越多的大型结构［7］也会采用拓扑优化方法来达到减重的目的。

在航空航天领域中，飞机机身、发动机吊架、机翼前缘翼肋等大尺寸且具有高精度要求的零件，动辄需要几十万个网格单元，直接求解带来的计算负担让人望而却步，因此寻找一种适用于具有大量网格或大型结构的拓扑优化方法具有重要研究意义。

拓扑优化过程需要反复迭代，每次迭代过程都需要对结构进行性能分析，有限元计算占比较大，其中网格数量直接影响最终的计算时间。PRAMOD等［8］提出了利用自适应相场方法模拟脆性断裂，将所需网格数量减少为原来的1/10； BAIGES等［9］采用自适应有限元网格和降阶模型使计算效率提高11倍；INNERBERGER等［10］通过对三角形网格的细化，更加快速地求解了非线性偏微分方程的迭代线性化方法引起的问题。许多研究人员通过网格自适应的方法来对有限元网格数量进行控制，从而减轻计算负担。

基于改进四叉树和比例边界有限元法的自适应设计域拓扑优化方法——王  浩  王江北  罗浩东等

中国机械工程 第35卷 第5期 2024年5月

网格自适应算法中，四叉树算法的灵活性强、计算效率高、精度高，最早被应用于计算机图像学领域，是一种图像压缩算法，后被迁移到工程计算中，主要应用在求解尖端裂纹、应力强度因子［11-14］等方面。孙立国等［15］将水平集函数与图像四叉树法相结合对网格进行自适应划分，求解了复合型应力强度因子； GRAVENKAMP等［16］将四叉树与比例边界有限元法（scaled boundary finite element method， SBFEM）相结合对复杂几何形状的静态、谐波、模态和瞬态分析等的数值算例进行了分析。CORNEJO等［17］利用自适应网格与有限元法和离散元法耦合的方法，实现了在变量场曲率较大处细化网格，并对裂纹的扩展进行了分析。WIJESINGHE等［18］将数字图像网格与SBFEM相结合，用于岩土边坡稳定性分析的数值技术研究。

在断裂力学中裂纹源位置是已知的，四叉树算法可以准确地在裂纹尖端处进行网格自适应划分。拓扑优化过程中生成的拓扑结构具有未知性，传统四叉树算法无法准确得到所需细化网格的位置，因此在拓扑优化过程中传统的四叉树细化策略失效。

改进四叉树算法将拓扑优化结构中的密度值作为细化策略的判定条件，从而应用到拓扑优化中。细化结构边界处的网格，可以大幅度提高拓扑优化结构边界处的精度。经改进四叉树算法细化后的拓扑优化结构具有不同等级的网格单元，首先通过SBFEM更新整体结构刚度矩阵，然后进行迭代优化。

本文将改进四叉树算法和SBFEM相结合自适应地更新设计域内网格，提出了一种基于固体各向同性材料惩罚模型（solid isotropic material with penalization， SIMP）方法的自适应设计域（adaptive design domain，ADD）拓扑优化算法（以下简称“ADD-SIMP方法”）。本文详细介绍了ADD-SIMP方法，给出了SBFEM的推导过程，并通过两个数值案例对所提方法进行了验证。ADD-SIMP方法可以通过更少的初始网格数量获得边界处相同精度的求解结果，大幅度降低了计算成本，可为大型结构的拓扑优化提供参考。

1  ADD-SIMP模型

ADD-SIMP方法的整体框架如图1所示，具体流程如下：

（1）给定结构几何参数、材料参数、优化目标和约束等条件。

（2）根据需求设置初始网格划分以及最小网格单元等级等参数。

（3）有限元计算得到每个单元的位移场、刚度矩阵等有限元参数，计算结构柔度和灵敏度，并对单元灵敏度进行更新。

（4）通过优化准则（optimality criteria，OC）进行初步优化，当结构初步稳定时执行步骤（5）。

（5）通过改进四叉树算法对结构边界网格进行自动加密，以达到精确结构边界位置的目的；迭代步驟（2）～步骤（5），直到结构最小网格单元等级满足步骤（1）中规定的最小网格单元等级时不再对网格进行细化。

（6）对该结构继续进行拓扑优化，直至最终结构满足收敛条件后输出该结构。

1.1  SIMP方法

固体各向同性材料惩罚模型（SIMP）［19］方法因其简单易懂、灵活性强、可扩展性好等优点成为目前主流的拓扑优化方法之一。但SIMP法生成的拓扑结构往往会存在大量的伪密度单元，如图2所示，在结构边界处伪密度单元数量约为实体单元数量的3倍，这会致结构边界处的求解精度受到限制，反映了较低的精度水平。精度要求高的大型结构就需要更多数量的网格来提高精度。

本文基于将结构的体积作为约束、最小柔度作为目标的拓扑优化问题进行了网格自适应的研究，其中优化问题的表达式如下：

其中优化问题的表达式为

find X=（x1，x2，…，xN）T∈D

min C=∑Ni=1xp0iuTikiui

s.t. V=φV0=∑Ni=1xiVi

0

式中，X为结构设计域内单元的相对密度；xi、ui、ki、Vi分别为第i个单元的相对密度、位移矢量、刚度矩阵、体积；N为设计域内单元的总数；D为结构总设计域；C为结构柔度；p0为惩罚因子，工程上一般取p0= 3；V为结构总体积；V0为结构初始体积；φ为目标体积分数；xmin为单元相对密度值的下限，为了计算稳定，防止出现奇异性，一般设置xmin=0.001。

1.2  改进四叉树算法

四叉树是一种数据结构，主要应用于计算机图像学领域，本文引入一种改进四叉树算法对拓扑优化过程中的网格进行自适应划分。通过采用新的细化策略、网格单元和节点存储方式，实现拓扑优化过程中的网格自适应控制，以提高计算效率并提高边界处的精度，其中网格单元储存于单元树（cellMap）中，用于储存结构网格单元的材料、位置等信息，具体参数含义详见表1；节点存储于节点树（nodeMap）中，用于存储节点位置和编号信息，具体参数含义详见表2。两个数据结构共同完成结构拓扑优化过程的网格实时更新。

1.2.1  网格单元与节点存储

图3所示为单元细化后新增子单元与新增节点情况，其中新增网格单元顺序如图3a所示，新增节点顺序如图3b所示，角节点的序号不变，将新增节点作为新节点插入节点树中。

每个网格单元有唯一的标识符与之对应，单元标识符的计算公式依靠单元结构中的单元等级lc和单元索引ic共同决定。lc=0表示该网格为初始网格，lc=h表示初始网格划分h次后的h级网格，网格等级每增加1级，网格单元尺寸则会变为原来的1/2，边界处的精度也随之提高。

每次划分每个单元时最多可以划分为4lc个单元，为防止重复单元标识符出现，假设每级单元完全细化，并在此基础上计算新网格单元编号，因此对于一个n行m列的结构中当前网格单元等级为h、单元索引为（x，y）的单元，其唯一标识符kc的计算表达式如下：

kc=2h（y-1）m+x+∑hlc=14（lc－1）mn （2）

细化后出现的单元索引可由下式计算得到：

ic，pa=（p，q）

ic，ch（1，1）=（2p－1，2q－1）

ic，ch（1，2）=（2p－1，2q）

ic，ch（2，1）=（2p，2q－1）

ic，ch（2，2）=（2p，2q）（3）

其中，p、q为父级单元的索引ic，pa。得到子网格的索引ic，ch后通过式（2）计算出子网格的编号。

需要说明的是，节点唯一标识符kn由 x、y方向上的标识符kn，x和kn，y构成，其计算由给定结构的尺寸参数、最小细化等级共同决定。最小细化等级为初始网格大小和设定的最小单元尺寸共同决定，单个网格边界上最多可以划分2lc，max+1个节点（其中lc，max为当前最高单元等级），单元真实最小尺寸smin为单元尺寸sc除以2lc，max，根据单元位置坐标Pn（Pn，x，Pn，y）可计算得到节点的kn，x和kn，y，则节点唯一标识符kn通过下式确定：

sc2lc，max=smin  lc，max∈Z

kn，x=floor（Pn，xsmin）

kn，y=floor（Pn，ysmin）

kn=kn，x_kn，y（4）

其中，floor（·）表示向下取整函数。

由于节点标识符不具有连续性，因此节点标识符的作用仅为确定设计域中的某个唯一节点，为后续计算的顺利进行需要引入连续的节点编号nn作进一步处理，具体流程如下：

（1）计算新增节点的唯一标识符kn。根据新增节点的位置信息，利用式（4）可以计算出该节点的kn。

（2）依据步骤（1）中计算出的kn，在节点树中搜索出对应的节点。

（3）若节点树中不存在相应kn的节点（即返回值为空），则新增节点，赋予其节点编号（即当前节点树中最大节点编号值加1），并加入节点树；若已存在，则保持原节点所有信息不变。

有限元计算需要将网格单元中的节点按逆时针排序，因此网格单元通过左下角节点的位置进行定位，如图4中网格单元3的节点顺序为｛1，2，5，3，6，4｝。

图4中，节点5、6（仅列举出了图中标出节点）为T型节点，SBFEM与FEM计算得到的位移场的误差最大为8%，最小为0.9%。同时为防止某个网格单元具有过多的T型节点而导致计算数值不稳定，因此新增网格单元等级控制策略，其具体步骤如下：

（1）读取当前单元的ic，并计算得到相邻所有单元索引的集合。

（2）根据式（2）计算出相邻单元的kc值。

（3）读取单元樹中对应kc值下的单元信息，若返回空数组则表示单元不存在，通过计算得出父级单元的ic，pa并通过式（2）计算出kc，pa值，读取单元树对应的网格单元信息，若仍返回空数组则表示单元等级差超过1，并对父级单元的上级单元进行网格划分。父级单元的索引ic，pa可表示为

ic，pa=k    ic∈2k    k=1，2，…

k+1ic∈2k+1k=1，2，…（5）

1.2.2  细化策略

当拓扑优化中两次输出结构的单元最大改变量小于设定阈值时对结构进行细分，该阈值为细化阈值，为网格是否需要细化的判定条件。细化阈值的选取影响了网格细化时主承力结构是否出现。为保证拓扑优化结构的稳定，将细化策略分为网格细化和网格检查两步。

网格细化的目标单元主要是生成的拓扑优化中所有伪密度单元和已经稳定的边界单元，其余单元不参与网格的细化。通过网格细化可以获得不同等级的网格单元，从而大幅度减少网格数量，提高网格的利用率和计算效率。网格细化后的效果如图5所示。

网格检查的目的是对细化后的整体网格单元进行检查，对网格细化后不满足要求的网格单元进行细化。网格检查的目标为当前最大单元等级的上两级单元。网格检查步骤可以有效缓解锯齿边界的出现，锯齿边界如图6所示。

锯齿边界主要是由网格细化后主承力结构在后续变化中移动至粗网格边界处所产生的。锯齿边界出现的原因在于有限元方法由传统有限元法替换为SBFEM后，灵敏度的改变引起优化准则出现不稳定现象。

细化策略中网格细化的具体步骤包括：

（1）从单元树中读取所有网格单元信息，提取的目标单元为所有未细化的网格单元。

（2）对未细化单元中所有伪密度单元进行细化。将惩罚因子p0设置为3，当单元相对密度x=0.1时，被惩罚后参与计算的实际相对密度值约为0.001；当x=0.8时，被惩罚后的相对密度值约为0.5，因此过程中参与细化的网格单元的单元相对密度值位于0.1～0.8。

（3）对于单元相对密度值大于0.8的网格单元，读取与该单元相邻的4个单元的密度值进行判断，若相邻单元中存在单元相对密度值小于0.8的网格单元，表明该单元处在结构边界处，则对单元进行网格细化。

网格检查的具体步骤包括：

（1）读取所有处于当前网格最高等级上一级单元的网格，图7中灰色单元为当前最高等级单元，红色单元为上一级网格单元。

（2）获取中心单元所在九宫格中所有的网格单元，计算9个单元的平均相对密度值。

（3）若平均相对密度值小于1且大于0.1，表明该单元为锯齿边界处单元，则对该单元进行细化。

图8为经过网格检查后，结构经第二次网格细化之后的效果图，可以看出锯齿边界得到明显改善。

1.3  灵敏度分析与网格无关滤波器

对于细化后的网格，灵敏度更新策略与传统方式不同。在传统有限元中以每个网格单元为单位，其刚度矩阵为常数，因此每个单元对整体结构的灵敏度可以直接求解。随着网格被细化后悬挂节点的出现，传统的有限元算法失效，因此需要引入SBFEM来计算刚度矩阵等有限元信息。

SBFEM中每个网格由若干个线单元构成，细化后会产生不同等级单元相互接触的情况，如图9所示，其中3号单元包含6个线单元，而4号单元只包含4个线单元，计算得到3号单元的灵敏度要高于4号单元的灵敏度，因此由多个单元组成的网格的灵敏度要高于其他四边形网格的灵敏度。

如图10所示，单元柔度C与单元长度无关，与线单元的数量nl有关。

单元灵敏度受到网格中单元个数的显著影响，为降低由于线单元数量不同所带来的影响，进行了归一化处理。具体操作是，将所有包含线单元的网格统一归一化为四边形单元。这样的处理旨在降低不同数量的线单元对单元灵敏度的影响。为此，引入了一个新的修正公式来更加准确地评估单元灵敏度，以减少评估结果因网格中线单元的数量差异而产生的偏差，修正公式如下：

d=p0xp0－1C［1－（1－4nl）2］=8p0xp0－1Cnl－2n2l（6）

（a）C=5.7847  （b）C=5.7847，nl=4  （c）C=6.0104，nl=5

（d）C=6.4386，nl=6  （e）C=7.8104，nl=7

results of different type of cells

其中，为修正后的单元柔度，修正后的值如表3所示。

为防止出现棋盘网格现象，需要对过滤半径rmin内所有网格单元进行灵敏度过滤处理。为节省计算资源，从中心单元j开始，搜索以单元j为中心、过滤半径rmin内所有未划分的网格单元，读取rmin内所有网格单元的灵敏度，并根据下式对中心網格单元的灵敏度进行更新：

dc^j=1xj∑Mm=1H^m∑Mm=1H^mxmcmxm（7）

H^m=max（0，rmin－dist（j，m））

{M∈N|dist（j，m）≤rmin}

式中，c^j为修正后中心单元j的柔度；xj为当前中心单元j的相对密度；M为过滤半径rmin内的单元个数；xm为过滤半径内第m个单元（以下简称m单元）的相对密度；cmxm为m单元的灵敏度；H^m为m单元对中心单元的影响系数；dist（j，m）表示m单元与中心单元j的距离，距离越大，m单元对中心单元的影响系数越小。

搜索网格单元时可分为如下两种情况。

（1）高等级网格单元为中心时的情况如图11a所示，当将一个高等级（更细分）的网格单元作为中心进行搜索时，过滤半径内可能会发现一些等级低于中心单元（尺寸较大）等级的网格单元。若搜索时与中心单元等级相同的邻近单元1不存在，则会返回空结构数组，此时需要通过式（5）来计算获取网格单元1的父级单元索引，并调用父级网格单元的灵敏度来参与计算，以确保考虑到高等级中心网格单元周围较大尺寸单元的影响。

（2）低等级网格单元为中心时的情况如图11b所示，当将一个低等级（即尺寸较大）的网格单元作为中心时，过滤半径内可能存在一些高等级的网格单元。与第一种情况不同，这时由于邻近单元（以灰色表示）实际存在，搜索不会返回空结构数组。在这种场景下，需要增加额外的判断条件来检查rmin半径内的网格单元是否已经被细分。如果这些网格单元已经被细分，则需要通过式（3）计算获取子网格的单元索引，并调用子网格单元的灵敏度来参与计算，以确保在以较大尺寸单元为中心时能够适当考虑到周围被细分的较小网格单元的影响。

1.4  网格自适应流程

ADD-SIMP方法中的网格自适应主要是通过改进四叉树算法来对边界处的网格进行自动细化，并进行拓扑变量的更新。具体流程如下：

（1）检查细化条件。判断是否满足网格细化的阈值条件。

（2）网格细化。若满足细化条件则将伪密度单元划分为四个小网格单元，并传递拓扑信息。

（3）网格检查。对细化后的网格进行检查，处理非正常网格。

（4）更新设计变量。将新生成的小网格与未细化的大网格组合，并参与后续拓扑优化。

（5）检查最小网格要求。判断是否满足最小网格要求，若不满足则返回步骤（1）。

（6）优化输出结构。若满足最小网格要求则完成优化并输出结构。

2  SBFEM方法

2.1  理论推导

由于网格局部划分，新增的悬挂节点使传统的全解析有限元计算方法失效。SBFEM是一种半解析算法［20］。传统有限元方法通过使用所有角节点进行网格单元的拼装，而SBFEM则通过将四边形单元离散为四个线单元，并将具有悬挂节点的线单元分割成两个线单元的组合，从而避免了悬挂节点处无法计算的问题，因此相较于传统的有限元方法，SBFEM更适用于改进四叉树算法对网格进行划分后的有限元计算。

对于以柔度为目标的拓扑优化流程的有限元分析，主要需要求解单元上的位移场和单元刚度矩阵，因此面临的主要问题是求解伪密度单元上的应变场和具有悬挂节点的单元刚度矩阵的问题。

笛卡儿坐标系中（x，y）与自然坐标系（ξ，η）的转换关系如图12所示，其中位于边界上的节点B（xB，yB）所对应的自然坐标系下的计算表达式如下：

xB=ξx（η）=ξN（η）x

yB=ξy（η）=ξN（η）y（8）

其中，N（η）表示2节点线单元形函数，η为边界插值坐标；ξ为放缩系数，当ξ=1时表示节点B位于边界上，ξ=0表示节点B位于缩放中心O处。

SBFEM中位移方程u在自然坐标系下的表达式如下：

u（ξ，η）=ξN（η）u（9）

应变位移关系如下：

ε（ξ，η）=Lu（ξ，η）（10）

其中，L为微分算子矩阵，ε为应变方程。则线单元上应力的表达式如下：

σ（ξ，η）=Dε（ξ，η）（11）

其中，D为刚度系数矩阵。

根据虚功原理可以得到控制方程，其中E0、E1、E2为与线单元几何有关的常系数矩阵。设式中自由度为n，通过引入变量式，将控制方程变为自由度为2n的一阶常微分方程，具体表达式如下：

E0ξ2u″（ξ）+（E0+ET1－E1）ξu′（ξ）－E2u（ξ）=0

X（ξ）=u（ξ）q（ξ）

ξX′（ξ）=－ZX（ξ）

Z=E－10ET1－E－10－E2+E1E－10ET1－E1E－10（12）

式中，u（ξ）为节点位移矢量；q（ξ）为节点力矢量；u′（ξ）、u″（ξ）分别为u对变量ξ的一阶偏导和二阶偏导;X′（ξ）为X（ξ）的一阶偏导;Z为哈密顿矩阵。

通过求解哈密顿矩阵的特征值λ和特征向量矩阵Φ，再结合边界条件中的特殊点，可以求解得到积分常数c，其表达式如下：

X（ξ）=Φuλuc

Φqλqc（13）

结合式（12）和式（13），可得等效节点力矢量方程为

q（ξ）=ΦqΦ－1uu（ξ）=Keu（ξ）（14）

其中，Φu、Φq和λu、λq分别为与u和q相关的特征向量矩阵和特征值；Ke为单元刚度矩阵。根据式（14），当ξ =1时可以得到边界处力与位移的关系，则宏观结构上力与位移的关系表达式为

F=KU（15）

其中，F为结构外载荷矢量；U为结构位移矢量；K为总体刚度矩阵。结合式（14）和式（15），可得总体刚度矩阵的计算公式为

K=ΦqΦ－1u（16）

则SIMP插值后的單元刚度矩阵可表示为

Ke=ρp0eΦqΦ－1u

式中，ρe为单元密度。

2.2  有限元信息更新

在SIMP方法中，所有网格都要通过有限元计算获得应力场等其他有限元信息，其中关键是得到每个网格单元的刚度矩阵。经过改进四叉树算法细化后的结构会增加大量具有悬挂节点的网格单元，因此，每次经过网格细化后的整体结构刚度矩阵需要得到实时更新。

SBFEM算法所需要的数据包括：①全部节点的坐标信息；②全部节点的唯一编号信息；③全部网格单元上的节点信息；④网格节点按逆时针排序信息。以上所有信息均由改进四叉树算法进行存储，以便后续调用。具体的操作如下：

（1）从单元树中读取所有未细化单元中的单元密度信息以及包含全部按逆序排序的节点编号信息。

（2）从节点树中读取所有节点位置信息与节点编号信息。

（3）通过式（15）计算获得单元刚度矩阵并将其保存。

（4）结合步骤（1）中获得的单元密度信息，依据式（16）计算获得经过惩罚后的单元刚度矩阵。

（5）先对结构整体刚度矩阵进行组装，后计算所有节点上的位移。

3  数值案例和分析

本节通过两个案例对ADD-SIMP与SIMP方法在计算时间、迭代次数、所用网格数、结构边界的清晰度几个方面进行了对比研究，并对其中主要参数的影响规律进行了分析，其中包含初始网格数量、细化阈值、过滤半径。所有计算均在同一台计算机上进行，计算机参数如表4所示。

3.1  悬臂梁案例分析

如图13所示，案例一为一个左端固支的悬臂梁，在右端最下角节点处施加外载荷为1 N的外力，悬臂梁的尺寸为120 mm×40 mm，弹性模量为1 MPa，泊松比为0.3。

SIMP方法的优化结果如图14所示，当该方法获得1 mm左右结构精度时，网格划分至少应为480×160，网格数量为76 800个，单个网格尺寸为0.25 mm，此时受伪密度单元影响，结构边界处的精度约为0.75 mm，优化后的结果如图14a所示，计算时间为60 642 s（约17 h）。

ADD-SIMP方法的结果由不同等级网格单元组成，所有网格线均显示，初始网格设置为120×40、初始网格尺寸为1 mm、最小网格单元等级为2、细化阈值设置为0.08的最终优化结果如图15所示。最终所用网格总数为26 205，计算时间为711 s（约12 min），最小网格的尺寸为0.25 mm。

在保持原始结构比例不变的前提下，SIMP方法需要将设计区域划分为270×90共24 300个网格，此时与ADD-SIMP方法所用网格数量相近。在270×90的结构中，每个网格的尺寸约为0.44 mm。SIMP方法最终的优化结果如图14b所示，整个计算过程耗时为6644 s，是采用ADD-SIMP方法所需时间的9倍。SIMP方法的最小网格尺寸为采用ADD-SIMP方法时的两倍左右，因此ADD-SIMP方法具有更优的边界精度。

结合图14和图15可以看出，利用ADD-SIMP方法得到的最终结构与SIMP方法得到的120×40主结构相似（图14c），随着网格数量的增加，ADD-SIMP方法兼顾了270×90与480×160网格中的主要分支结构，减少了细小分支数量，具有更好的效果。

从图16a中可以看出，SIMP方法优化结果中网格尺寸相同，结构中已稳定区域（实体单元和空白区域）充斥着大量的网格单元，但是对边界处的精度并没有贡献，因此网格利用率低；相比之下，在图16b所示的ADD-SIMP方法结果中，整体结构由三级尺寸不同的网格构成，已稳定区域的网格尺寸大且稀疏，结构边界处网格尺寸小且密集。在ADD-SIMP方法中，结构边界处的最小网格尺寸与SIMP方法中相同，因此两者在边界处的精度几乎没有差别。然而，ADD-SIMP方法在构建整体结构时所需的网格总数显著少于SIMP方法所需的网格总数，表明ADD-SIMP方法在提高计算效率的同时保持了边界精度，优化了网格的利用率。

从表5中可以看出，ADD-SIMP方法在达到相同精度水平时，最终所用的网格数量为26 205个，相较于SIMP方法在480×160结构所需的网格数量，ADD-SIMP方法仅用了其1/3的网格数量。两种方法在针对结构最终目标函数的优化结果上相似，但是在边界处精度相同的条件下，ADD-SIMP方法的计算效率显著提高，可达到SIMP方法计算效率的85倍。此外，当与具有相似网格数量的SIMP方法相比时，ADD-SIMP方法的计算效率也有显著提高，可达到原来的9倍。值得注意的是，在这种比较条件下，ADD-SIMP方法边界处的精度提高至原来的两倍。这些结果表明了ADD-SIMP方法在提高计算效率方面具有优势。

3.2  简支梁

如图17所示，案例二为一个两端固支的简支梁，在中点处施加外载荷为1N的外力，悬臂梁的尺寸为240 mm×40 mm，弹性模量为1 MPa，泊松比为0.3。

对于ADD-SIMP方法，初始网格划分为120×20共2400个初始网格，最小网格单元等级设置为2，即最小网格的尺寸为0.5 mm，细化阈值设置为0.06；为保证相同精度，传统SIMP的网格划分应设置为480×80共38 400个网格。两种方法的拓扑优化惩罚因子和过滤半径均选取为3。ADD-SIMP方法经过30次迭代后的优化最终效果如图18a所示。

ADD-SIMP方法计算所用时间为217.6 s（约0.06 h），SIMP方法经过201次迭代，所用时间为24 019.2 s（约6.67 h），ADD-SIMP方法的计算效率为原来的110倍，对于同样迭代次数的SIMP方法，计算所用时间仍需要3584.9 s；在目标函数的优化上，ADD-SIMP方法的目标函数值小于SIMP方法的目标函数值；ADD-SIMP方法最终所用网格数量为15 231个，仅用原方法的1/2网格数量即可得到结构与精度相似的最终结构。

此外从图18b中可以看出，SIMP方法得到的结果中细枝结构仍较多，ADD-SIMP方法得到的结构在一定程度上抑制了细小分支的出现。

3.3  参数分析

3.3.1  初始网格数量

对于ADD-SIMP方法，将初始网格划分为30×10、60×20、120×40三组进行分析对比，如表6和表7所示。为保证最小网格单元尺寸相同，均设置为0.25 mm，因此上述三组的最小网格单元等级设置分别为4、3、2。

从表6中可以看出，三种初始网格划分情况下所得到的最终网格数量均在26 000个左右，相对于原来的76 800个网格，所需网格数量仅为原来的1/3。

从表7中可以看出，初始网格数量为60×20的结构计算时间最短，为556.0 s（约9 min），时间仅为SIMP方法的1/109，计算效率提高了100倍以上；三组的目标函数值相近，SIMP方法最终目标函数值为182，因此使用ADD-SIMP方法可能导致目标函数值略微偏高；但是三种不同初始网格划分的情况下，最终的目标函数值相差不大，可知初始网格数量对结构轻量化程度的影响较小；然而，初始网格数量对循环次数的影响显著，计算时间会随着初始网格数量的增加先减少后增加。

表8～表10分别展示了不同初始网格数量下结构的优化过程，可以看出，随着初始网格数量由少变多，结构在第一次进行网格细化时所处状态不同。初始网格过少时，一次细化的结构并未出现主承力结构，而是处于乱序状态，直到后续网格数量增加才出现稳定的主承力结构，并趋向于稳定。当初始网格数量变多，结构首次细化时就

已经具有相对清晰的主承力结构，但是第一次细化所处的循环数也更大，进入第一次网格细化的时间会被延长。此外大量的初始网格数会直接增加计算负担，最终导致计算时间变长。

从表8中还可以发现，当初始网格数量较少时，结构中出现的承力结构数量少、尺寸大，但是随着后续网格数量的增加，结构有出现小分支的趋势，因此ADD-SIMP方法在一定程度上具有减少结构细小分支的作用。

综上可知，当初始网格数量选择合适时，不仅可以大幅度缩短计算时间，還可以达到减少结构细小分支的效果。

3.3.2  细化阈值

在ADD-SIMP方法中，细化阈值是作为网格细化策略的关键参数，用于决定何时对结构的网格进行细化。具体来说，该阈值设定了一个标准：只有当结构中网格单元的密度变化最大值不超过此阈值时，才会执行网格细化操作。由此可知，细化阈值的设定直接影响到拓扑优化过程中网格细化的时机。选择合适的细化阈值是确保优化过程既高效又精确的关键因素。

在60×20的初始网格和最小单元等级为3的前提下，选择细化阈值为0.06、0.08、0.1三个数值进行对比分析，如表11所示，可以看出，细化阈值设置取值越大，第一次进行网格细化时所处的循环数越小，则计算时间越短；但是细化阈值设置过大时可能导致结构主承力结构并未出现（如表12中0.1细化阈值下第一次细化结果），需要配合多次细化才能保证结果的准确性；若细化阈值取值偏小，则会延长结构进行细化的时间，从而降低计算效率；若细化阈值取值过小，甚至不会进入网格细化操作。

3.3.3  过滤半径

过滤半径为SIMP方法中的主要影响参数之一，分别设计2、2.5、3、3.5、4五组过滤半径，分析过滤半径对ADD-SIMP方法的影响规律。对案例一初始网格数量的研究中，可以发现初始网格为60×20、细化阈值为0.08、最小网格单元等级设置为3的一组计算时间最短。

从表13中可以看出，随着过滤半径的增大，当过滤半径为2.5～4之间时迭代次数先减少后增加，计算时间也是先增加后减少随后又快速增加，最终网格数量整体呈现增加趋势。

由于过滤半径影响每次迭代中每个单元上的计算成本，则过滤半径选择过大会导致每个单元上灵敏度计算过滤的成本增加，但当过滤半径为3左右时，迭代次数会减少，从而导致最终的计算时间缩短。此外由于过滤半径内所有单元都会参与中心单元的灵敏度计算，因此导致过滤半径内的每个单元都会获得部分周围实体单元的相对密度，从而导致边界处的伪密度单元数量增加，进而影响ADD-SIMP方法细化后的网格数量。

从表14中可以看出，过滤半径对最终结构会产生较大的影响，过滤半径选择较大时，靠近约束端的细小分支结构会被整合，同时伪密度单元数量会增加；当过滤半径过小时，会具有更多的细小分支，但是边界处的伪密度单元数量会减少。

4  结论

对于大型结构拓扑优化方法计算成本高和固体各向同性材料惩罚模型（SIMP）方法在优化后结构边界模糊的问题，通过改进四叉树算法、比例边界有限元法（SBFEM）和自定义细化网格的策略，本文提出了基于SIMP方法的自适应设计域（ADD）拓扑优化算法（以下简称“ADD-SIMP方法”），并对其中的参数设置进行了分析，在悬臂梁和简支梁上进行了应用。从本文应用的案例中可得出如下结论：

（1）改进四叉树算法通过定义细化策略可以实现在拓扑优化过程中对结构的边界自动进行划分网格操作，从而提高最终结构在边界处的计算精度。

（2）若设置最小网格单元等级为2时，仅需要1/3左右的初始网格数量即可得到同边界精度的优化结果。

（3）ADD-SIMP方法得到的结构分支较少，可制造性更高，在边界处的平滑程度和精度与传统方法相似。

（4）ADD-SIMP方法的在计算成本方面远优于传统方法，最终网格数量仅为原来的1/3，最快的计算速度效率是原始SIMP方法的100倍。

对于网格数量更多、结构更庞大、计算需要数天的大型结构，采用ADD-SIMP方法可以节省更多的计算资源、缩短设计周期。

参考文献：

［1］  郭西园， 卫钟可， 杜建国， 等. 基于3D打印的减速箱体拓扑优化及应用［J］. 航天制造技术， 2022， 235（5）：59-61.

GUO Xiyuan， WEI Zhongke， DU Jianguo， et al. Topology Optimization and Application of Gearbox Based on 3D Printing［J］. Aerospace Manufacturing Technology， 2022， 235（5）：59-61.

［2］  孟亮， 仲明哲， 李文彪， 等. 面向增材制造的航空发动机外部系统支架拓扑优化设计［J］. 中国机械工程， 2022， 33（23）：2822-2832.

MENG Liang， ZHONG Mingzhe， LI Wenbiao， et al. Topology Optimization Design of Aero-engine External System Brackets for Additive Manufacturing［J］. China Mechanical Engineering， 2022， 33（23）：2822-2832.

［3］  梁健， 李曉杰， 谢炯. 拓扑优化在工业机器人造型设计中的应用［J］. 机械设计， 2020， 37（9）：128-133.

LIANG Jian， LI Xiaojie， XIE Jiong. Application of Topology Optimization in Industrial Robot Modeling Design［J］. Journal of Machine Design， 2020， 37（9）：128-133.

［4］  朱鹏程， 杨志勋， 阎军， 等. 基于拓扑优化方法的海洋柔性管缆限弯器新型设计研究［J］. 船舶力学， 2023， 27（3）：415-426.

ZHU Pengcheng， YANG Zhixun， YAN Jun， et al. New Configuration Design Methodology of Bend Restrictor of Marine Flexible Pipes and Cables Based on Topology Optimization［J］. Journal of Ship Mechanics， 2023， 27（3）：415-426.

［5］  刘英杰， 胡强， 赵新明， 等. 汽车发动机连接支架拓扑优化及增材制造研究［J］. 中国机械工程， 2023， 34（18）：2238-2247

LIU Yingjie， HU Qiang， ZHAO Xinming， et al. Research on Topology Optimization and Additive Manufacturing of Automotive Engine Connection Brackets［J］. China Mechanical Engineering， 2023， 34（18）：2238-2247