连续体结构拓扑优化研究现状探讨

1.结构优化发展的回顾

自然界提供给人类可利用的资源是有限的，因此材料的有效利用一直是人类追求的目标，也是许多研究领域不变的话题，并伴随着结构优化理论和方法的产生而发展。结构优化的实质就是在满足各种约束条件（如应力、位移、频率、重量等）的情况下，寻找一个使得目标函数（如重量、造价等）最小的设计方案。结构优化的发展始终与结构分析、数学规划理论，以及计算速度的发展同步。

计算机出现以前，由于缺乏高速的计算工具来进行分析验证，以及系统的方法用于指导和改进结构设计，结构的优化是依靠人们世世代代累积起来的经验进行的。结构优化设计的理论，最早可以追溯到17世纪伽利略和伯努利对变截面梁形状优化问题的研究。那个时代的计算是靠手工，而且数学理论还不完备。结构优化领域中，第一个最重要的贡献是由Maxwell（1890）和Michell（1905）提出单载荷仅有应力约束条件下最小重量桁架结构布局的基本理论。

20世纪40、50年代，大量的研究工作集中于飞机结构的最小重量设计。

这期间提出了满应力设计（Full Stress Design，FSD）准则，它起始于“同步失效”概念，也就是结构的各个组成部分同时达到容许强度或失稳安全限度，由此得出一组联立方程，它的解析解提供了构件截面的优化尺寸。这样可以充分利用材料，得到最轻的结构设计方案。优化准则法的理论基础是由Prager和Sheu建立的，他们的研究结果表明用连续体问题导出的微分方程作为优化准则，微分方程的解就确定了结构的最优形状[200-201]。

20世纪60年代后，借助计算机的计算能力和数值优化方法的快速发展，结构优化领域才有了进展。有限元法和数学规划的日趋成熟确保了作为设计工具的结构优化方法的出现。1960年，Schmit在桁架优化方面的工作是现代结构优化诞生的标志[202]先在结构优化中引入数学规划理论，并与有限元方法结合，求解多种载荷工况下弹性结构的最小重量问题。

2.结构优化的分类

按照设计变量类型和优化层次，结构优化可以分为截面尺寸优化、形状优化和拓扑优化，如图1.6所示。从其发展的历程而言，优化技术也是按照这个顺序由低向高进行的。任何一个结构设计过程也都可分为3个阶段。首先是概念阶段，在这个阶段中选择结构系统和组成形式，还包括应用以往的经验来定性判断结构的行为；其次是初步的设计阶段，这个阶段决定结构的外形及几何形状；最后是详细的设计阶段，确定结构的具体尺寸。

结构优化的早期研究主要集中于截面尺寸优化问题，例如，寻找一个桁架各单元最佳的横截面尺寸，或者一个平板结构的最优厚度。截面尺寸优化是结构优化设计中最早开展研究的领域也是发展相对成熟的领域，它以截面尺寸为设计变量，通过寻求最优截面尺寸达到结构重量，体积或造价最小。约束条件可以是应力约束、位移约束、局部稳定、频率约束、动响应约束等。常用的截面尺寸优化方法有准则法和数学规划法两种。在这类问题中，设计区域是固定的，而且在优化过程中设计区域不发生改变。尺寸优化问题可以被看作是在详细的设计阶段开始执行优化计算的。

随着结构优化的进一步发展，人们开始考虑寻找结构的最优边界问题。

图1.6　结构优化的分类

（a）尺寸优化；（b）形状优化；（c）拓扑优化

形状优化是通过调整结构内外边界形状来改善结构性能和达到节省材料的目的。结构形状优化从对象上区分，主要有杆系结构和连续体结构。杆系结构的形状优化，一般选择节点坐标（位置）作为设计变量。对于连续体结构，结构的边界形状常采用适当的曲线/曲面方程、或一组基函数再附加可以自由变化的参数来描述，此时形状优化就可以选择这些自由参数作设计变量（特点是设计区域的形状可以变化），但是设计区域内的拓扑关系在形状优化过程中是固定不变的。这就要求在设计的初始阶段确定结构的最优拓扑形式，但是在没有任何先验知识的条件下，这往往很难做到。即使得到了最优的拓扑关系，在形状优化过程中还需要几次重新划分有限元网格，若要在设计过程中改变结构拓扑，则将使设计问题更加复杂。形状优化问题可以被看作是在初步的设计阶段开始执行优化计算的。

