李林 张雪彬 刘涛 章俊

(重庆大学航空航天学院,重庆 400044)

引言

自由地控制波的传播一直是波动调控领域中的热门话题,如利用透镜[1]、光纤[2]、消声器[3]等对光波/声波的调控.作为一种亚波长准二维超材料,超表面[4]的出现使得波的各种反常、新颖调控变得可能,最近十几年得到了众多科研工作者的关注.超表面研究始于电磁波领域[5],随后被引入到声波[6-7]和固体弹性波领域[8-9].作为一种经典的弹性波,弯曲波在工程中十分常见,是梁板类工程结构在横向动态载荷激励下产生的一种波动形式.利用超表面/超材料对弯曲波进行调控在梁板类结构隔振减振[10-11]、基于梁板类结构开发振动能量收集装置[12-15]和各种传感器以及结构健康检测等领域有重要应用.截至目前,科研工作者利用超表面/超材料已经实现了弯曲波的聚焦[16-21]、反常折射和反射[22-25]、隐身[26-35]、完美吸收[36]、分波[37-38]、非对称传输[39-40]和拓扑保护[41-42]等调控.

通过设计结构或材料某种参数可改变的微结构来构成超表面的子功能单元是目前常用的手段,利用这一参数的改变来调控弹性波在子功能单元中的等效相速度,从而实现超表面中每个子功能单元的透射/反射波相位满足一定的空间分布,达到波前调控目的.截止目前,基于上述思想设计的弯曲波超表面子功能单元主要有以下几类: 组合梁型[23]、附着台柱直梁型[24,43]、变厚度梁型[25]和锯齿形声子晶体梁型[21,44].组合梁型子功能单元由两部分不同材料的直梁轴向拼接而成,通过改变这两部分直梁的长度比来调节弯曲波在其中的传播速度.附着台柱直梁型子功能单元由直梁上附着不同几何尺寸的柱状体构成,附着柱状体几何尺寸的不同会改变弯曲波有效相速度.变厚度梁型子功能基于梁中弯曲波相速度与梁厚度相关这一事实而提出.锯齿形声子晶体梁型子功能单元为波纹型梁,通过改变波纹状轮廓曲线的高度来调节梁的等效弯曲刚度,从而调节弯曲波有效相速度.有关超表面的更多细节,读者可以参考北京交通大学陈阿丽教授课题组最近发表的综述文章[4].

随着对超表面研究的深入,研究人员发现当入射角大于某一临界角时,用来预测折射型超表面行为的广义斯涅耳定理(GSL)不再适用.因此,表达式如下的高阶衍射理论[45-47]被提了出来,即

其中kt,r和ki分别为透射/反射波和入射波波数,G=2π/Γ,Γ为超表面中单个超胞在其排布方向的宽度.一般情况下,单个超胞由多个子功能单元(subunit)组成,且这些子功能单元的透射/反射波相位覆盖[0,2π]区间且满足一定的空间分布规律.n表示衍射波阶次,为第n阶衍射波角度.在本文中入射波和衍射波具有相同的波数,将ki=2π/λi代入式(1)可以得到如下关系

其中 η=λi/Γ,λi为入射波波长.当n=1 时,高阶衍射定理即退化为GSL.与GSL 相比,高阶衍射定理除了可以预测第1 阶衍射波,还可以预测其他阶次的衍射波,因此高阶衍射定理可被用来设计分波超表面.基于高阶衍射定理设计的超表面又称为超栅(metagratings),因其具有类似光栅的周期结构.根据式(2)中 η值,基于高阶衍射定理设计的超栅可以分为三类: η >2.0,1.0<η<2.0和η <1.0.由式(2) 可知,第一类超栅中只存在0 阶衍射模式,而且与入射角无关,因此这一类超栅可以实现全角度反射[8,10].第二类超栅中存在0 阶和±1 阶模式,但由于0 阶在子功能单元中传播次数最多[36],很容易被阻尼耗散掉,因此0 阶可以忽略,+1 阶对应透射波模式,−1 阶模式为透射波还是反射波取决于超胞中子功能单元的数目,如果子功能单元的数目为奇数,为透射波,反之为反射波.因此,这类超栅可用来实现非对称传输调控[48].第三类超栅中的衍射模式比前两类超栅都多,因此这类超栅更容易实现分波功能[38,49].

