常见的禁带构造机理有3种，分别可以构造布拉格禁带、局域共振禁带和耦合禁带。布拉格禁带是周期性变化的材料或几何特性与弹性波相互作用形成的，结构中各组分的弹性常数和密度越大越容易产生禁带。这一禁带由于其宽频特性在减振方面有一些应用[12-14]，但其受制于产生的布拉格条件，导致该禁带往往位于高频，因此难以实现低频减振[15-18]。为了突破这一限制，Liu等[19]提出了局域共振禁带的构造机理。该禁带主要由结构中谐振子的共振产生，导致禁带波长可以随其共振频率变化，从而实现低频隔振。然而其禁带带宽很窄，虽然有研究对其进行拓宽。如朱兴一等[20]就研究了谐振子的排列方式、材料性质等对禁带调节的影响，以获得宽频禁带。但这种定性研究大多是为了研究禁带的演变规律，而非为了满足工程减振的应用需求[3, 10-21]。耦合禁带是结构中2种波由于耦合效应产生的[22-23]，调节一种波的传播特性产生的耦合禁带可以影响另一种波，因此具有丰富的设计空间和调节裕度，为突破布拉格禁带和局域共振禁带的局限提供可能。但机械结构内耦合禁带[24-25]随着结构确定而确定，因此无法针对特定频带进行设计。近年来，压电材料的出现为实现可调耦合禁带提供了可能。
压电材料作为一种成熟的商用智能材料，结合电路后可以对结构中的能量进行收集或耗散[8, 26-29]。因此，Li等[30]在压电周期板中施加了电感电路网络，形成压电网络板，利用其中的电波和结构中的弯曲波通过压电效应进行耦合，从而产生了可调且频带较宽的耦合禁带。此处的"电波"实际上是一种电磁波，只是由于存在电感等电学元件，导致其波长增加到可以与弯曲波在弹性波传导的低频段产生耦合禁带的程度。故分析压电结构时，一般认为该电磁波的磁场分量与电场分量相比可以忽略[31]。用术语"电波"来代替一般意义上的"电磁波"，并只用电场量(电势)对其进行描述。
考虑到文献[30]是针对可调耦合禁带生成的可行性研究，并未分析不同压电片形状对耦合禁带的影响。实际上，大部分基于压电周期结构的研究都集中在通过设计主动或被动电路的方式，以调整禁带的位置和带宽[32-33]，较少有研究探索压电片形状和尺寸对禁带特性的影响。然而为了在低频得到宽频禁带，压电片几何尺寸的设计也是一个十分关键的因素。同时在航空航天等领域，轻质化一直是设计的核心需求，对于附加质量和形状有着严格的限制。因此为了在最小质量比的条件下得到最佳性能的耦合禁带，需要从压电网络板中耦合禁带的影响规律入手，从而为选取最佳的压电片尺寸提供理论依据。
考虑到研究尺寸变化需要计算多种模型的频散曲线，为了提高计算效率，采用了减缩的波有限元算法[34]。波有限元算法的主要优势在于可以借助商业有限元软件建立模型，导出动力学矩阵后施加周期性边界条件便可求解系统的频散特性，确定禁带分布[35]。在减缩模型方面，Zhou等[36-37]将固定界面模态综合法[38]与波有限元算法结合，将元胞(结构最小可重复单元)内部自由度用模态坐标进行代替，从而进一步提高了计算效率。但这种针对机械系统的减缩波有限元算法无法直接用于机电耦合系统。因此，对该算法进行了改进，推导了适用于压电周期板的减缩波有限元算法。
考虑到结构中主要依赖弯曲波和电波的耦合，因此对保留模态阶数和仅保留弯曲模态的情况进行了误差分析。在此基础上，针对市场上常见的圆形和正方形压电片，研究了其尺寸对耦合禁带性能的影响。同时设计电学参数，对比了不同压电模型在同一位置产生耦合禁带时的性能差异。为工程应用中如何在特定频带设计合适的耦合禁带提供指导。最后利用强迫响应对耦合禁带的减振效果进行了验证。为压电周期结构的合理设计提供了计算工具和方法指导。
1 压电结构的减缩波有限元算法 用波有限元算法分析周期结构的波传导特性包含2个主要步骤。首先，通过商用软件或有限元子程序获得一个元胞的动力学方程组；其次，在元胞的边界自由度(与相邻元胞共有)上施加周期边界条件，构造求解自由波动特性的特征值问题。对于压电结构的波传导特性分析也是如此，只需在建立元胞动力学模型时考虑压电材料及其电路的影响即可。本节将按照上述步骤阐述具有网络电路的压电板的波动特性分析流程，并着重介绍提出的元胞动力学模型减缩算法。
1.1 压电网络板的元胞动力学模型 本文研究的压电网络板如图 1(a)所示。
 |
