區間模型下聲子晶體的帶隙優化研究

劉 堅, 陳俊煌, 夏百戰, 滿先鋒

(湖南大學 汽車車身先進設計制造國家重點實驗室, 長沙 410082)

聲子晶體是一種新型功能材料，其組分的材料參數呈周期性變化。聲子晶體帶隙的形成機理通常有兩種，分別是Bragg散射機理[1]及局域共振機理[2]。Bragg散射形成的彈性波帶隙所對應的波長一般和晶格常數相當，其第一帶隙的中心頻率處所對應的波長一般為晶格常數的兩倍。Bragg散射的這一特性有礙于低頻減振降噪。Liu等[3]在2000年提出了局域共振聲子晶體概念。他將硅橡膠包裹的鉛球按簡單立方晶格排列方式，周期性地嵌入環氧樹脂基體中，形成三維三組元聲子晶體。這種聲子晶體產生的帶隙所對應的波長比晶格常數大兩個數量級，突破了Bragg散射的波長限制。Goffaux等[4]理論上證實二維三組元聲子晶體同樣存在局域共振帶隙，并結合類Fano現象和近似的機械振動模型初步揭示了局域共振帶隙的形成機理。Ho等[5]在實驗室制得了可用于低頻隔音的局域共振材料。Zhang等[6-7]發現在局域共振型聲子晶體的高頻段存在較寬的帶隙。Hirsekorn等[8-9]分析了三種不同材料組成的二維聲子晶體的傳播性質。

幾何結構參數對于聲子晶體聲學特性的影響較大。改變聲子晶體的結構參數，能改變聲子晶體帶隙寬度及帶隙范圍等。基于這一性質，優化設計逐漸被應用于聲子晶體優化設計。Romero-Garcia等[10]通過對孔的分布形狀及數量進行優化，使得聲子晶體中的聲波聚焦和稀疏性能得到改善。通過對共振腔形狀尺寸的優化，Wang等[11]得到了具有完全帶隙的三維孔狀聲子晶體。這些優化都是基于確定性物理模型。在實際工程應用中，由于制造誤差和測量誤差，以及多變的環境因素和不可預測的外部激勵，聲子晶體的不確定性是難以避免的。聲子晶體帶隙對這些不確定參數比較敏感，XIA等[12]在最近的一篇文章研究發現，霍姆赫茲共振腔聲子晶體的帶隙及有效體積模量對水溫的變化比較敏感。如果不考慮這些不確定性因素，會導致得到的優化結構不可靠。因此，在對聲子晶體進行優化時，需要考慮不確定性參數的影響。

概率模型是常用的不確定數值分析模型。然而在樣本數據有限的情況下，通常難以獲得不確定參數精確概率密度函數。針對這種情況，Moore[13]提出了區間模型的方法。Monte-Carlo法[14]是種最簡單最穩健的區間分析方法。夏日戰等[15]將區間不確定模型引入聲學超材料，采用Monte-Carlo 法分析區間不確定性對聲學超材料性能的影響；接著在此基礎上，構造區間模型下聲學超材料的可靠性優化模型，采用優化算法對優化模型進行求解。根據Monte-Carlo 法的概率收斂性，其計算精度隨著樣本數據的增加逐漸提高。但是巨額的計算負擔使Monte-Carlo 法難以適用于復雜的工程實際問題[16]。Liu等[17]利用Chebyshev展開法分析不確定非線性系統動力學響應分析。Xia等[18]采用Chebyshev展開法分析時變區間模型下結構的動態響應。Wu等[19-20]將Chebyshev展開法應用于區間模型下汽車動力學分析和多體機械系統的不確定分析。

目前，國內外對于聲子晶體的研究基本都是在確定幾何尺寸、材料屬性、環境因素的情況下進行，而鮮有關于聲子晶體的不確定性研究。針對這一情況，本文將區間分析方法引入聲子晶體，對其進行不確定分析。本文將兩個不確定性變量與兩個設計變量作為Chebyshev多項式的四個變量，構建每條能帶的Chebyshev代理模型，并基于該代理模型用Monte-Carlo 法求解能帶的變化區間。以聲子晶體帶隙最大化為目標函數，以期望帶隙為約束條件，構建基于Chebyshev代理模型的區間聲子晶體可靠性優化模型。最后，采用遺傳算法對該區間優化模型進行求解。該方法的優點是只需一組初始的有限元計算數據構建Chebyshev代理模型，后續優化皆基于計算效率較高Chebyshev代理模型，避免了傳統優化方法中有限元模型的重復計算，極大地減小了計算負擔。