不论是超表面还是超栅,目前的设计大多基于理论直觉、实践经验和两者的结合,这样设计出的超表面和超栅调控性能往往不是最优.此外,此类方法往往只针对单一频率进行设计,导致设计出来的超表面和超栅工作频率范围非常窄,严重限制它们在实际工程中的应用.为了解决这个问题,目前提出了两种办法: 一种是引入主动控制技术[50];另一种为实施拓扑优化设计.主动调控技术通常需要引入额外的结构或装置,这将增大超表面的复杂性.作为一种被动宽频调控技术,拓扑优化设计[51]将数学优化算法引入到超表面的设计中,通过对子功能单元微结构拓扑优化设计来满足透射/反射波相位和幅值要求.例如,Rong 等[52-53]建立了基于多目标遗传算法NSGA-II 的弹性波超表面拓扑优化设计框架,通过在目标函数中设置多个目标相移,优化得到了不同频率下具有不同相位的子功能单元,提出了一种在不同频率有不同调控功能的弹性波超表面.Li等[38]基于遗传算法设计了一类超胞由两个子功能单元组成的结构简单超表面,实现了单一频率下弯曲波的分波和聚焦调控.目前,基于拓扑优化设计超表面得到了越来越多的关注,但宽频超表面的拓扑优化设计工作还较少.拓扑优化所得到的结构形状一般比较复杂,加工是个难题,庆幸的是随着3D打印技术的发展,这类结构的加工难度在逐渐降低.因此,结合拓扑优化和3D 打印技术开展弹性波宽频超表面研究变得可行.

本文基于拓扑优化方法提出一种弯曲波宽频分波超栅设计方法.具体地,利用多目标算法NSGAII 搭建一套二维平面应力稳态动力学问题拓扑优化框架,然后利用这套优化平台优化超栅子功能单元,最后对优化得到的超栅进行数值表征、利用3D 打印技术制备试件开展实验验证,并与其他设计方法进行比较,进一步验证本文提出的设计方法的优势.

1 基于高阶衍射定理的弯曲波宽频分波超栅拓扑优化设计和表征

在前面的引言中已经提到,当 η ≥2 时,n≡0,此时衍射波场中只有0 阶通道存在.随着η的减小,衍射阶次会增多.图1 所示为根据式(2)计算得到的两个不同 η值情况下衍射角与入射角之间的关系.从图中可知,当 η <1时,存在多个衍射阶次,且不同阶次的衍射波有不同的传播方向.因此,可以设计η值小于1.0 的超栅来实现分波功能.本文考虑设计一种弯曲波宽频分波超栅.

图1 不同衍射阶次的衍射角与入射角的关系Fig.1 Relationship between diffraction angles of different diffraction orders and incident angles

为了制备简单起见,本文设计的分波超栅每个超胞由2 个子功能单元组成,且第一个子功能单元选为与主板等厚度的直梁,其尺寸为长度L=50 mm,沿y方向的宽度W=34.36 mm,厚度h=2 mm.第二个子功能单元由优化得到,优化目标为在设定频率范围内由这两个子功能单元所组成的分波超栅功能保持不变.本文设定频率范围为: 中心频率为 ωc=2.4 kHz,上下限在中心频率基础上各偏移10%,即[2160,2640] Hz.图2 所示为第二个子功能单元的二维平面应力优化模型,其中左右两边灰色区域代表与超栅相连接的主板,中间黄色区域为第二个子功能单元的优化区域.此优化区域由5 个胞元周期排列构成,优化区域的总长度为L,总高度为H.考虑到拓扑优化可能得到的结构复杂性,灰色区域和黄色区域采用同一种材料.本文具体采用的是3D 打印材料VeroPureWhite,这样即使优化得到的结构较复杂,也可通过现有的3D 打印技术进行制备.VeroPureWhite 材料的实测杨氏弹性模量E=3.2394 GPa,密度ρ=1185 kg/m3,泊松比 ν=0.4185[21].在有限元模拟中,沿着z方向的单位位移加载在图中的黑色虚线上.为了减少反射波,模型的左右两端设定为完美匹配层(PML).本文考虑了三种不同的优化区域总高度H,目的是想研究不同H情况下能否得到工作机制不同的胞元结构.子功能单元确定之后,按照图3 排列即可形成超栅整体结构.超胞中相邻子功能单元间留有间隙,本文中间隙设为1 mm.

