浅谈拓扑优化在建筑结构领域中的应用（三）
这次先暂时不谈建筑与结构，来谈谈拓扑优化及其基本概念。
一直以来，结构优化设计（Optimal Structure Designing）都是力学领域的一个重要分支，它以力学为基础，集数学、计算机科学甚至图形学等众多不同学科于一体，是指为了满足某种目标（如材料消耗最少化、刚度最大化等），在某种约束下（如体积/重量约束、几何约束等），而对结构进行的优化设计。
目前，结构优化一般可分为三个不同的层次，包括：尺寸优化（Size Optimization）、形状优化（Shape Optimization）以及拓扑优化（Topology Optimization），它们分别对应着实际设计中的三个不同阶段。
其中，尺寸优化一般以壳体的厚度或者构件的截面尺寸等作为设计变量（设计变量即优化设计中需要改变的参数）；形状优化是在保证结构拓扑关系不变的前提下，以构件的截面形状、节点的空间位置以及连续体的形状等作为设计变量，即改变模型的几何特征（修改几何边界）；拓扑优化的概念源自于拓扑学，是指在某种优化准则和约束条件下，在模型的设计域（即模型需要优化的区域）确定实体域的数量及布局，对建筑结构整体而言就是确定各构件（梁、柱、墙等）的布局及其节点的连接形式与位置。简单的说，就是我就给你这么多材料，你需要确定这些材料在设计域内的空间分布（也就是说，以材料在设计域内的空间分布作为设计变量），使所得结构在这个条件下满足最大刚度或满应力等优化准则。
浅谈（一）中已经简单介绍了离散体结构（Discrete structure）拓扑优化方法和连续体结构（Continuum structure）拓扑优化方法，两种方法各有优劣。
以桁架为代表的离散体结构，其拓扑优化主要是研究在满足一定的边界条件下桁架间各杆件的最佳尺寸及其最佳连接方式，具有较少的设计变量，计算效率较高。然而在许多设计领域中，设计的往往不是类桁架结构，可能就是在一块板里开几个洞，那么对于这种结构，离散体是很难描述的，此时就需要采用连续体结构拓扑优化方法。然而，对于杆件数目成千上百的输电铁塔、体育馆网架屋盖等结构，若采用连续体模型去描述，那得需要多精细的三维有限元模型啊。在这种问题下，考虑到当前的计算能力，明显离散体拓扑优化是更好的选择。
对于离散体结构拓扑优化问题，Dorn等在1964年首先提出“基结构”（Ground Structure）的概念[1]。所谓的基结构，是先确定荷载作用点、约束作用点以及其他可能的结点等一系列点的分布，然后将各个点用杆件两两连接而成的初始结构。而基结构法（Ground Structure Method）就是对基结构的杆件截面尺寸采用优化算法进行优化的过程，即根据满应力原理，根据杆件内力逐步调整杆件截面，使结构所有杆件的应力都趋于指定的应力值。然后去掉结构中杆件截面面积小于消除基准值的杆件，从而消除低效杆件，最终得到满足目标约束下的最优结构形式。
除了基结构法，像粒子群算法、蚁群算法等一系列搜索类的算法也大量应用于离散体结构拓扑优化问题（看来通过人工智能来进行结构优化也是指日可待了）。其优化方法与前面所说的基结构法一致，只是应用了搜索类算法求解离散问题特有的优势。
看下面的一个基结构法算例[2]：图2中各杆的弹性模量E= 206. 9GPa,材料密度ρ= 7.8×103 kg /m³,许用应力[σ]= 172. 4MPa,荷载P1= P2= 0. 45MN,各杆的初始截面积均取为64. 5 cm²,截面积下限取为10-3cm²。在只考虑应力约束的条件下，通过PSO算法，得到了图3所示的优化结果。杆2、5、6被去除，而剩下的杆1、3、4、7、8、9的面积也分别优化至51.6、51.6、25.8、36.5、36.5、36.5cm²。
直至今天，国内外关于离散体结构拓扑优化的相关研究工作已经相当成熟，取得了不少的成果。然而，相比之下，连续体结构拓扑优化具有更大的应用范围，但其难度更大，挑战性也更强，还存在着许多问题亟需解决，其相关理论研究和实际应用工作是当下结构优化设计领域的热点。
相比于离散体结构拓扑优化，在过去的很长一段时间里，由于计算量庞大，算法复杂等原因，连续体结构拓扑优化的发展迟缓，甚至一度陷入瓶颈：在上世纪五十年代之前，对结构优化问题的解决方法只限于最经典的数学解析方法，但是由于对连续体结构拓扑优化问题的解答需要进行复杂的计算分析，而人工计算效率低下，使之不能在实际工程中大范围的推广应用。后来由于有限元法的推广、计算机计算性能的快速提升以及数学规划方法的出现，推动了结构优化的快速发展。近年来，不同的优化技术和方法层出不穷，主要包括：均匀化方法、变密度法、渐进结构优化法、水平集法等。
连续体结构拓扑优化一般是指在给定的载荷和约束条件下，确定连续体内有无孔洞以及孔洞的位置、数量和形状等，以得到使目标性能最佳的材料分布。由于理论上连续体的设计域中坐标点的数目总是无穷的，在这些点上都存在1或0的变量值来表示该点处材料的存在与否，但是，我们没办法通过数值方法将无穷多的点进行处理。为此，在有限元方法的启发下，我们将连续体进行离散，划分为若干个单元，基于特定算法判断各单元保留或去除。
从本质上来说，由于设计变量要么为0（单元去除），要么为1（单元保留），拓扑优化其实就是一种0-1离散变量的组合优化问题，若设计域中存在n个设计变量，就需要从2^n个方案中寻找最优解。当n足够大时，就会出现“组合爆炸”的问题，所需的分析次数远远超出计算机的计算能力之外，而导致无法计算。为了克服该困难，人们往往将设计变量松弛为0-1的连续变量，以利于使用成熟的基于连续变量的导数优化算法，从而避免离散变量问题的求解困难。但中间密度单元在现实中是不存在且难以制造的，因而在优化过程中，需要对中间密度进行惩罚，使之向0或1两极转变。而变密度法则是在此背景下诞生的。
明确拓扑优化中的拓扑描述方式和材料插值模型非常重要，这是一切后续优化方法的基础[3]。在变密度法中，对材料分布的描述可分为单元变量与节点变量两大类。基于单元变量的变密度法是在在对结构进行有限元网格划分后，为网格中的所有结构单元分别赋予0-1连续取值的变量相对密度ρ，其中以ρ=0表示该结构单元处不存在材料或者说由一种假想的“空白材料”所填充，以ρ=1表示该结构单元处充满了实体材料，而0<ρ<1则表示结构单元处所填充的是一种具有微结构孔洞的材料，其材料物理属性（一般指弹性模量）可根据特定的函数关系由实体材料与“空白材料”的属性值插值所得。根据材料物理属性插值函数的选取不同，变密度法又可分为固体各向同性材料惩罚法(Solid Isotropic Material with Penalization，SIMP )[4-7]和材料属性的有理近似法(Rational Approximation of Material Properties, RAMP)[8]等。基于节点变量的变密度法则是在结构网格中选取一定数目的空间点（可包括但不限于网格节点）赋予变量ρ，然后根据给定的映射关系得到各结构单元内的材料分布，即结构单元内的相对密度不一定是常量。
下图为分别基于单元变量与节点变量的变密度法中对结构材料分布的描述示意图，其中蓝线为网格单元的边线，黑色代表实体材料，白色代表空白材料，灰色代表中间密度材料。
渐进结构优化法（ESO）是Xie与Steven在达尔文进化学说的启发下所提出来的[9]，其基本思想是：按照一定的优化准则，逐步去除低效单元，保留高效单元，最终展现出结构中的最佳传力路径，从而得到更为高效的拓扑结果。某种意义上来说，ESO方法是变密度法的一种变形。在ESO方法中，也需要先将连续体结构用有限元网格进行划分，离散为若干个单元。但不同的是，在ESO方法中单元的相对密度只有0（空白单元）和1（实体单元）两种，分别对应被删除和被保留的单元。
下图为ESO/BESO法中对结构材料分布的描述示意图，其中蓝线为网格单元的边线，黑色代表实体材料，白色代表空白材料。
上面描述的变密度法和ESO方法可以说都属于材料本构模型的插值问题。
从另一个角度出发，在连续体结构拓扑优化中，材料分布的变化实际上一般也等同于材料交界面的变化，其中，二维平面结构与三维板壳结构的几何边界为闭合曲线，三维实体结构的几何边界为闭合曲面。因此，连续体结构的拓扑优化过程可转化为几何意义上闭合曲线或闭合曲面的演化过程。水平集方法正是跟踪界面移动的数值计算技术[10]，可以隐式地表达动态变化中的连续体结构几何边界。后来，Sethian 和 Wiegmann 将水平集方法引进了结构拓扑优化领域[11]。
采用水平集函数进行拓扑描述时，在任意时刻只需求出水平集函数场的零等值线或零等值面即可得到此时的结构几何边界，并不需要对几何边界本身进行显式的描述，因此可以自然地处理复杂的几何界面及其拓扑变化的情况（如孔洞的分离、融合、消失等）。
水平集法是一种以高维函数去描述低维目标的方法。举个简单的例子，如下方左图所示的圆，以原点为圆心，半径为r=4，则可通过定义方程x^2+y^2=4^2对其进行显式描述。水平集方法则是先选取覆盖该圆的足够大区域，然后在区域内定义水平集函数标量场，最后在区域内找到水平集函数值为0的空间点所构成的等高线，此等高线即为所需描述的圆。其中，水平集函数标量场的数值实施具体过程为，先对区域进行网格划分，然后在各网格节点上分别赋予特定的标量函数值，其取值可以但不限于设为该节点与原点距离的平方与r^2之差，所得函数场如下方中图所示。因此，如下方右图所示，黑线为网格线，红线为水平集函数的零等高线即所描述的圆。
那么，如何利用水平集函数描述结构拓扑情况呢？很简单，一般情况下，零水平集表示的就是结构的几何边界，一般将水平集函数场量大于0的地方视为结构内部，充满实体材料，反之水平集函数场量小于0的地方视为结构外部，充满空白材料。
类比变密度法基于节点变量的材料分布描述示意图，水平集法无需引入中间密度材料，如下图所示，其中蓝线为网格单元的边线，黑色代表实体材料，白色代表空白材料。
在二维平面应力结构的拓扑优化应用中，对于同一结构模型，上述三种拓扑描述方法所得结构如下图所示：
最后我们还得确定优化所需要的数值算法。变密度法多采用的是优化准则法（Optimality Criteria, OC）[12,13]，收敛速度快，但其一般适用于只单约束条件下的问题优化。ESO方法则是基于达尔文进化学说的思想，在优化迭代过程中不断删除低效的单元，以达到“进化”的目的。至于水平集法，其中一种方法就是在每次迭代时通过求解反应扩散方程确定边界的演化过程[14]。关于上述优化方法涉及的数值算法，内容复杂且篇幅过长，在此处就不展开来讲了，想要深入了解的可以阅读有关的文献。
确定了对连续体结构拓扑的描述，建立好结构的边界条件、优化目标和约束条件，选定合适的数值算法，便可以对结构进行优化了。下面给出了一个基于SIMP方法进行拓扑优化设计的基本流程图。
关于该流程图需要补充两点：一是在SIMP方法中，若以结构整体应变能最小（即刚度最大）为优化目标，则单元的应变能灵敏度将会是优化准则法的控制参数，根据该值及当前的设计变量（单元相对密度）来更新变量。二是过滤技术，这是为了处理优化过程中的数值不稳定现象而提出的一种方法，可以简单理解为对各材料单元的应变能灵敏度进行平均化修正。
本文只是简单地介绍了一下拓扑优化的基本概念以及目前流行的一些优化方法，并未涉及更深层次的优化过程等内容。若看了本文，有兴趣进一步入门学习拓扑优化的，可以看一下周平章博士在知乎写的一篇文章： 周平章：如何入门拓扑优化研究？
参考文献：
[1] Dron W.S., Gromory R.E., Greeberg H.J. Automatic design of optimal
structures[J]. http://J.de Mechanics, 1964(3):25~52.
[2] 袁代林,陈虬. 桁架结构拓扑优化的微粒群算法[J].西南交通大学学报,2007(01):94-98.
[3] 李祖宇. 面向建筑结构设计的水平集拓扑优化方法研究[D]. 华南理工大学, 2017.
[4] Sigmund O. Design of material structures using topology optimization[D]. Technical University of Denmark , 1994.
[5] Bendsoe MP and Sigmund O. Material interpolations in topology optimization[J]. Archive of Applied Mechanics, 1999, 69: 635~654.
[6] Theocaris PS and Stavroulakis GE. Optimal material design in composites: An iterative approach based on homogenizaed cells[J]. Computer Methods in Applied Mechanics and Engineering, 1999, 169: 31~42.
[7] Rietz A. Sufficiency of a finite exponent in SIMP (power law) methods[J]. Structural and Multidisciplinary Optimization, 2001, 21, 159-163.
[8] STOLPEM and SVANBERG K. An alternative interpolation scheme for minimum compliance topology optimization[J]. Structural and Multidisciplinary Optimization, 2001, 22: 116~124.
[9] Xie Y.M. and Steven G.P. A simple evolutionary procedure for structural
optimization[J]. Computers and Structures, 1993, 49: 885～896.
[10] Osher S，Sethian J A, Fronts propagating with curvature dependent speed：algorithms based on Hamilton-Jacobi formulations, Journal of Computational Physics, 1988, 79:12-49.
[11] Sethian J.A. and Wiegann A. Structural boundary design via level set and immersed interface methods[J]. Journal of Computational Physics, 2000, 163(2)：489—528．
[12] 钱令希．工程结构优化设计[M]．北京：水利电力出版社，1983.
[13] 程耿东．工程结构优化设计基础[M]．北京：水利电力出版社，1984.
[14] Yamada T, Izui K, Nishiwaki S, et al. A topology optimization method based on the level set method incorporating a fictitious interface energy[J]. Computer Methods inApplied Mechanics & Engineering, 2010, 199(45-48):2876-2891
文章被以下专栏收录
