基于相场描述的拉压不对称强度结构拓扑优化

张晓鹏， 康 柯， 杨东生

(1.大连理工大学 工程力学系，工业装备结构分析国家重点实验室，大连 116024；2.中国运载火箭技术研究院，北京 100076)

1 引 言

拓扑优化设计作为一种有效的概念设计工具，目前已应用于航空航天、高速列车以及船舶与海洋工程结构等重大装备的研发和设计中以及近年来，结构拓扑优化技术已获得长足发展，并产生了均匀化方法、各向同性实体材料惩罚法SIMP(solid isotropic material with penalization)、水平集法、渐进结构优化方法ESO(evolutionary structural optimization)、相场法和独立连续映射法ICM(independent continuous mapping method)等一系列有效的拓扑优化方法[1-3]。

在实际结构设计中，结构部件都需要满足一定的强度设计要求，即结构的应力水平需要满足指定的约束要求，在此基础之上建立应力约束下的结构拓扑优化方法。应力约束下的拓扑优化研究的主要难点之一是应力奇异性问题，随着相关研究的开展，目前常用的处理方法有ε放松法[4]和qp放松法[5]等。考虑应力约束的结构优化的另一个难点是约束过多，与结构柔顺性和体积全局约束优化问题相比，应力约束属于局部约束，对于大规模优化问题来说，面临着约束数目多及灵敏度分析和优化求解复杂等问题。近年来，在实际处理中，可以用聚合函数将众多的局部应力约束聚合为一个应力约束，相应的聚合函数包括K-S函数[6]和 P -norm 函数[7,8]等。为了处理应力约束问题，人们不断地改进相关技术。隋允康等[9]把局部应力约束问题转化为结构整体的应变能约束问题，减少了约束的数量。荣见华等[10]建立了一种改进的基于应力及灵敏度的结构拓扑双方向渐进优化算法，该方法能够较大程度减少解的振荡次数。刘宏亮等[11]引入了稳定转换法来修正P -norm 聚合后的应力约束，克服了拓扑优化过程中的迭代振荡和收敛困难问题。王选等[12]对双向渐进拓扑优化方法进行改进，实现了有效降低应力集中效应的结构优化设计。杨迪雄等[13]引入了稳定转化法STM(stability transformation method)和应力约束的违反约束集列式，实现了应力约束拓扑优化问题的稳定迭代和可接受的设计解。

大部分考虑材料失效的拓扑优化研究都是基于Mises应力屈服准则的，能够很好地预测金属材料的破坏行为。然而，有很多的非金属材料在受到拉伸和压缩时会表现出不同的应力极限，在这种情况下，von Mises应力准则是不合适的。以应力不变量定义的Drucker-Prager准则作为最简单的塑性屈服模型，可以反映材料不同的许用应力极限。此前，仅有少量的研究考虑了拉压不同强度约束的拓扑优化方法。如Luo等[14]首先提出了基于变密度法的考虑结构拉压强度不同的拓扑优化设计方法；Jeong等[15]等研究了基于多种不同失效准则下韧性和脆性结构应力拓扑优化方法;Bruggi等[16]在应力约束的拓扑优化中引入了柔顺性，探讨了不同的目标函数和约束对优化结果的影响。

相场(Phase field)方法能够很好地描述固体-液体过渡等复杂材料的扩散界面特性[17]，并已经成功用于结构拓扑优化设计中[18]。与变密度法相比，基于相场法可以通过在实体域和孔洞域之间使用显式曲线来减少中间设计变量，使后处理更加直观，减少后处理对优化结果的影响，在求解应力约束下的优化问题中具有一定的优势[19]。然而，在少量基于相场描述的结构应力拓扑优化研究中，尚未见基于相场模型的考虑拉压强度不同约束下的结构拓扑优化的相关研究。

因此，本文通过引入Drucker-Prager准则来描述塑性材料的屈服行为，提出了具有该屈服行为的连续体结构的拓扑优化设计。优化问题可以描述为，基于相场法和在Drucker-Prager应力约束下，找到设计域内体积最少的材料布局。利用qp放松法和 P -norm 凝聚函数使拓扑问题具有数值可解性，引入稳定转化法来降低迭代过程中的振荡问题，加速收敛。利用伴随变量法进行灵敏度分析，通过求解Allen-Cahn方程更新相场函数设计变量，实现拉压不同强度下的结构拓扑优化，并通过典型拓扑优化算例验证了本文方法的有效性。