图2 超栅子功能单元二维拓扑优化模型示意图Fig.2 Schematic of topology optimization 2D models for the subunits in a supercell

图3 薄板弯曲波超栅整体结构俯视图,本文中子功能单元数为2Fig.3 Top view of metagratings for flexural waves in thin plates.In this paper,the number of subunits is 2

1.1 定义适应度函数

优化胞元区域被离散成一定密度的像素块集合,每个像素块的质量密度都是一个独立的设计变量,并通过适应度函数最小化来确定像素块的密度分布.在本文中,适应度函数定义如下

其中,S为频率区间[2160,2640] Hz 内的采样频率点数.考虑到数值计算成本,本文中S=3,对应的频率分别为频率区间的下限、中心频率和频率区间的上限.其中和是第j个采样频率下的目标相移和目标透射率;ϕj和Tj是优化过程中第j个采样频率下的可行相移和可行透射率.可行相移 ϕj=∆s·k0,∆s为经过第一和第二两个子功能单元的透射波相应波峰或波谷的位置差,Tj等于透射波幅值之比.由于本文研究的是透射分波超栅,所以目标透射率=1(j=1,2,···,S).不同频率下,目标相移是不一样的,具体由以下公式确定

式(3)中的vios为罚函数.在本文中,它由两部分组成.第一部分为

其中nd,Ad分别为非连通区域的数量和非承载区域的像素数量.f1,f2为相应的罚参数.在本文中,f1=0.5,f2=0.1.罚参数满足f1>f2以便有效地促进结构的连通性.这一部分罚函数的作用是对不可行的拓扑结构进行惩罚,比如断开的拓扑形状、不能提供承载作用的像素[54-55].第二部分为

其中ft=0.1 是罚参数,nt为薄弱连接区域的个数.在本文中,通过判断沿x和y方向单列连续的拓扑结构的长度是否小于预设的最小加工尺寸来确定薄弱连接.这一部分的作用是为了提高拓扑结构的可制造性和可靠性.除了这些罚函数外,在计算适应度函数之前,还应对拓扑形状中的铰接情况进行处理以提高优化所得结构的承载能力和可加工性.

本文采用了“粗−细”策略的遗传算法[56],在第一阶段,优化区域x方向被划分为20 个网格,在第二阶段则被划分为 40 个网格.y方向的网格密度与x方向始终保持一致.在这两个阶段,染色体种群的个体数均为Np=100,染色体交叉概率为Pc=0.9,变异概率为Pm=0.02,总的迭代次数设定为1000 次.

具体的拓扑优化流程与已有文献报道中的流程基本一致[57],但是本文在子代的选择方面做了一些改变,并没有完全采用精英化策略[54],而是在轮盘赌和精英化策略之间采用概率选择.具体地,在本文中,轮盘赌和精英化选择的概率各为50%.这样处理的好处是可以避免遗传算法过早地收敛.

1.2 胞元优化结果

图4 所示为在中心频率2.4 kHz 下,优化区域高度为H=3h情况下的胞元拓扑优化结构进化历程,从此图可以看到,优化从随机生成的连接性较差的结构开始,经过一百代左右优化之后,拓扑结构的连通性明显改善.之后,遗传算法搜索到有效的拓扑特征并慢慢演化结构,结构连通性在整个演化过程中逐渐提高.

图4 优化区域高度H=3h,在2.4 kHz 中心频率激励下的胞元拓扑结构的进化历程Fig.4 Evolutionary history of topological structure of a unit cell as the optimization region height H=3h and under excitation of a 2.4 kHz loading

