基于可行域调整的多相材料结构拓扑优化设计

俞燎宏，荣见华，唐承铁，李方义

1.长沙理工大学 汽车与机械工程学院，长沙 410114 2.宜春学院 物理科学与工程技术学院，宜春 336000 3.长沙理工大学 工程车辆安全性设计与可靠性技术湖南省重点实验室，长沙 410114

结构和材料优化设计已从传统的单一材料扩展到了多相材料的结构拓扑优化设计。Thomsen[1]较早采用均匀化方法研究了两相材料拓扑优化问题。Yin和Ananthasuresh[2]将峰值函数引入到SIMP (Solid Isotropic Material with Penalty)方法，并开展了柔性机构多相材料拓扑优化设计。孙士平和张卫红[3-4]提出了多相材料微结构的多目标优化设计模型，并针对多相材料结构拓扑优化中出现的棋盘格现象，提出了新的周长约束控制方法。Gao和Zhang[5]提出了以质量为约束的多相材料优化模型，并对多相材料的插值模型进行了比较研究。Stegmann和 Lund[6]提出了一种用于复合材料铺层优化的离散材料优化(Discrete Material Optimization，DMO)方法。在该方法中，将离散的纤维方向特征化为一种材料，DMO插值方案也是SIMP的扩展。而后，Hvejsel和Lund 通过统一参数化将SIMP和RAMP (Rational Approximation of Material Properties)推广到多相材料结构优化[7]。Blasques等[8-9]利用改进的SIMP多相材料插值模型，开展了复合材料层合梁的多相材料拓扑优化设计。由于在该类基于SIMP的多相材料插值的模型[6-9]中，刚度矩阵明显存在同一单元的不同设计变量的乘积项，故采用可分离和凸展开的移动渐近线方法(MMA或GCMMA)很难求解含2个以上约束的多相材料拓扑优化问题。最近，DMO拓展到了层压复合材料的稳健屈曲优化和并发的多尺度优化，并形成了较完善的材料插值模型和求解框架[10-12]。Stegmann和 Lund[6]指出DMO模型的一个困难是优化问题的总体设计空间为非凸空间，并采用仅有2个单元的简单复合材料单层板设计例子展示了优化问题的多解性、优化解与初始解的关联性。对于拥有大量变量的多相材料结构拓扑优化问题，采用现有基于SIMP的方法进行多相材料结构拓扑优化求解难于获得全局优化解[6,12-13]。

Tavakoli和Mohseni[14]提出了交替主动相算法求解多相材料优化问题，先将多相材料优化问题转化为序列二元相材料优化子问题，然后采用传统优化准则法进行求解。之后，Tavakoli进一步改进该方法，并求解了弹性动载荷下的复合材料设计问题[15]。Lieu和Lee[16]基于序列二元相拓扑优化子问题求解思路和等几何分析措施，研究了多相材料多分辨率的结构柔顺度拓扑优化设计问题。Park和Sutradhar[17]研究了三维多相材料多分辨率的结构柔顺度拓扑优化问题。在上述方法的序列二元相材料优化子问题模型中，仅考虑了一相材料体积约束，且没有涉及到对不同初始优化拓扑的多样性设计研究和比较。Zuo和Saitou[18]提出了序列多相材料SIMP插值方法，考虑质量和成本为约束，开展了结构柔顺度最小化的多相材料拓扑优化设计。杜义贤等[19]将该模型推广到了多相材料柔性机构拓扑优化设计问题，然而该模型使得结构性能函数在多相材料惩罚模型曲线的材料变相点处的导数不存在。贾娇等[20]提出了一种基于多相材料的一体化模型，研究了稳态热传导优化设计问题。张宪民等[21]提出一种并行策略，研究了多相材料柔顺机构多目标拓扑优化设计问题。龙凯等[22]在结构轻量化设计要求下，提出了基于多相材料的动态拓扑优化方法。另外, 还有水平集方法[23-24]、渐进结构优化方法[25]和相场法[26-27]也用于多相材料的拓扑优化设计。

