弹性支撑条件下一维声子晶体结构振动特性研究

何东泽， 史冬岩， 王青山， 马春龙,4

(1.哈尔滨工程大学 机电工程学院，哈尔滨 150001；2.中南大学 机电工程学院，长沙 410083；3. 中南大学 高性能复杂制造国家重点实验室，长沙 410083;4.哈尔滨职业技术学院 汽车学院，哈尔滨 150001)

声子晶体的物理概念提出是基于凝聚态物理领域内的光子晶体的研究基础，是指由两种或两种以上介质组成的周期性复合材料结构，并具有弹性波的带隙特性。当弹性波在声子晶体中传播时，在某些频率范围内不可传播，相应的频率范围称为带隙。相反，弹性波可以传播的频率范围称为通带[1]。声子晶体的分类是按照声子晶体在空间中的维数进行划分，通常分为一维、二维和三维声子晶体。一维声子晶体结构较为简单，在减振降噪应用中较为广泛，许多装备可以将分析模型简化为一维结构，比如火箭、船舶等，是机械工程、声学、力学领域中较为重要的分析模型[2]。

目前，对于声子晶体振动特性研究的数值方法主要分为传递矩阵法、平面波展开法以及集中质量法。传递矩阵法在一维声子晶体振动特性的计算中应用较为广泛，通过有限个传递矩阵的相乘可以得到相应的频率响应，具有较小的计算量，但是在处理二维以及三维声子晶体振动特性问题较为复杂[3-4]；平面波展开法是声子晶体振动特性研究中较为常见的一种计算方法，可以对各种维度的声子晶体进行计算。同时，在处理结构较为复杂的声子晶体结构时较为困难，具有一定的局限性[5-6]；集中质量法是基于离散化的思想，将连续介质中的质量密度进行集合，转化到有限的节点上，通过将连续系统进行离散化处理，具有较强的收敛性。同时，该方法可以开展一维、二维、三维的声子晶体振动特性研究，具有较广的应用范围[7]。

对于一维声子晶体结构振动特性的研究，目前存在较多的科研成果：Witarto等[8]采用方差分解对一维声子晶体带隙进行了整体灵敏度研究，分析输入参数对一维声子晶体中第一频带间隙的影响特性；Nogaty等[9]开展外界温度对一维压电声子晶体可调谐振频率的影响，研究发现温度对其带隙以及局部共振频率存在着一定的影响；刘扭扭等[10]设一种具有低频共振带隙的细直梁声子晶体结构，采用传递矩阵法对其纵向波振动特性进行对比，并于有限元法进行比较，验证计算的正确性；舒海生等[11]设计一种新型的布拉格型声子晶体T型杆，采用传递矩阵法对其面内和面外振动带隙特性进行研究，结果证明该结构具有良好的减振特性。

目前对于一维声子晶体振动特性的研究模型主要集中于理想条件下，对于弹性支撑条件下一维声子晶体振动特性的研究较少。本文采用回传射线矩阵法(method of reverberation-ray matrix, MRRM)开展弹性支撑条件下一维声子晶体结构的振动特性研究。因其引入双坐标系统，令结构总体散射矩阵变为常数矩阵，避免数值奇异问题以及超越函数问题，对一维声子晶体结构振动特性的研究具有良好的适用性。

本文基于欧拉-伯努利梁理论，通过弯曲波控制微分方程，采用回传射线矩阵法对弹性支撑下一维声子晶体结构进行模型建立，研究并分析弯曲波在弹性支撑条件下一维声子晶体结构中的传输特性。分别开展几何参数以及多种条件下弹性支撑刚度对一维声子晶体结构振动特性的研究，比较分析频率响应函数曲线的变化规律，对带隙的起始频率、截止频率以及带隙宽度进行分析，总结几何参数以及多种条件下弹性支撑刚度的影响，为一维声子晶体的研究提供数值方面的理论数据。

1 模型的建立

本文所建立的弹性支撑条件下一维声子晶体结构，如图1所示。将两种材料属性以及几何尺寸不同的单元梁结构沿着轴向方向均匀分布，进而得到本文所建立的声子晶体结构。图1中：L1和L2分别为由材料A和材料B构成的单元梁的长度;E1,A1以及ρ1分别为材料A的弹性模量、横截面面积以及质量密度；E2,A2以及ρ2分别为材料B的材料参数以及尺寸参数。

图1 弹性支撑下一维声子晶体结构模型Fig.1 The model of one-dimensional phononic crystal structure with elastic supporting

图1中：K为弹性支撑底座的刚度参数；支反力R(x,t)可以表示为

R(x,t)=Kv(x,t)

(1)