| 图 1 压电网络板及元胞示意图 Fig. 1 Illustration of plate with piezoelectric network and its unit cell |
| 图选项 |
图 1(a)中压电片周期性贴在基板上，相邻压电片通过相同的电感相连，以构成电路网络。该周期结构可以由虚线方框内的元胞沿x方向和y方向延拓得到。考虑到机电耦合系统内含有电场和机械场的自由度，在后续处理时需要分别处理，故画出元胞所有机械场自由度向量a分布如图 1(b)所示。由图 1(a)可知，施加电路网络后，会在元胞边界处增加4个电学节点，其分布如图 1(c)所示。
考虑到结构中不只存在一种机械波，需要确保其他波(如剪切波等)不会对耦合禁带的减振效果产生影响[11]。因此如图 1(c)所示，将压电片对称贴在基板两面的相同位置，且极化方向相同，从而保证弯曲波上下表面刚度一致，结构内的剪切波不会导致元胞发生弯曲变形。图中，v表示单节点的自由度。另外设置压电片与基板相连面接地，另一面并联后与电感电路相连，可以使拉压波产生的电荷相互抵消，从而拉压波也可以与电波解耦[10]。
由于难以获得图 1(a)中压电网络板波动特性的理论解，故通过波有限元算法对图 1(c)中的元胞进行建模和求解，以获得系统内的禁带分布。
由周期结构理论[39]可知，系统波动特性求解时需要给元胞施加复数形式的边界条件，而大部分有限元软件难以完成这一设置。同时利用有限元软件建立元胞模型时且不支持导出含电路网络的动力学矩阵，故难以直接建立这一元胞的动力学方程。因此此处利用软件ANSYS建立该压电系统的机械场模型，导出其动力学矩阵后，编程实现与电场动力学矩阵的组集，以得到压电网络板元胞的动力学方程。对其施加边界条件后，可求解其波动特性。基板选取SHELL181单元，压电片由SOLID5单元构成。基板材料为环氧树脂，密度为1 180 kg/m3，弹性模量为4.35×109 N/m2，泊松比为0.37。压电片材料为PZT-5H，其材料常数如表 1所示，表中: ρp为压电片密度，s11E、s12E、s55E为弹性柔顺常数，d31为压电应变常数，ε33T为厚度方向的恒应力介电常数。固定基板尺寸为80 mm×80 mm×5 mm，选取最常用的正方形和圆形2种形状的压电片，分别研究其长度和直径对耦合禁带的影响。
表 1 压电片材料参数 Table 1 Material parameters of piezoelectric patches
| 参数 | 数值 |
| ρp/(kg·m-3) | 7 500 |
| s11E/(m2·N-1) | 13×10-12 |
| s12E/(m2·N-1) | -4.29×10-12 |
| s55E/(m2·N-1) | 22×10-12 |
| d31/(C·N-1) | -1.86×10-10 |
| ε33T/(F·m-1) | 3.009×10-8 |
表选项
同时为了纵向对比形状对禁带的影响，设置每一组正方形和圆形压电片的面积相等，即Lp2=πDp2/4，Lp为正方形压电片边长，Dp为圆形压电片直径。
建立不含电路的压电网络板元胞有限元模型后，可以得到如下的动力学方程：
 | (1) |
式中：a为节点位移向量；f为节点力向量；Qi为压电片表面电荷量；vi为压电片上表面电势；m为元胞的质量矩阵；k为元胞的刚度矩阵；kp为耦合矩阵；Cp为常应变下压电片的电容值，下标p表示压电片相关参数；下标i为内部节点；上标T为转置符号。
为了方便建立该电感电路网络的动力学方程，利用有限元算法的思路，建立单个电感的单元动力学矩阵后，依据节点位置进行组集，可得该电路网络的动力学方程如下：
 | (2) |
式中：Le为网络中的电感值；ω为谐响应下的频率；v=[vi v5 v6 v7 v8]T为电学节点的电势；Q=[Qi Q5 Q6 Q7 Q8]T为电学节点对应的电荷分布；下标5、6、7、8分别表示左边节点、底边节点、右边节点和上端节点，其位置如图 1(c)所示。
由于施加电路网络后，压电片上表面电势vi的电荷满足Qi=0。将方程和按照这一关系进行组集，得到压电网络板元胞的动力学方程如下：
 | (3) |
式中：E为电场的动刚度矩阵，其中的分块矩阵为
对方程施加周期性边界条件后即可求解未减缩模型的波动特性。为了确保计算结果的正确性，每个波长内至少需要有6~12个网格，导致模型自由度较大，计算费时。同时考虑到需要计算不同压电片尺寸和形状下的频散曲线，为了节约计算时间，在求解前对方程进行减缩处理。
1.2 压电网络板元胞减缩模型及波动特性求解 该减缩主要基于Craig-Bampton方法，通过借鉴固定界面模态综合法的思路[35]，将自由度分为保留自由度和减缩自由度，将保留自由度固定后，对减缩自由度进行模态分析，获得将减缩自由度转到模态坐标的转换矩阵，即可完成减缩。由于电场的动刚度矩阵E与频率有关，将电场自由度归为减缩自由度求解得到的固定界面主模态无法进行对角化，导致不同模态间存在耦合，从而截断高阶模态会导致计算出错。故只对机械场内部自由度ai进行减缩，保留自由度为qR=[abT vi vbT]T，下标R意为保留节点。包括机械场边界节点ab，电学内部节点vi，电学边界节点vb。按照这一关系对方程中的矩阵重新划分为
 | (4) |
式中：fR=[fbT 0 QbT]T为机械场边界自由度、电学内部自由度以及电学边界自由度对应的边界力，此时只有一个电学内部自由度。根据区分的机械场内部自由度ai和保留自由度qR，将式(4)的刚度矩阵K、质量矩阵M和电学矩阵H重新分块，有
 | (5) |
 | (6) |
 | (7) |
由固定界面模态综合法可知，将内部自由度ai转化到模态坐标y的转换矩阵B为
 | (8) |
式中：I为单位矩阵；Ψ为固定界面主模态，通过固定保留自由度qR=0，求解如下特征值方程获得
 | (9) |
考虑到高模态对描述内部自由度的贡献较小，故仅保留前n阶模态，最终得到的固定界面主模态为Ψ=[ψ1T ψ2T … ψnT]T。保留模态阶数选取的原则为
 | (10) |
即保留模态的固有频率ωi不超过分析频带的上界ωN的若干倍。通过选取不同的系数αf来控制保留的模态阶数。将式(8)代入方程(4)，并左乘转换矩阵B的转置，即可得到减缩后的动力学矩阵：
 | (11) |
式中：
 | (12) |
式中：PC=[y vi]T，包括模态坐标y和电学内部自由度vi；qb=[abT vbT]T为广义坐标的边界自由度，包括边界的位移自由度ab和电势自由度vb; Fb=[fbT QbT]T为广义坐标边界节点的力向量。考虑到求解前需要对边界节点施加周期性边界条件，故需要按照节点位置进行划分。另外由于内部节点不受外力，因此可进一步减缩求解矩阵的规模：
 | (13) |
式中:
由Bloch边界理论可知[39]，节点自由度应当满足如下关系：
 | (14) |
式中：
 | (15) |
同样，根据节点力平衡可知：
 | (16) |
