基于COMSOL局域共振聲子晶體薄板振動帶隙研究

王翌偉，徐曉美，林 萍

（南京林業(yè)大學(xué) 汽車與交通工程學(xué)院，南京 210037）

局域共振聲子晶體板是一種周期結(jié)構(gòu)功能材料，具有負等效參數(shù)[1]，通過調(diào)節(jié)單胞結(jié)構(gòu)參數(shù)可獲得低頻段內(nèi)的禁帶[2]。在禁帶內(nèi)彈性波的激勵下，聲子晶體板的散射體間產(chǎn)生共振，并與彈性波相互作用，從而抑制彈性波的傳播。

近些年來，局域共振聲子晶體板已成為振動控制領(lǐng)域的研究熱點。Oudich等[3]研究了貼附一層或雙層柱的板結(jié)構(gòu)的帶隙機制，通過調(diào)節(jié)結(jié)構(gòu)參數(shù)，可以改變局域模態(tài)與板中Lamb模態(tài)之間的耦合方式，從而獲得低頻帶隙。Jin等[4]提出了具有空心柱的聲子晶體板，通過調(diào)整空心柱內(nèi)部的直徑，能使布拉格帶隙和共振帶隙相互作用，從而拓寬了帶隙。廖濤等[5]將壓電材料加入聲子晶體板，基于有限元法證明了該結(jié)構(gòu)的可調(diào)帶隙特性?？堤P等[6]也基于有限元法對所提出的聲子晶體板帶隙特性及其影響因素進行了研究，證明了所提結(jié)構(gòu)具有優(yōu)越的低頻寬帶隙特性，并提出拓寬該結(jié)構(gòu)低頻帶隙的方法。孫向洋等[7]基于有限元法分析了經(jīng)典單面柱聲子晶體的帶隙特性及影響因素，并依此進行了結(jié)構(gòu)優(yōu)化。Zhu等[8]利用有限元法研究了一種新型二維聲子晶體板，該結(jié)構(gòu)能通過調(diào)節(jié)預(yù)應(yīng)變實現(xiàn)低頻帶隙的可調(diào)諧性。Zhang 等[9]針對車輛噪聲、振動與聲振粗糙度（Noise、Vibration and Harshness，NVH）問題，提出一種具有填充圓柱體的聲子晶體板，并通過色散模擬研究了其帶隙特性及影響因素。

鑒于此，本文在驗證COMSOL有限元模擬方法用于聲子晶體振動帶隙計算可靠性的基礎(chǔ)上，針對低頻振動控制問題，提出一種局域共振聲子晶體薄板結(jié)構(gòu)，應(yīng)用有限元模擬方法研究其帶隙結(jié)構(gòu)和元胞結(jié)構(gòu)參數(shù)對振動帶隙的影響，并考察該聲子晶體薄板內(nèi)部振動波的傳輸特性，以為類似結(jié)構(gòu)薄板的設(shè)計與應(yīng)用提供理論支撐。

1 彈性波方程及Bloch定理

在線彈性、各向同性、體積無限大的均勻介質(zhì)中，忽略阻尼的影響，彈性波動方程的矢量形式可寫作[10]：

式中：r為位置矢量，▽為哈密頓算子，u為質(zhì)點位移矢量，ρ為介質(zhì)的密度，μ和λ為介質(zhì)材料的拉梅常數(shù)，ω為特征圓頻率。

由于聲子晶體的周期性和對稱性，使得晶格中的場也具有一定的對稱性。根據(jù)Bloch 定理，平移周期性的線性系統(tǒng)的本征場可表示為[11]：

式中：k為Bloch波矢，其幅值函數(shù)具有與晶格相同的平移周期性，即：

結(jié)合式(3)和Fourier級數(shù)數(shù)學(xué)定理，式(2)表述本征場的波函數(shù)可以寫為：

式中：G為倒格矢。由于聲子晶體點群的對稱性，波矢k的取值范圍壓縮在第一不可約Brillouin 區(qū)內(nèi)。因此，只需沿著第一不可約Brillouin 區(qū)邊界上的值遍歷一次，就能求得聲子晶體的本征頻率。

2 帶隙計算有限元法

目前，COMSOL等有限元軟件具有設(shè)置Floquet周期性邊界的功能，本節(jié)基于彈性波方程和Bloch定理，通過算例闡述應(yīng)用COMSOL軟件計算聲子晶體振動帶隙的方法，并驗證其計算結(jié)果的可靠性。