1 能帶結構的有限元計算方法

對于各向同性的線彈性體材料，角頻率為ω的彈性波傳播方程為

(1)

(2)

(3)

圖1(a)為聲子晶體單胞。根據Bloch定理[21]，位移向量u(r)可以用如下公式表示

u(r)=ei(k·r)uk(r)

(4)

其中，uk(r)表示周期性矢量函數，k=(kx,ky)表示第一布里淵區的波矢量，如圖1(b)所示。

用有限元方法數值求解方程(4)，得到離散化的廣義特征值方程如下

Ku=ω2Mu

(5)

式中:K和M分別表示整體的剛度矩陣及質量矩陣;u表示節點位移。

剛度矩陣和質量矩陣的計算公式為

(6)

(7)

式中:B表示應變矩陣;N表示形函數矩陣;Ve表示整個單胞區域。

假設波入射方向上有N個聲子晶體單胞，邊界則滿足Born-von Karman條件[22]

f(r)=f(r+Na)

(8)

式中:a表示晶格基矢;N表示整數;f(r)表示在波矢方向的周期函數，其周期性為a。

(a) 聲子晶體單胞(b) 不可約布里淵區域

圖1 聲子晶體單胞有限元模型及其不可約布里淵區域

Fig.1 Finite element model of the phononic crystal unit cell and the corresponding irreducible Brillouin zone

將式(5)和(8)聯立，即可求解給定波矢k下的特征頻率。將求得的特征頻率代到控制方程(1)中，即可得到該頻率下的本證模態u(r)。將波矢k對單胞結構的不可約Brillouin區進行掃掠，即可得到該聲子晶體的能帶結構。

2 區間模型的構造及可靠性優化理論

減振降噪是聲子晶體的一個重要性能。局域共振聲子晶體能夠用小尺寸產生低頻帶隙，突破Bragg散射的波長限制，從而為聲子晶體在低頻減振降噪方面的應用提供了可能。本文將不確定區間模型引入局域共振聲子晶體，以帶隙最大化為目標函數，以期望頻帶為約束，構建區間模型下聲子晶體的可靠性優化模型，并采用Chebyshev展開法對局域共振帶隙進行優化。

2.1 Chebyshev展開法

Monte-Carlo法可用于區間模型下聲子晶體帶隙分析。但其高昂的計算成本嚴重限制其工程應用價值。本節討論將Chebyshev多項式展開引入到區間聲子晶體不確定性分析，構建區間聲子晶體的Chebyshev代理模型。

對于，第n階Chebyshev多項式定義如下[23]

Cn(x)=cos (nθ)

(9)

其中n表示非負整數。在區間上，Chebyshev多項式的遞推關系可表示如下

(10)

L個變量的Chebyshev多項式定義為

Cn1,n2,…,nx(x1,x2,…,xL)=

cos (n1θ1)cos (n2θ2)…cos(nLθL)

(11)

其中，θi=arccos(xi), (i=1, 2, …,L)

聲子晶體的能帶曲線的Chebyshev多項式表示如下[24]

fs(x,[WTHX〗k)=

(12)

其中,p表示下標i1,i2,…,iL等于零的數量。下標s(s=1,2,3…,10)表示第s條能帶。fi1,…,iL(k)表示多項式系數。向量x包含L個變量。系數fi1,…,iL(k)為波矢k的函數。

系數fi1,…,iL(k)可由下式推得

fi1,i2,…,iL(k)=

cosi1θ1…cosiLθL)dθ1…dθL

(13)

其中,L表示變量數，下標i1,…,iL=0,1,2,…,n。

用梅勒積分公式[25]進行轉換后可得到

Ci1,…,iL(xj1,…,xjL) ≈

cosi1θj1…cosiLθjL

(14)

直接基于有限元模型，采用Monte-Carlo法分析聲子晶體帶隙的變化范圍，其計算精度隨著樣本量的增加而趨于精確解。但是有限元模型的反復計算會帶來巨額的計算負擔。通過Chebyshev代理模型替換有限元模型，僅需采用有限的樣本數據構造Chebyshev多項式，極大地提高了聲子晶體帶隙的分析效率。

2.2 區間可靠性優化模型

