彈性支撐條件下一維聲子晶體結構振動特性研究

何東澤， 史冬巖， 王青山， 馬春龍,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)通過對多種條件下彈性剛度對一維聲子晶體結構的頻率響應函數曲線進行計算，比較分析不同條件下彈性剛度對帶隙的影響規律。可以發現，各種條件下彈性支撐剛度對第一帶隙的影響程度較為明顯，對第三、第四帶隙的影響程度較小。