根據(jù)Bloch 理論，對彈性波在整個晶格中傳播行為的研究可以轉(zhuǎn)換到單個元胞及其第一不可約Brillouin區(qū)中進行。元胞的廣義本征值方程為[12]：

式中：Ks為元胞剛度矩陣，Ms為元胞質(zhì)量矩陣，u(v)為位移特征向量。

要使元胞中每個點都滿足Bloch 定理，只要元胞邊界上的所有點滿足Bloch 邊界條件即可，即要滿足式(7)。

式中：a為晶格常數(shù)。

COMSOL 軟件中的Floquet 周期性邊界設(shè)置可約束周期性結(jié)構(gòu)的邊界位移，從而實現(xiàn)Bloch 邊界條件，其滿足的方程為：

式中：kF為波矢，在軟件中通過設(shè)置kx、ky和kz表示。usource和udestination分別表示Floquet 周期性邊界條件的源端和目標(biāo)端，rdestination-rsource的結(jié)果等于晶格常數(shù)。

在COMSOL軟件中，對于周期元胞與相鄰元胞的非接觸面，無須設(shè)置邊界條件；對于周期元胞與相鄰元胞的接觸面，則需在元胞邊界的對應(yīng)節(jié)點上設(shè)置Floquet周期性邊界條件。以二維方形元胞為例，其Brillouin 區(qū)如圖1 所示，陰影部分是第一不可約Brillouin區(qū)。首先，取x和y方向的兩條邊或邊界作為兩個源項，選擇相對的邊或邊界作為目標(biāo)項，即在元胞的x方向和y方向各應(yīng)用一次周期性邊界條件。然后，添加波矢k的設(shè)置，k∈[0，3]，0到1區(qū)間定義為M-Γ邊的波數(shù)，1到2區(qū)間定義為Γ-X邊的波數(shù)，2到3區(qū)間則定義為X-M邊的波數(shù)。為了與Floquet周期性邊界條件中波矢的設(shè)置相一致，將k分解為kx和ky表示。接著，對波矢設(shè)置參數(shù)化掃描，掃描范圍0～3，掃描步長0.01，如此波矢能沿著二維元胞的第一不可約Brillouin 區(qū)M-Γ-X-M的邊界遍歷一次。最后，選擇相應(yīng)的特征頻率求解器求解，即每給定一個k值，就能根據(jù)式(8)，求解式(6)所示的特征值方程。

圖1 二維方形元胞的Brillouin區(qū)

為驗證有限元法帶隙計算的可靠性，以文獻[13]所示的典型單面柱狀聲子晶體板為例進行帶隙計算。圖2 給出了該聲子晶體板的元胞結(jié)構(gòu)及其Floquet 周期性邊界設(shè)置，元胞的結(jié)構(gòu)參數(shù)為：晶格常數(shù)a=10 mm，橡膠板厚度d=2 mm，鋼柱半徑r=4.8 mm、高度h=10 mm，其色散曲線及帶隙如圖3所示，圖中陰影部分是結(jié)構(gòu)的第一條完全帶隙。雖然從圖3 可以看出，該算例中的元胞結(jié)構(gòu)具有較寬的低頻局域共振帶隙，但由于鋼柱密度和尺寸比橡膠均質(zhì)板大得多，工程實際中橡膠板很難提供鋼柱所需的支撐強度，故而此結(jié)構(gòu)在工程中很難實現(xiàn)應(yīng)用。

圖2 算例元胞結(jié)構(gòu)及其周期性邊界設(shè)置

圖3 算例的色散曲線及帶隙

表1 比較了有限元法和文獻[13]的數(shù)值方法計算得到的帶隙相關(guān)結(jié)果，從其誤差可見，運用有限元法計算帶隙是可靠的。

表1 帶隙計算結(jié)果對比

3 聲子晶體薄板帶隙特性

3.1 聲子晶體薄板結(jié)構(gòu)