式中，v(x,t)为弹性支撑下一维声子晶体结构的弯曲位移。

2 控制方程

根据欧拉伯—努利梁理论，弹性支撑条件下一维声子晶体结构的控制方程可以表示为[12]

(2)

式中：E为弹性模量；I为惯性系数；ρ为质量密度。根据傅里叶变换公式，将一维声子晶体结构的控制方程进行频域范围内的转化，表示为

(3)

式中：ω为频率参数；v(x,ω) 为频域范围内的弯曲位移。根据上述频域范围内的控制方程，可以得到一维声子晶体结构的弯曲位移齐次解，具体可以表示为

v(x;ω)=a1eik1x+a2eik2x+d1e-ik1x+d2e-ik2x

(4)

式中：a1,2和d1,2分别为入射波参数与出射波参数；k1,k2为一维声子晶体结构中的轴向波数，具体表达式为

(5)

相应的，弹性支撑条件下一维声子晶体结构的其余物理量与轴向位移的微分关系可以表示为

(6)

式中：φ(x,ω) 为旋转位移；M(x,ω) 为弯矩；Q(x,ω) 为剪切力。根据式(6)中表示的弯曲位移与其余变量之间的关系，可以得到旋转位移、弯矩以及剪切力的简谐波表达形式，为

φ(x,ω)=α1(a1eik1x-d1e-ik1x)+α2(a2eik1x-d2e-ik2x)，

Q(x,ω)=β1(a1eik1x-d1e-ik1x)+β2(a2eik2x-d2e-ik2x)，

M(x,ω)=γ1(a1eik1x+d1e-ik1x)+γ2(a2eik2x+d2e-ik2x)

(7)

式中：α1,2为旋转位移影响参数；β1,2为剪切力影响参数；γ1,2为弯曲影响参数，定义为

(8)

根据弯曲位移，旋转位移，剪切力以及弯矩的简谐波表达形式，对位移矩阵δ以及力矩阵f进行定义，具体表达为

δ=Aδa+Dδd，
f=Afa+Dfd

(9)

式中：δ为位移向量，δ=[v(x,ω),w(x,ω)]T;f为力向量,f=[Q(x,ω),M(x,ω)]T；a=[a1,a2]T,d=[d1,d2]T分别为入射波向量以及出射波向量；Aδ,Dδ分别为入射波、出射波位移矩阵；Af,Df分别为入射波、出射波力矩阵。具体表达式为

(10)

3 回传射线矩阵法

3.1 单元划分以及节点定义

基于回传射线矩阵法单元划分的思想，对弹性支撑条件下一维声子晶体结构振动特性进行研究时，需要对整体结构进行单元划分以及节点定义[13]。一维声子晶体结构是由不同种类的单元梁在轴向方向上进行顺序排列，这与回传射线矩阵法的计算思想相吻合。因此，对弹性支撑条件下一维声子晶体结构进行单元划分以及节点定义，如图2所示。其中：N为梁单元的总体数量；n为节点总体数量。经过单元划分以及节点定义之后，采用回传射线矩阵法对一维声子晶体结构进行频率响应函数计算。

图2 单元划分以及节点定义Fig.2 Unit division and node definition

3.2 相位关系

通过上述对弹性支撑条件下一维声子晶体结构进行单元划分以及节点定义之后，对每个单元梁单元进行局部坐标的建立以便对回传射线传播方式进行描述，单元梁的局部坐标如图3所示。

图3 单元梁局部坐标示意图Fig.3 The local coordinates in the unit beam

从图3中可以看出，分别在单元梁k的节点i和j处建立了局部坐标(xij,yij)与(xji,yji)。其中：xij为由节点i向节点j方向映射；xji为由节点j向节点i映射。通过对弹性支撑条件下一维声子晶体结构局部坐标的建立，为后续的单元节点位移协调条件以及力平衡关系建立基础。同时，梁单元k的位移与力关系可以表示为

vij(xij)=-vji(lji-xji)，Qij(xij)=Qji(lji-xji)，

φij(xij)=φji(lji-xji)，Mij(xij)=-Mji(lji-xji)

(11)

根据梁单元k的位移变量与力变量之间的关系，可以得到单元梁内的相位关系表达式，为

[vij,φij]T=[-vji,φji]T，δij=Tδδji，
[Qij,Mij]T=[Qji,-Mji]T，fij=Tffji

(12)

式中：Tδ和Tf分别为相位关系中力向量与位移向量的转换矩阵，具体可以表示为

(13)

通过上述分析，得到了任意梁单元k中位移变量与力变量之间的转换关系。下一步，将对相邻梁单元k-1和梁单元k之间节点i处的力变量以及位移变量之间的关系进行表述，具体表示图4所示。