聲子晶體的帶隙特性是聲子晶體的重要特性之一。帶隙的最大化可以使得聲子晶體的禁帶變寬，在工程實際中降噪效果更為明顯。因此，本文的目標函數是聲子晶體帶隙的最大化。可以表達為

max{Fmax-Fmin}

(15)

式中帶隙的下邊界頻率Fmin和上邊界頻率Fmax可以表示為

Fmin=inf{Δf}

Fmax=sup{Δf}

(16)

其中，Δf表示帶隙。由于不確定參數的存在，Fmin和Fmax均不是一個確定的數值，而是一個變化的區間。圖2為聲子晶體能帶的局部放大圖，圖中Fmin,1和Fmin,2分別表示帶隙下邊界頻率Fmin變化區間的下限及上限。Fmax,1和Fmax,2分別表示帶隙上邊界頻率Fmax變化區間的下限及上限。以區間[Fmin,2,Fmax,1]的最大化為優化目標，優化后的帶隙是保守的。若以區間[Fmin,1,Fmax,1]、[Fmin,2,Fmax,2]或[Fmin,1,Fmax,2]中的任意一個為優化的帶隙時，由于不確定參數導致的帶隙區間波動，實際帶隙不一定完全處于目標區間中，優化結果存在一定的風險。

圖2 聲子晶體能帶局部放大圖Fig.2 Local enlarged drawing of phononic crystal energy band

因此，本文以區間[Fmin,2,Fmax,1]變化范圍的最大化為目標函數，式(15)的更為具體的表達形式應為

max{Fmax,1-Fmin,2}

(17)

其中，Fmax,1和Fmin,2可用如下形式表示

?

?

(18)

在實際工程應用中，某個頻率段的噪聲或振動對于生產是有害的，因此往往需要屏蔽特定頻率段的振動或噪聲。本文將期望的頻帶作為約束條件，若[f1,f2]為所期望的頻帶，則目標函數中的帶隙[Fmin,2,Fmax,1]的上邊界頻率Fmax,1要不小于f2，下邊界頻率Fmin,2要不大于f1。該約束條件可表示為

Fmin,2≤f1≤f2≤Fmax,1

(19)

綜上可知，聲子晶體的區間可靠性優化模型為

max {Fmax,1-Fmin,2}

(20a)

s.t.Fmin2≤f1≤f2≤Fmax 1

?

?

(20b)

2.3 區間可靠性優化模型的優化流程

遺傳算法模擬達爾文的進化論機理來尋找最優解，其過程簡單，又具有很好的收斂性與魯棒性。因此，本文選用遺傳算法作為優化算法。以Chebyshev代理模型為基礎，聲子晶體區間優化模型的求解流程如下：

步驟1將設計變量R1…Rn作為個體進行編碼，并生成一組初始群體；

步驟2基于Chebyshev代理模型，用Monte-Carlo 法求解能帶的變化范圍，計算個體適應度值；

步驟3計算約束條件的值；

步驟4判斷可靠性約束條件是否成立？若成立，則進行下一步；若不成立，則更新上一個群體，并回到第二步；

步驟5判斷是否達到演化循環次數(給定的演化循環次數為10 000次)，若達到則結束并輸出結果；若未達到則對群體P(t)進行一輪選擇、交叉、變異運算之后可得到新一代的群體P(t+1)，并回到第二步。

3 數值算例

3.1 含梳狀夾層二維局域共振聲子晶體模型

本文所采用的聲子晶體為含梳狀夾層的聲子晶體[26]，如圖3所示。它的內核為金屬，基體為聚合物，中間的夾層呈梳狀，由橡膠和空氣交替組成。聲子晶體的晶格常數為a，中間梳狀橡膠夾層由16個均勻分布的扇形單元組成，每個扇形單元的中心角為π/8，夾層的內外半徑分別是r1和r2。聚合物基體的材料參數分別是ρ=1 200 kg/m3，E=3.5×107Pa，γ=0.49。金屬內核的材料參數分別是ρ=8 950 kg/m3，E=2.1×1011Pa，γ=0.29。幾何參數初始值分別為a=20 mm，r1=5.4 mm，r2=8.0 mm。

圖3 含梳狀夾層的聲子晶體單胞Fig.3 Phononic crystal unit with a comb like coating

3.1.1 幾何結構參數對帶隙的影響