本文所研究的局域共振聲子晶體薄板結(jié)構(gòu)如圖4 所示，圖5 是其晶格常數(shù)為a的元胞結(jié)構(gòu)，包括三層：上、下層是厚度分別為h1和h2的EVA 板，中間層是厚度為h0的橡膠板；上、下層內(nèi)嵌的橡膠圓柱半徑為r，高度與EVA層厚度相等。元胞初始尺寸參數(shù)選為：a=10 mm、h0=h1=h2=1 mm、r=4 mm，材料參數(shù)見表2。在該結(jié)構(gòu)中，EVA材料是基體部分，主要起結(jié)構(gòu)支撐作用，橡膠材料充當(dāng)散射體。由于未貼附金屬材料作散射體，所以該結(jié)構(gòu)的總體平均密度較小，并且避免了常規(guī)結(jié)構(gòu)中貼附的質(zhì)量塊在振動環(huán)境下易掉落的問題，使板的性能更加穩(wěn)定。此外，該聲子晶體薄板結(jié)構(gòu)僅由EVA材料和橡膠材料組成，板的柔韌性與彈性都較好，更易被應(yīng)用于復(fù)雜的振動環(huán)境中。

圖4 聲子晶體薄板結(jié)構(gòu)

圖5 聲子晶體薄板元胞結(jié)構(gòu)

表2 聲子晶體薄板材料參數(shù)

3.2 帶隙特性

建立上述聲子晶體薄板元胞的有限元模型，采用有限元法計算其振動帶隙，計算結(jié)果如圖6所示。

根據(jù)彈性波傳播理論，由圖6 可知，板中有Lamb 波和水平剪切波，其中，Lamb 波分為反對稱Lamb波（彎曲波）和對稱Lamb波。圖6中有三條能帶從G點出發(fā)，沿三個方向分開，其中，兩條能帶在G點附近呈線性上升，分別對應(yīng)著板中的水平剪切波和對稱Lamb波，由于在低頻范圍內(nèi)這兩種波的傳播特性與在均勻介質(zhì)中傳播特性相近，因此總體表現(xiàn)為線性。第一條能帶對應(yīng)著板中的彎曲波，從圖6可以看到，在第一條能帶和第二條能帶間沒有任何彎曲波穿過，形成了一條完全帶隙，即彎曲波帶隙，帶寬為88 Hz，帶寬頻率范圍為448 Hz～536 Hz。

圖6 振動帶隙結(jié)構(gòu)圖

為進一步說明板的振動帶隙形成機理，圖7 給出了圖6 中幾個關(guān)鍵點的振動模態(tài)圖，圖中用散布在結(jié)構(gòu)中的箭頭表示區(qū)域位移的相對大小和方向。第一條能帶是形成彎曲波帶隙起始頻率的關(guān)鍵，從圖7 中A 點的振動模態(tài)可以看出，元胞中間的橡膠圓柱以及與橡膠圓柱相連的橡膠薄板均向z軸正方向運動，而EVA板基本不動，即振動集中在中間，此時A 點對應(yīng)于結(jié)構(gòu)發(fā)生了局域共振，使得彎曲波與之發(fā)生強烈耦合，致使彎曲波的能量局限在橡膠材料中而不能繼續(xù)傳播，從而使第一條能帶被截斷，繼而開始出現(xiàn)帶隙。帶隙在B點所對應(yīng)的平直帶處截止，觀察B點的振動模態(tài)可知，橡膠圓柱和橡膠板依舊向z軸正方向運動，而EVA 基體板開始向與橡膠部分相反的方向運動，該運動模式標(biāo)志著局域共振的消失，此后，彎曲波便能夠在元胞內(nèi)通過。

圖7 A～F點的振動模態(tài)

從圖7還可以看出，在帶隙上方的頻段內(nèi)，能帶圖中還存有一些近似水平的平直帶，觀察相應(yīng)C、D、E、F 點的振動模態(tài)，同樣發(fā)現(xiàn)振動位移基本都在橡膠圓柱和橡膠板內(nèi)，同時EVA 基體基本不動，可見相應(yīng)的結(jié)構(gòu)均發(fā)生了局域共振。雖然這些局域共振模式不能像A 模式一樣形成完全帶隙，但是它們的存在有利于方向帶隙的產(chǎn)生，如在圖6 的C 點附近有一條Γ-X方向的彎曲波帶隙（圖中陰影）。

3.3 帶隙特性影響分析

為探究元胞彎曲波帶隙的調(diào)控規(guī)律，指導(dǎo)帶隙特性優(yōu)化，分析了結(jié)構(gòu)參數(shù)（包括圓柱半徑r，中間橡膠層厚度h0，上、下EVA層厚度h1、h2，晶格常數(shù)a）對振動帶隙的影響規(guī)律。圖8至圖12的左側(cè)縱坐標(biāo)對應(yīng)帶隙的起始頻率和截止頻率，右側(cè)縱坐標(biāo)對應(yīng)帶隙寬度。