图4 相邻单元位移与力变量示意图Fig.4 The displacement and force variables of adjacent unit beams

图4表示相邻梁单元k-1与k之间节点i的位移变量以及力变量关系，相应的表达式为

(14)

式中：δij和fij分别为梁单元ij的位移向量以及力向量；δi(i-1)和fi(i-1)分别为梁单元i(i-1)的位移向量以及力向量。根据式(9)以及式(12)，可以得到任意梁单元k内入射波参数以及出射波参数之间的表达式，为

(15)

aij=Pijdji
aji=Pijdij→ak=PKdK′

(16)

a=Pd′=PUd

(17)

式中：P为整体相位矩阵；a，d′，d分别为整体入射波、出射波向量以及修正后的出射波向量，具体表达形式为

a=[a12,a21,…,aij,aji,…,a(n-1)n,an(n-1)]T，

d′=[d21,d12,…,dji,dij,…,dn(n-1),d(n-1)n]T，

d=[d12,d21,…,dij,dji,…,d(n-1)n,dn(n-1)]T

(18)

为了保证整体出射波向量与整体入射波向量中的元素保持相同的排列顺序，采用整体转换矩阵U=diag[U1,U2,…,Uk,…,UN]T进行顺序的重新排列组合，Uk=anti-diag[U0,U0]为局部单元转换矩阵；U0=diag[1,1]为基础单元转换矩阵。根据一维声子晶体结构总体单元划分以及节点定义，可以得到一维声子晶体结构的总体相位矩阵为

P=diag[P1,P2,…,Pk,…,PN-1,PN]4N×4N

(19)

3.3 散射关系

根据式(14)中的节点i处的位移变量以及力变量之间的关系式，可以得到节点i处的散射关系表达式，为

Aiai+Didi=0

(20)

式中：ai={ai(i-1),aij}T和di={di(i-1),dij}T为节点i处的波传播向量；Ai和Di为对应的散射参数矩阵，为

(21)

di=Siai,i=2,…,n-1

(22)

根据广义位移协调条件以及广义力平衡条件，结合外力源向量，可以得到弹性支撑条件下一维声子晶体结构节点1和节点n的散射关系表达式，为

A1a1+D1d1=q1,
Anan+Dndn=qn

(23)

式中：A1,n和D1,n分别为两端节点的散射参数矩阵；a1,n和d1,n分别为两端节点1和节点n的入射波向量以及出射波向量；q1,n为节点1和节点n的外力向量。通过式(23)中节点1和节点n的散射关系表达式进行化简，具体可以表示为

d1=S1a1+s1,
dn=Snan+sn

(24)

式中：S1和Sn为节点1和n处的散射矩阵；s1和sn分别为节点1和n的外力源向量。根据弹性支撑条件下一维声子晶体结构的节点定义，对各个节点的散射矩阵按照节点顺序进行组装，可以得到一维声子晶体结构的总体散射关系式，为

d=Sa+s

(25)

式中：S为一维声子晶体结构的总体散射矩阵；s为总体外力源向量，具体可以表示为

S=diag[S1,S2,…,Sn-1,Sn],
s=[s1,0,…,0,sn]T

(26)

根据式(17)以及式(25)中弹性支撑条件下的一维声子晶体的整体散射矩阵以及整体相位矩阵，对总体入射波向量a进行消除，可以得到整体出射波向量d的表达式，为

d=(I-R)-1s

(27)

式中，R=SPU为一维声子晶体结构的回传射线矩阵，由总体散射矩阵，总体相位矩阵以及转换矩阵组成。

3.4 振动特性与频响函数

当弹性支撑条件下一维声子晶体结构受到外界周期载荷的激励，需要将外界激励转化为稳态响应的形式[14]。同时，需要将力矩阵以及位移向量转化为稳态响应的表达形式，为

δ(x,t,ω)=δ(x,t)eiωt,
f(x,t,ω)=f(x,t)eiωt

(28)

将相应的频率代入力向量与位移向量的表达式，可以得出时域范围内的广义力与广义位移向量，具体可以表示为

f(x,t,ω)=[AFOPU+DFO](I-R)-1seiωt，
δ(x,t,ω)=[AVOPU+DVO](I-R)-1seiωt

(29)

式中：AFO与DFO为整体结构的力向量传播矩阵；AVO与DVO为整体结构的位移向量传播矩阵。进而可以求出频率范围内的稳态响应，得出对应的频率响应函数曲线。

4 数值计算与特性研究

