《1.引言》
1.引言
随着计算机科学技术的飞速发展，结构优化设计已成为获得轻量化和高性能结构的最重要手段之一。通常，根据设计变量的类型，结构优化分为尺寸优化、形状优化和拓扑优化三大类。其中，拓扑优化因无需初始构型，可获得意想不到的创新设计结果，被认为是最通用的优化设计方法。本质上，拓扑优化使用优化技术以发现材料在设计域中的最优位置，使得结构在满足一定约束条件下特定性能指标最优。经过40多年的发展，拓扑优化已成功应用于重大工程结构的设计中，如航空航天、汽车和生物医学等领域[1,2]。现阶段，主要的几种拓扑优化方法包括变密度方法[3,4]、水平集方法[5,6]、进化方法[7]等。有兴趣的读者可以参阅参考文献[8]以了解更多细节。在这些方法中，变密度法采用单元伪密度来描述结构拓扑，具有较高的计算效率和稳定性，是目前发展最成熟且应用最为广泛的技术。
增材制造技术（AM），俗称为3D打印，通过层层累积材料制造结构零件[9–11]。这种制造工艺推翻了传统的“减材”“等材”的制造理念，使得具备复杂的几何特征的创新结构的制备成为可能[12–16]。因此，通过使用增材制造工艺对拓扑优化设计结果进行快速试制，然后进行实验来评估结构性能，可有效地缩短设计周期[17–19]。然而，使用变密度方法获得的拓扑优化结果是以单元密度来描述的，缺乏基本的几何特征，如点、线、面和体。这种描述方法通常无法直接制造，且难以用于如形状优化的后续研究中。因此，需要引入人工干预来解释拓扑优化结果以获得计算机辅助设计（CAD）模型，但是此过程往往耗时耗力[2]。因此，迫切需要研究如何将拓扑优化结果快速有效地转换成可进一步修正和制造的CAD模型。
许多文件格式，如STL格式、初始图形交换规范（IGES格式）等被开发出来，用于描述适用于不同制造技术的CAD模型。其中，STL格式由3D Systems公司[20]首先提出，以其简单的格式广泛应用于快速原型制造领域，这种文件格式采用三角面片表达实体表面数据，每个三角形面片用三角形法向量的三个分量及三角形的三个顶点坐标来表示，不包含模型颜色、纹理或其他常见的CAD模型属性等。现阶段，许多其他软件包（如AutoCAD、Solid-Designer和Unigraphics）都可以读取和导入此文件格式，以进行快速的原型制作、3D打印和计算机辅助制造。IGES是用于描述参数化CAD模型的文件格式之一[21]，是不同CAD软件程序（如Maya、Pro/ENGINEER、Softimage和CATIA）之间信息交换的标准格式。现阶段，大多数商用CAD软件都可以将IGES文件转换为用于增材制造的STL文件，然而无法从STL文件获得IGES文件。
文献中已经提出了许多方法来从拓扑优化结果获得CAD模型。其中，密度等值线法是现有文献中使用最广泛的方法之一，该方法通过使用提取密度场的等值线提取拓扑优化结果边界[22–24]。但是，这种方法需要通过多次测试，才能确定适当的密度阈值以获得有效的结果。此外，对于几何模型较为复杂的拓扑优化结果[25]，采用该方法获得的CAD模型中可能会存在非常薄的组件或者不连续的结构部分（孤岛），并且表面较为粗糙，因此后续需要多次重复的模型修订和人为干预来获得理想的CAD模型。除此方法外，其他一些研究人员尝试使用基本形状模板[26,27]，如具有特定参数的球、圆柱体和矩形，通过拟合手段来获得拓扑优化结果的近似模型。由于该方法中只能使用有限数量的形状模板，因此研究者通常将研究重点放在如何提高拟合精度上[28,29]。拓扑优化结果的其他几何重建方法主要基于插值函数，如B样条曲线、二次曲面样条曲线、T样条曲线和非均匀有理B样条[30–35]。这些方法适用于任意复杂的表面，并且可以无缝地连接到IGES文件。然而，使用这种方法需要事先准备稠密的点云，这使得边界拟合的计算量过大。Yi和im[36]最近提出了使用基本的参数化特征，如直线、圆弧、圆和圆角等，来确定拓扑优化结果的边界。与其他方法相比，这种方法处理具有复杂边界形状的优化结果时，具有相对适中的拟合变量数量，但是此方法需要事先得到一个不包括非常薄的部分和孤岛结构且可用0和1表示的拓扑优化结果。
本文提出了一种可自动将拓扑优化结果转换为能直接3D打印的STL文件和以及方便进一步精细化设计的参数化CAD模型（IGES文件）的方法。该方法的流程图如图1所示。首先，为了避免产生非常薄的组件和孤岛现象，该方法首先提出一种识别拓扑优化结果微小特征的策略。其次，为了获得一组密集的边界点，在保证结构基本特征的前提下，提出一种边界加密化过程。基于此，可获得拓扑优化结果的STL文件。最后，本文提出一种自适应B样条拟合方法来获得一个边界光滑且可用于形状优化的参数化CAD模型。这种将拓扑优化设计转换为STL或IGES格式的新实现方法可以弥合设计和制造之间的间隙，为应用AM技术实现拓扑优化结果的快速试生产以及使用形状优化进行进一步精细化设计创造了条件，有望缩短生产周期所需的时间。
《图1》
图1. 所示方法的流程图
本文的组织结构如下：第2节提供了基于变密度法的拓扑优化的基本概述，在第3节和第4节中，以悬臂梁为例介绍了新方法的具体实现过程，在第5节中，展示了其他两个数值案例来说明所提出方法的鲁棒性。
《2.基于变密度法的拓扑优化》
2.基于变密度法的拓扑优化
《2.1.概览》
2.1.概览
拓扑优化试图回答在规定的域内材料的最佳分布是什么的问题，最早可追溯到Bendsøe和Kikuchi[37]提出的基于均匀化方法的结构拓扑优化。在此基础上， Bendsøe [38]与Zhou和Rozvany [3]提出了固体各向同性材料惩罚（SIMP）理论，这被视为拓扑优化发展历史上的里程碑，也是变密度法的基本雏形。SIMP方法的基本思想是通过有限元网格对设计域进行离散化，通过优化单元密度ρi来判断该单元有无材料，有材料定义其单元密度为1，无材料则定义单元密度为0。图2显示了一个结构的拓扑模型，其中黑色元素（ρi = 1）代表材料， 白色元素（ρi = 0）代表空白。
《图2》
图2. 基于单元密度的拓扑优化示意图。
单元密度的数学表达可表示为如下：
式中，Ωi表示第i个单元的域；Ω为设计域；Ωmat表示优化后的结构拓扑。为了解决这个离散优化问题，第i个元素的弹性模量Ei定义如下：
式中，p为惩罚因子。当p=1时，这个问题与Cheng和Olhoff[39]首次研究的经典“变厚度板优化”问题相对应。研究表明，此时会出现较多的灰色单元，即0<ρi<1，这些单元是没有物理意义的。因此，引入了p>1以产生离散的“0-1”设计。应该指出的是，选择太低或太高的p值分别会导致优化结果中存在太多的灰度单元或者优化模型过快收敛到局部最小值，因此，一个适当的p值是有必要的。大量的数值算例表明，当p值取为3时可以确保收敛到接近0-1的解。
本文以文献中广泛研究的最小柔顺性问题为例作为算法说明对象，其数学公式可以写成如下所示的有限元形式：
式中，ρ为设计变量向量，它由单元密度ρe组成；m为设计变量的数量；Rm表示m个实数的集合；c为结构柔顺性，并用作目标函数；K为全局刚度矩阵；U为位移向量；F为外部载荷向量。变量g表示体积约束，其表达形式如下：
式中，pe为单元密度；ve为单元体积；γ为体积分数；V为设计域的总体积。在本文中设定γ=0.3，结构的材料属性是各向同性的，杨氏模量E0=1 MPa，泊松比ν=0.3。并使用基于梯度的优化方法即移动渐近线（MMA）方法[40]求解拓扑优化模型式（3）。Sigmund[41]给出了一个99行的MATLAB代码，用于实现拓扑优化并使结构的柔顺性最小化，这为拓扑优化提供了一个基本的流程图。在快速88行拓扑优化MATLAB代码[42]的基础上，文献[43]提供了一个高计算效率的MATLAB代码来解决三维拓扑优化问题。
《2.2.悬臂梁的拓扑优化》
2.2.悬臂梁的拓扑优化
在本节中，给出了图3中所示的悬臂梁结构的拓扑优化结果。设计域尺寸为120 mm×40 mm，设计域左侧固定，并在右下角点处受沿y轴负方向上的集中载荷。 采用四节点平面应力单元将设计域分割为30×10个网格。图4为拓扑优化结果，其中黑色元素表示该单元内填充材料，白色元素表示该单元为孔洞。需要说明的是， 拓扑优化结果中往往存在大量的灰色单元，即单元密度低于1且大于0，此类单元并不具备真实的物理意义，难以制备。因此，将拓扑优化结果中的灰色单元解释为清 晰的0-1可读且可用于3D打印的方法非常重要。
《图3》
图3. 设计模型的示意图。
《图4》
图4. 拓扑优化结果。
Liu和Tovar[43]尝试应用密度阈值法或密度等值线法将拓扑优化结果转换为清晰的0或1表示的结果，并编写了MATLAB函数Top3dSTL_v3来获得可直接3D打印的STL文件。图5分别显示了使用这两种方法获得的STL模型。从图中可以看出，这两种方法获得的结果中包含不受欢迎的薄部件（仅依靠点与点连接）或孤岛。本文的目的是解决如上两个问题。
《3.将拓扑优化结果转换为 STL 格式》
3.将拓扑优化结果转换为 STL 格式
在本节中将提出一种新的后处理流程，主要通过以下两个主要步骤对拓扑优化结果进行处理，以获得没有仅靠点与点连接的薄部件和孤岛现象的清晰0-1优化结果。在此之后，采用MATLAB函数Top3dSTL_v3将优化结果转为STL文件，以适用于3D打印机直接制备。
《3.1.提取结果的骨架以识别细小特征》
3.1.提取结果的骨架以识别细小特征
如图5所示，点状连接或孤岛现象通常出现在尺寸较小的结构特征中。因此，为了避免这两种现象的发生， 首先应识别出细小的结构特征。本文采用骨架提取及其内切圆半径来识别这些细小特征。具体步骤如下。
《图5》
图5.使用Liu[43]的方法获得的STL模型。（a）密度阈值法（cube.stl）；（b）密度等值线法（iso.stl）。
步骤1.1：截取单元密度并重新定义网格。首先确定密度阈值ρcutoff的值，当第i个元素的密度ρi≥ρcutoff时，将其密度定义为ρi=1，而当ρi≤ρcutoff时将其密度定义为ρi=0，至此可以获得清晰的0-1优化结果。图6给出了采用不同阈值时获得的拓扑优化结果，工程师可以根据他们的经验选择其中之一。接下来，网格被细分以便提取结果的骨架。图7展示了网格细分的过程，每个网格被分成nc×nc个较小的细小网格，其中nc>0是正值。应该注意的是，若不进行网格细化，可能会导致结构骨架提取失败。
《图6》
图6.针对不同阈值的拓扑优化结果。（a）ρcutoff = 0.4；（b）ρcutoff = 0.5；（c）ρcutoff = 0.6；（d）ρcutoff = 0.7。
《图7》
图7. 网格细化。
步骤1.2：提取结果的骨架并去除毛刺。定义所有至少两点与基本拓扑边界轮廓相切的内切圆圆心单元为结构骨骼单元，如图8所示。在实际使用中，使用MATLAB函数skel和spur分别提取骨架和去除毛刺。内切圆的半径被定义为骨架的特征尺寸。
《图8》
图8. 结构骨架的示意图。
步骤1.3：识别细小特征。根据所使用的3D打印机的制造精度，确定最小制造尺寸llimit，定义特征尺寸（Rweak）小于制造精度一半（llimit）的骨架单元为细小特征。
图9提供了用于识别拓扑优化结果中的细小特征的流程图。
《图9》
图9.用于识别拓扑优化结果中细小特征的流程图。（a）拓扑优化结果；（b）步骤1.1 ：截取单元密度并重新定义网格；（c）步骤1.2 ：提取结果的骨架并去除毛刺；（d）步骤1.3 ：识别细小特征。
《3.2.保特征的边界细化》
3.2.保特征的边界细化
本节提出一种可以保留原结构特征的边界细化过程，以避免在结果中出现节点与节点连接、孤岛现象以及锯齿状边界。图10提供了该过程的流程图，其中包括以下步骤
《图10》
图10. 保存特征的边界细化流程图。
步骤2.1 ：网格重细化。重新细化原始拓扑优化结果的网格。此过程与步骤1.1类似，但不包括提前截断拓扑优化结果的密度。
步骤2.2 ：处理细小特征。将细小特征的单元元素密度增加0.1 ：
步骤2.3 ：密度过滤。应用密度过滤技术避免锯齿边界：
式中，
为密度过滤后的单元密度；h(i, j) = max(rfilter–||xi – xe||)为权函数，其中xi和xe分别表示第i个和第e个元素中心的坐标；rfilter为过滤半径，在本方法中取rfilter为1。 在完成此过程后，可以获得更平滑的边界。
步骤2.4 ：在保证体积不变的前提下重新截断单元密度。在本步骤中，通过阈值来截断过滤后获得的单元密度，该阈值是使用二分法计算的，根据体积守恒 而确定。在此之后，按步骤1.2中所述计算骨架的特征尺寸，并以此决定是否存在细小特征。如果细小特征仍然存在，回到步骤2.2，如果不再存在细小特征，进行步骤2.5。
步骤2.5 ：获得STL文件。在以上步骤的基础上， 应用密度阈值或密度等值线法获得清晰0-1的拓扑优化结果。进而使用MATLAB函数Top3dSTL_v3将结果转换为STL文件，如图11所示。与图5所示的结果相比，采用本文所提出方法得到的结果更适合于实际制造。
《图11》
图11.使用本文所提出的方法获得的STL模型。（a）密度阈值法（cube.stl）；（b）密度等值线法（iso.stl）。
《4.参数化 CAD 模型（IGES 文件）》
4.参数化 CAD 模型（IGES 文件）
在前一节中，我们提供了一个用于将拓扑优化结果转换为可用0或1表示的清晰结果的流程，其边界相对平滑并且边界点云密度更高。STL文件可以直接由3D打印机制造，以实现快速试生产。但是，如前文所述，STL模型不适用于模型校正。工程师无法修改此模型以根据制造约束条件或应力分布进行修正，因此参数化CAD模型是更好的选择。Chacón等[35]应用三次均匀B样条曲线来获得二维（2D）拓扑优化结果的拟合边界。但是，该方法所需控制点数目较多，进行形状优化时设计变量数目较多，计算量较大。在本节中，为了获得具有尽可能少控制点的参数化CAD模型，本文提出一种适用于形状优化的自适应B样条拟合。其具体过程如下。
步骤3.1：提取边界点并使用三次均匀B样条曲线进行拟合。提取边界点并使用三次均匀B样条曲线来获得拓扑优化结果的拟合曲线。图12展示了具有97个控制点的拟合结果，其中，红点是边界节点，黑点是控制点，蓝线是拟合曲线。三次均匀B样条曲线的函数表达为s(u)。
《图12》
图12.三次均匀B样条曲线拟合结果。（a）边界节点；（b）拟合结果（97个控制点）。
步骤3.2：减少控制点的数量。为了减少控制点的数量，我们提出了一种用于确定在非均匀B样条曲线拟合中使用的控制点位置的方法。其基本思想是减少直线和二次曲线中控制点的数量，均匀B样条曲线的曲率[K(ui)]和曲率的导数[Vi(Ui)]可以通过如下方程计算：
式中，||·||为2-范数；s′(ui)和s′′(ui)分别为拟合曲线s(u)的一阶和二阶导数。图13（a）表示用于边界拟合的三次均匀B样条曲线，图13（b）表示其中一条均匀B样条曲线的曲率和曲率导数。
《图13》
图13.自适应B样条拟合的示意图。（a）使用三次均匀的B样条曲线进行边界拟合；（b）均匀B样条曲线的曲率和曲率导数；（c）自适应B样条拟合（43个控制点）。
步骤3.3：获得标准的IGES文件。根据几何理论，一个大的s′(ui)值表示曲线的转向，而一个大的s′′(ui)值表示转弯与直线之间的相邻点。在本文中，选择阈值kthreshold，将s′′(ui)>kthreshold的点作为节点[即图13（a）中所示的绿点]，将两个节点之间的中点设置为控制点[即图13（a）中所示的黑点]。通过此步骤可以获得非均匀的拟合曲线。得到的自适应B样条拟合如图13（c）所示。该样条曲线只使用了43个控制点，其中红色点是边界节点，黑色点是控制点，蓝色线是拟合曲线。
通过使用MATLAB函数igesout，可以从该B样条曲线中获得标准的IGES文件。图14（a）、（b）分别示出了CAD软件和计算机辅助工程（CAE）软件中的拓扑优化结果的参数模型。图14（c）展示了使用熔融沉积技术（FDM）制造的样品。
《图14》
图14.拓扑优化结果的参数化模型。（a） 犀牛（CAD） 软件；（b）ANSYS（CAE）软件；（c）使用FDM制造的样品。
《5.数值模拟算例》
5.数值模拟算例
在本节中，我们提供了两个数值模拟算例来展示所提出方法的有效性和鲁棒性。此外，在第二个例子中，我们对结构进行形状优化以减少应力集中。
《5.1.热传导问题》
5.1.热传导问题
在第一个例子中，我们考虑了热传导问题。设计域是如图15（a）所示的一个尺寸为50 mm×50 mm的方形板。所有边都被设定为绝热边界，红点的温度被设定为T=0，设计域的中心有一个热源Q。树形的拓扑优化结果如图15（b）所示。
《图15》
图15.（a）热传导问题的设计域；（b）拓扑优化结果。
接下来，我们将所提出的方法应用于该拓扑优化结果以获得CAD模型。按照之前所述过程，我们首先识别细小特征，然后在保留特征的同时对边界进行细化（即步骤1和步骤2）。图16显示了每个子步骤后的结果。在这两个步骤中所使用的参数为ρcutoff=0.2和nc=8。
《图16》
图16. 进行步骤1与步骤2之后获得的热传导问题的结果。（a）步骤1.1和1.2 ；（b）步骤1.3 ；（c）步骤2。
接下来，分别使用均匀B样条拟合和自适应拟合获得两个参数化CAD模型。对于均匀B样条拟合结果，需要439个控制点才可以使结果平滑，而自适应拟合只需要201个控制点。图17展示了使用这两种方法获得的CAD模型以及使用FDM制备的样品。其中，红点是边界节点，绿点是控制点，蓝线是拟合曲线。
《图17》
图17. 拟合结果的CAD模型和使用FDM制造的相应产品。（a）均匀拟合结果；（b）自适应拟合结果。
《5.2.L 型梁》
5.2.L 型梁
在本节中，我们使用了L型梁作为案例来验证所提出方法的有效性，并且对获得的参数化模型进行形状优化以减少拐角处的应力集中。设计域及其拓扑优化结果如图18所示。套用与第一个案例相同的过程。图19显示了步骤1和步骤2之后的结果。在此示例中使用参数ρcutoff=0.5和nc=8。
《图18》
图18.（a）L型梁的设计域（单位：mm）；（b）拓扑优化结果。
《图19》
图19. 步骤1和步骤2之后获得的L型梁的结果。（a）步骤1.1和1.2 ；（b）步骤1.3 ；（c）步骤2。
然后分别使用均匀拟合和自适应拟合获得两个参数化CAD模型。其中，均匀拟合需要113个控制点才能获得平滑的统一的拟合结果，而自适应拟合只需要57个控制点。我们使用本文提出的新方法，获得的参数化CAD模型如图20所示。其中，红色点是边界节点，黑色点是控制点，蓝色线是拟合曲线。
《图20》
图20. 拟合结果的CAD模型。（a）均匀拟合结果；（b）自适应拟合结果
通过使用MATLAB函数igesout，可以获得一个基于自适应拟合结果的标准IGES文件，并建立一个如图21所示的有限元模型。该模型顶侧的所有自由度都受到约束，并在右侧中部（图21中的红色箭头处）施加P=100 N的力。有限元模型的网格尺寸为1，材料的弹性模量为236 MPa，泊松比为0.3。
《图21》
图21. 参数化模型的有限元模型。
为了减少拐角处的应力集中，我们对此模型进行形状优化，优化模型如下所示：
式中，σmax为载荷P作用下的结构最大应力；xi，yi为第i个控制点的坐标，并且被设置为设计变量；Ncp为控制点的数量；VA为优化结构的体积。其他符号的含义与式（3）中的符号含义相同。本文考虑了体积约束和应力约束。
本文使用多岛遗传算法（MIGA）求解该优化模型，优化后的结果如图22（b）所示。基于有限元方法同时分析这两个结果，得到的von Mises应力云图如图23所示。这两个结果的比较表明，形状优化后，结构中承受的最大von Mises应力从110 MPa减小到77.3 MPa。
《图22》
图22.形状优化前后的优化结果。（a）形状优化前的结果；（b）形状优化后的结果。
《图23》
图23. 两个结果的Von Mises应力云图。（a）形状优化前的Von Mises应力云图；（b）形状优化后的Von Mises应力云图。
《6.结论》
6.结论
本文提出了一种将拓扑优化结果转换为适用于增材制造的STL文件的自动实现方法。在所提出的方法中， 首先提取了拓扑结果的骨架并识别出结构细小细节，并且使用滤波方法来确保细小特征得到保存。在此过程之后，获得的拓扑优化结果具有更密集的边界点云和相对平滑的边界。在所获得的结果中引入了自适应B样条拟合，以获得具有更少控制点的光滑参数化CAD模型，这 适用于可提高结果性能的形状优化。这种解释拓扑优化 结果的方法减少了制备优化结果所需的时间，并有助于 标准化拓扑优化的解释过程。
《致谢》
致谢
感谢国家自然科学基金（11332004，1372073， 11402046）、111引智项目（B14013）以及中央高校基本科研业务费（DUT15ZD101）的支持。
《Compliance with ethics guidelines》
Compliance with ethics guidelines
Shutian Liu, Quhao Li, Junhuan Liu, Wenjiong Chen, and Yongcun Zhang declare that they have no conflict of interest or financial conflicts to disclose.
京公网安备 11010502051620号