图5 所示为三种不同优化区域高度下胞元的最终拓扑优化结构.图6 所示为以图5 中三种拓扑结构为胞元分别形成的三种一维声子晶体梁的能带图,其中最低的黑色扩散曲线对应弯曲波模式.扩散曲线结果显示这三种声子晶体梁的工作机理相同,均为通过改变等效弯曲刚度从而实现有效相速度的改变.图7 所示为由图5 中三种优化结构分别组成的子功能单元在2.4 kHz 中心频率激励下的透射模拟结果.图8 所示为由这三种拓扑结构所构成的子功能单元在三个采样点频率下的透射率和相移与目标值的比较.在中心频率下,每个胞元准确地实现了所需的相移;在其他两个频率点,实际相移与目标相移存在小的偏差,这样的结果符合预期;透射率均大于0.9,表明所设计的透射超栅将具有较高的透射率.总的来说,H=2h情况下优化得到的胞元(见图5(a))性能最优.

图5 三种不同优化区域高度下所得到的胞元最终拓扑优化结构Fig.5 The optimized topological structures of unit cells under three different optimization region heights

图6 以图5 中三种拓扑结构为胞元分别形成的三种一维声子晶体梁的能带图Fig.6 Band structures of the three 1D phononic crystals constituted by the three units presented in Fig.5

图6 以图5 中三种拓扑结构为胞元分别形成的三种一维声子晶体梁的能带图(续)Fig.6 Band structures of the three 1D phononic crystals constituted by the three units presented in Fig.5 (continued)

图7 三种拓扑优化结构分别组成的子功能单元在中心频率激励下透射波行为数值模拟表征结果Fig.7 Numerical characterization of wave transmission behavior through three subunits composed by each of the three optimized structures under the central frequency excitation

图8 分别由三种拓扑优化结构构成的三种子功能单元在三个采样点频率下的(a)相移和(b)透射率Fig.8 (a) Phase shifts and (b) transmittances of subunits composed by the three optimized structures at three different frequencies

1.3 超栅整体结构数值模拟表征

基于图5 中优化得到的胞元,形成第二个子功能单元(见图7),然后将超胞(第一和第二个子功能单元)沿y方向周期排列即可形成超栅(见图3).考虑到数值模拟和实际加工成本以及加工技术限制,本文中的超栅由四个超胞组成,胞元沿y方向的宽度为34.36 mm.因此,中心频率下超栅的 η=0.707.图9 所示为由图5 中三种胞元分别形成的三种超栅整体结构.为了后面表述方便,这三种超栅分别被称为超栅1、超栅2 和超栅3.为了验证这三种分波超栅的性能,建立如图10 所示有限元数值模型,弯曲波由线源激发,单位载荷沿z方向加载.图11 展示了这三种分波超栅在三个采样点频率激励下的响应结果.从图中可以清晰地看出,正入射的弯曲波在经过超栅之后被分成了两束,传播角度与高阶衍射理论预期角度(图中标注值)非常吻合,且在三个频率下的角度基本一致,说明超栅在本文关注的频率范围内功能保持不变.

图9 由图5 中三种胞元分别形成的三种超栅整体结构Fig.9 Three metagratings constituted by each of the three unit cells presented in Fig.5

图10 超栅结构有限元数值模拟模型Fig.10 Numerical simulation models of metagratings with the finite element method

图11 三种超栅在三个采样点频率激励下响应的数值模拟结果Fig.11 Numerical simulation results for the three metagratings under excitation of three different frequencies

1.4 实验验证

为了进一步验证本文基于优化设计所得到的宽频超栅性能,本文利用3D 打印技术制备了两种试件.一种是与图9(b)所对应的超栅整体结构,另外一种是基于优化所得到的单个子功能单元,如图12 所示.实验测试平台如图13 所示,主要设备包含一套激光扫描测振仪系统(Polytec PSV-600)、一个功率放大器、一个相位转换装置(自制)和一台信号发生器(ATA-308,Agitek).为了产生弯曲波,在试件正反面相同位置各粘一排压电片(无锡海鹰).由信号发生器产生的电压信号经过功率放大器放大,然后经过相位转换装置出来电压相反的两路信号.随后,这两路信号分别输入到正反面的压电片上,在相反的两个电压作用下,上下表面的压电片将产生反向的面内变形,相当于对试件施加一谐变弯矩,从而产生弯曲波.离面质点振动速度由激光测振仪测量.为了减少边界的反射波,试件的边界被粘上一层阻尼材料(Blu-Tack).整个试件使用弹性绳悬挂在框架上.