本文主要针对弹性支撑下一维声子晶体结构振动特性进行研究。通过对一维声子晶体结构端点施加位移载荷，求出频率响应函数曲线并于有限元结果进行对比，验证本文计算方法的正确性。其次，开展相应的参数化研究，探寻各个参数对弹性支撑下一维声子晶体结构振动特性的影响规律，为一维声子晶体振动特性的研究拓展新的途径。本文构成一维声子晶体结构的材料名称以及属性，如表1所示。

表1 材料名称及参数

本文节采用有限元分析法中通过应用较为广泛的仿真模拟软件ANSYS Workbench进行分析，结合谐响应分析模块进行仿真计算，一维结构的位移激励幅值设置为1 mm; 采用Full算法进行网格自由划分。

本文中，采用由有机玻璃与铝组成的三周期一维声子晶体结构作为数值计算分析模型。其中，横截面尺寸设置为1 mm×1 mm，各个单元梁结构的长度设置为L1=L2=50 mm，晶格常数L=L1+L2=100 mm。图5表示弹性单元支撑条件下一维声子晶体结构模型，不同单元梁结构A和B分别设置不同的单元弹性支撑条件，单元支撑刚度分别设置为KA和KB。

图5 弹性单元支撑条件下一维声子晶体结构模型Fig.5 The model of one-dimensional phononic crystal structure with elastic unit supporting

4.1 理论计算与有限元仿真对比

本节对弹性支撑条件下一维声子晶体结构振动特性进行研究。通过由回传射线矩阵法所计算得到的频率响应曲线与有限元法进行对比，如图6所示。可以看出，二者吻合情况较为良好。因此，本文所提出的方法对于弹性支撑条件下一维声子晶体结构的频率响应计算较为准确。其中，弹性支撑参数KA=KB=1 000 Pa。

图6 数值计算与有限元仿真对比Fig.6 The comparison of the results by numerical calculation and finite element analysis

4.2 几何参数的影响

为了研究几何参数对一维声子晶体结构振动特性的影响规律，分别开展横截面宽度、厚宽比以及晶格长度的影响规律研究。通过改变参数的数值大小，根据所得到的频率响应曲线，对带隙起始频率以及截止频率进行数值提取并进行绘制。在本节中，一维声子晶体结构承受线性弹性支撑条件，其长度参数以及材料构成参数与4.1节中算例相同。

4.2.1 横截面宽度以及厚宽比的影响

如图7(a)所示，从图7(a)中可以看出，随着横截面宽度的逐渐增大，第一带隙以及第二带隙的起始频率以及截止频率呈现增长趋势。同时，第一带隙与第二带隙对应的带隙宽度同样随着横截面宽度的增加而增大。图7(b)表示厚宽比对弹性支撑条件下一维声子晶体结构第一带隙以及第二带隙的影响规律。可以发现，随着厚宽比的增大，带隙起始频率、截止频率以及带隙宽度均存在增长趋势。

图7 横截面宽度及厚宽比的影响Fig.7 The effect of section width and thickness to width radius on the first and second band gap

4.2.2 晶格长度的影响

图8表示不同晶格长度所对应的频率响应函数曲线前四带隙所对应的带隙起始频率以及截止频率的变化情况。可以看出，随着晶格长度的不断变化，前四带隙的起始频率以及截止频率均存在减小的趋势。同时，第一带隙的起始频率以及截止频率衰减程度较小，第四带隙的较大。对于带隙宽度而言，不难发现，随着晶格长度的逐渐增大，带隙宽度逐渐减小，但是对于第四带隙而言，带隙宽度的缩减程度较小。

图8 晶格长度对前四带隙的影响Fig.8 The effect of lattice constant on top four band gaps

4.3 弹性支撑刚度的影响

4.3.1 单元支撑刚度的影响

为了比较分析单元支撑刚度对一维声子晶体结构振动特性的影响，分别设置不同位置下的单元弹性支撑条件，分析各个情况下弹性支撑刚度对带隙变化的影响规律。图9表示仅单元梁A承受单元支撑条件下，单元支撑刚度KA对前四带隙的起始频率以及截止频率的变化情况。其中，单元B刚度设置为0，即KB=0，表示仅单元梁A承受弹性支撑。可以看出，随着弹性支撑刚度的增加，第一、第二带隙的起始、截止频率均存在增长趋势，但是第一带隙宽度逐渐较小。对于第三、四带隙而言，起始频率、截止频率以及带隙宽度变化程度较小，基本稳定。

图9 单元支撑刚度KA对前四带隙的影响(KB=0)Fig.9 The effect of the elastic unit supporting stiffness KA on top four band gaps(KB=0)

