基于图像的拓扑优化和后处理技术在2D结构问题中的应用

软件X 14（2021）100701原始软件出版物用于拓扑优化的IbIPP-用MATLAB编写的基于图像的建模和后处理代码Osezua Ibhadodea，Zhao，Zhidong Zhanga，Ali Bonakdarb，Ehsan Toyserkania多尺度增材制造实验室，滑铁卢大学，200 University Avenue West，Waterloo，Ontario，N2L 3G1，Canadab西门子能源加拿大有限公司，9505 Liesse-de-Liesse，蒙特利尔，魁北克，H9 P 1A 5，加拿大ar t i cl e i nf o文章历史记录：2021年2月12日收到2021年4月12日收到修订版，2021年保留字：拓扑优化MATLAB基于密度的方法SIMP斜坡BESO水平集方法后处理3D打印a b st ra ct这项工作提出了一个框架，用于初始化2D拓扑优化的结构问题，同时提供将优化结果后处理为可通过增材制造技术直接生产的STL模型的能力。方便地命名为基于图像的拓扑优化和后处理（IbIPP），它利用准备好的输入图像来表示设计利用拓扑优化的优势为某些物理学创建清晰的结果，可以通过挤出或旋转无缝地获得STL文件。几个数值例子来证明这个工具的有效性版权所有2021作者。由爱思唯尔公司出版这是CC BY许可下的开放获取文章（http://creativecommons.org/licenses/by/4.0/）中找到。代码元数据当前代码版本V1.0此代码版本使用的代码/存储库的永久链接https://github.com/ElsevierSoftwareX/SOFTX-D-21-00033代码海洋计算胶囊https://codeocean.com/capsule/6741412/tree法律代码许可证GPL-3.0使用git的代码版本控制系统使用MATLAB的软件代码语言、工具和服务编译要求，操作环境依赖性MATLAB 2017Rb如果可用，链接到开发人员文档/手册https://github.com/CADmaniac/IbIPP问题支持电子邮件ooibhadode@gmail.com1. 动机和意义拓扑优化是结构优化的一种形式，它解决了结构设计过程中的材料分布问题。它通过基于目标和一个或多个约束[1，2]获得最佳结构布局来实现这一点。拓扑优化生成具有复杂特征的优化设计，这些复杂特征通常难以用成形和减材制造技术来生产。它变得特别有吸引力，因为这个制造瓶颈现在可以通过增材制造或3D打印技术解决。在这项工作中，*通讯作者。电子邮件地址：osezua. uwaterloo.ca（Osezua Ibhadode）.https://doi.org/10.1016/j.softx.2021.100701开发了一个软件平台，以帮助部署用于研究和教育目的的拓扑优化。多年来，有几个令人印象深刻的研究和教育开源代码的拓扑优化。他们使得这一领域的教学、学习和研究相对容易，尤其是对于那些新手来说。从Sigmund的经典99行、基于密度的SIMP代码[3]到更新的代码，如新一代99行代码[4]、使用BESO方法的ANSYS APDL代码[5]、对序列整数规划使用的详细说明[6]、使用几何投影方法[7]等，研究和学术界已经获得了大量的开源拓扑优化代码。然而，仍然有几个进步，可以提高灵活性和易用性，2352-7110/©2021作者。 由Elsevier B.V.出版。这是一篇开放获取的文章，使用CC BY许可证（http://creativecommons.org/licenses/by/4.0/）。可在ScienceDirect上获得目录列表SoftwareX期刊主页：www.elsevier.com/locate/softxOsezua Ibhadode，Zhidong Zhang，Ali Bonakdar等.软件X 14（2021）1007012∑∑nVρe≤∑SIMP⎩ee将多个封装无缝地集成为一个封装，或者使制造具有实际效果在Wei等人[8]报告的介绍中，可以推断出与3D（例如参见[ 4，18 - 20 ]）设计问题相比，更多数量的拓扑优化代码处理2D（例如参见[3，4，9-由于软件和硬件的改进，3D情况更加实用，并且在实现中越来越流行，许多实际设计问题可以简化为2D形式，而不会损失显着的精度。不幸的是，许多可用的开源代码的一个共同属性是在定义初始自由形式域、载荷和支撑条件时的刚性此外，用户无法轻松地将从优化中获得的结果转换为可重新配置的模型。随着增材制造技术的进步，拥有一个开源代码，允许用户几乎无缝地从结构问题到使用拓扑的可打印模型，将带来巨大的好处对于基于密度的方法的具有体积约束的最小柔度问题nminC=FT U=UT KU=Ee（ρe）uT k0ue=1服从：K（ρ）U（ρ）=Fe=11（1）fV00≤ρe≤ 1其中C、F、 U、 K、 V、 f和ρ分别是设计问题的柔度、力、带下标e的变量与有限元有关。两个主要的材料插值函数，在这项工作中定义为γ，是（2）和（3）中提出的具有惩罚的固体各向同性材料（SIMP）[25优化，无需与太多软件接口工具，同时节省计算成本。有些作品γ（ρ）=ρP（二）在[12，21]中可以看到考虑的自由形式设计其他EF-为网络拓扑优化提供软件平台，ρeγRAMP（ρe）=1+q（ 1−ρe）（三）在[20，22，23]中可以看到。开发PolyMesher [12]的一个原因是帮助任意设计域的初始化，然而，由于大多数开源拓扑优化代码使用均匀的正方形网格元素，适应性可能是一个问题。PETSc工具[20]是一个强大的开源拓扑优化代码，但是，它是用C++编写的，不支持任意域初始化。[15，22-24 ]中的其他人总的来说，这些问题加上难以适应其他流行的开源代码，使IbIPP因此，IbIPP是一个Matlab代码，它全面涵盖了自由形式的初始化、几种拓扑优化方法的部署和2D拓扑优化的后处理此外，如果所采取的方法是简单的和成本有效的，则将有利于研究和学术界，因为已经存在一些完善的、相对昂贵的和复杂的拓扑优化工具。由于大多数研究工作都集中在2D问题上，因此将结构设计问题表示为图像可以实现上述目标。这是真的，因为几乎所有能够计算的设备（例如台式机、笔记本电脑、电话等）有基本和免费的图像软件或应用程序。因此，根据一些颜色代码规则，可以开发用作优化工具的输入的图像可以实现用于图像处理、网格化和特征识别的子例程，以将图像转换为对拓扑优化有用的数据集。此外，拓扑优化的结果大多是二进制值元素的矩阵，使用很少的函数程序，可以通过拉伸或旋转将结果处理为可用的STL文件。图像初始化和后处理技术在以下章节和补充材料中有充分的描述。2. 理论背景其中P和q分别是SIMP和RAMP的惩罚因子。因此，每个元件的杨氏模量E e（以及实际上解决设计问题的物理学所对于SIMP，Ee（ρe）=Emin+γe（ρe）（E0−Emin）−Ee（ρe）=γe（ρe）E0−对于RAMP（4）（4）在[3]和[11]中可以看到对基于密度的方法的广泛讨论和实现与（1）类似的（5）中的柔度最小化问题也适用于硬杀伤双进化结构优化（BESO），但没有材料函数。最小C（单位e，u）埃莱服从：K（εe）U（εe）=FnVe=fV0（5）e=1e=0或其中，Δε是每个元件的设计变量，并且被限制为二进制值。总体积计算如下：Vi=Vi−1（1±cer）（6）其中cer是确定在当前迭代中从先前迭代的体积中添加或移除多少材料的演变比率。对进化方法的详细回顾可以在[30]中看到水平集方法在优化过程中通过求解一个方程来隐式定义结构边界：最常见的是Hamilton-Jacob方程。水平集函数φ可以采用（7）中的三个等值面值中的任何一个：φ1≥φ（x）>0φx∈φ\φ2.1. 拓扑优化方法在几种拓扑优化方法中，φ（x）=00> φ（x）≥ −1 x∈D\n（七）多年来，三种方法被广泛应用：基于密度的方法，进化方法和水平集方法，没有流行的顺序。其中x是一个域，而x是x中的一个位置。D是一个固定的设计域，其中将为优化过程确定结构布局。使用边界Γu处的固定位置Osezua Ibhadode，Zhidong Zhang，Ali Bonakdar等.软件X 14（2021）1007013（）下一页D=×× ×⎤[客户端]⎥Γ⎡⎢⎤⎥⎡⎢⎤⎥和在边界Γt处的牵引力ti，连续场中的柔度最小化问题定义如下[14]：（11）可以改写为infF=0φti uidrg（i，j，k）=d1，1d1，2···d1，nelxd2，1d2，2···d2，nelx...（十二）好吧Ω.···⎦S. t.Vmax≤ 0其中ddnely，1dnely，2· · ·dnely，nelx获得作为divCijkluk，l=0 onn （8）ui=uiinΓui、j1I j∑cxy（x，y）∈<$ij（13）ti=tiinrt，npi、、、，ji， j其中Vmax是最大体积约束，Cijkl是弹性张量，ui和ti是位移和载荷变量，其中ui和ti分别是恒定位移和载荷值。水平集方法在[14]和[8，31]中进行了充分讨论，而MATLAB代码在[14]和[8]中实现。2.2. 数字图像为了初始化在前一节中讨论的拓扑优化方法的2D设计问题，可以开发并处理准确地描绘问题的特征（诸如负载、边界条件、保留和非设计区域）的图像作为输入。使用基本的图像矢量化和处理技术，可以分析分配给设计问题的特征的颜色，并且此后提取初始化所需的变量。数字图像包括几个像素，类似于由有限元素。图1（a），由函数定义的图像f（x， y， z）有m n个像素[32]，每个像素的强度为cx，y。对于真彩色表示，图像函数可以在（9）中表示。f（x，y，z）=f（0， 0，c0，0）f（0， 1，c0，1）···f（0，n− 1，c0，n−1）f（1，0，c1，0）f（1，1，c1，1）···f（1，n−1，c1，n−1）...其中npi，j是区域内元素（i，j）吉尔吉在（11）中，nely和nelx分别是均匀分割图像中y和x方向上的元素数，而di，j是由i， j表示的区域的平均颜色强度。读者可以参考[32，33]和[34]了解更多关于这个主题的细节和相关信息。3. 软件描述IbIPP是几个MATLAB函数的集合，可实现设计域初始化、轻松部署和执行几个开放源代码MATLAB代码，用于拓扑优化和生成3D打印模型。使用颜色代码，设计问题表示为图像格式，并应用几个子程序将数据转换为有用的输入拓扑优化。除了用于网格化的单位长度的均匀正方形有限元外，参数是可定制的。IbIPP中所有可能的输入参数的详细信息在补充材料或GitHub存储库中的表S1中列出。3.1. 软件构架在Matlab命令提示符下输入以下代码初始化程序：.好吧..···..⎥⎦ibipp（domain，nelx，volfrac，. ）f（m−1， 0，cm−1，0）f（m−1， 1，cm−1，1）f（m−1，n− 1，cm−1，n−1）（九）其中x和y表示像素位置，而c表示强度，对于RGB图像，c是大小为3的向量，并且总数据大小为mn3.如果这个图像被假定为一个设计问题的域，如图所示。1（b）中的每个元素该域将表示像素区域λ。因此，（9）可以被解释为（10）由有限元组成的新矩阵可以从图像导出为：g（i，j，：）=fαi， j，βi， j，：（10）其中g（i，j，：）包含f的子矩阵，并且αi，j和βi，j分别是f中像素区域 i ， j的行索引和列索引，如图1所示。 1（b）. 域函数g（ i， j， k）可以写为：g（i，j，k）g（1， 1，d1，1）g（1， 2，d1，2）···g（1，nelx，d1，nelx）g（2，1，d2，1）g（2，2，d2，2）···g（2，nelx，d2，nelx）=...其中domain是字符串格式的图像文件的名称（推荐使用Portable Network Graphics-png），volfrac是对问题的体积分数约束，而''...“”表示用户可以通过可选输入和名称-值输入声明其他几个参数。在ibipp.m（主文件）中，主要部分是默认参数声明、必需、可选和名称-值输入的输入解析器、优化参数声明、负载和边界条件参数声明、后处理参数声明、域初始化、负载解析、边界条件、保留和非设计域、拓扑优化和后处理。代码的关键方面是为域初始化、载荷解析、边界条件、保留和非设计域以及后处理编写的子例程。3.1.1. 域解释（imageprocessor.m）在ibipp.m运行之前，必须根据某些设定的要求准备代表设计问题的图像。任何允许用户绘制2D草图的软件或工具并指定用于着色的RGB值。用户鼓励使用工具，使素描，和颜色应用程序-.好吧..···..⎥⎦不允许渲染功能，如“环境”设置。有8个独特的特征可识别g（nely，1，dnely，1）g（nely，2，dnely，2）···g（nely，nelx，dnely，nelx）（十一）通过程序：设计和非设计元素，点荷载和压力节点，仅在x中的固定节点，仅在y中的固定节点，Osezua Ibhadode，Zhidong Zhang，Ali Bonakdar等.软件X 14（2021）1007014=-Fig. 1. 图像的坐标表示（a）通过普通图像处理包，（b）通过有限元的类似表示，其中m和n分别是图像的垂直和水平轴上的像素数，nely和nelx是图像的垂直和水平轴上的元素数，代表性的有限元网格，α和β是（10）中的子矩阵f的列和行索引y，保存的元素。设计域根据颜色强度分为这些特征这些特征被分配了一些颜色代码，如表1所示。要启用多个负载，表1中的点负载和压力的R和G值必须是5的倍数，并且必须从点负载的R 200和G开始压力负荷为100。这有助于保持将荷载工况与其相应幅值相关联的顺序和方法。随后如何提取载荷节点以及施加的量值在补充材料的第1图2，半MBB和另一个稍微复杂的设计问题说明了如何为设计问题生成图像。据观察，对于图1中的载荷和边界条件，在图2（a，f）中，直接受影响的设计域中的区域由根据表1的适当颜色代码描绘。这些区域的大小取决于用于优化的分辨率。分辨率越精细，图中这些彩色负载和边界区域应覆盖的面积就越小。2（b，g）. 值得注意的是，用户对nelx的选择应该谨慎选择。图2（c，h到e，j），使用粗略的分辨率可能导致对数字图像的错误解释，而更精细的分辨率确保所有特征都能很好地复制。为了在正确的图像解释和由于应用更精细的分辨率而产生的计算成本之间取得平衡，nelx应该大致为较小边界框尺寸与最小特征的最小尺寸之比的两倍的最小值。图像判读由子程序实现[dom，nely]=imageprocessor（domain，nelx）在该子例程的第2行到第11行中，图像被读取、离散化并提取RGB值。根据表1，元素值从第12行应用到第25行，如图1所示。3 .第三章。imageprocessor.m的输出是dom和nely，它们是矩阵f（i， j）和沿y轴的元素数。准备图像文件时应注意的几点：鼓励用户使用表1中提供的RGB范围的最小值（白色除外）（255，255，255），彩色区域不应以另一种颜色的线为边界（例如，绿色区域不应以黑线为边界）。最后鼓励用户浏览提供的示例，以便更好地了解如何为设计问题准备图像。3.1.2. 荷 载 、 支 撑 、 非 设 计 和 保 留 范 围 的 解 析 度（loadandsupport.m）在将数字图像解释为具有代表设计特征的元素编号的矩阵这是在 ibipp.m 的 第 59 行 中 完 成 的 ， 其 中 启 动 了 子 例 程loadandsupport.m。子例程由以下函数命令[F，fixeddofs，NonD，MusD，volfrac，edofMat1]=loadandsupport（nelx，...nely、Fmag、Fang、Pmag、dom、pre_support、pre_load、volfrac）除了nelx、nely、dom和volfrac之外，该函数的其他输入是Fmag，在补充材料中的图S1中被描绘为F的力大小的矢量，Fang是被描绘为φ的力角度的对应矢量，Pmag是压力大小的矢量，pre_support 是保留承载固定节点的元素的选项，并且pre_load是保留定义载荷节点的元素的选项。在这个子程序中，第5行到第9行声明了域的节点数和自由度的连接矩阵补充材料第1节中详细说明了荷载、支座、保留和非设计域的实施。3.1.3. 后处理（datatostl.m）一些实际的设计问题可以公式化或近似为2D问题。通过此后处理功能，用户可以获得可重新打印的模型，特别是3D打印技术。使用IbIPP，用户可以初始化任意设计域，使用三种主要方法中的任何一种运行拓扑优化，并具有足够的参数控制，并几乎无缝地生成可打印的STL模型IbIPP中有两种建模类型：挤出和旋转，而且，对称可以Osezua Ibhadode，Zhidong Zhang，Ali Bonakdar等.软件X 14（2021）1007015图二. 单元尺寸对ibipp.m.图像判读的影响。这示出了（a，f）要优化的设计问题，（b，g）表示设计问题的输入图像，以及通过（c，h）10× 5（6）（d，i）60× 31（36）和（e，j）150× 77（90）网格元素的域离散化被应用。为了启动子例程datatostl.m，名称字符串中的.stl必须作为函数输入。当以下函数调用运行时，子例程将运行datatostl（nelx，nely，xPhys，tx，form，ht，lr，theta，symm）tx是STL文件名，必须以.stl结尾，并且是字符串格式，form是建模类型，ht是基于挤出的建模的因素，并且当它乘以最小2D尺寸时，它定义了挤出长度，如图所示。在图4（a）中，lr是用于基于革命的建模的因子，其定义优化拓扑与轴的距离Osezua Ibhadode，Zhidong Zhang，Ali Bonakdar等.软件X 14（2021）1007016图3.第三章。 通过颜色代码和相应元素值进 行 图 像 判 读 ， 如表1所示。图四、Ib IP P 中 提 供 的 建模选项。（a）挤压（b）旋转（c）对称。中心这个距离是当lr乘以nelx时得到的。 旋转轴总是平行于优化的2D模型的左侧，如图所示。4（b），θ是旋转角。最后，对于对称设计情况，symm允许用户生成优化拓扑的镜像模型。当symm采用字符串输入left、bottom、right或top时，指定对称建模和镜像轴，如图 4（c）.在datatostl.m中的第2行中，在伪密度矩阵的每个元素上放置0.5的阈值，使得下面的任何值第4至21行仅在symm具有值时获得伪密度矩阵第23至35行根据建模类型生成优化拓扑的3D矩阵。虽然第26到28行中的简单指令处理挤出，但Treeby和Cox在他们的k-waveMatlab中开发的子例程revolve 2D. m[35]在第34行中用于将2D矩阵旋转为3D。在revolve2D.m中，通过添加第60至75行修改了子例程，以适应任何旋转角度。最后，Sven [36]开发的子例程stlwrite.m将IbIPP的工作流程如图所示。图5示出了在前几节中详细描述的主要子例程之间的交互。调整流行的开源拓扑优化代码Osezua Ibhadode，Zhidong Zhang，Ali Bonakdar等.软件X 14（2021）1007017−5−−表1设计特征及其颜色代码表示。设计特征颜色代码元素值f（i，j）指定域名RGB200 1<非指定结构域RGB>200 0点荷载200≤Ra≤ 255，G= 0，B= 0 20+R2005压力200≤R≤230，100≤Ga≤150，B=0 30+G−100保留区R=0， 200≤G≤ 255，B=0 4固定区域R=0，G= 0， 200≤B≤255 51沿x固定的区域100≤R≤ 150，G= 0， 200≤B≤255 52沿y固定的区域=0，200≤G≤255，200≤B≤255 53注意：ibipp将任何颜色代码超出所列规则的元素解释为0或白色。a必须是5的倍数，以允许多个荷载工况;此外，R和G必须分别从200和100开始图五. IbIPP的软件体系结构显示了主要的子程序。(such如top88.m、esoL.m和levelset88.m），以适应IbIPP，详见补充材料第2部分第1节。4. 说明性实例4.1. 实施例1提出了两个简单的问题;第一个是图1中的半MBB/MBB问题。6（a），这是多年来文献中拓扑优化的基准设计之一[4，37[38]在图。6（d）.对于半MBB，图1B中的输入图像。图6（b）具有代表力位置的具有颜色代码200、0、0的红色区域、代表仅在y轴上固定的支撑元件的具有颜色代码0、200、200的青色区域以及代表在x轴上固定的支撑元件的具有颜色代码100、0、200的紫色区域只有轴。为了锤头码头图如图6（d）所示，有四个力，其区域从左到右具有颜色代码（200， 0， 0）、（205，0， 0）、（210， 0， 0）、（215， 0， 0），在IbIPP中，当力位置位于载荷区域的中心时，从支撑元件获得固定的自由如图6（b，e）所示，为了更好地表示设计问题，载荷和支撑区域尽可能小，以减少这些区域中的元件数量，从而精确定位力的位置并减少固定自由度的数量然而，请记住，载荷和支撑区域越小，优化所需的分辨率就越高，因此用户应该在模型精度和计算成本之间取得平衡。图中的优化拓扑结构。6（c，f）通过输入以下命令获得：ibipp（'half_mbb.png'，500，0.4）→半MBBibipp对于半MBB问题，只有三个输入参数，IbIPP调用几个默认值。因此，用户可以使用三个参数运行优化， 对于多个荷载工况，如锤头墩所示，输入参数中的力大小和角度应采用矢量形式，同时确保右矢量元素对应于右荷载工况，如图所示。 6（e）.4.2. 实施例2为了实现IbIPP的后处理子程序，本文提出了一种带有作者实验室缩写（MSAM）的跨接器设计，如图1所示。第7（a）项进行调查。一半的设计被认为是优化，以节省计算成本和时间，因为可以在中间画一条对称线。可以观察到，在图1中，附加区域被添加为固定区域。7（b）只是以确保更多的材料被分配到头部的。 在图 中 对模型进行优化后 ，第 7（c）条，Osezua Ibhadode，Zhidong Zhang，Ali Bonakdar等.软件X 14（2021）1007018见图6。 （a）半MBB和（d）锤头墩的拓扑优化，使用（b）和（e）中所示IbIPP的数字图像输入，以及优化结构分别载于（c）及（f）见图7。使用IbIPP生成优化的XML模型的工作流。(a)设计问题与负载和边界条件，（b）输入图像（c）优化拓扑结构（d）通过对称和挤出生成模型（e）打印的图形。后处理子程序DataSt1.m接收密度矩阵和挤出所需的其它数据。通过连接密度矩阵和它的镜像矩阵，并应用扩展长度，得到图1中的优化模型。7（d）.出于制造目的，例如3D打印零件，可能需要根据优化的分辨率进行平滑处理MeshLab [40]等开源代码可用于平滑和其 他网格相关操作 。在本 例中，对 生成的模型执行Laplacian [41]和Taubin [42]平滑操作。如图所示，通过熔融沉积成型（FDM）使用PLA打印出了该模型。7（e）.用户应该记得，只有当名称-值输入“modelName”作为输入参数包含在函数声明中时，才会启动后处理，如下所ibipp5. 影响IbIPP提供了一个平台，可以在MATLAB中执行从任意域初始化到数据后处理的2D拓扑优化。还有其他几个开源和商业的拓扑优化工具具有相当令人印象深刻的功能，但是，存在一些缺点。正如本文介绍中所强调的，大多数可用MATLAB代码的一个最大缺点是无法初始化自由形式的域Osezua Ibhadode，Zhidong Zhang，Ali Bonakdar等.软件X 14（2021）1007019并进行数据后处理。因此，该平台不被认为是更强大的开源和商业代码的替代品，而是一种简单，具有成本效益和整体性的方法，用于在MATLAB中实现2D拓扑优化问题，用于研究和教育目的。考虑到这一点，IbIPP提供了各种优势：(a) 它打破了初始化任意设计域的障碍，并在表示载荷和边界条件方面提供了灵活性。(b) 它使研究人员能够测试各种设计问题，而无需频繁修改或编写新的代码行。(c) 它部署了几种拓扑优化方法，确保研究人员，教师和学生可以使用各种这些方法具有相对的灵活性。它还有助于增加新的方法（补充材料对此作了(d) 在许多情况下，对优化结果进行实验测试。IbIPP提供了一个简单易用的平台来生成拉伸、镜像或旋转模型.对于低分辨率优化，生成的模型可能由具有许多锯齿状边缘的特征组成;可以采用基于网格的开源工具来平滑特征。目前，IbIPP用于多尺度增材制造（MSAM）实验室，以协助开发拓扑优化的新模型，因为它与多物理场应用，设计相关负载和多孔特征考虑因素有关。IbIPP有望使拓扑优化领域的研究人员、教师和学生6. 结论开发了一个基于图像的二维拓扑优化初始化和后处理程序。它接受根据某些颜色代码规则开发的输入图像，这些颜色代码规则正确表示拓扑优化设计问题中的特征。根据用户的选择，它还可以获得基于优化拓扑的可打印模型。这些能力已经得到了证明，进一步的工作，以改善该工具的部分，并将代码翻译成其他语言将完成。竞合利益作者声明，他们没有已知的竞争性财务利益或个人关系，可能会影响本文报告的工作致谢这项工作得到了加拿大自然科学和工程研究委员会（NSERC）、安大略省南部联邦经济发展署（Fed-Dev Ontario）、西门子能源加拿大有限公司和石油技术开发署（PTDF）的资助。附录A. 补充数据与本文相关的补充材料可以在https://doi.org/10.1016/j.softx.2021.100701上找到。引用[1]Sigmund O.拓扑优化：理论，方法和应用。七二四第2版2003年。[2]西 格 蒙 德 ·O ， 穆 特 ·K 。 拓 扑 优 化 方 法 。 比 较 审 查 。2013;1031http://dx.doi.org/10.1007/s00158-013-0978-6[3] 西格蒙德·奥用mat-lab编写的99行拓扑优化代码. Struct Multiplier Optim2001;21（2）：120-7. http://dx.doi.org/10.1007/s001580050176.[4]Ferrari F，Sigmund O.新一代99行Matlab代码，用于柔性拓扑优化及其扩展到3D。StructMultiplexpOptim2020;62：2211-28.http://dx.doi.org/10.1007/s00158-020-02629-w网站。[5]林华，徐A，Misra A，赵R. 本文介绍了一种基于动态演化速率的BESO方法的 多约 束结 构拓 扑优 化 设计 程序 。Struct Multiplexp Optim 2020;2229-54.http://dx.doi.org/10.1007/s00158-020-02588-2网站。[6]梁毅，程国.进一步阐述了用序列整数规划和正则松弛算法进行拓扑优化的方法和128行MATLAB代码。Struct Multifunctional Optim 2020;61（1）：411-31.网址：//dx.doi.org/10.1007/s00158-019-02396-3网站。[7]Smith H，Norato JA.用几何投影法进行拓扑优化的MATLAB代码。StructMultiplexp Optim 2020;1579-94. http://dx.doi.org/10.1007/s00158-020-02552-0网站。[8]魏平，李智，李新，王美云.基于径向基函数的参数化水平集方法拓扑优化的88行MATLAB代码。Struct Multiplexp Optim 2018;831-49.[9]苏雷什·K一个199行的Matlab代码，用于拓扑优化中的Pareto最优跟踪。StructMultiplexpOptim2010;42 （ 5 ） ： 665-79.http://dx.doi 的 网 站 。org/10.1007/s00158-010-0534-6。[10]查利斯VJ。一个用Matlab写的离散水平集拓扑优化代码. Struct MultiplexpOptim 2010;41（3）：453-64. http://dx.doi.org/10的网站。1007/s00158-009-0430-0。[11][10]张晓刚，张晓刚. MATLAB中的高效拓扑优化，使用88行代码。StructMultiplexp Optim 2011;43（1）：1-16. http://dx.doi.org/10.1007/s00158-010-0594-7网站。[12]Talischi C，Paulino GH，Pereira A，Menezes IFM. PolyMesher：一个用Matlab 编 写 的 通 用 多 边 形 网 格 生 成 器 . Struct Multiplexp Optim 2012;45（3）：309-28. http://dx.doi.org/10.1007/s00158-011-0706-z.[13]Xia L，Breitkopf P.在Matlab中使用拓扑优化和基于能量的均匀化方法设计材料 。 结 构 多 线 程 优 化2015;52 （ 6 ） ： 1229-41 。http://dx.doi.org/10.1007/s00158-015-1294-0网站。[14]张文忠，张文忠.基于水平集的反应扩散方程拓扑优化方法的Matlab代码。Struct Multiplexp Optim 2015;1159-72. http://dx.doi.org/10.1007/s00158-014-1190-z网站。[15]Victoria M，Querin OM，Díaz C，Martí P. liteITD一个用于连续体结构拓扑设计的MATLAB图形用户界面（GUI）程序。Adv Eng Softw 2016;100：126-47. http://dx.doi.org/10.1016/j.advengsoft.2016的网站。07.012.[16]Sanders ED，Aguiló MA，Paulino GH.具有任意体积和质量约束的多材料连续 体 拓 扑 优 化 。 计 算 方 法 应 用 机 械 工 程 2018;340 ： 798-823 。http://dx.doi.org/10.1016/j.cma.2018的网站。01.032[17]陈强，张X，朱B.几何非线性结构的213行拓扑优化程序。Struct MultiplexpOptim 2019;59（5）：1863-79. http://dx.doi.org/10.1007/s00158-018-2138-5网站。[18]Liu K ， Tovar A. 一 个 用 Matlab 编 写 的 高 效 三 维 拓 扑 优 化 代 码 。 StructMultiplexp Optim 2014;50 （ 6 ） ： 1175-96. http://dx.doi.org/10 的 网 站 。1007/s00158-014-1107-x。[19]Amir O，Aage N，Lazarov BS.基于多重网格CG的高效拓扑优化。StructMultiplexpOptim2014;49 （ 5 ） ： 815-29.http://dx.doi 的 网 站 。org/10.1007/s00158-013-1015-5。[20]Aage N，Andreassen E，Lazarov BS.使用PETSc进行拓扑优化：一个易于使用 、 完 全 并 行 、 开 源 的 拓 扑 优 化 框 架 。 Struct Multiplexp Optim 2015;51（3）：565-72. http://dx.doi.org/10.1007/s00158-014-1157-0.[21] 蔡S，张伟.自由设计域的应力约束拓扑优化。计算方法应用机械工程2015;289（2018）：267 http://dx.doi.org/10.1016/j.cma.2015.02.012网站。[22]周永华，林熙。自动结构最佳化系统之改良影像判读与模型化技术StructMultiplexp Optim 2010;40（1-6）：215-26. http://dx.doi.org/10.1007/s00158-008-0352-2网站。[23] Gamache JF，Vadean A，Noirot-Nérin JF，Beaini D，Achiche S.基于图像 识 别 的 桁 架 密 度 拓 扑 优 化 方 法 。 Struct Multiplexp Optim 2018;58（6）：2697-709. http://dx.doi.org/10.1007/s00158-018-2028-x.[24]作者：Lin CY，Chou YH.集成拓扑和形状优化的自动结构优化系统。中文研究 所 工 程 师 / 中 华 民 国 。 Cheng Hsuch K'an 2008;31 （ 5 ） ： 745-56.http://dx.doi的网站。org/10.1080/02533839.2008.9671429。Osezua Ibhadode，Zhidong Zhang，Ali Bonakdar等.软件X 14（2021）10070110[25]放大图片作者： Rozvany GIN，Zhou M. 无均匀化广义形状优化. StructOptim1992;4（3-4）：250-2.http://dx.doi.org/10的网站。1007/BF01742754。[26]Zhou M，Rozvany GIN. COC算法，第二部分：拓扑，几何和广义形状优化 。 计 算 方 法 应 用 机 械 工 程 1991;89 （ 1-3 ） ： 309-36. http ：//dx.doi.org/10.1016/0045-7825（91）90046-9.[27]Sigmund O. 拓 扑 优 化 中 的 材 料 插 值 方 法 。 Arch Appl. Mech 1999;69 （ 9-10）：635-54. http://dx.doi.org/10的网站。1007/s004190050248。[28][10]张志华，张志华.用边界识别-载荷演化（BILE）模型对设计相关压力载荷下 的 结 构 进 行 拓 扑 优 化 。 结 构 多 线 程 优 化2020;62 （ 4 ） ： 1865-83 。8.http://dx.doi.org/10.1007/s00158-020-02582-[29]Stolpe M ，Svanberg K.拓 扑最 佳 化惩 罚 法之 轨迹 研 究。 Struct MultiplierOptim 2001;21（2）：128网址：//dx.doi.org/10.1007/s001580050177网站。[30]夏丽，夏勤，黄新，谢玉明。先进结构和材料的双向渐进结构优化：全面综述。Arch Comput Methods Eng 2016;（November）. http://dx.doi.org/10的网站。1007/s11831-016-9203-2。[31]杨文，李文.基于水平集方法引入虚拟界面能的拓扑优化方法。计算方法应用机械工程2010;199（2876）：45-8-2891.http://dx.doi.org/10.1016/j.c