幾何結構參數是影響聲子晶體帶隙的一個重要因素。本節將探討夾層內半徑r1和外半徑r2對聲子晶體帶隙的影響。表1為幾何參數r1和r2的變化情況。其中，Case1表示外半徑r2和晶格常數a保持不變，而內徑r1從5 mm逐漸變化到7 mm；Case2表示內徑r1和晶格常數a保持不變，而外徑r2從7.5 mm逐漸變化到9.5 mm。Case1和Case2所對應的帶隙變化如圖4所示。

表1 聲子晶體幾何結構參數Tab.1 Geometrical structure parameter of phononic crystal

由圖4(a)可知，當聲子晶體內徑r1從5 mm逐漸變化到7 mm時，帶隙的上下界逐漸向高頻帶移動，且帶隙的絕對寬度逐漸變大。由圖4(b)可知，當聲子晶體外徑r2從7.5 mm逐漸變化到9.5 mm時，帶隙上下界逐漸向低頻帶移動。由此可知，幾何結構參數對聲子晶體的局域共振帶隙有顯著影響。通過對幾何結構參數的合理調整，可有效優化聲子晶體局域共振帶隙。因此，本文將幾何參數r1和r2作為設計變量對聲子晶體進行結構優化。

(a) r1

(b) r2圖4 結構參數對聲子晶體帶隙的影響Fig.4 Effects of structural parameters on the band gap of phononic crystals

3.1.2 聲子晶體帶隙的區間不確定分析

在工程實踐中，不確定性廣泛而必然地存在著。如果在忽略這些不確定因素的情況下對聲子晶體進行優化設計，那么得到的結果往往不可靠。

在區間模型中，不確定參數被定義為變化范圍已知的區間變量。梳狀夾層聲子晶體的橡膠夾層在工程實踐中受溫度等影響，其材料屬性易發生變化。故本文以橡膠夾層的密度ρ和楊氏模量E取為區間變量，如表2所示。

表2 區間變量的變化范圍及不確定度Tab.2 The range and uncertainty of interval variables

在區間變量E和ρ的變化范圍內等間距地各取10個樣本點，共獲得10×10=100組樣本數據。由圖5可知，第三條能帶與第四條能帶之間形成帶隙，帶隙下邊界頻率Fmin和上邊界頻率Fmax都不是確定的值而是變化的區間。用Fmin2和Fmin1分別表示帶隙下邊界頻率Fmin變化區間的上下限。用Fmax2和Fmax1分別表示帶隙上邊界頻率Fmax變化區間的上下限。表3為Fmin、Fmax和帶隙Δf的區間變化范圍及不確定度。由表3可以看出，Fmin、Fmax和帶隙Δf的不確定度分別為2.6%、2.7%和7.6%。而不確定性參數的不確定度為5%，聲子晶體帶隙的不確定度比不確定參數的不確定度略大。

由此可知，在受到不確定性參數影響時，聲子晶體的帶隙不是確定的值，而是變化的區間。而且，聲子晶體的帶隙的不確定度比不確定性參數中的最大不確定度還大。

圖5 不確定參數對帶隙的影響Fig.5 The influence of uncertain parameters on band gap表3 Fmin、Fmax和Δf的變化范圍及不確定度Tab.3 The range and uncertainty of Fmin、Fmax and Δf

區間變量變化范圍中間值不確定度Fmin /Hz[136.28, 143.41]139.852.6%Fmax /Hz 3[969, 1 071]286.952.7%Δf/Hz[135.92, 158.28]147.107.6%

3.1.3 Chebyshev代理模型精度驗證

考慮到聲子晶體橡膠夾層易受周圍環境因素的影響，橡膠夾層的密度ρ和楊氏模量E均為區間變量，其變化范圍分別為[969, 1 071]kg/m3、[0.95×105,1.05×105]Pa。橡膠夾層的剪切模量為γ=0.47。選定內徑r1和外徑r2為設計變量，其初始為r1=5.4 mm，r2=8.0 mm，其設計范圍為[5, 7]mm、[7.5, 9.5]mm。因此，Chebyshev多項式的變量向量x包含設計變量r1和r2以及區間變量ρ和E。采用三階Chebyshev多項式擬合能帶曲線。按照Chebyshev多項式抽樣法則，每個變量抽取4個樣本點，故在每個波矢k下共抽取44=256個樣本點。即采用有限元模型計算256組樣本數據，并基于這些樣本數據構建每條能帶的Chebyshev代理模型。