圖8 圓柱半徑r對振動帶隙的影響

圖8 為元胞中橡膠圓柱半徑對振動帶隙的影響。由圖可知，隨著圓柱半徑r的增大，帶隙的起始頻率和截止頻率都向低頻區(qū)域大幅度移動，當(dāng)半徑r從2 mm 增至4.5 mm，起始頻率下降了近60 %；同時，帶隙的寬度從63 Hz拓寬到92 Hz。

圖9 為元胞中間橡膠層厚度對振動帶隙的影響。由圖可見，隨著橡膠層厚度增加，帶隙的起始頻率和截止頻率都向高頻區(qū)域移動，帶隙先變寬再變窄。

圖9 中間層厚度h0對振動帶隙的影響

圖10為元胞的上、下層EVA板厚度對振動帶隙的影響。由圖可見，隨著上、下層厚度增加，振動帶隙的起止頻率和截止頻率都向高頻區(qū)域移動，且?guī)峨S之變寬。

圖10 上、下層厚度h1、h2對振動帶隙的影響(h1=h2)

由于晶格常數(shù)a受限于圓柱半徑r，因此其對振動帶隙的影響可分兩種情況討論。

一種是圓柱半徑r不變，即仍為4 mm，考察振動帶隙隨晶格常數(shù)a的變化規(guī)律，如圖11 所示。從圖中可以看出，隨著晶格常數(shù)增大，帶隙起始頻率略有下降，截止頻率也隨之下降，同時帶隙寬度逐漸變小。

圖11 晶格常數(shù)a對振動帶隙的影響(r=4 mm)

另一種是圓柱半徑與晶格常數(shù)的比值r/a不變，即r/a=2/5，考察晶格常數(shù)a對振動帶隙的影響。通過對不同晶格常數(shù)a開展多組仿真試驗，得到能較好反映晶格常數(shù)對振動帶隙影響規(guī)律的對比數(shù)據(jù)，如圖12 所示，圖中晶格常數(shù)的變化間隔為5 mm。由圖可以看出，隨著晶格常數(shù)增大，帶隙的起始頻率和截止頻率都明顯向低頻區(qū)域移動，帶寬也隨之迅速減小。因此，總體來說，大的晶格常數(shù)會產(chǎn)生更低頻段和更窄的振動帶隙。

圖12 晶格常數(shù)a對振動帶隙的影響(r/a=2/5)

總之，對于此聲子晶體薄板，要使振動帶隙產(chǎn)生于低頻區(qū)域，可以增加圓柱半徑或增大晶格常數(shù)，但增大晶格常數(shù)會使帶隙寬度變窄。增加薄板中間層厚度或上、下層厚度都會使帶隙往高頻區(qū)域移動。因此，為使聲子晶體薄板能在低頻區(qū)域產(chǎn)生較寬的振動帶隙，對上述結(jié)構(gòu)參數(shù)進行了優(yōu)化組合研究，其中晶格常數(shù)a的3個參考值分別為10 mm、15 mm和20 mm，相應(yīng)的圓柱半徑r的參考值分別為4 mm、6 mm和8 mm，以維持r/a的比值為定比2/5；上、中、下三層薄板的厚度參考值分別為0.5 mm、1 mm 和1.5 mm，且使上、下層板的厚度相等，即h1=h2。

基于上述結(jié)構(gòu)參數(shù)，對枚舉法所得的27組參數(shù)組合方案開展了0到1 000 Hz范圍內(nèi)的振動帶隙仿真研究。以汽車上200 Hz～400 Hz 的中、低頻為目標(biāo)頻段，選擇帶隙寬度在目標(biāo)頻段內(nèi)占比最大的參數(shù)組合為最優(yōu)方案。由仿真結(jié)果可知，當(dāng)晶格常數(shù)為15 mm、圓柱半徑為6 mm、三層薄板厚度均為1.5 mm 時，所得局域共振聲子晶體薄板的振動帶隙在目標(biāo)頻段內(nèi)的占比最大（29%）。此時，振動帶隙的起始頻率為299 Hz、截止頻率為357 Hz。圖13即是通過枚舉優(yōu)選出的聲子晶體薄板的振動帶隙結(jié)構(gòu)圖。可見，通過結(jié)構(gòu)參數(shù)調(diào)整，已使聲子晶體薄板的振動帶隙移至低頻段，但帶隙寬度有所減小。

