论文部分内容阅读
摘要：为实现更加先进的拓扑优化算法，研究采用反应扩散方程的水平集结构拓扑优化方法，通过理论推导给出算法中的参数选择建议。该方法允许在拓扑优化过程中生成新的孔洞，初始结构无须包含孔洞，不需要重新初始化步骤，从而可提高算法的收敛性。针对传统拓扑优化中主要采用体积约束、以柔度最小为目标和体积保留率设定存在一定主观性的问题，探究不同体积保留率下的结构应力水平的变化规律，结果显示可以依据结构最大应力水平与体积保留率的变化规律确定最优体积保留率。
关键词：
反应扩散方程; 水平集; 拓扑优化; 体积保留率; 收敛
中图分类号：TB115.2
文献标志码：B
Topology optimization method and its implement of
level set structure using reaction diffusion equation
LAI Song， WANG Kun
（School of Energy and Power Engineering， Huazhong University of Science and Technology， Wuhan 430064， China）
Abstract：
To realize a more advanced topology optimization algorithm， the level set structure topology optimization method using reaction diffusion equation is studied. The parameter selections in the algorithm are suggested by theoretical deduction. In this method， the generation of new holes in the process of topology optimization is allowed， and the initial structure does not need to set holes， and the initialization step need not restart， and then the convergence of the algorithm can be improved. In traditional topology optimization， the constraint is volume， and the objective is minimum flexibility， and the setting of volume retention rate is subjective. As to this issue， the variation law of structural stress level under different volume retention rate is explored. The results show that the optimal volume retention rate can be determined according to the variation law of the maximum stress level and volume retention rate of the structure.
Key words：
reaction diffusion equation; level set; topology optimization; volume retention rate; convergence
0 引 言
拓扑优化方法是在给定设计域中确定材料的分布、形状和连通性以达到某一最优目标的方法。[1]与尺寸优化和形状优化相比，拓撲优化有更大的设计自由度，因此可用于在概念设计阶段产生最佳的拓扑结构，进而得到可观的经济收益。在最近几十年里，结构拓扑优化取得长足发展，以变密度法为代表的拓扑优化方法广泛应用于结构、传热、流动和多物理场等领域。拓扑优化模块已经嵌入到某些商业有限元软件中，并应用于实际工程中。[2]在众多的拓扑优化方法中，最常见的是变密度法和水平集法。变密度法引入一个在0～1范围内连续变化的虚拟密度，并以此作为设计变量，表示材料的有无。由于虚拟密度是连续的，可以视作普通的场变量，所以变密度法与现有成熟的有限元法耦合良好。变密度法可以利用惩罚机制减少中间密度单元，但是无法得到真正的“0-1”分布，所以可能会在考虑材料非线性、应力约束等情况下出现问题。水平集法采用零水平集函数显式地表示结构边界，因此可以得到“0-1”分布，特别适用于考虑边界的拓扑优化。水平集方法的缺陷是过度依赖初始结构、迭代过程受限于时间步长、优化过程更加复杂等。
以柔度最小为目标的优化问题一般都具有良好的数学性质，因此大多数结构拓扑优化算法均以柔度最小为目标。但是，在实际工程结构设计中，有时工程师更看重的设计目标是强度。目前，与应力相关的优化目标在结构拓扑优化设计中还存在一定的困难。[3]一方面，应力的局部性质要求约束必须作用在每一个材料点上，导致约束的量非常大;另一方面，在基于密度的方法中，密度非常低的单元反而可能会出现较高的应力值，即应力奇异现象。
针对第一个问题，常见的解决方案是通过一个全局应力函数将大量局部应力约束拟合成一个全局约束，常见的全局函数有p-范数函数和K-S函数。LE等[4]将分区域的全局应力约束函数和自适应正则化方案相结合，在一定程度上可缓解应力集中现象。XIA等[5]提出改进的p-范数全局应力函数以控制最大应力，该全局函数可节约大量计算资源，但有时不能完全估计结构的局部应力特征，无法起到控制最大应力的作用。CHEN等[6]发现某些实际工程中使用最大主应力失效准则而不是VON MISES应力准则，因此将最大主应力作为约束，采用BESO作为优化求解器进行计算。EMMENDIERFER等[7]使用与SIMP法一样的应力松弛、p-范数聚集等技术，研究柔性机构的应力约束优化问题。PICELLI等[8]通过在每个迭代步中求解一个子问题得到优化的边界速度场，算例显示优化后的结构有非常平滑的边界。LIAN等[9]基于DSC方法更新网格，结合形状和拓扑优化解决最大VON MISES应力问题。 针对第二个问题，利用应力松弛方法，如qp-松弛法[3]、插值函数惩罚法[4，10]，可以消除应力奇异现象。但是，由于应力与结构的强非线性相关性，在尖角、孔洞和边界处会出现应力集中现象[11]，极易导致拓扑优化出现不稳定、收敛困难甚至无法收敛等问题。尽管如此，在工程应用中基于应力约束的拓扑优化方法只能应用于一些简单的结构，如板、壳、杆和混凝土结构等，而数学性质更好的柔度最小化目标仍然是更广泛使用的模型。
鉴于应力问题的复杂性，本文以结构柔度最小为优化目标，研究基于反应扩散方程的水平集结构拓扑优化方法，并考察方法中各参数对优化迭代过程与拓扑结构的影响。在不同体积约束下对下承式桥梁进行优化，分析得到结构应力与体积保留率的变化规律及其适用条件。
1 拓扑优化问题
1.1 基于水平集的拓扑结构描述
在水平集法中，结构拓扑的数学表达式是通过引入水平集函数实现的，见图1。当结构的边界由零水平集{x（x）=0}表示时，结构的拓扑可确定为
（x）<0， x∈D/Ω
（x）=0， x∈Ω
（x）>0， x∈Ω/Ω
（1）
式中：D为固定的设计域;Ω为实体材料域，Ω是D的子集;Ω是实体材料域边界。
定义Heaviside函数为
H（）=1， ≥0
0， <0
（2）
实体材料域Ω可以表示为
Ω={x|H（（x））=1}
（3）
2.2 结构场控制方程
假设实体域材料为各向同性线弹性材料，线弹性结构的微分控制方程为
-div（σ（u））=0，u∈Ω
u=0，u∈ΓD
σ（u）n=t，u∈ΓN
（4）
式中：u为位移场变量;ΓD为狄利克雷边界;ΓN为纽曼边界;t为作用在纽曼边界上的载荷，在本文中是一个常量，不随优化进程变化;σ（u）为应力张量，其定义为
σ（u）=Cε（u）
（5）
ε（u）=12（SymbolQC@u+（SymbolQC@u）T）
（6）
式中：C为弹性张量;ε（u）为应变张量。
2.3 柔度最小优化问题
优化目标为柔度最小，约束为体积约束，即
min J（u）=∫Ω（ε（u））TCε（u）dV
s.t. G（Ω）=∫ΩdV-V∫DdV=0
（7）
式中：V为给定的体积保留率。
2.5 增广拉格朗日乘子法
增广拉格朗日乘子法是大多数水平集法使用的优化算法，能将优化问题（式（7））转化为无约束问题并求解。增广拉格朗日函数定义为
L（u，λ，η）=J（u）+λG（Ω）+η2（G（Ω））2
（8）
式中：η为惩罚因子，在本文中，惩罚因子在整个迭代过程保持不变;λ为拉格朗日乘子，其更新方案为
λk+1=λk+ηG（Ω）
（9）
3 水平集函数更新
3.1 反应扩散方程
传统的水平集函数更新Hamilton-Jacobi方程无法生成孔洞，需要重新初始化步骤。为避免这一问题，可使用基于反应扩散方程的水平集更新方式[12]，即
t=αSymbolQC@2+Vn，∈D
（t=0）=0，∈D∪D
n=0，∈D/DN
=1，∈DN
（10）
式中：α為扩散系数;Vn为法向速度场;αSymbolQC@2为扩散项;DN为非设计域的边界。αSymbolQC@2可增加拓扑优化问题的正则性，促进水平集函数更新过程的稳定性，所以上式可以用标准伽辽金有限元法离散求解，以简化有限元法求解水平集方程的过程。此外，扩散项使形状灵敏度更加平滑，因此调整扩散系数α可以控制拓扑结构的复杂度，但扩散项有时也令灵敏度过于平坦，从而导致收敛出现问题，应综合考虑。
3.2 法向速度场的计算
法向速度场Vn可为水平集函数更新提供最高速下降方向，由增广拉格朗日函数的形状灵敏度求得。形状灵敏度反映边界移动对目标函数的影响，可以由物质导数公式和伴随方法[13]推导得到。
柔度目标的形状灵敏度为
δJδ=∫Ω（-（ε（u））TCε（u））VndΓ
（11）
体积约束的形状灵敏度为
δGδ=∫ΩVndΓ
（12）
因此
δLδ=∫Ω（-（ε（u））TCε（u）+λ+ηG（Ω））VndΓ
（13）
法向速度场为
Vn=V0-λ-ηG（Ω）
（14）
其中
V0=（ε（u））TCε（u）
（15）
4 数值实现
为能够在固定网格中运用有限元法求解位移场，位移场必须在整个设计域中都有定义，而不只是存在于实体材料域中。为将位移场拓展到整个设计域，使用弹性模量非常小的弱材料填充无材料区域D-Ω，避免整体刚度矩阵奇异。参照变密度法，采用虚拟密度ρ表示有限元单元的材料属性（有或无），因此材料的弹性模量Eρ可表示为
Eρ=Emin+ρ（E-Emin）
（16）
式中：E为实体材料的弹性模量;Emin为弱材料的弹性模量，一般取Emin=E×10-3。 离散后有限元单元Ωe的虚拟密度ρe定义为
ρe=∫ΩeH（）dV∫ΩedV
（17）
式（17）表示实体材料区域体积与整体单元体积的比，在数值计算时分3种情况：
（1）当单元全部是实体材料时，单元节点水平集函数值都是非负的，ρe=1。
（2）当单元全部是弱材料时，单元节点水平集函数值都是负值，ρe=0。
（3）当实体材料边界穿过单元时，ρe是0～1范围内的中间值，此值可用高斯积分法求得。需要注意的是，与变密度法相比，虚拟密度有严格的物理意义，而且只有边界穿过的单元是中间密度，不会出现中间密度单元过多的情况，因此不需要对中间密度进行惩罚。
求解位移场的有限元方程为
KU=F
（18）
式中：K为整体刚度矩阵;U为位移向量矩阵;F为载荷向量矩阵。
水平集方程采用有限元法求解，即
1ΔtK1+K2t+Δt=K11Δtt+∫ΩeVndV
（19）
K1=∫ΩeNTNdV
（20）
K2=∫Ωe（ΔN）T（ΔN）dV
（21）
式中：Δt为时间步长;t+Δt为Δt时间步长后的水平集函数，是待求量;t为当前时间的水平集函数值，是已知量;N为水平集函数的插值函数矩阵。
实际上，式（14）定义的速度场只存在于实体材料边界，因此为求解水平集方程，速度场必须扩展到整个设计域。采用自然扩展法，当位移场在整个设计域都有定义时，Vn自然扩展到整个设计域，
V0=（ε（u））TCε（u）
（22）
另外，在计算中，结构场的狄利克雷边界（位移量已知的边界）与纽曼边界（边界力不为0的边界）的速度值应设置为0，表示这2种边界不随拓扑变化。
为让各种不同问题的参数相差不大，可以对目标函数的速度场进行归一化处理，即
V0=V0max（V0）
（23）
迭代收敛判据有3个，只有同时满足3个条件才算收敛：（1）迭代次数大于20;（2）Jk+1-JkJk≤ε;（3）vol-V≤0.005，即当前体积分数vol与体积保留率V的差值小于或等于0.005。
完整的优化计算流程见图2。
5 数值算例与讨论
所有模型实体材料弹性模量均取E=1，弱材料弹性模量Emin=1×10-3，泊松比v=0.3。有限元求解结构场和水平集方程的网格单元均为双线性单元。水平集函数初始值为1，拉格朗日乘子λ初始值为0。
5.1 各参数对迭代过程和拓扑结构的影响
以二维悬臂梁为例，见图3。悬臂梁的尺寸为1.0 m×0.8 m，左端固定，在右端中间施加方向向下、大小为6 N的集中载荷，离散后悬臂梁网格数量为5 120个。
在有限元求解水平集方程的过程中，时间步长Δt是一个非常重要的参数。当水平集求解次数m=10、η=0.01、α=8×10-5、V=0.4时，取3个不同时间步长Δt时的优化结果对比见图4。3个时间步长下都可得到稳定的结果。
随着Δt增大，结构分支越来越小直至被完全删除，这说明时间步长越大复杂性越小，得到的拓扑结构越简单。这是因为时间步长越大，节点上的水平集函数值变化越大，结构变化更加剧烈，结构中分支的稳定性减弱，导致迭代过程中这些分支被删除。然而，过小的时间步长会让迭代次数增多，增加计算时间，因此也不一定可取。
求解位移场比求解水平集方程耗费更多的时间，因此为节约计算成本，传统水平集方法往往在一次迭代（求解一次位移场）中多次求解水平集方程，即m≥2。但是，在使用反应扩散方程更新水平集函数的方法中，时间步长并不受CFL条件的限制，因此m≥2的必要性值得商榷。研究发现，优化结果更多地取决于时间步长Δt与次数m的乘积（即总时间步长），而不仅仅是单个参数。如果两者乘积一致，那么得到的拓扑优化结构大致相同。固定总时间步长为3的同时取不同的m，需要的求解时间见表1。由此可知，求解位移场花费的时间远大于求解水平集方程花费的时间，增大m可以减少求解位移场的次数、减少计算时间，当网格数量较多且单元大小不一致时，这种效果会更加显著，这是因为刚度矩阵组装会花费更多的时间。但是，过大的m也不可取，例如当m=20时的迭代次数与m=10时相比并没有减少。综上所述，可以取较大的时间步长，若将总时间步长分解成m个小时间步长求解的方法可以减少计算时间，还是适用的，然而m不宜过大，应根据实际迭代情况适当调整。
惩罚因子η的取值也很重要，不同η时的优化结构对比见图5，不同η时的迭代曲线见图6。
过大的惩罚因子会使优化陷入局部最优甚至无法收敛的情况，当η=0.020时算法收敛后得到一个无用的结构;过小的惩罚因子导致迭代次数增加，η=0.002时算法迭代197步收敛，而η=0.010时只需要迭代67步。同时，小惩罚因子的迭代曲线更加平稳。综合考虑稳定性和效率，在一般情况下，η=0.010比较合适。当存在较严重的应力集中现象时，为确保迭代的稳定性，惩罚因子应适当减小。
扩散系数α可以调整优化结构的复杂度，α越大得到的结构越简单，实际应用时可根据需要调整，具体可参见文献[12]。此外，根据单元中心点水平集函数值的正负情况判断单元材料有无的几何映射方案可以不出现中间密度单元，但依然可以得到稳定的结果，只是需要比有中间密度的方案花费更多的迭代次数。
5.2 基于反应扩散方程的水平集法与变密度法对比
以悬臂梁为例，网格数量为160×80=12 800个，左端固定，右下角施加1 N竖直向下的载荷，优化参数取Δt=0.10、m=10、η=0.010、α=8×10-5、V=0.4。变密度法采用ANDREASSEN等[14]的方案，使用PDE过滤器，过滤半径r=3。2种方法的拓扑优化结果对比见图7。水平集法得到的结构边界更清晰，拓扑更复杂，结构的柔度也更小。水平集法的结构复杂度可以通过扩散系数调节，变密度法除过滤半径r外没有其他参数可以影响结构拓扑复杂性，雖然改变r可以得到不同的拓扑结构，但中间密度单元的数量和比例也会改变，对计算结果不利。 2種方法的迭代曲线见图8。变密度法收敛只迭代64次，而水平集法迭代168次才收敛，说明变密度法的收敛过程在稳定性和速度方面要优于水平集法。因此，水平集法仍然需要在算法方面进一步研究以提高其稳定性和效率。
5.3 不同体积保留率V下最优结构的最大等效应力变化规律
以某桥梁模型拓扑优化为例，其初始结构见图9（a），尺寸为10 m×4 m×4 m，桥梁下底面两端固定，桥面施加10 N/m2方向向下的面载荷。桥面以下的区域为非设计域，占总体积的6.25%。不同体积保留率下的优化结果分别见图9（b）～（d）。
不同体积保留率V下优化结构的最大应力见图10。由此可知：随着V增大，优化结果的最大等效应力一直减小;当V<0.12时，最大等效应力随V的减小迅速增大;当V>0.12时，最大等效应力减小幅度不大。这是因为随着体积分数的减小，一些加强性结构被删除，导致结构拓扑发生较大的改变，从而造成最大等效应力迅速增大。V越小优化越困难，所花费的时间越多，因此V不应过小，在本文模型中最佳的V为0.12。因此，针对具体的优化问题，体积保留率应该存在一个最佳值，可以使材料消耗与结构应力处于最优搭配。
6 结束语
针对以柔度最小为目标的结构优化问题，实现基于反应扩散方程的水平集法，重点考察时间步长和水平集方程求解次数这2个影响参数，结果认为水平集法的时间步长不受CFL条件限制，无论时间步长为何值，都能得到稳定的优化结果。然而，总时间步长（即时间步长与求解次数的乘积）过大会使优化结构更加简单。当总时间步长一定时，适当增加水平集求解次数可以减少迭代次数、节约时间。将基于反应扩散方程的水平集法与变密度法进行对比，水平集法能够得到更清晰的边界结构，能够调节结构复杂度但收敛缓慢。根据不同体积保留率下最优结构的最大等效应力变化规律可知，针对具体的优化问题，体积保留率应该存在一个最佳值，可以使材料消耗与结构应力处于最优搭配，该值可以通过考察体积保留率与最大应力或结构柔度之间的变化规律获得。
参考文献：
[1] VAN DIJK N P， MAUTE K， LANGELAAR M， et al. Level-set methods for structural topology optimization： A review[J]. Structural and Multidisciplinary Optimization， 2013， 48（3）： 437-472. DOI： 10.1007/s00158-013-0912-y.
[2] 张召颖， 张帆， 邹洵， 等. 基于Ansys Workbench的T形结构优化设计[J]. 计算机辅助工程， 2019， 28（3）： 35-38. DOI： 10.13340/j.cae.2019.03.007.
[3] DUYSINX P， BENDSE M P. Topology optimization of continuum structures with local stress constraints[J]. International Journal for Numerical Methods in Engineering， 1998， 43（8）： 1453-1478. DOI： 10.1002/（SICI）1097-0207（19981230）43：8〈1453：：AID-NME480〉3.0.CO;2-2.
[4] LE C， NORATO J， BRUNS T， et al. Stress-based topology optimization for continua[J]. Structural and Multidisciplinary Optimization， 2010， 41（4）： 605-620. DOI： 10.1007/s00158-009-0440-y.
[5] XIA Q， SHI T L， LIU S Y， et al. A level set solution to stress-based structural shape and topology optimization[J]. Computers & Structures， 2012， 90/91： 55-64. DOI： 10.1016/j.compstruc.2011.10.009.
[6] CHEN A B， CAI K， ZHAO Z L， et al. Controlling maximum first principal stress in topology optimization[J]. Structural and Multidisciplinary Optimization， 2020， 63： 327-339. DOI： 10.1007/s00158-020-02701-5.
[7] EMMENDIERFER H， FANCELLO E A， SILVA E C N. Stress-constrained level set topology optimization for compliant mechanisms[J]. Computer Methods in Applied Mechanics and Engineering， 2020， 362： 112777. DOI： 10.1016/j.cma.2019.112777.
[8] PICELLI R， TOWNSEND S， BRAMPTON C， et al. Stress-based shape and topology optimization with level set method[J]. Computer Methods in Applied Mechanics and Engineering， 2018， 329： 1-23. DOI： 10.1016/j.cma.2017.09.001. [9] LIAN H J， CHRISTIANSEN A N， TORTORELLI D A， et al. Combined shape and topology optimization for minimization of maximal VON MISES stress[J]. Structural and Multidisciplinary Optimization， 2017， 55（5）： 1541-1557. DOI： 10.1007/s00158-017-1656-x.
[10] 占金青， 龍良明， 刘敏. 热弹性结构全局应力约束拓扑优化设计[J]. 机械科学与技术， 2019， 38（9）： 1386-1392. DOI： 10.13433/j.cnki.1003-8728.20190005.
[11] 张维生. 基于水平集方法的应力相关拓扑优化问题研究[D]. 大连： 大连理工大学， 2013.
[12] YAMADA T， IZUI K， NISHIWAKI S， et al. A topology optimization method based on level set method incorporating a fictitious interfaceenergy[J]. Computer Methods in Applied Mechanics and Engineering， 2010， 199（45-48）： 2876-2891. DOI： 10.1016/j.cma.2010.05.013.
[13] ALLAIRE G， JOUVE F， TOADER A M. Structural optimization using sensitivity analysis and a level-set method[J]. Journal of Computational Physics， 2004， 194（1）： 363-393. DOI： 10.1016/j.jcp.2003.09.032.
[14] ANDREASSEN E， CLAUSEN A， SCHEVENELS M， et al. Efficient topology optimization in MATLAB using 88 lines of code[J]. Structural and Multidisciplinary Optimization， 2011， 43（1）： 1-16. DOI： 10.1007/s00158-010-0594-7.
（编辑 武晓英）
关键词：
反应扩散方程; 水平集; 拓扑优化; 体积保留率; 收敛
中图分类号：TB115.2
文献标志码：B
Topology optimization method and its implement of
level set structure using reaction diffusion equation
LAI Song， WANG Kun
（School of Energy and Power Engineering， Huazhong University of Science and Technology， Wuhan 430064， China）
Abstract：
To realize a more advanced topology optimization algorithm， the level set structure topology optimization method using reaction diffusion equation is studied. The parameter selections in the algorithm are suggested by theoretical deduction. In this method， the generation of new holes in the process of topology optimization is allowed， and the initial structure does not need to set holes， and the initialization step need not restart， and then the convergence of the algorithm can be improved. In traditional topology optimization， the constraint is volume， and the objective is minimum flexibility， and the setting of volume retention rate is subjective. As to this issue， the variation law of structural stress level under different volume retention rate is explored. The results show that the optimal volume retention rate can be determined according to the variation law of the maximum stress level and volume retention rate of the structure.
Key words：
reaction diffusion equation; level set; topology optimization; volume retention rate; convergence
0 引 言
拓扑优化方法是在给定设计域中确定材料的分布、形状和连通性以达到某一最优目标的方法。[1]与尺寸优化和形状优化相比，拓撲优化有更大的设计自由度，因此可用于在概念设计阶段产生最佳的拓扑结构，进而得到可观的经济收益。在最近几十年里，结构拓扑优化取得长足发展，以变密度法为代表的拓扑优化方法广泛应用于结构、传热、流动和多物理场等领域。拓扑优化模块已经嵌入到某些商业有限元软件中，并应用于实际工程中。[2]在众多的拓扑优化方法中，最常见的是变密度法和水平集法。变密度法引入一个在0～1范围内连续变化的虚拟密度，并以此作为设计变量，表示材料的有无。由于虚拟密度是连续的，可以视作普通的场变量，所以变密度法与现有成熟的有限元法耦合良好。变密度法可以利用惩罚机制减少中间密度单元，但是无法得到真正的“0-1”分布，所以可能会在考虑材料非线性、应力约束等情况下出现问题。水平集法采用零水平集函数显式地表示结构边界，因此可以得到“0-1”分布，特别适用于考虑边界的拓扑优化。水平集方法的缺陷是过度依赖初始结构、迭代过程受限于时间步长、优化过程更加复杂等。
以柔度最小为目标的优化问题一般都具有良好的数学性质，因此大多数结构拓扑优化算法均以柔度最小为目标。但是，在实际工程结构设计中，有时工程师更看重的设计目标是强度。目前，与应力相关的优化目标在结构拓扑优化设计中还存在一定的困难。[3]一方面，应力的局部性质要求约束必须作用在每一个材料点上，导致约束的量非常大;另一方面，在基于密度的方法中，密度非常低的单元反而可能会出现较高的应力值，即应力奇异现象。
针对第一个问题，常见的解决方案是通过一个全局应力函数将大量局部应力约束拟合成一个全局约束，常见的全局函数有p-范数函数和K-S函数。LE等[4]将分区域的全局应力约束函数和自适应正则化方案相结合，在一定程度上可缓解应力集中现象。XIA等[5]提出改进的p-范数全局应力函数以控制最大应力，该全局函数可节约大量计算资源，但有时不能完全估计结构的局部应力特征，无法起到控制最大应力的作用。CHEN等[6]发现某些实际工程中使用最大主应力失效准则而不是VON MISES应力准则，因此将最大主应力作为约束，采用BESO作为优化求解器进行计算。EMMENDIERFER等[7]使用与SIMP法一样的应力松弛、p-范数聚集等技术，研究柔性机构的应力约束优化问题。PICELLI等[8]通过在每个迭代步中求解一个子问题得到优化的边界速度场，算例显示优化后的结构有非常平滑的边界。LIAN等[9]基于DSC方法更新网格，结合形状和拓扑优化解决最大VON MISES应力问题。 针对第二个问题，利用应力松弛方法，如qp-松弛法[3]、插值函数惩罚法[4，10]，可以消除应力奇异现象。但是，由于应力与结构的强非线性相关性，在尖角、孔洞和边界处会出现应力集中现象[11]，极易导致拓扑优化出现不稳定、收敛困难甚至无法收敛等问题。尽管如此，在工程应用中基于应力约束的拓扑优化方法只能应用于一些简单的结构，如板、壳、杆和混凝土结构等，而数学性质更好的柔度最小化目标仍然是更广泛使用的模型。
鉴于应力问题的复杂性，本文以结构柔度最小为优化目标，研究基于反应扩散方程的水平集结构拓扑优化方法，并考察方法中各参数对优化迭代过程与拓扑结构的影响。在不同体积约束下对下承式桥梁进行优化，分析得到结构应力与体积保留率的变化规律及其适用条件。
1 拓扑优化问题
1.1 基于水平集的拓扑结构描述
在水平集法中，结构拓扑的数学表达式是通过引入水平集函数实现的，见图1。当结构的边界由零水平集{x（x）=0}表示时，结构的拓扑可确定为
（x）<0， x∈D/Ω
（x）=0， x∈Ω
（x）>0， x∈Ω/Ω
（1）
式中：D为固定的设计域;Ω为实体材料域，Ω是D的子集;Ω是实体材料域边界。
定义Heaviside函数为
H（）=1， ≥0
0， <0
（2）
实体材料域Ω可以表示为
Ω={x|H（（x））=1}
（3）
2.2 结构场控制方程
假设实体域材料为各向同性线弹性材料，线弹性结构的微分控制方程为
-div（σ（u））=0，u∈Ω
u=0，u∈ΓD
σ（u）n=t，u∈ΓN
（4）
式中：u为位移场变量;ΓD为狄利克雷边界;ΓN为纽曼边界;t为作用在纽曼边界上的载荷，在本文中是一个常量，不随优化进程变化;σ（u）为应力张量，其定义为
σ（u）=Cε（u）
（5）
ε（u）=12（SymbolQC@u+（SymbolQC@u）T）
（6）
式中：C为弹性张量;ε（u）为应变张量。
2.3 柔度最小优化问题
优化目标为柔度最小，约束为体积约束，即
min J（u）=∫Ω（ε（u））TCε（u）dV
s.t. G（Ω）=∫ΩdV-V∫DdV=0
（7）
式中：V为给定的体积保留率。
2.5 增广拉格朗日乘子法
增广拉格朗日乘子法是大多数水平集法使用的优化算法，能将优化问题（式（7））转化为无约束问题并求解。增广拉格朗日函数定义为
L（u，λ，η）=J（u）+λG（Ω）+η2（G（Ω））2
（8）
式中：η为惩罚因子，在本文中，惩罚因子在整个迭代过程保持不变;λ为拉格朗日乘子，其更新方案为
λk+1=λk+ηG（Ω）
（9）
3 水平集函数更新
3.1 反应扩散方程
传统的水平集函数更新Hamilton-Jacobi方程无法生成孔洞，需要重新初始化步骤。为避免这一问题，可使用基于反应扩散方程的水平集更新方式[12]，即
t=αSymbolQC@2+Vn，∈D
（t=0）=0，∈D∪D
n=0，∈D/DN
=1，∈DN
（10）
式中：α為扩散系数;Vn为法向速度场;αSymbolQC@2为扩散项;DN为非设计域的边界。αSymbolQC@2可增加拓扑优化问题的正则性，促进水平集函数更新过程的稳定性，所以上式可以用标准伽辽金有限元法离散求解，以简化有限元法求解水平集方程的过程。此外，扩散项使形状灵敏度更加平滑，因此调整扩散系数α可以控制拓扑结构的复杂度，但扩散项有时也令灵敏度过于平坦，从而导致收敛出现问题，应综合考虑。
3.2 法向速度场的计算
法向速度场Vn可为水平集函数更新提供最高速下降方向，由增广拉格朗日函数的形状灵敏度求得。形状灵敏度反映边界移动对目标函数的影响，可以由物质导数公式和伴随方法[13]推导得到。
柔度目标的形状灵敏度为
δJδ=∫Ω（-（ε（u））TCε（u））VndΓ
（11）
体积约束的形状灵敏度为
δGδ=∫ΩVndΓ
（12）
因此
δLδ=∫Ω（-（ε（u））TCε（u）+λ+ηG（Ω））VndΓ
（13）
法向速度场为
Vn=V0-λ-ηG（Ω）
（14）
其中
V0=（ε（u））TCε（u）
（15）
4 数值实现
为能够在固定网格中运用有限元法求解位移场，位移场必须在整个设计域中都有定义，而不只是存在于实体材料域中。为将位移场拓展到整个设计域，使用弹性模量非常小的弱材料填充无材料区域D-Ω，避免整体刚度矩阵奇异。参照变密度法，采用虚拟密度ρ表示有限元单元的材料属性（有或无），因此材料的弹性模量Eρ可表示为
Eρ=Emin+ρ（E-Emin）
（16）
式中：E为实体材料的弹性模量;Emin为弱材料的弹性模量，一般取Emin=E×10-3。 离散后有限元单元Ωe的虚拟密度ρe定义为
ρe=∫ΩeH（）dV∫ΩedV
（17）
式（17）表示实体材料区域体积与整体单元体积的比，在数值计算时分3种情况：
（1）当单元全部是实体材料时，单元节点水平集函数值都是非负的，ρe=1。
（2）当单元全部是弱材料时，单元节点水平集函数值都是负值，ρe=0。
（3）当实体材料边界穿过单元时，ρe是0～1范围内的中间值，此值可用高斯积分法求得。需要注意的是，与变密度法相比，虚拟密度有严格的物理意义，而且只有边界穿过的单元是中间密度，不会出现中间密度单元过多的情况，因此不需要对中间密度进行惩罚。
求解位移场的有限元方程为
KU=F
（18）
式中：K为整体刚度矩阵;U为位移向量矩阵;F为载荷向量矩阵。
水平集方程采用有限元法求解，即
1ΔtK1+K2t+Δt=K11Δtt+∫ΩeVndV
（19）
K1=∫ΩeNTNdV
（20）
K2=∫Ωe（ΔN）T（ΔN）dV
（21）
式中：Δt为时间步长;t+Δt为Δt时间步长后的水平集函数，是待求量;t为当前时间的水平集函数值，是已知量;N为水平集函数的插值函数矩阵。
实际上，式（14）定义的速度场只存在于实体材料边界，因此为求解水平集方程，速度场必须扩展到整个设计域。采用自然扩展法，当位移场在整个设计域都有定义时，Vn自然扩展到整个设计域，
V0=（ε（u））TCε（u）
（22）
另外，在计算中，结构场的狄利克雷边界（位移量已知的边界）与纽曼边界（边界力不为0的边界）的速度值应设置为0，表示这2种边界不随拓扑变化。
为让各种不同问题的参数相差不大，可以对目标函数的速度场进行归一化处理，即
V0=V0max（V0）
（23）
迭代收敛判据有3个，只有同时满足3个条件才算收敛：（1）迭代次数大于20;（2）Jk+1-JkJk≤ε;（3）vol-V≤0.005，即当前体积分数vol与体积保留率V的差值小于或等于0.005。
完整的优化计算流程见图2。
5 数值算例与讨论
所有模型实体材料弹性模量均取E=1，弱材料弹性模量Emin=1×10-3，泊松比v=0.3。有限元求解结构场和水平集方程的网格单元均为双线性单元。水平集函数初始值为1，拉格朗日乘子λ初始值为0。
5.1 各参数对迭代过程和拓扑结构的影响
以二维悬臂梁为例，见图3。悬臂梁的尺寸为1.0 m×0.8 m，左端固定，在右端中间施加方向向下、大小为6 N的集中载荷，离散后悬臂梁网格数量为5 120个。
在有限元求解水平集方程的过程中，时间步长Δt是一个非常重要的参数。当水平集求解次数m=10、η=0.01、α=8×10-5、V=0.4时，取3个不同时间步长Δt时的优化结果对比见图4。3个时间步长下都可得到稳定的结果。
随着Δt增大，结构分支越来越小直至被完全删除，这说明时间步长越大复杂性越小，得到的拓扑结构越简单。这是因为时间步长越大，节点上的水平集函数值变化越大，结构变化更加剧烈，结构中分支的稳定性减弱，导致迭代过程中这些分支被删除。然而，过小的时间步长会让迭代次数增多，增加计算时间，因此也不一定可取。
求解位移场比求解水平集方程耗费更多的时间，因此为节约计算成本，传统水平集方法往往在一次迭代（求解一次位移场）中多次求解水平集方程，即m≥2。但是，在使用反应扩散方程更新水平集函数的方法中，时间步长并不受CFL条件的限制，因此m≥2的必要性值得商榷。研究发现，优化结果更多地取决于时间步长Δt与次数m的乘积（即总时间步长），而不仅仅是单个参数。如果两者乘积一致，那么得到的拓扑优化结构大致相同。固定总时间步长为3的同时取不同的m，需要的求解时间见表1。由此可知，求解位移场花费的时间远大于求解水平集方程花费的时间，增大m可以减少求解位移场的次数、减少计算时间，当网格数量较多且单元大小不一致时，这种效果会更加显著，这是因为刚度矩阵组装会花费更多的时间。但是，过大的m也不可取，例如当m=20时的迭代次数与m=10时相比并没有减少。综上所述，可以取较大的时间步长，若将总时间步长分解成m个小时间步长求解的方法可以减少计算时间，还是适用的，然而m不宜过大，应根据实际迭代情况适当调整。
惩罚因子η的取值也很重要，不同η时的优化结构对比见图5，不同η时的迭代曲线见图6。
过大的惩罚因子会使优化陷入局部最优甚至无法收敛的情况，当η=0.020时算法收敛后得到一个无用的结构;过小的惩罚因子导致迭代次数增加，η=0.002时算法迭代197步收敛，而η=0.010时只需要迭代67步。同时，小惩罚因子的迭代曲线更加平稳。综合考虑稳定性和效率，在一般情况下，η=0.010比较合适。当存在较严重的应力集中现象时，为确保迭代的稳定性，惩罚因子应适当减小。
扩散系数α可以调整优化结构的复杂度，α越大得到的结构越简单，实际应用时可根据需要调整，具体可参见文献[12]。此外，根据单元中心点水平集函数值的正负情况判断单元材料有无的几何映射方案可以不出现中间密度单元，但依然可以得到稳定的结果，只是需要比有中间密度的方案花费更多的迭代次数。
5.2 基于反应扩散方程的水平集法与变密度法对比
以悬臂梁为例，网格数量为160×80=12 800个，左端固定，右下角施加1 N竖直向下的载荷，优化参数取Δt=0.10、m=10、η=0.010、α=8×10-5、V=0.4。变密度法采用ANDREASSEN等[14]的方案，使用PDE过滤器，过滤半径r=3。2种方法的拓扑优化结果对比见图7。水平集法得到的结构边界更清晰，拓扑更复杂，结构的柔度也更小。水平集法的结构复杂度可以通过扩散系数调节，变密度法除过滤半径r外没有其他参数可以影响结构拓扑复杂性，雖然改变r可以得到不同的拓扑结构，但中间密度单元的数量和比例也会改变，对计算结果不利。 2種方法的迭代曲线见图8。变密度法收敛只迭代64次，而水平集法迭代168次才收敛，说明变密度法的收敛过程在稳定性和速度方面要优于水平集法。因此，水平集法仍然需要在算法方面进一步研究以提高其稳定性和效率。
5.3 不同体积保留率V下最优结构的最大等效应力变化规律
以某桥梁模型拓扑优化为例，其初始结构见图9（a），尺寸为10 m×4 m×4 m，桥梁下底面两端固定，桥面施加10 N/m2方向向下的面载荷。桥面以下的区域为非设计域，占总体积的6.25%。不同体积保留率下的优化结果分别见图9（b）～（d）。
不同体积保留率V下优化结构的最大应力见图10。由此可知：随着V增大，优化结果的最大等效应力一直减小;当V<0.12时，最大等效应力随V的减小迅速增大;当V>0.12时，最大等效应力减小幅度不大。这是因为随着体积分数的减小，一些加强性结构被删除，导致结构拓扑发生较大的改变，从而造成最大等效应力迅速增大。V越小优化越困难，所花费的时间越多，因此V不应过小，在本文模型中最佳的V为0.12。因此，针对具体的优化问题，体积保留率应该存在一个最佳值，可以使材料消耗与结构应力处于最优搭配。
6 结束语
针对以柔度最小为目标的结构优化问题，实现基于反应扩散方程的水平集法，重点考察时间步长和水平集方程求解次数这2个影响参数，结果认为水平集法的时间步长不受CFL条件限制，无论时间步长为何值，都能得到稳定的优化结果。然而，总时间步长（即时间步长与求解次数的乘积）过大会使优化结构更加简单。当总时间步长一定时，适当增加水平集求解次数可以减少迭代次数、节约时间。将基于反应扩散方程的水平集法与变密度法进行对比，水平集法能够得到更清晰的边界结构，能够调节结构复杂度但收敛缓慢。根据不同体积保留率下最优结构的最大等效应力变化规律可知，针对具体的优化问题，体积保留率应该存在一个最佳值，可以使材料消耗与结构应力处于最优搭配，该值可以通过考察体积保留率与最大应力或结构柔度之间的变化规律获得。
参考文献：
[1] VAN DIJK N P， MAUTE K， LANGELAAR M， et al. Level-set methods for structural topology optimization： A review[J]. Structural and Multidisciplinary Optimization， 2013， 48（3）： 437-472. DOI： 10.1007/s00158-013-0912-y.
[2] 张召颖， 张帆， 邹洵， 等. 基于Ansys Workbench的T形结构优化设计[J]. 计算机辅助工程， 2019， 28（3）： 35-38. DOI： 10.13340/j.cae.2019.03.007.
[3] DUYSINX P， BENDSE M P. Topology optimization of continuum structures with local stress constraints[J]. International Journal for Numerical Methods in Engineering， 1998， 43（8）： 1453-1478. DOI： 10.1002/（SICI）1097-0207（19981230）43：8〈1453：：AID-NME480〉3.0.CO;2-2.
[4] LE C， NORATO J， BRUNS T， et al. Stress-based topology optimization for continua[J]. Structural and Multidisciplinary Optimization， 2010， 41（4）： 605-620. DOI： 10.1007/s00158-009-0440-y.
[5] XIA Q， SHI T L， LIU S Y， et al. A level set solution to stress-based structural shape and topology optimization[J]. Computers & Structures， 2012， 90/91： 55-64. DOI： 10.1016/j.compstruc.2011.10.009.
[6] CHEN A B， CAI K， ZHAO Z L， et al. Controlling maximum first principal stress in topology optimization[J]. Structural and Multidisciplinary Optimization， 2020， 63： 327-339. DOI： 10.1007/s00158-020-02701-5.
[7] EMMENDIERFER H， FANCELLO E A， SILVA E C N. Stress-constrained level set topology optimization for compliant mechanisms[J]. Computer Methods in Applied Mechanics and Engineering， 2020， 362： 112777. DOI： 10.1016/j.cma.2019.112777.
[8] PICELLI R， TOWNSEND S， BRAMPTON C， et al. Stress-based shape and topology optimization with level set method[J]. Computer Methods in Applied Mechanics and Engineering， 2018， 329： 1-23. DOI： 10.1016/j.cma.2017.09.001. [9] LIAN H J， CHRISTIANSEN A N， TORTORELLI D A， et al. Combined shape and topology optimization for minimization of maximal VON MISES stress[J]. Structural and Multidisciplinary Optimization， 2017， 55（5）： 1541-1557. DOI： 10.1007/s00158-017-1656-x.
[10] 占金青， 龍良明， 刘敏. 热弹性结构全局应力约束拓扑优化设计[J]. 机械科学与技术， 2019， 38（9）： 1386-1392. DOI： 10.13433/j.cnki.1003-8728.20190005.
[11] 张维生. 基于水平集方法的应力相关拓扑优化问题研究[D]. 大连： 大连理工大学， 2013.
[12] YAMADA T， IZUI K， NISHIWAKI S， et al. A topology optimization method based on level set method incorporating a fictitious interfaceenergy[J]. Computer Methods in Applied Mechanics and Engineering， 2010， 199（45-48）： 2876-2891. DOI： 10.1016/j.cma.2010.05.013.
[13] ALLAIRE G， JOUVE F， TOADER A M. Structural optimization using sensitivity analysis and a level-set method[J]. Journal of Computational Physics， 2004， 194（1）： 363-393. DOI： 10.1016/j.jcp.2003.09.032.
[14] ANDREASSEN E， CLAUSEN A， SCHEVENELS M， et al. Efficient topology optimization in MATLAB using 88 lines of code[J]. Structural and Multidisciplinary Optimization， 2011， 43（1）： 1-16. DOI： 10.1007/s00158-010-0594-7.
（编辑 武晓英）