共求解聲子晶體的十條能帶。由于只有第三條能帶與第四條能帶之間產生帶隙，故只分析第三條、第四條能帶的精度以及其產生的帶隙精度。圖6(a)和圖6(b)分別表示r1=6.5 mm和r2=8.5 mm時，基于有限元模型和Chebyshev代理模型的能帶圖。由圖6可以看出，基于Chebyshev代理模型得到的能帶與基于有限元模型得到的能帶匹配的非常好，Chebyshev代理模型得到的能帶上下界與有限元模型得到的能帶上下界的誤差非常小。

(a)

(b)圖6 r1=6.5 mm，r2=8.5 mm時有限元方法的能帶圖及 Chebyshev方法的能帶圖

Fig.6 The energy band structure based on the finite element method，and the Chebyshev method (r1=6.5 mm，r2=8.5 mm)

有限元模型每分析一次的時間為10 h。若Monte-Carlo法的樣本數據為100時，采用有限元模型直接分析聲子晶體帶隙的變化范圍時，其計算成本將達到1 000 h。采用三階Chebyshev代理模型分析聲子晶體帶隙的變化范圍時，構建Chebyshev代理模型的時間為72.6 h，基于100組樣本點分析帶隙變化范圍的時間為4.3 h。由此可見，若直接采用有限元模型來進行優化，由于極其昂貴的計算成本而將導致優化過程難以實現。本文通過引入Chebyshev代理模型，在保證精度的前提下，極大減少優化時間，提高優化效率。

3.1.4 區間可靠性優化

該聲子晶體的區間可靠性優化模型為

max {Fmax,1-Fmin,2}

(21a)

s.t.Fmin2≤200 Hz≤450 Hz≤Fmax 1

E∈[0.95×105, 1.05×105] Pa

ρ∈[969, 1 071] kg/m3

5 mm≤r1≤7 mm

7.5 mm≤r2≤9.5 mm

(21b)

采用遺傳算法求解基于Chebyshev代理模型的區間聲子晶體的可靠性優化模型，得到了聲子晶體最大帶隙及其對應的設計變量。種群規模為N=20，交叉概率為Pc=0.7，變異概率為Pm=0.3，最大迭代代數為Gene=60。圖7給出了聲子晶體最大帶隙的收斂過程，圖中為最優的個體對應的適應度值(即最大帶隙)隨進化次數變化的趨勢。由圖7可知，聲子晶體最優帶隙在進化代數達到10代時實現收斂。收斂后得到的最優帶隙為290.69 Hz。

圖7 聲子晶體的帶隙優化收斂圖Fig.7 Optimal convergence graph of phononic crystal bandgap

圖8為優化前的能帶圖，設計變量為初始設計變量r1=5.4 mm，r2=8.0 mm。帶隙頻帶為[143.69, 280.26] Hz，帶隙寬度為136.57 Hz，且僅頻帶[200, 280.26] Hz落在期望帶隙[200, 450] Hz內。圖9為優化后的能帶圖，優化后的設計變量為r1=6.99 mm，r2=8.13 mm。優化后的帶隙頻帶為[199.85, 490.54] Hz，帶隙寬度為290.69 Hz，完全包含期望頻帶[200, 450] Hz。相對于優化前的帶隙，優化后的帶隙完全滿足約束條件，且帶隙寬度從136.57 Hz拓寬到290.69 Hz。

圖8 優化前能帶圖Fig.8 The energy band structure before optimization

圖9 優化后能帶圖Fig.9 The energy band structure after optimization

4 結 論

本文針對不確定性廣泛存在于聲子晶體，并嚴重影響其物理性質這一現狀，將區間模型引入聲子晶體，描述其模型參數的不確定性。數值分析結果表明，Chebyshev代理模型能高效且較精確地預測區間模型下聲子晶體的帶隙變化范圍。在考慮區間不確定性的條件下，本文以Chebyshev 代理模型為基礎構建聲子晶體優化模型。優化結果表明，優化后的帶隙相對于優化前有大幅度拓寬，且滿足期望頻帶這一約束條件，聲子晶體的聲音屏蔽性能得到極大地改善。本文所提出的方法在聲子晶體的工程實際應用中有廣泛的應用前景，不僅極大地提高了不確定條件下聲子晶體帶隙特性的可靠性，而且通過構建代理模型的方式避免了聲子晶體優化設計的巨額計算成本。