2 基于应力约束的连续体结构的拓扑优化

2.1 基于相场描述的拓扑优化方法

设计域的状态可以通过最小化范德华自由能方程来确定：

图1 由相场模型确定的扩散界面

(1)

(2)

(3)

(4)

(5)

在基于相场描述的拓扑优化方法具体实现中，为了使拓扑构型不断演变，需要对两相材料的过渡区域材料性能(如弹性模量)进行插值，与SIMP法处理方式类似，在过渡区域相场点通常引入插值函数p对中间设计变量进行幂指数插值，其表达式为

图2 双势阱函数图解

2.2 Drucker-Prager准则

在材料拉压强度不同的研究中， Drucker-Prager准则因其简单和实用性广泛应用于非金属材料设计中。通过对许用拉压应力不同材料的实验，可以验证Drucker-Prager数值预测结果的准确性。Drucker-Prager准则的屈服方程可表示为

(7)

式中α和H为材料常数，I1为应力张量的第一不变量，J2为应力偏量的第二不变量，其函数表达式为

I1(σ) =σx x+σy y+σz z

(8)

式中σ= [σx xσy yσz zσx yσy zσz x]为质点应力。

对于Drucker-Prager准则而言，F(σ) = 0定义的屈服面是主应力空间中的一个圆锥体(图3(a))。根据这一屈服准则，当F(σ) ≤ 0时，材料不会发生失效。当α=0时，该屈服准则简化为von Mises失效准则，屈服面变为主空间中的一个圆柱体(图3(b))。因此，von Mises失效准则可以看做Drucker-Prager准则的一个特例,Drucker-Prager准则通过增加一个修改项来施加静水压力分量对破坏效果的影响。材料常数α和H可以表示为

(9)

式中σc和σt分别为材料的许用压缩强度和许用拉伸强度。

图3 主应力空间中Drucker-Prager准则屈服面和Mises准则屈服面

2.3 拓扑优化列式

基于相场描述的结构拓扑优化模型，考虑满足Drucker-Prager强度准则下最小化结构的体积优化问题。其优化列式为

s. t.KU=F

式中V为设计域的总体积，K，U和F分别为总体刚度矩阵、结构位移向量和载荷向量，KU=F为线性平衡方程。

基于相场描述的优化方法，选取观察相场点的相场函数值为设计变量(数值实现中相场点位置与单元中心位置重合)。考虑应力的拓扑优化问题，常面临当设计变量相对密度趋于0时，单元的局部应力并非为0的应力奇异性问题。为克服这一问题，本文采用了qp放松法[5](在低密度区域放宽对局部应力的约束)，则式(10)应力约束公式改写为

(11)

式中取p=3和q=2.5。当设计变量为1时，式(11)与式(10)是等效的;当设计变量趋于0时，式(11)得到的是放松后的单元应力值，可以保证局部应力的连续性。

2.4 应力约束的处理方法

为处理应力约束拓扑优化中的大量局部约束问题，采用 P -norm 函数对局部应力约束进行凝聚，以减少约束的数量，其表达式为

(12)

式中Fe为每个单元的应力约束,表示对局部应力pn约束的聚合参数。当参数pn趋近于无穷大时，聚合函数的值等于局部约束的最大值，即

(e= 1,2,…,N)(13)

将所有局部约束转换为一个全局约束可以显著减少计算工作量。当pn值相对较小时，全局约束函数不能充分反映局部约束下的可行域;但当pn值过大时，全局应力约束函数呈现高度非线性，将导致最终结果的不稳定性和不收敛。因此，在具体数值算例中，本文选择pn=6。为克服在连续体拓扑优化中局部约束的数量过多导致的聚合函数精度不够的问题，本文采用稳定转换法对全局应力聚合函数进行修正[13]。修正后应力能聚约束可以用P -norm 函数重新表示为

(14)

式中ce为修正参数。对于第n次优化迭代，其修正参数可以计算为

(15)

值得注意的是，应力修正方案是基于稳态转化法得来的，特别适用于克服应力约束拓扑优化中由于应力的高度非线性所造成的迭代振荡问题[13]。

2.5 灵敏度分析