图12 实验试件Fig.12 Experiment specimens

图12 实验试件(续)Fig.12 Experiment specimens (continued)

图13 本研究实验测试平台Fig.13 Experimental setup in this study

图14～ 图16 所示为超栅整体结构在三个采样点频率激励下的实验结果.结果表明,在三个频率下,入射波经过超栅之后均被分成两束,这两束透射波相对入射波的折射角度为 ±45°,与理论值和前面的数值模拟结果均一致,表明本文基于拓扑优化所设计的超栅确实在设定的宽频范围内具有所设计的功能.

图14 激励频率为2160 Hz 时,测得的三个时刻的超栅离面速度场Fig.14 Snapshots of measured out-of-plane velocity fields under the excitation of 2160 Hz at three instants

图15 激励频率为2400 Hz 时,测得的三个时刻的超栅离面速度场Fig.15 Snapshots of measured out-of-plane velocity fields under the excitation of 2400 Hz at three instants

图16 激励频率为2640 Hz 时,测得的三个时刻的超栅离面速度场Fig.16 Snapshots of measured out-of-plane velocity fields under the excitation of 2640 Hz at three instants

除了上面的超栅整体结构,本文还测试了单个子功能单元的结果.图17 所示为第一个子功能单元和由图5(b)胞元所构成的第二个子功能单元在三个不同频率激励下的实验结果.从图中可以清楚地看到经过它们之后的透射波相移与理论预期非常吻合,进一步验证了本文基于拓扑优化方法设计宽频超栅的可行性.

图17 三个采样点频率激励下,拓扑优化梁与直梁实测相位比较Fig.17 Comparison of phases of transmitted waves between optimized beams and straight beams at three different frequencies

2 与其他方法设计的超栅的宽频优势比较

为了进一步展示本文拓扑优化设计超栅的宽频特性,把本文的方法与已有的两种弯曲波超栅设计方法进行了对比.第一种方法为前面引言中提及的锯齿形声子晶体梁方法;另一种方法为变厚度梁方法[10,25].基于这两种方法本文设计了相应的超栅结构,分别称为锯齿型超栅和直梁型超栅.然后对超栅整体结构进行数值表征,并与本文的超栅进行对比.两种超栅的超胞均由两个不同高度的子功能单元构成,根据相位要求,两直梁高度分别为2 mm 和6.4 mm,两个锯齿形声子晶体梁的高度分别选为1h和3.9h,其他参数与本文的超栅一致.图18 和图19 分别展示了直梁型和锯齿型分波超栅在[2160,2640] Hz 范围内的数值模拟结果.可以看出,直梁型分波超栅的性能较差,锯齿型分波超栅在前两个频率激励下的性能也不太理想,两者均无法在本文所考虑的宽频范围内实现稳定的分波功能.本节内容表明本文提出的基于拓扑优化的宽频超栅设计方法非常有必要,为宽频超栅设计提供了一种可行路径.

图18 直梁型分波超栅数值模拟结果Fig.18 Numerical simulation results of straight-beam-type wave-splitting metagratings

图19 锯齿型分波超栅数值模拟结果Fig.19 Numerical simulation results of zigzag-type wave-splitting metagratings

3 结论

本文提出了一种基于拓扑优化的弯曲波宽频超栅设计方法.利用遗传优化算法设计一系列满足高阶衍射理论的亚波长胞元结构,这些胞元组成的超胞产生的相位能够覆盖[0,2π]范围,并在目标宽频范围内具有高透射率.然后将这些超胞组合以形成具有特定相位梯度的超栅,并通过数值模拟和实验验证了超栅的分波功能和宽频特性.此外,通过与直梁型分波超栅和锯齿型分波超栅进行对比,发现本文的超栅宽频性能明显优于这两种超栅,显示出本文提出的基于拓扑优化的宽频超栅设计方法的必要性和有效性.本文提出的方法同样也适用于其他宽频超栅的设计.