劉 堅, 陳俊煌, 夏百戰(zhàn), 滿先鋒

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

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

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

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

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

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

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

(1)

(2)

(3)

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

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

(4)

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

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

Ku=ω2Mu

(5)

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

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

(6)

(7)

式中:B表示應變矩陣;N表示形函數(shù)矩陣;Ve表示整個單胞區(qū)域。

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

f(r)=f(r+Na)

(8)

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

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

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

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

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

2 區(qū)間模型的構造及可靠性優(yōu)化理論

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

2.1 Chebyshev展開法

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

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

Cn(x)=cos (nθ)

(9)

其中n表示非負整數(shù)。在區(qū)間上，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等于零的數(shù)量。下標s(s=1,2,3…,10)表示第s條能帶。fi1,…,iL(k)表示多項式系數(shù)。向量x包含L個變量。系數(shù)fi1,…,iL(k)為波矢k的函數(shù)。

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

fi1,i2,…,iL(k)=

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

(13)

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

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

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

cosi1θj1…cosiLθjL

(14)

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

2.2 區(qū)間可靠性優(yōu)化模型

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

max{Fmax-Fmin}

(15)

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

Fmin=inf{Δf}

Fmax=sup{Δf}

(16)

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

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

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

max{Fmax,1-Fmin,2}

(17)

其中，F(xiàn)max,1和Fmin,2可用如下形式表示

?

?

(18)

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

Fmin,2≤f1≤f2≤Fmax,1

(19)

綜上可知，聲子晶體的區(qū)間可靠性優(yōu)化模型為

max {Fmax,1-Fmin,2}

(20a)

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

?

?

(20b)

2.3 區(qū)間可靠性優(yōu)化模型的優(yōu)化流程

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

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

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

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

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

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

3 數(shù)值算例

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

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

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

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

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

(a) r1

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

3.1.2 聲子晶體帶隙的區(qū)間不確定分析

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

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

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

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

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

圖5 不確定參數(shù)對帶隙的影響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

區(qū)間變量變化范圍中間值不確定度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代理模型精度驗證

考慮到聲子晶體橡膠夾層易受周圍環(huán)境因素的影響，橡膠夾層的密度ρ和楊氏模量E均為區(qū)間變量，其變化范圍分別為[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以及區(qū)間變量ρ和E。采用三階Chebyshev多項式擬合能帶曲線。按照Chebyshev多項式抽樣法則，每個變量抽取4個樣本點，故在每個波矢k下共抽取44=256個樣本點。即采用有限元模型計算256組樣本數(shù)據(jù)，并基于這些樣本數(shù)據(jù)構建每條能帶的Chebyshev代理模型。

共求解聲子晶體的十條能帶。由于只有第三條能帶與第四條能帶之間產(chǎn)生帶隙，故只分析第三條、第四條能帶的精度以及其產(chǎn)生的帶隙精度。圖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法的樣本數(shù)據(jù)為100時，采用有限元模型直接分析聲子晶體帶隙的變化范圍時，其計算成本將達到1 000 h。采用三階Chebyshev代理模型分析聲子晶體帶隙的變化范圍時，構建Chebyshev代理模型的時間為72.6 h，基于100組樣本點分析帶隙變化范圍的時間為4.3 h。由此可見，若直接采用有限元模型來進行優(yōu)化，由于極其昂貴的計算成本而將導致優(yōu)化過程難以實現(xiàn)。本文通過引入Chebyshev代理模型，在保證精度的前提下，極大減少優(yōu)化時間，提高優(yōu)化效率。

3.1.4 區(qū)間可靠性優(yōu)化

該聲子晶體的區(qū)間可靠性優(yōu)化模型為

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代理模型的區(qū)間聲子晶體的可靠性優(yōu)化模型，得到了聲子晶體最大帶隙及其對應的設計變量。種群規(guī)模為N=20，交叉概率為Pc=0.7，變異概率為Pm=0.3，最大迭代代數(shù)為Gene=60。圖7給出了聲子晶體最大帶隙的收斂過程，圖中為最優(yōu)的個體對應的適應度值(即最大帶隙)隨進化次數(shù)變化的趨勢。由圖7可知，聲子晶體最優(yōu)帶隙在進化代數(shù)達到10代時實現(xiàn)收斂。收斂后得到的最優(yōu)帶隙為290.69 Hz。

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

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

圖8 優(yōu)化前能帶圖Fig.8 The energy band structure before optimization

圖9 優(yōu)化后能帶圖Fig.9 The energy band structure after optimization

4 結 論

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