考虑可制造性约束的声子晶体多目标拓扑优化1)

曹蕾蕾 武建华 樊 浩 张传增 孙霖霖

* (长安大学道路施工技术与装备教育部重点实验室,西安 710064)

† (德国锡根大学土木工程系,德国锡根D-57068)

引言

声子晶体是一种新型周期性人工超材料,它可以使得某些频率范围内弹性波的传播衰减或禁止,其特有的带隙特性使其在负折射、声聚焦、声吸收、声隐身、声成像领域具有广阔的应用前景[1-2].带隙的产生和分布是声子晶体应用的关键问题,因此,利用拓扑优化方法进行声子晶体带隙设计与优化引起了国内外学者的广泛关注.

关于声子晶体的拓扑优化方法主要有两类:梯度类算法和非梯度类算法.梯度类算法通过处理目标函数和约束条件的灵敏度信息,使可行解沿问题的梯度方向移动更新,可通过较少迭代次数完成问题的求解.目前,该方面的研究主要采用移动渐近线算法(method of moving asymptotes,MMA)[3-6]和双向进化结构优化(bi-directional evolutionary structural optimization,BESO)方法[7-8]对声子晶体进行带隙优化.虽然梯度类算法求解效率高,但存在容易陷入局部最优解的弊端.非梯度类算法是通过对适应度函数的不断评估和多次迭代,在设计域内找到满足要求的近似全局最优解.其中,遗传算法在声子晶体拓扑优化中应用最为广泛,如标准遗传算法(standard genetic algorithm,SGA)[9-15]、自适应遗传算法(adaptive genetic algorithm,AGA)[16-18]、多精英遗传算法(multiple elitist genetic algorithm,MEGA)[19]、快速非支配排序遗传算法(non-dominated sorting genetic algorithm II,NSGA-II)[20-24]与不同的数值计算方法相结合皆可实现对声子晶体的带隙优化.

虽然声子晶体的拓扑优化近年来已经取得了一系列的重大进展,但优化中存在的“棋盘格效应”使得结构中往往出现某种孤立材料单元、超细而导致刚度较弱的链接单元、微结构组元几何形貌十分复杂、微结构边界和各组元边界或界面粗糙等异常情况,这对于实际加工制造十分不利,因而大大限制了声子晶体的应用.针对优化结果常常出现孤立材料单元这一问题,本文开展考虑可制造性的二维多相声子晶体拓扑优化研究,通过对优化过程中拓扑构型的连通性分析以限制最小材料连通域的占比,避免出现孤立材料单元,进而建立提高其可制造性的优化模型.在此基础上,应用NSGA-II 对模型进行求解,得到了多目标帕累托(Pareto)解集.具有代表性的数值算例结果表明本文方法可在满足带隙性能等指标的前提下同时兼顾到制造可行性的要求,有利于声子晶体优化结果的实际加工制造和推广应用.

1 考虑可制造性的声子晶体拓扑优化方法

1.1 带隙计算理论模型

主要考虑二维声子晶体可视为非均质各向同性弹性材料的情况,忽略体力影响时,其弹性波动方程为

式中,ρ,λ 和 µ 分别是质量密度和拉梅常数,u是位移向量.∇ 和 ω 分别为梯度算子和角频率.假定弹性波仅在xy平面内传播,则有 ∂u/∂z=0,于是方程(1)可分解为以下两个解耦的二维波动方程

仅考虑由式(2)描述的面内波动情形而不考虑由式(3)描述的面外或反平面波的传播.根据周期性结构中的Bloch 定理,位移向量u(r) 可用如下公式表示

式中,uk(r) 为周期性函数,k=k(x,y) 为波矢.利用有限元方法求解,式(2)可转化为以下特征方程

式中,K和M分别为整体刚度矩阵和整体质量矩阵,U是节点位移向量.给定波矢求解式(5)中的特征频率可得到能带结构或频散曲线.为了方便起见,能带结构图的纵坐标为带隙归一化频率 Ω=ωa/(2πct),其中a为晶格常数,ct为基体的剪切波波速.

1.2 优化模型

实施声子晶体的拓扑优化,主要思路是采用有限元网格对单胞结构进行离散化,通过优化算法在网格中填充合适的材料,得到满足优化目标要求的拓扑构型.由于遗传算法采用二进制变量,特别适合建立两相材料与拓扑构型之间的关系[12,19,24],即采用0 和1 分别代表两种不同的材料,故声子晶体的拓扑优化多采用遗传算法.对于多相声子晶体,可通过引入“二进制向量”的方式对离散单元的材料进行编码,即

式中,φi(i=1,2,···,NE) 为任意离散单元的二进制编码(NE为单元数),bs为二进制编码值,t为向量长度,NM为材料数.以4 相材料为例,有