多相材料连续体拓扑优化是典型的非线性问题，存在大量的局部解[6,28]。从而无法确保基于凸规划的梯度方法优化获得其全局最优解。依据工程实际需求，王博等[28]提出了一种连续体拓扑优化多样性设计方法。但不同初始拓扑导致优化解差异性的研究文献较少。如何获得多相材料拓扑优化设计问题的多个不同的局部解，并最终得到较好的优化解，值得深入研究。本文基于近年来荣见华等提出的可行域调整方法[29-31]，构建与原优化问题等效的多个含2个主动材料相的拓扑优化子问题，改进交替主动相算法，并采用光滑化对偶算法求解。最后，通过算例验证本文方法的可行性和稳健性, 并说明本文方法解决多相材料结构拓扑优化问题的工程价值。

1 多相材料的结构拓扑优化模型

1.1 多相材料插值模型

参考文献[29-31]，采用RAMP模型，构建单元体积和单元刚度矩阵的表达式为

(1a)

(1b)

(2)

(3)

(4)

1.2 多相材料拓扑优化模型及其近似模型

考虑多工况作用，则第l组工况载荷作用下的结构柔顺度Cl(ρ)表达式为

l=1,2,…,L

(5)

式中：K(ρ)为结构整体刚度矩阵；Fl为第l组工况载荷；lu为第l组工况载荷作用下所对应的结构位移向量；L为工况数。

借鉴文献[32]，采用η范数凝聚函数将多个工况的结构柔顺度优化问题转化为一个单目标的优化问题。以结构柔顺度最小化为目标、实体材料体积为约束的多相材料拓扑优化模型表示为

(6)

为了增强优化算法的寻优能力和稳健性，参考文献[31]，引进变体积限可行域调整技术，将优化模型式(6)转化为近似优化模型，即

(7)

(8a)

(8b)

式中：ρs,(k)为第k外循环迭代步求解获得的结构中第s相材料的拓扑变量向量；βs为经验参数，其取值范围为[0.002, 0.02]。

(9)

(10)

式(9)中，结构柔顺度的2阶导数采用MMA近似式[33-34]近似获得。为了处理数值不稳定性问题，采用文献[25,30]中的灵敏度过滤方法对上述目标函数的1阶、2阶导数进行过滤。

2 改进的交替主动相算法及求解优化

2.1 改进的交替主动相算法

(11a)

(11b)

(12)

将优化模型式(7)的求解转化为式(13)二元相拓扑优化子模型的求解，

(13)

如果a和b两相材料均为实体材料。根据式(9)，结构柔顺度对主动设计变量的1阶导数和2阶导数可分别表示为

(14)

(15)

根据式(10)，材料体积对主动设计变量的1阶导数和2阶导数可分别表示为

(16)

(18)

(19)

如果第a相材料为实体材料，第b相材料为空洞材料。根据式(9)，结构柔顺度对主动设计变量的1阶导数和2阶导数分别表示为

(20)

(21)

材料体积对主动设计变量的1阶导数和2阶导数，可由式(16)和式(17)得到。

因此，二元相拓扑优化子模型式(13)中近似目标函数和材料体积可分别表示为

(22)

(23)

如果第b相材料为空洞相材料，则只有与第a相材料相关的体积约束和其近似函数。

2.2 拓扑优化子模型的求解

参考文献[29-31]，将人工变量w=[w1w2…wmλ]T引入近似优化子模型式(13)。在y0邻域范围内可形成与优化子模型式(13)近似等效的m(m-1)/2个二次规划优化模型，即

(24)

式中：r+1表示第k+1外循环迭代步两相材料优化设计迭代的组号；第k+1外循环迭代步的前(m-1)(m-2)/2组两相材料优化设计时，mλ=2；第k+1外循环迭代步的最后的m-1组两相材料优化设计时，mλ=1。其中，

(25)

s=a,b

(26)

为了减少计算量，采用光滑化对偶方法求解优化模型式(24)，其求解详细步骤见文献[29]。

当满足式(27a)和式(27b)时，优化迭代收敛。

(27a)

(27b)

式中：εs和ε为收敛参数，取εs=ε=0.001。

2.3 交替主动相优化设计流程

假设结构由m相材料(含空洞材料)组成，按材料弹性模量值从大到小的顺序，依次用数字1, 2, …,m表示材料相。在交替主动相求解过程中，首先，选取第1相和第2相材料，记a=1,b=2，形成二元相拓扑优化子模型，采用光滑化对偶算法求解该子模型，更新设计变量。然后，判断b=m是否成立，如果不成立，则设置b=b+1，形成新的二元相拓扑优化子模型，求解并更新设计变量；如果b=m成立，则判断a=m-1是否成立，如果不成立，则设置a=a+1及b=a+1，形成另一个新的二元相拓扑优化子模型，求解并更新设计变量；以此类推。最后，如果a=m-1成立，则判断优化的收敛性。结构中所有材料按照上述方式两两组合，依次求解优化，最后达到多相材料的合理布局。交替主动相优化设计流程图如图1所示。