根据链式法则，式(14)对设计变量的灵敏度可以表示为

(16)

式中P -norm 聚合函数对Drucker-Prager应力的灵敏度表达式为

(17)

式(8)中表示Drucker-Prager准则的应力张量第一不变量和应力偏量的第二不变量可以用矩阵表示为

(18)

(e= 1,2,…,N)(19)

(20)

本文采用qp放松而将所有应力聚合为一个全局应力约束[13]，因此伴随变量λ可得

(21)

(22)

3 数值算例

以应力约束拓扑优化为例，验证所提方法的有效性。在所有的算例中，材料的杨氏模量和泊松比分别为E=1 MPa和ν=0.3。为了便于优化迭代的振荡控制，在进行基于稳定转化法的应力校正时，设定qn=0.5，ce0= 1。在优化过程中，当目标函数在最近的50次迭代中相对波动小于1%或迭代次数超过300次时，迭代终止。

3.1 底端固定梁拓扑优化

算例针对具有不同拉压强度材料组成结构进行分析和优化。在拉压强度相同时，即σc=σt=5.8 MPa 时，其拓扑优化结果如图5所示。图6 显示最大应力以及材料体积比在拓扑优化中的迭代过程，优化过程中目标函数能够稳定减小并趋于稳定，最后优化在迭代228次时收敛。图6的迭代历史给出了演化的拓扑构型，分别是第20，60，100，140，180和220迭代步的构型。由拓扑构型可以看出，该优化问题收敛较快，在第20步时，就能够得出大致的拓扑构型，在第130步左右迭代趋于稳定。

图4 底端固定梁的设计域和矩形设计域的初始布局

图5 拉压强度相同时的拓扑结果

当拉压强度不同时，分别考虑σc= 1.5σt=8.7 MPa和σc= 2σt=11.6 MPa两种情况，其拓扑优化结果如图7所示。优化结果显示，随着许用拉压应力比的增加，受拉杆与水平线的夹角增加，受压杆与水平线的夹角减小。同时，由于许用压应力变大，优化后结构整体的体积减小，受压杆变细。文献[14]通过理论推导得出了在不同拉压比下的最优构型，与本文结果基本保持一致。该结果也说明了在许用拉压强度不同时，考虑Drucker-Prager准则的结构拓扑优化具有必要性。表1给出不同拉压强度下，拓扑优化收敛时的材料体积、迭代次数和最大应力约束违反率，可以看出，不同拉压强度下，基于相场描述的拓扑优化方法都能稳定收敛并达到较高的应力约束精度。

图6 最大应力和材料体积的迭代历史

图7 不同拉压比下拓扑优化结果

表1 不同拉压应力比下拓扑优化的体积比、迭代次数和最大应力约束违反率

3.2 L型梁拓扑优化算例

值得注意的是，图9(a)的优化结果呈现一个直角，这一设计虽然可能导致应力集中，但对于拉压不同强度来说，因材料更抗压所以在直角区域仍能满足应力约束(图10(a))，鲜有的考虑拉压不同强度的应力优化研究[14]也印证了这一结果。采用拉压不同强度材料时，两种优化设计的Drucker-Prager屈服函数云图如图10所示，可以看出,考虑拉压不同强度的优化设计能够完全满足Drucker-Prager屈服约束约束(最大值为0)，如图10(a)所示。然而将考虑von Mises应力约束下的优化结果采用式(7)进行D -P准则校核时，其应力最大处的F(σ)max=5.67 MPa,违反了Drucker-Prager屈服约束，如图10(b)所示，说明利用拉压强度相同的拓扑优化结果并不能满足拉压非对称强度约束。

图8 L型梁设计域

图9 L型梁的拓扑优化结果

图10 L型梁的Drucker-Prager屈服函数云图

4 结 语

本文基于相场描述的拓扑优化方法，通过使用P -norm 函数聚合应力约束，大量减少了约束数目，提高了求解效率，同时引入稳定转换法来减少聚合函数导致的聚合误差，实现了应力较为精确的约束。同时，引入Drucker-Prager准则，考虑了非金属材料的优化问题，数值算例表明，与Mises应力准则相比，在不同拉压比下，结构的拓扑构型有很大的不同。因此，在对非金属材料进行优化设计时，考虑Drucker-Prager准则是很有必要的。