因为结构的尺寸和形状优化问题不能保证从一个最优的初始拓扑形式开始进行优化计算，所以尺寸和形状优化可能得到是一个次最优的结果而不是全局最优结果。为了弥补尺寸和形状优化固定拓扑的局限性，拓扑优化问题被提了出来[204-205]。这是继尺寸优化和形状优化后，在结构优化领域出现的一种富有挑战性的研究方向。

对于给定的设计问题，在没有任何尺寸、形状和拓扑等先验信息的条件下，拓扑优化必须在所谓的设计区域内确定材料的最优布局，这种布局能极小化给定的目标函数，并满足强加的约束。因此，拓扑优化问题被看作是材料分配问题，需要确定设计区域中每一点材料的特征[206]。其优化模型描述的困难和数值算法的巨大计算量使得连续体结构拓扑优化的发展较慢。自从Bends-e等提出基于均匀化理论的拓扑优化方法[207]后，优化这个领域才有了快速的发展，各种拓扑优化方法相继出现。典型的二维和三维连续体拓扑优化的目的是确定设计区域内空洞的数量及其位置等特征。拓扑优化的理想情况是同时得到结构的最优几何结构（即尺寸和形状）和拓扑，有时这类问题也被叫做布局优化。拓扑优化可以被看作是从概念阶段开始进行优化的。

众所周知，这种材料分配问题一般是病态的，通常没有解，这种病态性典型表现在网格的依赖性方面，也就是说，一个细的网格将产生好的设计，当然此设计不是最优的，因为更细的网格将产生更好的设计[208]。它的复杂性，拓扑优化是结构优化中更具挑战性的研究课题[209]。

3.结构拓扑优化

结构拓扑优化可分为离散体结构拓扑优化和连续体结构拓扑优化两大类[204]。离散体结构拓扑优化设计是从桁架结构研究开始的，可追溯到1904年Michell提出的桁架结构的设计理论，其设计结果中二力杆两两正交，杆的拉应力和压应力都达到最大。在一般的桁架拓扑优化问题中，通常假定外力、支承和节点己经给定，要求确定节点之间杆件的最优连接情况及杆件横截面积，使结构的重量或造价最小，同时满足应力、节点位移和结构柔顺性等功能要求。20世纪60年代，随着有限元方法的发展，人们开始数值求解较复杂结构的优化问题。1964年，Dorn等人提出了基结构法（Ground Structure Approach），将数值计算方法引入到拓扑优化领域中，克服了Michell桁架理论的局限性，从而使拓扑优化领域重新得到了发展。基结构法是目前桁架结构拓扑优化的主要方法，其原理就是把给定的初始设计区域离散成足够多的设计单元，对于由所有设计单元的节点组成的集合，在每两个节点之间用杆件连接起来，形成所谓的基结构，然后按照某种优化策略和准则从这个基结构中将某些不必要的单元删除掉，用保留下来的单元组成结构的最优拓扑形式。1974年，Prager应用解析方法求解了拓扑优化问题，奠定了这一研究领域的理论基础[210]。1990年，Rozvany和Kirsch等人建立了桁架结构拓扑优化的最优准则，使用这个准则可以从包含所有可行杆单元的基本结构中确定出最优结构[211-213]。在桁架结构拓扑优化中的奇异最优解的问题上，程耿东等人的研究指出杆件的应力函数在其零截面积处的不连续性是造成奇异最优解存在的根本原因，并提出了零截面积和杆件极限应力的概念。他认为奇异最优解所对应的设计点并非是设计空间中的孤立点，而是位于设计空间中某个低维退化子域的端点。整个设计空间仍然是连通的，其形状就像一个水母，由若干不同维数的可行子域组成。正是由于杆件极限应力的存在，才会在某些情况下影响杆件的删除，从而使得通常的数值优化算法无法获得问题的真正最优解。针对这一问题，人们提出了结构拓扑优化的ε-放松算法和延拓算法[214-215]。这一研究工作得到人们广泛的关注。(www.xing528.com)