式中：λx=eμx；λy=eμy; μx=-jkxLx; μy=-jkyLy分别为沿着x和y方向的传播常数，其幅值和相位分别反映波在经过一个周期单元的传播后产生的幅值衰减和相位变化。kx和ky为沿着x和y方向的传播常数的波数，其实部为传播常数，虚部为衰减常数。Lx和Ly为元胞沿着x和y方向的尺寸，此处由于选取了正方形元胞，有Lx=Ly=80 mm。
将式(14)和式(16)代入方程(13)可得
 | (17) |
方程中包含3个未知量ω、λx和λy。其中(λx, λy)表征了板结构内波的方向性分布，可以将其转化为波矢(kx, ky)。由周期结构理论可知，只需考虑分布在布里渊区([(-π/Lx, π/Lx), (-π/Ly, π/Ly)])中波矢(kx, ky)的分布，就可以刻画系统的传播特性。然而根据结构的对称性，往往只需计算得到位于简约布里渊区([(0, π/Lx), (0, π/Ly)])的分布。因此在不同频率ω下，固定波矢中的一个值以求解另一个。从而可以得到该频率下波空间内所有传播波的分布，如果存在某些角度范围内不存在传播波，则为禁带。
固定频率ω和(λx, λy)中的一个(譬如λy)，就可以将其转化为二次特征值问题进行λx的求解[11]，从而获得减缩后压电网络板沿某一方向的波动特性。另外对于减缩模型求解出的特征向量
| $$ | (18) |
对方程(4)直接施加周期性边界条件求解，则可得到未减缩模型的频散曲线。
1.3 保留模态对计算精度的影响分析 选取尺寸为Lp=40 mm的正方形压电片，设定电感Le=0.08 H。选取0~2 000 Hz为目标频带，保留模态系数αf=5。以验证减缩效果。模型1有限元模型如图 2所示。
 |
| 图 2 含正方形压电片元胞有限元模型 Fig. 2 Finite element model of unit cell with square lead zirconate titanate patch |
| 图选项 |
压电片贴于基板中央。计算其λy=1且未加电路的波动特性，得到系统沿x方向的频散曲线，结果如图 3所示，图中横坐标为频率，纵坐标为传播常数，表明弹性波相位的变化。可见其中存在3种波，其中由于压电片在基板中的周期分布，在755~951 Hz内产生了波1的布拉格禁带。
 |
| 图 3 未加电路时系统传播常数分布 Fig. 3 Distribution of system propagation constants without circuits |
| 图选项 |
施加电路网络后依然令λy=1，即求解系统沿x方向的频散曲线，结果如图 4所示，图中横坐标都为频率，图 4(a)中纵坐标为传播常数，图 4(b)中纵坐标为衰减常数。图中虚线为未减缩计算结果，实线为减缩计算结果。
 |
| 图 4 αf=5时减缩前后频散曲线对比 Fig. 4 Comparison of frequency dispersion curves with/without reduction when αf=5 |
| 图选项 |
由图 4可知，减缩前后频散曲线吻合。其中图 4中波2为施加电路网络后产生的电波，该网络形式下，波数随频率增加而增大，经过耦合禁带后进入永久禁带。波1和波2在目标频带各存在一个布拉格禁带，即满足Real(kx)=π/Lx。同时波1和波2在1 300 Hz附近，由于耦合作用产生了耦合禁带，该禁带的衰减常数和带宽都优于波1的布拉格禁带。由于压电效应只影响耦合区附近的频散曲线[22]，因此随着耦合结束，波1的波数依然随频率增加而减小，波2的波数随频率增加而增大。画出700 Hz下各个波的波形如图 5所示。
 |
| 图 5 700 Hz下各行波波形 Fig. 5 Wave shapes of propagation waves at 700 Hz |
| 图选项 |
可知波1和波2有弯曲变形，故为机电耦合波，波3为纵波，波4为剪切波。
减缩矩阵规模和效率提升如表 2所示。
表 2 减缩模型对比 Table 2 Comparison between full model and reduced model
| 模型对比 | 计算时间/s | 矩阵规模 | 自由度数 |
| 完整模型 | 2 843.17 | 8 365×8 365 | 1 175 |
| 减缩模型 | 227.98 | 313×313 | 297 |
| 增益/% | 91.98 | 96.25 | 74.72 |
表选项
由表 2可知，该减缩算法可以有效提升计算效率。另外选取弯机电耦合波1的传播常数做误差分析，结果如图 6所示。由图 6可知，其误差基本在0.01%附近，可见该算法在提升效率的同时还能保证结果的准确性。
 |
| 图 6 各频率下模型1内波1减缩前后误差 Fig. 6 Errors of wave 1 with/without reduction in model 1 at different frequency |
| 图选项 |
值得注意的是，之所以选取保留模态系数αf=5，是因为求解特征值方程时发现其特征频率较高，如果αf过低则会导致保留模态数过少而导致计算结果误差较大。图 7给出了αf=3时传播常数的对比，图中虚线为减缩计算结果，实线为未减缩计算结果。可见虽然减缩后纵波3和剪切波4的吻合较好，波1远离耦合区的频散曲线也在一定程度上吻合。但波1布拉格禁带后的波数以及波2的波数计算的误差都较大，因此在进行减缩时需要酌情选取αf的值。
 |
| 图 7 αf=3时减缩前后传播常数对比 Fig. 7 Comparison of propagation constants with/without reduction when αf=3 |
| 图选项 |
考虑到系数αf对剪切波和纵波影响较小，同时耦合禁带主要是弯曲波和电波为主。因此为了进一步降低保留的模态阶数，在αf=5时只考虑其中的弯曲模态。其对比如图 8所示。
 |