对于更多相声子晶体,可通过增加向量 φi的长度t来进行相应的编码.由离散单元的“二进制向量”编码组合可得到整个设计域的材料分布，即

式中，ϕ 为设计域的材料分布.在优化问题求解中,采用N×N的网格对单胞进行离散,由于结构的平移周期性和点群对称性,通过波矢遍历不可约布里渊区边界求解特征频率可得到能带结构,则设计域内单元数目减少至 (N/2)(1+N/2)/2,对应遗传算法染色体长度为t(N/2)(1+N/2)/2 .

为了避免出现孤立材料单元而难以或无法实际制造的问题,需要引入合适的可制造性约束条件.本文试图通过控制连通单元的占比达到限制孤立单元的目的.引入图像处理中四连通区域的定义,即从区域内一点出发,可通过上下左右4 个方向的移动组合,在不越出区域的前提下,能达到区域内任一像素的区域定义为四连通区域.将离散后的网格看作像素点,则可将此定义扩展到材料连通单元中,进一步由连通单元组成材料连通域,如图1 中所示.红、绿、黑不同标记的网格分别由3 种不同材料填充,其中,红色、绿色区域分别为独立的连通域,而黑色区域虽在角点处相连接,但属于两个连通域.

图1 连通域示意图Fig.1 The connectivity domain

按照编码方式为元胞网格填充材料后,可进行连通域的识别和其相对于元胞面积的占比计算,即

式中Ajk(ϕ) 和Ijk(ϕ) 分别为第k种材料的第j个连通域的面积和其相对于元胞的面积占比;A为元胞面积,NM为材料数,NQ为每种材料所包含的四连通域数目.Ijk(ϕ) 过小会使得声子晶体的可制造性下降,故通过限制Ijk(ϕ) 值保障其可制造性,相应的附加约束条件具体如下

式中,I(ϕ) 为Ijk(ϕ) 的最小值,r*为阈值,可根据实际对可制造性的要求和数值试验确定,r*越大则表示对结构可制造性要求越高.

对于声子晶体的带隙优化研究,常以最大化相对带隙宽度、最大化绝对带隙宽度、特定带隙中心频率等为优化目标.本文主要以轻质结构的减振隔振为应用背景,考虑到实际应用中通常是针对特定场合对特定频段的减振需求设计具有预期带隙的声子晶体,同时,往往还需兼顾结构轻量化设计的需求,故本文以在特定频率段带隙最宽和结构质量最小为优化目标.

综上所述,考虑可制造性约束的二维多相声子晶体多目标优化优化模型表述如下

式中,ϕ 为设计域材料分布;Δ ωn为第n带隙的带隙宽度,为了方便描述,指定第n条能带和第n+1 能带之间存在的带隙为第n带隙,即当第n+1 条能带的最小值大于第n条能带的最大值时,带隙宽度Δωn＞0,反之带隙宽度 Δ ωn=0,带隙不存在;F1(ϕ)为特定频率范围段内带隙宽度之和的占比,其值越大表示带隙对特定频率段的覆盖率越高,ωlow和ωupp分别为特定频率范围的上、下界,F2(ϕ) 为优化结构质量,N为网格的离散规模,ρs为单元质量密度,V为单元面积.

1.3 优化算法

对于式(9)的多目标优化问题,可采用快速非支配排序遗传算法 N SGA-II[25]进行求解,该方法能够有效生成 P areto 解集,使决策者能够根据需要从Pareto解集中选择出最适合的解.算法流程如图2 所示.

图2 基于NSGA-II 算法的声子晶体多目标优化流程图Fig.2 Flow chart of the multi-objective optimization of phononic crystals based on the NSGA-II algorithm

在计算种群内个体的适应度值时,将种群中不满足可制造性约束和虽满足可制造性约束但无带隙个体的带隙特性目标函数值赋予极小值0.000 1,可以增加满足可制造约束且存在带隙个体在非支配排序时处于更优的非支配排序层的可能性,便于后续选择操作时对此类优秀个体的保留.

2 数值算例

以二维正方晶格声子晶体面内模态为例,晶格常数为0.01 m,声子晶体单胞由4 种材料组成,即金(绿)、环氧树脂(灰)、硅橡胶(黄)和硫化橡胶(蓝),材料名称、编码 φi、质量密度 ρ、杨氏模量E和泊松比 ν 如表1 所示.以在给定归一化频段0～0.05 带隙最宽为优化目标,有 ωlow=0,ωupp=0.05 ;考虑可制造性约束要求,取r*=0.005 ;对单胞进行网格划分,取N=16,则设计域内的单元数目和染色体长度分别为36 和72;遗传算法的参数选取为:种群规模20,交叉率0.9,变异率0.02,采用二元锦标赛选择方式,单点交叉方式,设定迭代次数1000 为终止条件.