相应的，图10表示仅梁单元B承受弹簧支撑条件下(KA=0)，单元支撑刚度KB对前四带隙的起始频率以及截止频率的变化情况。其中，第一、第二带隙起始频率以及截至频率增长明显，增长趋势相似。同时带隙宽度基本保持不变；第三、第四带隙基本保持不变。

图10 单元支撑刚度KB对前四带隙的影响(KA=0)Fig.10 The effect of the elastic unit supporting stiffness KB on top four band gaps(KA=0)

通过上述分析可以发现，不同位置的单元弹性支撑条件对一维声子晶体结构第一带隙的影响较为明显，对第三、第四带隙影响程度较小。

4.3.2 线性弹性支撑刚度的影响

为了研究线性弹性支撑刚度对一维声子晶体结构振动特性的影响规律，设置单元支撑刚度KA=KB=K0，即梁单元A和B分别承受相同弹性支撑。通过改变K0的数值变化，比较分析线性弹性支撑刚度对一维声子晶体结构振动特性的影响。

通过图11中的带隙特性变化规律可以看出，线性弹性支撑刚度的增加使第一、第二带隙的起始频率以及截止频率逐渐增大，但是，第二带隙所对应的带隙起始、截止频率的增长趋势较第一带隙较小。对于第三、第四带隙而言，弹性支撑刚度的影响程度较小，带隙的起始、截止频率以及带隙宽度基本保持不变。

图11 线性弹性支撑刚度对前四带隙的影响Fig.11 The effect of linear stiffness constant on top four band gaps

4.3.3 非线性弹性支撑刚度的影响

为了开展非线性弹性支撑刚度对一维声子晶体结构振动特性的影响，分别开展非线性单元支撑刚度对一维声子晶体结构振动特性的影响，对前四带隙的起始、截止频率进行提取并绘制，总结非线性弹性支撑刚度的影响规律。

首先，设置单元支撑刚度KA=K0，KB=0.5K0，即梁单元A和B分别承受弹性支撑，但刚度大小不同。通过改变K0的数值变化，比较分析非线性单元支撑刚度对一维声子晶体结构振动特性的影响。如图12所示，随着单元支撑刚度K0的逐渐增大，第一、二带隙起始、截止频率均逐渐增大，但是第一带隙增加程度较为明显。同时，第一带隙的宽度逐渐减小。第二、第三、第四带隙宽度基本保持不变。

图12 单元支撑刚度K0对前四带隙的影响(KA=K0, KB=0.5K0)Fig.12 The effect of the elastic unit supporting stiffness K0on top four band gaps(KA=K0, KB=0.5K0)

相应的，设置梁单元A和B的单元支撑刚度为KA=0.5K0,KB=K0，对前四带隙的变化情况进行比较分析。通过图13可以看出，随着弹性支撑刚度K0的逐渐增大，第一、第二带隙起始频率以及截至频率逐渐增大，并且第一带隙增长趋势较为明显。对于带隙宽度而言，前四带隙的带隙宽度基本保持不变。

图13 单元支撑刚度K0对前四带隙的影响(KA=0.5K0, KB=K0)Fig.13 The effect of the elastic unit supporting stiffness K0on top four band gaps(KA=0.5K0, KB=K0)

通过上述分析可以看出，非线性单元支撑条件对一维声子晶体结构第一、第二带隙影响程度较为明显，对于第三、第四带隙在带隙起始频率、截止频率以及带隙宽度方面影响程度较小。

5 结 论

本文基于欧拉—伯努利理论，结合回传射线矩阵法对弯曲波在弹性支撑条件下一维声子晶体结构进行分析模型的建立，计算其频率响应函数曲线，对振动特性进行分析研究。同时，开展几何参数与弹性支撑条件对弹性支撑条件下一维声子晶体结构振动特性的影响规律研究，分析其变化规律，得到相应的结论。

(1)为了验证本文计算方法的正确性，采用有限元法对弹性支撑条件下一维声子晶体结构频率响应函数曲线进行计算，并于回传射线矩阵法进行对比，验证本文提出方法计算的正确性，为后续的研究奠定基础。

(2)通过对弹性支撑条件下一维声子晶体结构振动特性的计算，比较分析不同几何参数下频率响应函数曲线。可以发现，横截面宽度、厚宽比以及晶格长度对一维声子晶体结构的带隙影响较为明显，存在不同的影响规律。

(3)通过对多种条件下弹性刚度对一维声子晶体结构的频率响应函数曲线进行计算，比较分析不同条件下弹性刚度对带隙的影响规律。可以发现，各种条件下弹性支撑刚度对第一带隙的影响程度较为明显，对第三、第四带隙的影响程度较小。