连续体结构拓扑优化的重要发展源于1981年程耿东和Olhoff的工作，他们在研究最大刚度变厚度板最优设计时，发现最优解中包含许多各种尺寸的加强筋，具有非光滑的特征，这意味着在最优设计中必须引入复合材料以拓展优化设计空间[216，217]。这导致了随后的一系列研究进展，包括1984年Lurie等用G-收敛理论解释拓扑优化过程中的非光滑现象[218]；Kohn和Strang等引入松弛概念来处理拓扑优化中的病态变分问题，并讨论了这种松弛与均匀化之间的关系[219-221]；Rozvany等研究了在设计加强筋板中引进松弛的含义[222，223]。这些研究工作直接启发了1988年Bends-e和Kikuchi提出结构拓扑优化的均匀化方法，标志着对连续体结构拓扑优化的研究进入蓬勃发展的阶段[207]。

连续体结构拓扑优化可同时优化设计对象的拓扑、形状和尺寸。相对来说，它是快速发展的、年轻的研究领域，是一种创新性的设计方法。近期随着计算技术的快速发展，不同的技术和方法层出不穷。从描述拓扑优化问题方法的角度分析，连续体结构拓扑优化有两大类：基于微观结构的材料方法（Material Method）和基于宏观结构的几何法（Geometry Method），如图1.7所示。前者主要包括均匀化方法[224]、密度惩罚法[206]、渐进结构优化法[225]、独立连续映射模型方法[226-227]等；后者主要包括冒泡法[228]、拓扑函数方法[229]和水平集法[230-231]等。

在连续体结构拓扑优化问题的初始化阶段，基于微观结构的材料方法用具有微结构的孔状材料均匀填充整个设计区域（见图1.7左侧部分）。它通常使用固定的有限元网格来描述结构的几何形状，以及计算设计区域的位移场。通常设计区域都是长方形，有限元网格则是正方形，设计变量在每个单元内被认为是恒定不变的。有限元分析时，所用的材料属性与所在单元的密度是相关的，这个关系通常由微结构材料的物理模型确定，优化的目的在于确定连续体内每个单元内材料的有无状态。基于宏观结构的几何法使用各向同性材料填充整个设计区域（见图1.7右侧部分）。这类方法是在形状优化基础上进行的，设计区域内有限元网格不能固定，必须随着设计区域几何形状的变化不断重新划分网格。它是通过材料的增减和生成新孔洞来改变结构的拓扑关系。

图1.7　连续体结构拓扑优化分类

均匀化数学理论的发展，可以追溯到20世纪70年代，它主要用来预测复合材料的宏观材料属性，已经应用于很多物理和工程领域，如多孔介质中的流体流动、复合材料中的电磁场等[232]。事实上，力学和物理中关于连续介质的基本假设就可以看作是一种均匀化，因为材料是由原子或分子组成的。目前以均匀化理论为基础的优化方法已成为结构优化的主要途径[233]。它借助具有周期微结构的复合材料，将拓扑优化问题转化为复合材料微结构参数的尺寸设计问题，然后采用最优准则法或数学规划方法来寻找多孔介质的最优配置[234]。通过引入微结构，解决了材料分配问题的病态，它依靠微结构参数来松弛材料密度，使其可以在区间[0，1]内连续取值，代替了取值的离散集合{0，1}，这种设计空间的拓展，保证了拓扑优化最优解的存在。结构拓扑优化的均匀化方法已被许多研究者广泛应用，如最小柔度问题[235]、特征频率问题[236]和屈曲问题[237]等。但这种方法的不足之处在于，优化过程中需要确定最优的微结构和微结构的方向，有时显得过于繁琐；同时计算结果中常包含中间密度的微结构，这在实际工程应用中很难加工制造。然而这并不影响均匀化方法在结构拓扑优化中的重要位置，因为它能在理论上预测最佳的性能结构[238]。目前均匀化方法一般用于逆问题的求解和拓扑优化的理论研究方面，如解决材料微观胞元的优化设计问题。吴长春等用均匀化方法对蜂窝材料、复合材料的有效弹性模量等力学性能参数做了数值模拟计算，程耿东、刘书田、顾元宪等用均匀化理论对复合材料的性能进行了研究。