| 图 8 αf=5时只考虑弯曲模态减缩前后传播常数对比 Fig. 8 Comparison of propagation constants with/without reduction when αf=5 with only bending modes considered |
| 图选项 |
此时虽然波1和波2布拉格禁带前的波数计算结果都较好，但布拉格禁带后的波数依然存在可见误差。因此无法对保留的模态阶数进行进一步的选择。
2 结构尺寸对耦合禁带影响分析 由图 4可知，当电波和弯曲波在耦合禁带内，剪切波和纵波依然可以传播，因此只考虑结构中的传播波无法确定禁带的分布。故采用单波形占比的指标RSW以消除纵波和剪切波的影响。这一指标的详细定义参见文献[11]。
该指标主要通过设定某种特定变形的参考波模态Φref，计算其与完全波模态Φfull的模态置信因子(MAC)：
 | (19) |
式中：上标H为共轭转置。RSW是位于[0, 1]的实数，0表示完全波模态Φfull中不含有变形Φref，1则为完全波模态Φfull与参考波模态Φref变形完全一致。每种波的波模态包含7类自由度，Φfull=[uxT, uyT, uzT, θxT, θyT, θzT, λeT]T，即3个方向的平动自由度(ux, uy, uz)和3个方向的转动自由度(θx, θy, θz)以及电学自由度λe(磁通量)。对于波模态中弯曲变形的参考波模态有
 | (20) |
下标f代表弯曲波(Flexural Wave)。画出电波和弯曲波的占比如图 9所示。电波对应的参考波模态为只保留λe。由图 9(a)可知，弯曲波的RSWf除了耦合禁带附近都接近于1，表明该波完全由弯曲变形主导；图 9(b)中电波的RSWe也有类似的结论。在耦合禁带附近，弯曲波和电波的单波形占比变化剧烈(由接近1下降为0.5左右)，说明二者存在强烈的能量交换，同时验证了元胞压电片布置形式和极化方向的有效性。设置弯曲波的占比RSWf>0.2，以排除剪切波和纵波的影响。
 |
| 图 9 电波和弯曲波占比分布 Fig. 9 Distribution for RSWf and RSWe |
| 图选项 |
同时考虑到板为二维结构，需要考虑不同方向上耦合禁带性能的差异。所以计算不同λy下λx的值，以得到禁带的方向性分布。故使用减缩算法可以节约大量计算时间。
2.1 正方形压电片尺寸变化 参考表 1，计算5种边长不同的正方形压电片的机电系统，具体尺寸为Lp={20, 35, 40, 45, 50} mm，设定电感Le=0.1 H，使得在0~2 000 Hz产生耦合禁带和布拉格禁带，利用减缩算法研究其波动特性。
当LP=20 mm时，得到x方向频散曲线如图 10所示。其中布拉格禁带范围为736~774 Hz，耦合禁带范围为1 562~1 710 Hz。可知当压电片面积较小时，基板和散射体刚度相差较小，导致布拉格禁带较小，同时将机械能转化为电能的能力有限，因此耦合禁带效果也有限。但耦合禁带依然比布拉格禁带宽110 Hz，衰减常数是布拉格禁带的3倍。
 |
| 图 10 LP=20 mm时频散曲线 Fig. 10 Frequency dispersion curves when LP=20 mm |
| 图选项 |
调整λy的值，反复计算可以得到该结构参数下禁带的方向性分布，如图 11所示，图中横坐标为频率，纵坐标为禁带存在的角度。图中左侧为布拉格禁带分布，右侧区域为耦合禁带分布，考虑到选取的元胞为正方形，因此根据对称性只用画出[0，π/4]的区域即可表示其方向性。
 |
| 图 11 LP=20 mm时禁带方向性分布 Fig. 11 Directional distribution of band gaps when LP=20 mm |
| 图选项 |
由图 11可知，随着角度增加，布拉格禁带往高频移动，且带宽迅速降低。而耦合禁带虽然也向高频移动，但带宽近似不变。所以当压电片尺寸过小时，难以利用结构中的布拉格禁带，而耦合禁带依然有较好的性能。考虑到研究对象为二维板，因此主要研究其尺寸对禁带方向性的影响，x方向的波动特性可以从角度分布图中看出，因此后续仅展示角度分布图。
保持电感Le=0.1 H不变，进一步增大其尺寸到35 mm，可得禁带的方向性分布如图 12所示。增大压电片尺寸可以明显拓宽布拉格禁带和耦合禁带。在相同电感值的情况下，耦合禁带向低频移动，原因在于增大压电片尺寸使得其本征电容增大，导致电波向低频移动从而耦合禁带出现在较低频率范围。
 |
| 图 12 LP=35 mm时禁带方向性分布 Fig. 12 Directional distribution of band gaps when LP=35 mm |
| 图选项 |
当增大压电片尺寸后，布拉格禁带虽然得到了拓宽，但其带宽随着传播角度增加依然明显衰减。相比于20 mm的情况，耦合禁带虽然整体在向低频移动，同时其带宽也随着角度增加而略微降低，即沿着x方向带宽最宽。
当LP=40 mm时，画出其方向性分布如图 13所示。可见当尺寸增大到一定程度后，对布拉格禁带影响近似不变，然而压电片增大导致其电容增加，使得耦合禁带不断向低频移动，与布拉格禁带靠近后，二者相互作用，使得在较高的传播角度上具有较宽的禁带。因此在利用耦合禁带进行方向性隔振时，考虑到最佳角度附近禁带带宽变化剧烈，需要有针对性地设计角度范围，以达到最佳的禁带宽度。
 |
| 图 13 LP=40 mm时禁带方向性分布 Fig. 13 Directional distribution of band gaps when LP=40 mm |
| 图选项 |
当压电片尺寸增大到LP=45 mm时，在x方向上耦合禁带和布拉格禁带也开始相互影响。其方向性分布如图 14所示。与40 mm的情况相比，虽然其带宽较宽且频率较低，但随着角度增加其带宽迅速降低。此时耦合禁带更适合用于沿x方向的隔振或减振，以达到最佳减振效果。
 |
| 图 14 LP=45 mm时禁带方向性分布 Fig. 14 Directional distribution of band gaps when LP=45 mm |
| 图选项 |
当压电片尺寸增大到LP=50 mm时，布拉格禁带和耦合禁带的相互作用越发强烈，此时其方向性分布如图 15所示。可见二者几乎合二为一，耦合禁带具有了布拉格禁带的特征，即随着角度增加带宽迅速降低，然而原布拉格禁带在角度较高时带宽却近似不变。可见在调整压电片尺寸的
 |