表1 材料属性Table 1 Properties of the materials

种群规模Np取值对优化结果不会产生影响,但是会对优化收敛过程产生影响,即Np取值越大,所需得到Pareto 解集的迭代次数越少.对于本文考虑的问题而言,因在每一步迭代中要对种群中所有个体的带隙进行有限元计算,增大种群规模会导致单次迭代过程中有限元带隙计算次数大大增加,加大单次迭代的计算成本.因此,综合考虑以上因素并通过数值试验,本算例选取Np=20 .为了提高算法的收敛速度,在初始种群中引入具有带隙的“种子”结构,其余个体随机产生,“种子”单胞构型如图3 所示.

图3 “种子”结构Fig.3 The seed unit-cell

图4 为进化 1 000 代的Pareto 解集,横、纵坐标分别为目标函数F1和F2的值,可以看出两个目标之间的制约关系,即特定频段的带隙宽度之和的增大将使得结构质量也在同时变大.

从图4 中选择4 个非支配解A,B,C和D进行分析,其拓扑构型、能带结构和透射谱如图5 所示.由图5(a)可以发现,非支配解A,B,C和D的拓扑构型均由金、环氧树脂和硫化橡胶组成且具有相似的特征,即由金和环氧树脂组成的夹杂体被硫化橡胶包覆层包裹后嵌入环氧树脂基体中.图5(b)与图5(c)的对比表明能带结构的带隙范围与透射谱频率衰减范围吻合.由图5(b)的能带结构可知,解A,B,C和D的第3 带隙充分打开,其归一化带宽分别为0.017 1,0.022 0,0.025 2 和0.030 0,计算得到的对应目标函数值如表2 所示.

图4 附加可制造性约束的多目标Pareto 解集Fig.4 Pareto optimal solutions from the MOOP considering the manufacturing constraint

由表2 可知,结构A,B,C的带隙特定频段的覆盖率比D分别小42%,26% 和15%,结构质量比D分别小52%,44%和35%,因此,决策者可根据实际应用情形的目标权重从Pareto 解集中挑选自己所需要的声子晶体结构.

表2 附加可制造性约束优化结构目标函数值,S 为单目标优化结果Table 2 The objective function values of the optimized structures from the MOOP and the SOOP considering the manufacturing constraint

为进一步理解优化结果的带隙形成机理,对图5(b1～ b4)带隙的上、下边界所对应的振动模态进行分析.现以结构D为例进行说明(构型A,B,C的情况与之类似).构型D的第3 个带隙上、下边界处对应的振动模态如图6 所示.

图5 附加可制造性约束的多目标优化结果Fig.5 Pareto optimal solutions from the MOOP considering the manufacturing constraint

由图6 可以看出优化结果的带隙机理为局域共振型.其振动特性可近似的通过简化的等效“弹簧-质量”系统进行描述,相应的带隙上、下边界的频率可由下式进行估算[26-27]

图6 结构D 的带隙边界所对应的振动模态Fig.6 Vibration modes corresponding to the band-gap edges of the optimized structure D

其中,flow和fupp分别为带隙的下边界和上边界频率,K为包覆层(弹簧)的动态等效刚度,M1和M2分别为金散射体和环氧树脂基体的动态等效质量.

由图5 中的优化结果构型可以看出,优化构型中仅剩下3 种材料,材料初始分布中的硅橡胶在优化过程中被淘汰掉,包覆层全部由弹性模量更大的硫化橡胶构成,使得弹性包覆层动态等效刚度增加.由式(12)可知,当其他参数不变时,K的增大将导致带隙起始(下边界)和截止(上边界)频率均会出现不同程度的上升.由于截止(上边界)频率上升幅度更大,因而带隙范围变宽,带隙对特定频率段的覆盖率随之提高.