近年来，各国学者纷纷致力于对结构拓扑优化的均匀化方法进行改进，以进一步提高其计算效率。其中最具代表性的是固体各向同性材料惩罚法（Solid Isotropic Material with Penalization，SIMP），简称密度惩罚法。该方法以其理论简单明了，算法实现简单而著称，是目前主流的拓扑优化方法，也是工程上最有应用前景的一种拓扑优化方法[205]。该方法假设材料密度在单元内是常数，而材料属性用密度的指数函数来模拟。该方法提出后曾遭到质疑，因为没有相应的物理材料，它满足材料特性是其密度的幂函数，然而1999年Bends-e和Sigmund指出，只要满足简单的条件，例如指数大于等于3，泊松比等于1/3，很多情况下这个所谓的人工材料插值模型能够用微结构来模拟。为了保证解的存在性，在算法实施时密度惩罚法通常和边界结构长度约束、密度总变分约束或滤波技术相结合[208]。值得注意的是，结构拓扑优化的均匀化方法和密度惩罚法之间的细微差别：密度惩罚法接受中间密度值仅作为优化的工具，而均匀化方法接受中间密度值作为有效的设计特征，可以用微结构来实现[233]。目前密度惩罚法己被广泛应用于多约束、多物理和多材料的拓扑优化问题，因为这种问题中材料的各种宏观属性很难用均匀化方法求得。

结构拓扑优化的渐进结构优化法（Evolutionary Structural Optimization，ESO）包括硬杀（Hard-kill）和软杀（Soft-kill）两类方法，它源于等应力结构的设计。其基本思想是：设计区域中的非主动材料，也就是具有低应力或低应变能密度的材料，在结构中没有被充分利用，可以被删除，这个删除过程或者通过改变弹性模量把它作为应力的函数，或者直接从结构中删除低应力区的材料来实现[225]。人们已经设计了几种不同的演化算法，通常计算过程是反复进行有限元分析，按计算结果根据给定的材料演化策略，在低应力区删除一定比例的材料，或在高应力区添加一定比例的材料，逐渐逼近等应力结构，直至得到最优解，例如Xie和Steven设计的进化方法[244]，Baumgartner设计的软杀方法[245]，Hinton和Sienz设计的硬杀方法[246]，Querin和Reynolds设计的逆向自适应方法[247-248]等。设计这类优化方法的直接动机是为了避免在连续结构拓扑优化方法中使用过多的数学理论，而直接利用物理概念和计算机的计算能力来进行拓扑优化。同时，它具有与商业有限元程序衔接方便的优点，被广泛地应用于应力约束、频率约束[249]和屈曲约束[250-251]的连续体结构拓扑优化问题。但有的研究者指出，将这类方法拓展到非最小柔度设计，或多约束和多物理优化设计似乎存在问题，即使对于最大刚度结构设计问题，也通常得到非最优解[56]。

国内学者隋允康等[226，252-253]提出了一种独立连续映射模型方法（Independent Continuous Mapping，ICM），通过引入过滤函数、磨光函数及光滑映射变换，将拓扑变量从依附于截面、厚度等低层次变量中抽象出来，使之成为独立的层次，用[0，1]之间的连续变量代替{0，1}离散变量，从而避免了离散优化的困难。它采用磨光函数对拓扑变量进行磨光映射，再用过滤函数对拓扑变量进行过滤映射，通过映射反演完成拓扑变量“离散-连续-离散”的变化过程，磨光函数及过滤函数可取幂函数。它将桁架结构中的离散变量转换为连续设计变量，优化求解后再将连续设计变量转换为离散设计变量，这样在离散设计变量和连续设计变量之间建立了一一对应关系，从而建立完善了桁架结构拓扑优化模型，成功解决了多工况应力与位移约束下的桁架结构拓扑优化问题。在此基础上，将ICM方法推广到连续体结构拓扑优化中，研究了位移和应力约束下的连续体结构拓扑优化问题。建立以重量为目标的优化模型，克服了以柔顺度为目标难以处理多工况的不足，将桁架类结构和连续体结构拓扑优化建模方法进行统一[254-255]。