| 图 15 LP=50 mm时禁带方向性分布 Fig. 15 Directional distribution of band gaps with LP=50 mm |
| 图选项 |
同时，需要配合设计电学参数，使得布拉格禁带和耦合禁带都具有最佳性能。
2.2 圆形压电片尺寸变化 选取圆形的压电片，计算5种直径的机电系统，有限元模型如图 16所示。
 |
| 图 16 含圆形压电片元胞有限元模型 Fig. 16 Finite element model of unit cell with circular lead zirconate titanate patch |
| 图选项 |
设定电感不变Le=0.1 H，选取圆形尺寸为Dp={22.57, 39.49, 45.14, 50.78, 56.42}mm，计算0~2 000 Hz内禁带的波动特性和方向性分布。当DP=22.57 mm时，得到x方向频散曲线如图 17所示。
 |
| 图 17 DP=22.57 mm时频散曲线 Fig. 17 Frequency dispersion curves with DP=22.57 mm |
| 图选项 |
与同一面积的方形压电片相比，参考图 10可知，耦合禁带和布拉格禁带都相差不大，对比耦合禁带边界频率发现圆形压电片的耦合禁带略窄，但二者的中心频率不变。原因在于二者压电片电容近似相等，导致未耦合前电波具有一致的频散特性，因此中心频率不变。然而正方形压电片将机械能转化为电能的能力略优，导致耦合禁带略宽。
对比同一面积的2种压电片，得到其方向性分布如图 18所示。由图可知，同一面积下，不同性状对耦合禁带影响有限，正方形情况下耦合禁带带宽略宽。由于布拉格禁带变化不大，因此并未在图中列出。可见在不同传播方向上，正方形情况依然略优于圆形情况，当传播角度高到一定程度时，耦合禁带带宽才趋于一致。
 |
| 图 18 DP=22.57 mm时耦合禁带方向性对比 Fig. 18 Comparison of directional distribution for coupled band gaps with DP=22.57 mm |
| 图选项 |
进一步增加直径至DP=39.49 mm，计算其方向性分布与同等面积的正方形压电片对比如图 19所示。二者依然具有相似的方向性分布，只是随着面积的增大，在x方向的带宽差异扩大，这一差异随着角度增加而减小。
 |
| 图 19 DP=39.49 mm时禁带方向性分布 Fig. 19 Directional distribution of band gaps with DP=39.49 mm |
| 图选项 |
进一步增大压电片直径至45.14 mm(见图 20)，由于此时耦合禁带和布拉格禁带相互作用，导致圆形和方形的耦合禁带在高频时发生微小差异。但二者的中心频率一直保持不变。类似的现象在DP=50.78 mm时都可以看到，如图 21所示。
 |
| 图 20 DP=45.14 mm时耦合禁带方向性对比 Fig. 20 Comparison of directional distribution for coupled band gaps with DP=45.14 mm |
| 图选项 |
 |
| 图 21 DP=50.78 mm时耦合禁带方向性对比 Fig. 21 Comparison of directional distribution for coupled band gaps with DP=50.78 mm |
| 图选项 |
2.3 不同尺寸压电片电学参数设计 由2.2节分析可知，随着压电片面积的增大，不改变电感时，耦合禁带会向低频移动。当与布拉格禁带发生耦合时会导致禁带带宽随着传播角度增加而迅速降低。因此为了避免这一情况需要针对压电片对电学参数进行设计。
根据耦合禁带形成机理，耦合禁带产生于未耦合2种波相交的频率附近[23]。因此设计原则为：通过选取电学参数，使得未耦合电波具有相同的频散曲线。
考虑到压电片相当于电容，与纯电感电路相连后，可以得到耦合前电波的理论波动曲线：
 | (21) |
式中: ω为系统角频率；Lall为基板长度；k为波数。值得注意的是，此时应当选取压电片开路时的电容
 | (22) |
由2.2节可知，LP=35 mm，Le=0.1 H的参数下耦合禁带效果较好，因此选为基准。同时由于正方形和圆形压电片对耦合禁带中心频率影响较小，故同一面积的压电片的电学参数保持一致。由图 4可知，电波具有永久布拉格禁带，该禁带满足波数Real(k)=π/Lall。代入式(22)，可得禁带初始频率有
计算不同面积的开路电容后，得到不同尺寸对应的电感值如表 3所示。
表 3 压电片结构参数 Table 3 Structural parameters of piezoelectric patches
| 情况编号 | 长度/mm | 直径/mm | 电感/H | 覆盖率/% |
| 1 | 20 | 22.567 6 | 0.319 | 6.25 |
| 2 | 35 | 39.493 3 | 0.1 | 19.14 |
| 3 | 40 | 45.135 2 | 0.078 | 25.00 |
| 4 | 45 | 50.777 1 | 0.060 6 | 31.64 |
| 5 | 50 | 56.419 0 | 0.051 2 | 39.06 |
| 注：表中覆盖率为压电片面积与基板总面积之比。 | ||||
表选项
由表 3可知，当压电片面积增大时，在同一位置产生耦合禁带所需的电感值迅速下降，对比情况1和情况3可知，当面积增大为原来的4倍时，电感变为24.45%。因此在难以实现大电感的情况时，可以通过增加压电片的面积使得构造的耦合禁带满足预期效果。
画出情况1的方向性分布如图 22所示。由图 22可知，方形压电片性能依然略优于圆形压电片，但这2种形状对耦合禁带的方向性影响较小。
 |
| 图 22 情况1时耦合禁带方向性对比 Fig. 22 Comparison of directional distribution for coupled band gaps with configuration 1 |
| 图选项 |
同时对比图 18发现，在同一尺寸下，降低电感会使得耦合禁带向低频移动，且带宽随角度增加而降低，在靠近π/4时，原来不存在耦合禁带的角度下出现了禁带分布，即该电学参数下禁带的方向性更好。
由于选取了图 19的结构和电学参数作为对照，且表 3中情况3~情况5耦合禁带方向性分布类似，因此为了节约篇幅，仅展示情况4的结果，如图 23所示。
 |