圖13 優(yōu)化后元胞的帶隙結(jié)構(gòu)圖

4 振動傳輸損失

為考察振動帶隙內(nèi)聲子晶體薄板對波傳播的阻礙作用，討論了頻域上波的傳播特性?；贑OMSOL軟件固體力學(xué)模塊的頻域研究，建立了有限周期結(jié)構(gòu)在頻域上的傳輸特性分析模型，如圖14所示。圖中，x方向和y方向均有20 個元胞，元胞的結(jié)構(gòu)參數(shù)為優(yōu)選后的尺寸，所構(gòu)成的有限周期結(jié)構(gòu)尺寸為300 mm×300 mm×4.5 mm。

如圖14所示，在薄板的左側(cè)邊界施加垂直于板面的位移激勵w0，在另一側(cè)邊界拾取位移響應(yīng)w1，利用積分算子計算左右側(cè)z方向上的位移幅值，并用傳輸損失表征波的衰減，傳輸損失（Transmission loss，TL）計算式為[14]：

圖14 有限周期薄板的載荷示意圖

式中：win和wout分別為激勵邊界和響應(yīng)邊界z方向上的位移幅值，A1和A2為對應(yīng)的邊界面積。

仿真計算的頻率范圍設(shè)定為10 Hz～1 000 Hz，計算間隔為10 Hz，計算得到的聲子晶體薄板在頻域上的傳輸特性如圖15所示。

從圖15可以看出，在300 Hz～350 Hz頻率范圍內(nèi)，傳輸損失明顯增大，最大衰減量為71.6 dB，出現(xiàn)在頻率為310 Hz 處，衰減頻率范圍與圖13 所示的299 Hz～357 Hz振動帶隙基本吻合。

圖15 聲子晶體薄板在頻域上的傳輸特性

為進一步呈現(xiàn)振動帶隙內(nèi)板對波傳播的阻礙作用，圖16給出了帶隙內(nèi)頻率為310 Hz和帶隙外頻率為500 Hz時聲子晶體薄板z方向的振動位移。

圖16 聲子晶體薄板的振動位移圖

由圖16可以看出，對于振動帶隙內(nèi)的310 Hz頻率，板左側(cè)的位移激勵產(chǎn)生的彎曲波被局限在左側(cè)的呈局域共振模式的元胞中，板右側(cè)的z向振動位移幾乎為0；對于振動帶隙外的500 Hz頻率，聲子晶體薄板呈現(xiàn)出明顯波動，即彎曲波能通過薄板傳播從板的左側(cè)到達板的右側(cè)?？梢?，在振動帶隙內(nèi)，聲子晶體薄板能很好地控制板內(nèi)的彎曲波在相應(yīng)頻段內(nèi)的傳播，從而控制薄板的振動。

5 結(jié)語

探討了應(yīng)用COMSOL 有限元模擬方法開展聲子晶體振動帶隙計算的可靠性，模擬計算了所設(shè)計的局域共振聲子晶體薄板的振動帶隙，分析了其帶隙結(jié)構(gòu)和元胞結(jié)構(gòu)參數(shù)對振動帶隙的影響，考察了聲子晶體薄板內(nèi)波的傳輸特性。主要研究結(jié)論如下：

（1）利用COMSOL 有限元模擬方法開展聲子晶體振動帶隙計算是可靠的，與數(shù)值計算方法相比，帶隙起始頻率、截止頻率和帶寬的計算誤差都很小。

（2）對于所設(shè)計的局域共振聲子晶體薄板，元胞的結(jié)構(gòu)參數(shù)對帶隙特性具有明顯的影響。增加薄板上、中、下層厚度，會使帶隙的起止頻率和截止頻率都向高頻區(qū)域移動；增加元胞的橡膠圓柱半徑和晶格常數(shù)，會使帶隙向低頻區(qū)域移動，但圓柱半徑增大使帶隙變寬，而晶格常數(shù)增大使帶隙變窄。

（3）基于枚舉優(yōu)化參數(shù)組合得到的局域共振聲子晶體薄板可在設(shè)計的目標(biāo)頻段內(nèi)形成振動帶隙，薄板內(nèi)波的傳輸特性和薄板的振動位移，進一步證實了聲子晶體薄板在帶隙內(nèi)對波傳播的阻礙作用。