采用文献[14]中的方法进行优化时，也需按材料弹性模量值从大到小的顺序，依次用数字1, 2, …,m表示不同材料相。如果m=3，第1相材料必须设为弹性模量最大的材料，第2相材料设为弹性模量次之的材料，第3相材料设为弹性模量最小的材料。如果调换上述材料相顺序，则文献[14]的程序将无法运行。采用本文方法进行优化时，按照上述主动相方案和材料相顺序进行优化设计，获得的优化结果更佳。如果调换上述材料相顺序，获得的优化结果相对较差。

采用改进的交替主动相方案，将2个材料体积约束与式(8)变体积约束限配合使用，使得该方法适用于不同的初始拓扑变量，增强了该方法对多相材料初始拓扑选择的灵活性，而且优化过程更为稳健。与单相材料的拓扑优化方法以及没有两相主动相材料子模型分解的多相材料优化问题的求解方法相比，该方法只需要增加额外的内循环迭代步。

3 优化算例

3.1 单工况多相材料结构拓扑优化设计

采用文献[14]中的梁式结构设计算例，用来验证本文方法的可行性和有效性，并研究本文方法获得多个局部优化解及寻找较好的优化解的能力。图2为梁式结构设计域，采用无量纲表示，长度L=96，高度H=48，厚度为1。在结构底边的点A、点B和点C铅垂方向分别作用载荷F、2F和F，其中F=1。将图2所示的结构设计域划分为96×48个等尺寸四节点平面应力单元。

本算例设置2个优化问题，在优化问题1中，假设结构由两相实体材料和一相空洞材料组成。第1相和第2相实体材料的弹性模量分别为2.0和1.0，第3相材料(空洞材料)的弹性模量设为1.0×10-9。所有材料的泊松比为0.3。第1相和第2相实体材料目标体积比分别预定为35%和25%。

对于优化问题2，假设结构由三相实体材料和一相空洞材料组成。第1相、第2相和第3相实体材料的弹性模量分别为9.0、3.0和1.0，第4相材料(空洞材料)的弹性模量设为1.0×10-9。所有材料的泊松比为0.3。第1相、第2相和第3相实体材料目标体积比分别预定为20%、10%和10%。

文献[28]提出了一种连续体拓扑优化多样性设计方法，在该方法中，将目标函数定义为多个设计构型的目标性能加权之和，加入多样性度量的约束条件，并形成了基本的优化列式和求解方法。考虑到多相材料连续体拓扑优化问题的多解性[6,28]，本算例设置了不同的初始拓扑变量，构建多个不同的优化初始拓扑，研究不同初始拓扑对优化结果的影响。采用本文和文献[14]的方法求解本算例时，过滤长度尺寸rmin=3.0Δmin(Δmin为结构所有单元边长的最小尺寸)。

由表1可见，基于8个不同的初始多相材料结构拓扑，采用文献[14]方法程序获得的8个局部最优多相材料结构构型相似度很大，仅能找出3种有细微差别的局部最优多相材料结构拓扑(见表1中(a)、(b)和(d))，且其结构柔顺度大小相差不大。即文献[14]方法程序分别从8个不同的初始多相材料结构拓扑开始优化设计，最后仅仅得到3种有细微差异的局部最优多相材料结构构型。究其原因，采用基于准则法的传统拓扑优化求解过程中，设计变量的变化范围仅由移动限参数move控制[35]。当移动限参数move取小值，比如0.01, 传统拓扑优化方法[35]可以确保两外循环迭代步的拓扑变量矢量在一小的邻域内变化，从而可确保结构柔顺度的1阶近似式有较好的精度，即传统准则迭代式所用的结构柔顺度导数有较好的精度[35]。然而，当移动限参数move取小值时，基于传统准则迭代式获得的优化拓扑往往存在不期望的少量灰度小构件，即该构件的伪密度变量处于0和1之间，很难从最后拓扑中提取满足约束要求的清晰优化拓扑[35-36]。故Sigmund等建议移动限参数move取较大的值[35-36]，比如0.05或者更大的值。由于一般初始拓扑远离最优拓扑，当移动限参数move取较大的值时，传统准则迭代式使得初始迭代的两外循环迭代步的拓扑变量矢量常不处于一小的邻域内，其结构柔顺度的1阶近似式误差很大，即传统准则迭代式所用的结构柔顺度导数误差也很大。因此，初始的一些迭代步获得的解在一定程度上丧失了优化特性。实质上，文献[14]方法继承了该特点。尽管文献[14]方法的程序是一个教育版程序，但它基本反映了文献[14]方法的特征。

