赖松 王坤
摘要：为实现更加先进的拓扑优化算法，研究采用反应扩散方程的水平集结构拓扑优化方法，通过理论推导给出算法中的参数选择建议。该方法允许在拓扑优化过程中生成新的孔洞，初始结构无须包含孔洞，不需要重新初始化步骤，从而可提高算法的收敛性。针对传统拓扑优化中主要采用体积约束、以柔度最小为目标和体积保留率设定存在一定主观性的问题，探究不同体积保留率下的结构应力水平的变化规律，结果显示可以依据结构最大应力水平与体积保留率的变化规律确定最优体积保留率。
关键词：
反应扩散方程; 水平集; 拓扑优化; 体积保留率; 收敛
中图分类号：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。
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.
（编辑 武晓英）
猜你喜欢 约束密度体积 排开液体体积知多少中学生数理化·八年级物理人教版(2020年5期)2020-10-29谈拟柱体的体积高中生学习·高二版(2017年9期)2017-10-25马和骑师小学阅读指南·低年级版(2017年1期)2017-03-13谁的体积大幼儿智力世界(2016年11期)2017-02-21“密度”练习中学生数理化·八年级物理人教版(2015年12期)2016-01-25密度的应用趣谈中学生数理化·八年级物理人教版(2015年12期)2016-01-25密度的不变性与可变性中学生数理化·八年级物理人教版(2015年12期)2016-01-25CAE软件操作小百科（11）计算机辅助工程(2012年5期)2012-11-21化学问答中学生数理化·高二版(2008年10期)2008-06-17人类性行为要受到约束吗祝您健康(1993年3期)1993-12-30