另外,通过对比图5(a),4 种构型发现,金散射体的占比不断提高,即M1逐步增大,而K与M2基本保持不变.由式(12)可知,M1的增大将导致带隙起始(下边界)和截止(上边界)频率降低,而起始(下边界)频率相比截止(上边界频率降幅更大,因而带隙范围变宽,带隙对特定频率段的覆盖率随之提高.即F1与F2同时增大,这种制约关系与图4 十分吻合.

为了说明多目标和单目标优化结果的差异,本文对只考虑了目标函数F1的单目标带隙优化问题也进行了求解,图7 给出了单目标进化曲线、近似最优解S的拓扑构型、能带结构和透射谱.

从图7(a)中可明显看出拓扑构型随迭代进程的演化过程,即随着进化代数的增加,包覆层的硅橡胶(黄)分布减小而硫化橡胶(蓝)分布增加,金散射体分布变化不大;进化800 代左右后种群的最佳适应度收敛到稳定值,表明种群进化趋于完善,对应的个体即为近似最优解S,其拓扑构型(图7(b))由硫化橡胶包覆着金嵌入环氧树脂基体中,能带结构(见图7(c))显示可在第3 和第4 能带间打开归一化带宽为0.030 2 的带隙,由此计算得到结构S的目标函数值同样在表2 中给出以作对比.对比D和S的结果发现,两者对特定频段的覆盖率十分接近,虽然S比D的特定频段覆盖率高2%,但其质量却比D大18%,所以从结构整体的角度来讲,结构D比S更优,这说明多目标优化不仅可以得到和单目标优化结果接近的非支配解,还可以得到其他的多目标优化非支配解集,因此较单目标优化更具有优越性.

图7 单目标优化结果Fig.7 The optimized solutions from the SOOP

此外,为说明考虑可制造性优化的有效性和优势,图8 给出了未附加可制造性约束的多目标优化Pareto 解集.从图8 内选择与考虑可制造性约束的4 个非支配解的目标函数值具有可比拟性的解A′,B′,C′和D′分析,其拓扑构型分别如图9 所示.

图8 未附加可制造性约束多目标Pareto 解集Fig.8 Pareto optimal solutions from the MOOP without the manufacturing constraint

从图9 可以看出,未考虑可制造性约束的结果构型中A′和B′存在多个金(绿)、硅橡胶(黄)和环氧树脂(灰)孤立单元,C′和D′对角线位置存在多个环氧树脂(灰)孤立单元.对4 种结构单胞内的孤立材料单元进行统计,结果如表3 所示.

表3 单胞的孤立材料单元数Table 3 The number of the isolated elements in the unit-cells

图9 未附加可制造性约束多目标优化结果Fig.9 The optimized solutions from the MOOP without the manufacturing constraint

由以上结果分析可见,未考虑可制造性约束的多目标优化结果均存在孤立材料单元,而图5 中附加可制造性约束的优化结果则有效地避免了这种情况的发生,表明本文所提出的考虑可制造性约束的声子晶体多目标拓扑优化方法可在满足带隙等性能指标的前提下同时兼顾到实际制造加工的可行性,有效提高了优化结果的实用性.

3 结论

本文针对声子晶体拓扑优化结果常会出现制造性能欠佳的问题,通过引入可制造性附加约束,建立了考虑可制造性约束的二维多相声子晶体多目标拓扑优化模型,基于NSGA-II 和有限元法对二维四相正方晶格声子晶体面内模态进行多目标拓扑优化设计,得到了多目标Pareto 解集,数值算例表明了本文模型的有效性,所得主要结论如下.

(1) 通过对比引入可制造性附加约束的多目标优化与相应的单目标拓扑优化结果发现,Pareto 解集内可以找到与单目标优化结果接近的非支配解;但多目标的优化解在一定程度上可优于单目标优化近似最优解,表明多目标优化不仅可以得到和单目标优化结果近似非支配解,还可以得到其他的多目标非支配解集.

(2) 通过对比引入和不考虑可制造性附加约束的多目标优化结果发现,引入可制造性附加约束可有效避免通常优化结果中存在孤立材料单元的情况,其优化结果在满足带隙性能指标的同时可兼顾到制造加工的可行性,更具有实用价值.

本文算例针对二维正方晶格,验证了协同考虑带隙性能和可制造性约束的多目标优化方法的可行性和有效性,但该方法对其他类型的晶格,比如二维三角晶格和六角蜂窝晶格、甚至三维体心立方晶格等等,也同样适用.其中,单胞离散及材料编码方案的选择、可制造性约束的添加和优化算法的选择可根据具体问题进行确定.对于三维晶格而言,材料初始分布构型的选取、复杂三维网格的划分、计算规模大、优化效率低等仍是其拓扑优化问题的主要挑战[28],同时也是该方面研究工作未来需要加强的一个方向.

最后需要强调指出的是,影响声子晶体优化结果可制造性的因素很多,比如孤立材料单元问题,超细而导致刚度超弱链接单元问题,不同材料连通域的个数问题,微结构奇形怪状问题和微结构边界或组元界面非光滑问题等.本文提出的多目标优化模型主要以避免出现孤立材料单元而提高多相声子晶体优化结构的可制造性能,其他方法比如避免出现超细超弱链接单元[20,29-30],力求微结构边界或组元边界具有一定的光滑度[31]等也是提高其可制造性的有效措施,值得引起科研和实用工作者的高度重视.