| 图 23 情况4时耦合禁带方向性对比 Fig. 23 Comparison of directional distribution for coupled band gaps with configuration 4 |
| 图选项 |
由图 23可知，虽然设计的电学参数可以使得未耦合电波频散曲线一致，但由于各尺寸耦合刚度不同，使得随尺寸增加，禁带起始频率近似不变，中心频率往高频移动，使得耦合禁带整体变宽。
针对方形压电片，选取3个角度，可以画出禁带和终止频率随尺寸增大的影响，如图 24所示，图中横坐标为情况编号，纵坐标为频率。
 |
| 图 24 耦合禁带边界频率对比 Fig. 24 Comparison of boundary frequencies for coupled band gaps |
| 图选项 |
对比情况4和情况5可知，当压电片覆盖率从25%提升到39.06%时，耦合禁带起始频率近似不变，终止频率略有上升，可见其性能增加并不明显。从而在根据实际情况设计压电片时，并不是覆盖率越高越好，需要从实际情况出发，综合考虑覆盖面积和目标性能的影响。
3 耦合禁带强迫响应仿真验证 为了验证耦合禁带的减振效果，选用耦合禁带带宽较窄，衰减常数较低的20 mm正方形压电片建立机电系统，设计验证模型如图 25所示。该模型参考了文献[40-41]中的禁带强迫响应验证的思路，建立含有均质部分和周期部分的模型，相互对照以体现禁带对弹性波的控制作用。其中模型右半部分为压电网络板，相邻压电片施加电感Le=0.319 H以构成电路网络，左半部分为环氧树脂基板作为对照。考虑到其方向性分布如图 22所示，因此得到特定频率下该结构的动力学响应后，与该频率下耦合禁带的方向性分布对比，即可验证耦合禁带的方向性分布。同时与均质板部分对比，可以研究有限周期结构内耦合禁带的减振性能。
模型采用了四边自由边界条件，在中心施加沿板厚度方向的瞬时力，即图中F点。为了降低频散效应，激励信号选为施加汉宁窗的正弦信号。其中心频率设定为1 140 Hz，信号长度为3个周期，最大振幅为0.5 N，采样频率为中心频率的30倍，以保证波包的信号质量。初始位移、速度和加速度都是0。结构中不含任何形式的阻尼，强迫响应利用ANSYS17中的完全法进行求解，从而保证结果的准确性。
为了监测各个方向的传播特性，选取距离激励点6个元胞长度(即0.48 m)圆弧上的单元能量作为指标，考虑到结构的对称性，只选取了上半部分的单元，如图 25所示。对于压电网络板部分，传播角θ=0即为沿着x轴正方向传播，对于均质板部分，传播角θ=0即为沿着x轴负方向传播。计算结果如图 26所示，图 26(a)为弹性波传播角度为44°和89°时，动能随时间变化的分布。由于1 140 Hz位于耦合禁带内，导致压电网络板内不同角度动能随时间变化分布区分明显，1 140 Hz时位于θ=89°的耦合禁带内，位于θ=44°的通带内。因此可知禁带内的能量得到了明显控制。选取图 26(a)中虚线时刻(t=1.87×10-3 s)，研究监测点不同方向的能量分布可得图 26(b), 图中横坐标为弹性波传播的角度，纵坐标为动能分布。
 |
| 图 25 强迫响应计算模型 Fig. 25 Calculation model for forced response |
| 图选项 |
 |
| 图 26 1 140 Hz下强迫响应计算结果 Fig. 26 Forced response calculation results at 1 140 Hz |
| 图选项 |
图 26(b)阴影部分为存在耦合禁带的角度，可见在禁带内，相比于均质板部分，动能大幅下降，通带内二者相差不大，这验证了禁带对弹性波传播的控制能力。
为了进一步比较二者的差异，计算了均质板部分和压电板部分所有动能之和随时间变化的关系曲线，如图 27所示，图中横坐标为时间，纵坐标为动能分布。图 27(a)为情况1对应正方形压电片尺寸下模型的能量分布，可见虽然刚开始时二者相差不大，但随着时间增加，弹性波每经过一次元胞总动能都有所衰减。
 |
| 图 27 不同情况下强迫响应能量结果对比 Fig. 27 Comparison of forced response energy results in different situations |
| 图选项 |
图 27(b)为情况3对应正方形压电片尺寸下模型的能量分布，可见在耦合禁带频率内(1 300 Hz)依然计算强迫响应, 得到能量对比如图 27所示，可见随着压电片面积增大，能量转化效率提高，对振动控制效果越明显。
图 28进一步对比了这2种情况下，均质板部分和网络板部分同一时刻的能量场分布，图中横坐标为结构x方向的坐标，纵坐标为结构y方向的坐标。考虑到结构的对称性，仅画出上半部分的能量分布。图 28中由原点出发的2条直线间的角度范围为2.3节计算的耦合禁带分布，在0°~45°内，情况1耦合禁带的涵盖角度为0°~9.4°，情况3耦合禁带的涵盖角度为7.01°~26.7°。
 |