图3给出了与表1中(a)～(e)相对应的起始几个外循环迭代步的结构柔顺度变化曲线图。图4给出了采用相应于表1的不同初始多相材料拓扑和文献[14]方法求解优化问题1获得的第8外循环迭代步的多相材料结构拓扑及其柔顺度值。

表1采用本文方法和文献[14]方法求解优化问题1获得的局部最优结构拓扑及其柔顺度值

Table1LocaloptimalstructuraltopologiesandtheircompliancesobtainedbyproposedmethodandmethodinRef.[14]tosolveoptimizationproblem1

序号实体材料的初始拓扑变量本文方法文献[14]方法局部最优结构柔顺度局部最优结构柔顺度(a)ρ1,(0)i=0.35ρ2,(0)i=0.2541.899 057.585 2(b)ρ1,(0)i=0.40ρ2,(0)i=0.3039.629 957.362 1(c)ρ1,(0)i=0.99ρ2,(0)i=0.0139.900 557.620 9(d)ρ1,(0)i=0.80ρ2,(0)i=0.2039.521 257.638 2(e)ρ1,(0)i=0.60ρ2,(0)i=0.4039.551 157.6412(f)ρ1,(0)i=0.50ρ2,(0)i=0.5039.672 357.630 8(g)ρ1,(0)i=0.30ρ2,(0)i=0.7057.967 257.611 7(h)ρ1,(0)i=0.10ρ2,(0)i=0.9056.870 257.630 1

由图3和图4可见，初始几个迭代步的结构柔顺度变化很大，多相材料结构拓扑也变化很大，这些都进一步表明文献[14]方法具有上述特点。从而，多相材料结构拓扑优化特性的部分丧失导致迭代后期的优化解与初始多相材料拓扑的关联特性的部分丧失。从图4可知，仅需8个迭代步，原来5个差异较大的初始多相材料拓扑变为5个相似度很大的多相材料结构拓扑。既使后继的优化迭代比较稳健, 最终获得的5个局部最优多相材料结构构型的相似度也大。因为多相材料连续体结构拓扑优化问题存在大量的局部解，以及无法确保基于凸规划的梯度方法优化获得其全局最优解。因此，当寻求多相材料连续体拓扑优化问题的较好局部解时，上述特点成为文献[14]方法的弱点。

表2给出了结构上下部分的初始拓扑变量取值不同时，采用本文方法求解优化问题1获得的局部最优结构拓扑及其柔顺度值。从表1和表2可见，从不同的初始多相材料拓扑开始优化设计，采用本文方法获得的局部最优多相材料结构构型基本不同。

表2结构上下部分的初始拓扑变量取值不同时采用本文方法求解优化问题1获得的局部最优结构拓扑及其柔顺度值

Table2Localoptimalstructuraltopologiesandtheircompliancesobtainedbyproposedmethodtosolveoptimizationproblem1whentheupperandlowerpartsofstructurehavedifferentinitialtopologicalvariables

梁式结构设计域实体材料的初始拓扑变量局部最优结构柔顺度(a)上半部分初始拓扑变量ρ1,(0)i=0.60, ρ2,(0)i=0.40(b)下半部分初始拓扑变量ρ1,(0)i=0.30, ρ2,(0)i=0.2039.857 9(a)上半部分初始拓扑变量ρ1,(0)i=0.30, ρ2,(0)i=0.20(b)下半部分初始拓扑变量ρ1,(0)i=0.60, ρ2,(0)i=0.4039.893 9

表3采用本文方法和文献[14]方法求解优化问题2获得的局部最优结构拓扑及其柔顺度值