结构拓扑优化的冒泡法（Bubble）是Eschenauer等人于1994年提出的，目的是为了解决形状优化方法不能改变结构拓扑的问题[228]。当一个设计问题应用形状优化方法不能再进一步改进时，在结构中插入一个小圆孔，然后再应用形状优化方法来改变孔的形状和大小，以便更进一步地改进优化结果。

拓扑函数法采用了一种独立于有限元网格的拓扑结构表达方法，在设计区域内定义一个拓扑函数，此函数将设计区域内的每一点映射为一个数值，并以此值为基础选择一个阀值来明确决定单元材料有无。1999年，Sokolowski和Zochowski从数学角度结合Laplace方程和弹性平衡方程研究了二维开圆孔和三维开球形孔的拓扑导数理论[256]。2001年，Garreau等也开展了在弹性结构中开任意形状孔的拓扑导数理论的研究[257]。拓扑函数法的缺点是：虽然设计描述独立于网格，但有限元计算结果仍然依赖于网格，网格的变化将会导致计算结果的变化，因此设计目标函数值依赖于网格，特别是在网格细化的场合下，将导致优化算法的不稳定。另一个缺点是拓扑表达依赖于初始选定的几何体。

水平集方法（Level Set Method）求解结构拓扑优化问题的基本思想是引入一个水平集数φ=φ（x，t），它描述了结构的拓扑形式，即

式中，Ω代表结构实体材料当前所占有的区域，S是其边界，D是一个设计区域。

水平集函数包含了结构在优化过程中可能占有的一切构型，因此可以用水平集函数隐式地描述了结构拓扑，结构拓扑优化的过程是通过求解Hamilton-Jacobi方程来实现。

水平集方法主要是从界面传播等研究领域中逐步发展起来的，它是处理封闭运动界面随时间演化过程中几何拓扑变化的有效计算工具。Osher和Sethian[258]首先提出依赖时间的运动界面的水平集描述，其主要思想是将移动的界面作为零水平集嵌入高维的水平集标量函数中，这样由闭超曲面的演化方程可以得到水平集函数的演化方程，而嵌入的闭超曲面总是其零水平集，最终只要确定零水平集即可确定移动界面演化的结果。水平集方法是由Sethian和Wiegmann首先引入结构优化领域的[230]。采用水平集方法进行拓扑优化并使用响应泛函的灵敏度信息构造水平集发展所需要的速度场则是由Allaire[259]和Wang[260]等人的工作开始的。

梅玉林用水平集方法求解了结构拓扑优化、柔性机构设计和特定功能材料的设计，并在水平集拓扑优化的基础上，提出了违反约束向量水平集修正的返回映射算法，以便处理多材料、多约束和含有非设计区域的拓扑优化问题，并且将黎曼测度下的平均曲率流正则化技术引入到结构的拓扑优化中[231]。

贾海朋提出了水平集方法同进化算法相结合的改进的拓扑优化方法。该方法基于水平集方法实现边界演化，演化过程中采用进化算法中删除内部低应力区材料的策略生成新孔，消除了水平集方法不能生成孔的缺点，克服了常规水平集方法的拓扑形式对初始设计构型的依赖。同时，演化过程中根据需要自动生成新孔，提高了拓扑演化速度，改进了计算效率[261]。采用水平集方法求解连续体拓扑优化问题最大的优点在于它是以一种隐含的方式灵活地描述结构的拓扑变化。所有有关结构拓扑和结构边界的信息都体现在了水平集函数之中，在结构优化过程之中，无需显式地提取出结构的边界。