| 图 28 不同情况下能量云图对比 Fig. 28 Distribution of kinetic energy in different situations |
| 图选项 |
可见禁带内都没有能量分布。另外可以看出，在均质板部分，由于没有禁带，因此能量沿四周均匀传播。在网络板部分，由于在特定角度范围存在禁带，导致能量主要集中在45°方向。
虽然结构内的电学参数并非最优值，导致情况1的能量对比差异并不十分明显，但也足以验证耦合禁带对有限周期结构的振动控制能力突出，在低频宽频振动控制领域具有广阔的应用前景。
4 结论 在压电周期结构的设计中，需要尽可能提高禁带的性能，此时结构参数的设计至关重要。开发了适用于压电周期网络系统的减缩波有限元算法，在此基础上，结合单波形占比的概念研究了二维结构内耦合禁带的方向性分布。同时分析了不同压电片形状和尺寸对耦合禁带的影响，并提出了如何针对特定频率设计结构的电学参数，分析结果由强迫响应验证。得到如下结论：
1) 利用该减缩波有限元算法可以节约90%以上的计算时间。但为了保证计算精度，至少需要保留5倍计算目标频带内的所有模态阶数，仅选择特定变形的模态会使得高频的结果出现较大误差。
2) 圆形压电片和方形压电片随着尺寸增大都使得耦合禁带向低频移动。相同面积下，二者具有相似的禁带中心频率和方向性分布，方形压电片禁带带宽略宽于圆形压电片，因此在后续设计中，可以把压电片面积作为主要设计指标。
3) 提出了针对不同尺寸和形状的压电片如何在同一频率附近设计耦合禁带的准则。随着压电片面积的增大，所需的电学参数降低，耦合禁带宽度增强，然而压电覆盖率超过一定程度后，对禁带的宽度提升并不明显。因此在无法实现大电感等电学元件的情况下，可以通过设计压电片的尺寸，使得耦合禁带涵盖目标频带。
4) 对于带宽较弱的耦合禁带，利用强迫响应分析可知，即使在不含任何阻尼的情况下，对有限结构内的振动依然有明显的控制效果。
参考文献
| [1] | HUSSEIN M I, LEAMY M J, RUZZENE M. Dynamics of phononic materials and structures: Historical origins, recent progress, and future outlook[J]. Applied Mechanics Reviews, 2014, 66(4): 1-38. |
| [2] | 张荣英, 姜根山, 王璋奇, 等. 声子晶体的研究进展及应用前景[J]. 声学技术, 2006, 25(1): 35-42. ZHANG R Y, JIANG G S, WANGE Z Q, et al. Progress in researches of phononic crystal and the application perspectives[J]. Technical Acoustics, 2006, 25(1): 35-42. DOI:10.3969/j.issn.1000-3630.2006.01.009 (in Chinese) |
| [3] | ROCA D, YAGO D, CANTE J, et al. Computational design of locally resonant acoustic metamaterials[J]. Computer Methods in Applied Mechanics and Engineering, 2019, 345(1): 161-182. |
| [4] | NING L, WANG Y Z, WANG Y S, et al. Active control of elastic metamaterials consisting of symmetric double Helmholtz resonator cavities[J]. International Journal of Mechanical Sciences, 2019, 153-154: 287-298. DOI:10.1016/j.ijmecsci.2019.02.007 |
| [5] | PELAT A, GALLOT T, GAUTIER F. On the control of the first Bragg band gap in periodic continuously corrugated beam for flexural vibration[J]. Journal of Sound and Vibration, 2019, 446: 249-262. DOI:10.1016/j.jsv.2019.01.029 |
| [6] | FAN Y, COLLET M, ICHCHOU M, et al. Energy flow prediction in built-up structures through a hybrid finite element/wave and finite element approach[J]. Mechanical Systems and Signal Processing, 2016, 66-67(1): 137-158. |
| [7] | QI S, OUDICH M, LI Y, et al. Acoustic energy harvesting based on a planar acoustic metamaterial[J]. Applied Physics Letters, 2016, 108(26): 1-4. |
| [8] | WU Y G, LI L, FAN Y, et al. Design of dry friction and piezoelectric hybrid ring dampers for integrally bladed disks based on complex nonlinear modes[J]. Computers and Structures, 2020, 233: 1-19. |
| [9] | WEN S R, XIONG Y H, HAO S M, et al. Enhanced band-gap properties of an acoustic metamaterial beam with periodically variable cross-sections[J]. International Journal of Mechanical Sciences, 2020, 166: 1-10. |
| [10] | MOSCATELLI M, ARDITO R, DRIEMEIER L, et al. Band-gap structure in two- and three-dimensional cellular locally resonant materials[J]. Journal of Sound and Vibration, Academic Press, 2019, 454(18): 73-84. |
| [11] | LI L, JIANG Z, FAN Y, et al. Creating the coupled band gaps in piezoelectric composite plates by interconnected electric impedance[J]. Materials, 2018, 11(9): 1656. DOI:10.3390/ma11091656 |
| [12] | 舒海生, 张法, 刘少刚, 等. 一种特殊的布拉格型声子晶体杆振动带隙研究[J]. 振动与冲击, 2014, 33(19): 147-151. SHU H S, ZHANG F, LIU S G, et al. Vibration band gap of a special rod of phononic crystals[J]. Journal of Vibration and Shock, 2014, 33(19): 147-151. (in Chinese) |
| [13] | 舒海生, 董立强, 李世丹, 等. 布拉格型声子晶体弦的横向振动带隙研究[J]. 人工晶体学报, 2013, 42(9): 1918-1923. SHU H S, DONG L Q, LI S D, et al. Study on the lateral band gap of string of phononic crystals[J]. Journal of Synthetic Crystals, 2013, 42(9): 1918-1923. DOI:10.3969/j.issn.1000-985X.2013.09.036 (in Chinese) |
| [14] | 张若军, 肖勇, 温激鸿, 等. 四边固支局域共振型板的低频隔声特性研究[J]. 振动工程学报, 2016, 29(5): 905-912. ZHANG R J, XIAO Y, WEN J H, et al. Analysis of sound transmission through clamped locally resonant plate in low frequency[J]. Journal of Vibration Engineering, 2016, 29(5): 905-912. (in Chinese) |
| [15] | LIU L, HUSSEIN M I. Wave motion in periodic flexural beams and characterization of the transition between bragg scattering and local resonance[J]. Journal of Applied Mechanics, 2012, 79(1): 1-17. |
| [16] | XIAO Y, WEN J, WANG G, et al. Theoretical and experimental study of locally resonant and bragg band gaps in flexural beams carrying periodic arrays of beam-like resonators[J]. Journal of Vibration and Acoustics, 2013, 135(4): 1-17. |
| [17] | DONG Y, YAO H, DU J, et al. Research on local resonance and Bragg scattering coexistence in phononic crystal[J]. Modern Physics Letters B, 2017, 31(11): 1-10. |
| [18] | CAI C, WANG Z, CHU Y, et al. The phononic band gaps of Bragg scattering and locally resonant pentamode metamaterials[J]. Journal of Physics D: Applied Physics, 2017, 50(41): 1-7. |
| [19] | LIU Z, ZHANG X, MAO Y, et al. Locally resonant sonic materials[J]. Science, 2000, 289(5485): 1734-1736. DOI:10.1126/science.289.5485.1734 |
| [20] | 朱兴一, 钟盛, 叶安珂, 等. 声子晶体禁带特性及局域共振现象的试验研究[J]. 人工晶体学报, 2014(11): 2852-2859. ZHU X Y, ZHONG S, YE A K, et al. Experimental study on band gap properties and local resonance phenomenon of phononic crystals[J]. Journal of Synthetic Crystals, 2014(11): 2852-2859. DOI:10.3969/j.issn.1000-985X.2014.11.016 (in Chinese) |
| [21] | LEE T, IIZUKA H. Bragg scattering based acoustic topological transition controlled by local resonance[J]. Physical Review B, 2019, 99(6): 1-11. |
| [22] | MANCONI E, MACE B. Veering and strong coupling effects in structural dynamics[J]. Journal of Vibration and Acoustics, 2017, 139(2): 1-10. |
| [23] | MACE B R, MANCONI E. Wave motion and dispersion phenomena: Veering, locking and strong coupling effects[J]. The Journal of the Acoustical Society of America, 2012, 131(2): 1015-1028. DOI:10.1121/1.3672647 |
| [24] | THOMPSON D J, FERGUSON N S, YOO J W, et al. Structural waveguide behaviour of a beam-plate system[J]. Journal of Sound and Vibration, 2008, 318(1): 206-226. |
| [25] | YU D, WEN J, ZHAO H, et al. Vibration reduction by using the idea of phononic crystals in a pipe-conveying fluid[J]. Journal of Sound and Vibration, 2008, 318(1): 193-205. |
| [26] | LIU J, LI L, FAN Y. A comparison between the friction and piezoelectric synchronized switch dampers for blisks[J]. Journal of Intelligent Material Systems and Structures, 2018, 29(12): 2693-2705. DOI:10.1177/1045389X18778360 |
| [27] | LIU J, LI L, HUANG X, et al. Dynamic characteristics of the blisk with synchronized switch damping based on negative capacitor[J]. Mechanical Systems and Signal Processing, 2017, 95: 425-445. DOI:10.1016/j.ymssp.2017.03.049 |
| [28] | YE G, YAN J, WONG Z J, et al. Optimisation of a piezoelectric system for energy harvesting from traffic vibrations[C]//Proceedings of the IEEE Ultrasonics Symposium. Piscataway: IEEE Press, 2009: 759-762. |
| [29] | WANG W, LI L, FAN Y, et al. Piezoelectric transducers for structural health monitoring of joint structures in cylinders: A wave-based design approach[J]. Sensors, 2020, 20(3): 1-25. DOI:10.1109/JSEN.2020.2964488 |
| [30] | LI L, JIANG Z, FAN Y, et al. Coupled band gaps in the piezoelectric composite plate with interconnected electric impedance[C]//ASME 2018 Conference on Smart Materials, Adaptive Structures and Intelligent Systems, 2018: 1-27. |
| [31] | MEITZLER A H, BERLINCOURT D, WELSH F S, et al. IEEE standard on piezoelectricity: ANSI/IEEE Std 176-1987[S]. Piscataway: IEEE Press, 1988. |
| [32] | CHEN S, WANG G, SONG Y. Low-frequency vibration isolation in sandwich plates by piezoelectric shunting arrays[J]. Smart Materials and Structures, 2016, 25(12): 1-8. |
| [33] | WANG G, CHEN S. Large low-frequency vibration attenuation induced by arrays of piezoelectric patches shunted with amplifier-resonator feedback circuits[J]. Smart Materials and Structures, 2015, 25(1): 1-15. |
| [34] | DUHAMEL D, MACE B R, BRENNAN M J. Finite element analysis of the vibrations of waveguides and periodic structures[J]. Journal of Sound and Vibration, 2006, 294(1): 205-220. |
| [35] | FAN Y, COLLET M, ICHCHOU M, et al. A wave-based design of semi-active piezoelectric composites for broadband vibration control[J]. Smart Materials and Structures, 2016, 25(5): 1-15. |
| [36] | ZHOU C W, LAINé J P, ICHCHOU M N, et al. Numerical and experimental investigation on broadband wave propagation features in perforated plates[J]. Mechanical Systems and Signal Processing, 2016, 75: 556-575. DOI:10.1016/j.ymssp.2015.12.006 |
| [37] | ZHOU C W, LAINé J P, ICHCHOU M N, et al. Multi-scale modelling for two-dimensional periodic structures using a combined mode/wave based approach[J]. Computers and Structures, 2015, 154: 145-162. DOI:10.1016/j.compstruc.2015.03.006 |
| [38] | BAMPTON M C C, CRAIG J R R. Coupling of substructures for dynamic analyses[J]. AIAA Journal, 1968, 6(7): 1313-1319. DOI:10.2514/3.4741 |
| [39] | WILCOX C H. Theory of Bloch waves[J]. Journal d'Analyse Mathématique, 1978, 33(1): 146-167. DOI:10.1007/BF02790171 |
| [40] | MING-HUI L U, ZHANG C, FENG L, et al. Negative birefraction of acoustic waves in a sonic crystal[J]. Nature Materials, 2007, 6(10): 744-748. DOI:10.1038/nmat1987 |
| [41] | YANG W P, WU L Y, CHEN L W. Refractive and focusing behaviours of tunable sonic crystals with dielectric elastomer cylindrical actuators[J]. Journal of Physics D Applied Physics, 2008, 41(13): 135408. DOI:10.1088/0022-3727/41/13/135408 |