Table3LocaloptimalstructuraltopologiesandtheircompliancesobtainedbyproposedmethodandmethodinRef.[14]tosolveoptimizationproblem2

序号实体材料的初始拓扑变量本文方法文献[14]方法局部最优结构柔顺度局部最优结构柔顺度(a)ρ1,(0)i=0.50ρ2,(0)i=0.30ρ3,(0)i=0.2016.857 522.375 4(b)ρ1,(0)i=0.20ρ2,(0)i=0.30ρ3,(0)i=0.5020.884 622.398 8(c)ρ1,(0)i=0.40ρ2,(0)i=0.30ρ3,(0)i=0.3017.768 822.379 2(d)ρ1,(0)i=0.30ρ2,(0)i=0.30ρ3,(0)i=0.4019.444 322.635 8

因此，无论是三相或者四相材料结构拓扑优化，本文方法获得多个局部优化解及寻找较好的优化解的能力比文献[14]方法强得多。最近，Sanders等[12]提出了一种新的多相材料连续体结构拓扑优化方法, 在该方法中独立地更新与每一相材料的体积约束相关的拉格朗日乘子，并采用ZPR(Zhang-Paulino-Ramos)设计变量更新方案求解优化问题，实质上该方法是传统的基于准则的迭代法[35]在多相材料连续体结构拓扑优化中的推广。在本文方法中，结构柔顺度灵敏度计算量是主要的工作量，与现有方法的计算量相同[12]，多个子模型优化求解所需工作量不多，故对三相或者四相材料结构拓扑优化问题来说，本文方法的计算效率与Sanders等[12]方法的计算效率是相近的。对于更多相材料结构拓扑优化问题，本文方法的计算效率略有欠缺。然而，在现有的一些多相材料连续体结构拓扑优化方法中，如果采用自动满足与式(4)类似要求的多相材料插值方案[9], 则获得没有实体材料重叠的优化结构拓扑需很多的外循环迭代步；如果采用不能自动满足式(4)要求的多相材料插值方案[12]，则需在迭代算法中引入一些启发式处理措施以确保与式(4)类似的要求和实体材料不重叠的要求[12]，从而降低了获得的优化解的优化特性。这些情况和上述算例结果表明，本文方法具有更重要的工程应用价值。

3.2 多工况载荷下多相材料结构拓扑优化设计

图9为矩形结构设计域，长L=3.0 m，高H=1.0 m，厚度为0.01 m。假设第1相实体材料的弹性模量为200 GPa，第2相实体材料的弹性模量为100 GPa，两相实体材料的泊松比均为0.3。将初始设计域划分为120×40个等尺寸的平面应力单元，并指定两相实体材料目标体积比均为20%。

两工况载荷下结构优化问题(见图9(a))：矩形结构左右两侧固定，受2个工况载荷作用。工况1：集中载荷F1=10 kN沿铅垂方向作用在结构下边中间点A1；工况2：集中载荷F2=10 kN沿铅垂方向作用在结构上边中间点A2。

三工况载荷下结构优化问题(见图9(b))：矩形结构左右两侧中点固定，受3个工况载荷作用。工况1：集中载荷F1=10 kN沿铅垂方向作用在结构下边中间点B1；工况2：集中载荷F2=10 kN沿铅垂方向作用在结构上边左端点B2；工况3：集中载荷F3=10 kN沿铅垂方向作用在结构上边右端点B3。

由图10和图13可见，采用本文方法也能较好地解决多工况多相材料结构拓扑优化问题，可以找到结构柔顺度更小的优化解，再次验证了本文方法的可行性和有效性。

4 结 论

针对多相材料结构柔顺度拓扑优化问题，引进可行域调整技术，改进交替主动相算法，提出了一种新的寻优能力强的拓扑优化设计方法。通过理论推导和数值算例验证，得出以下结论：

1) 多相材料连续体拓扑优化是典型的非线性问题，存在大量的局部解。从而，无法确保基于凸规划的梯度方法优化获得其全局最优解。

2) 设置不同的优化初始拓扑，采用本文方法能够获得不同的优化结构，从中可以找到一个更优的多相材料结构拓扑。与文献[14]方法相比，本文方法具有强的寻找较好优化解的能力。

3) 引入变体积约束限调整技术和改进的交替主动相算法，使得采用本文方法获得的结构拓扑具有较好的互不重叠的实体材料分布特征。