基于SiPESC 平臺(tái)的聲子晶體能帶結(jié)構(gòu)分析算法研究及軟件實(shí)現(xiàn)

姜殿恒，陳飆松，張 盛，李云鵬

(大連理工大學(xué)工業(yè)裝備結(jié)構(gòu)分析國(guó)家重點(diǎn)實(shí)驗(yàn)室，遼寧，大連 116024)

聲子晶體是一類(lèi)擁有帶隙和缺陷態(tài)的特點(diǎn)的多介質(zhì)周期性結(jié)構(gòu)[1]。因其擁有帶隙和缺陷態(tài)的特性，聲子晶體可以阻隔與調(diào)控聲波與彈性波在介質(zhì)中的傳播[2?3]，具有豐富的工程應(yīng)用前景[4? 7]，可應(yīng)用在各類(lèi)減振降噪系統(tǒng)中[8?10]。聲子晶體的能帶結(jié)構(gòu)分析可使用多種數(shù)值分析方法，諸如傳遞矩陣法、平面波展開(kāi)法[11]、時(shí)域有限差分法[12]、多重散射理論[13?14]、有限元方法[15]等。各類(lèi)方法的離散方式不同，數(shù)值計(jì)算手段不同，具有不同優(yōu)勢(shì)和劣勢(shì)，從而適用于不同類(lèi)型的分析問(wèn)題。有限元方法具有不受問(wèn)題維度及模型的幾何特征所約束的特點(diǎn)，適用于各類(lèi)能帶結(jié)構(gòu)分析問(wèn)題，廣泛應(yīng)用于不同介質(zhì)、不同周期性的聲子晶體的能帶結(jié)構(gòu)分析。有限元方法計(jì)算量大、計(jì)算過(guò)程復(fù)雜，對(duì)于有限元方法的算法研究與軟件開(kāi)發(fā)是聲子晶體仿真方法研究的核心問(wèn)題。

聲子晶體的有限元方法能帶結(jié)構(gòu)分析通常借助于諸如COMSOL、ABAQUS 等成熟的商用有限元軟件[16]。此類(lèi)軟件沒(méi)有直接的能帶分析模塊，研究人員通過(guò)設(shè)定多個(gè)模態(tài)分析步，分別調(diào)用商業(yè)軟件中的廣義特征值計(jì)算模塊以完成聲子晶體仿真功能。COMSOL 軟件支持Hermitian 矩陣的廣義特征值求解，只能定義網(wǎng)格的疏密性，無(wú)法進(jìn)行獨(dú)立單元的修改操作，因而無(wú)法使用該軟件進(jìn)行模型拓?fù)鋬?yōu)化計(jì)算。ABAQUS 等結(jié)構(gòu)有限元軟件可對(duì)有限元網(wǎng)格單獨(dú)修改，卻不支持Hermitian矩陣的廣義特征值求解。

本文軟件基于工程與科學(xué)計(jì)算平臺(tái)SiPESC 的有限元分析模塊SiPESC.FEMS[17]開(kāi)發(fā)。SiPESC平臺(tái)為“微核心+插件”的設(shè)計(jì)模式，所有有限元分析功能被封裝為名為插件的軟件模塊(在后文中通稱軟件模塊)。本文依托于SiPESC.FEMS 軟件模塊中的模型文件導(dǎo)入功能、多類(lèi)型單元庫(kù)及后處理功能，實(shí)現(xiàn)了有限元方法的聲子晶體能帶結(jié)構(gòu)分析功能，完成了聲子晶體能帶結(jié)構(gòu)分析軟件。該功能可同樣用于各類(lèi)周期隔振結(jié)構(gòu)的仿真計(jì)算。通過(guò)將該軟件與COMSOL 對(duì)比，驗(yàn)證了其實(shí)體單元能帶結(jié)構(gòu)計(jì)算能力、結(jié)構(gòu)單元能帶結(jié)構(gòu)計(jì)算能力、大規(guī)模結(jié)構(gòu)能帶結(jié)構(gòu)計(jì)算能力。

1 有限元能帶結(jié)構(gòu)分析的基本理論

聲子晶體能帶結(jié)構(gòu)體現(xiàn)了彈性波在聲子晶體介質(zhì)中的傳播能力。通過(guò)將波矢遍歷不可約布里淵區(qū)的邊界，計(jì)算得到對(duì)應(yīng)波長(zhǎng)的波矢在聲子晶體中的傳播頻率，將其繪制為波矢頻率關(guān)系曲線即為聲子晶體的能帶結(jié)構(gòu)，又稱色散關(guān)系曲線。

波在理想無(wú)限周期聲子晶體中的傳播可以看作是一種平穩(wěn)狀態(tài)，Bloch 證明了當(dāng)簡(jiǎn)諧波在介質(zhì)中傳播時(shí)介質(zhì)位移符合簡(jiǎn)諧波位移場(chǎng)u(r,t)。

則單胞對(duì)應(yīng)周期維度的維度正方向邊界位移u+與維度負(fù)方向邊界位移u?符合Bloch 邊界條件是式(1)成立的必要條件。

因此聲子晶體的能帶分析問(wèn)題可以歸結(jié)為含Bloch 邊界條件約束的廣義特征值問(wèn)題：

顯然，用于模態(tài)分析的結(jié)構(gòu)線性有限元分析中的剛度矩陣K與質(zhì)量矩陣M為實(shí)對(duì)稱矩陣。Bloch 邊界條件中顯含復(fù)系數(shù)，將原實(shí)廣義特征值問(wèn)題變?yōu)閺?fù)數(shù)矩陣廣義特征值問(wèn)題。同時(shí)，由于約束的施加方法為酉合同變換，可通過(guò)Hermitian矩陣得定義公式證明約束后的矩陣K? 、M?為Hermitian矩陣，符合Hermitian 矩陣的定義式，即：

因此，聲子晶體能帶結(jié)構(gòu)分析，經(jīng)過(guò)上述公式推導(dǎo)，可整理為Hermitian 矩陣的廣義特征值問(wèn)題。能帶結(jié)構(gòu)分析往往只關(guān)心低階的前幾十階特征對(duì)，而非全部特征對(duì)，需通過(guò)稀疏Hermitian 矩陣得廣義特征值求解算法對(duì)其求解。

2 Hermitian 矩陣子空間迭代算法

Hermitian 矩陣是一類(lèi)復(fù)共軛對(duì)稱矩陣。結(jié)構(gòu)有限元分析中常用的實(shí)廣義特征值算法無(wú)法直接用于復(fù)Hermitian 矩陣的廣義特征值問(wèn)題的求解。本節(jié)在實(shí)數(shù)子空間迭代法的基礎(chǔ)上發(fā)展了用于Hermitian 矩陣子空間迭代算法。

子空間迭代法是一類(lèi)典型的特征值迭代求解方法，由BATHE[18]為求解結(jié)構(gòu)動(dòng)力學(xué)問(wèn)題引入力學(xué)領(lǐng)域。子空間迭代法代表了求解大型稀疏矩陣特征值的經(jīng)典思路，即通過(guò)多向量的逆迭代方法構(gòu)造模態(tài)變換基底，使用正交的變換基底構(gòu)建一個(gè)隨著子空間迭代步數(shù)增加逐步逼近原問(wèn)題特征值的特征子空間。

Hermitian 矩陣子空間迭代法的算法流程如圖1所示。

圖1 Hermitian 子空間迭代法算法流程Fig. 1 Hermitian subspace iterative algorithm

算法的計(jì)算細(xì)節(jié)描述如下：

1) 首先，構(gòu)造數(shù)值隨機(jī)、各列互相正交的實(shí)數(shù)試探矩陣。該矩陣可描述為由多列正交的試探向量組成：X0={x1,x2, ··· ,xs}。

2) 迭代求解該廣義特征值問(wèn)題直至問(wèn)題收斂。

3) 將該試探矩陣同Hermitian 質(zhì)量矩陣M相乘，得到用于逆迭代的復(fù)特征向量Un：

使用歸一化后的試探矩陣重新計(jì)算式(9)～式(16)，重復(fù)迭代直至該廣義特征值問(wèn)題收斂，得到指定個(gè)數(shù)的特征對(duì)。

并且，由于子空間變換矩陣Qn、特征向量Rn塊內(nèi)正交性，一組正交向量經(jīng)過(guò)正交變換所產(chǎn)生向量是必然正交的。因此，在Hermitian 矩陣的廣義特征值求解中，同樣無(wú)需做塊內(nèi)正交化。將Hermitian矩陣的共軛轉(zhuǎn)置退化為轉(zhuǎn)置，則可以推導(dǎo)得到用于實(shí)數(shù)子空間迭代法的歸一化公式。

3 定位格節(jié)點(diǎn)匹配方案

在聲子晶體計(jì)算中施加Bloch 邊界條件是其與普通結(jié)構(gòu)振動(dòng)問(wèn)題計(jì)算的最大區(qū)別。理想聲子晶體需要對(duì)周期維度上的所有節(jié)點(diǎn)施加Bloch 邊界條件，方法是，對(duì)空間中的主從邊界節(jié)點(diǎn)進(jìn)行一一對(duì)應(yīng)的暴力搜索節(jié)點(diǎn)匹配。經(jīng)典的三維空間暴力搜索點(diǎn)-點(diǎn)匹配算法的時(shí)間復(fù)雜度為O(n2)。當(dāng)聲子晶體的模型自由度增加時(shí)，節(jié)點(diǎn)匹配流程的計(jì)算時(shí)間占比大幅增加，匹配時(shí)間與廣義特征值分析處于同一量級(jí)，大幅增加了聲子晶體仿真的計(jì)算時(shí)間。因此，高效的聲子晶體能帶結(jié)構(gòu)分析，必須降低節(jié)點(diǎn)匹配階段的運(yùn)算時(shí)間。本文通過(guò)引入碰撞檢測(cè)領(lǐng)域中的定位格節(jié)點(diǎn)匹配方案，發(fā)展了用于聲子晶體能帶結(jié)構(gòu)分析的主從節(jié)點(diǎn)匹配策略。

3.1 定位格方案描述

圖2 定位格節(jié)點(diǎn)關(guān)系Fig. 2 sub-patches-node relationship

[]int為取整操作。將主節(jié)點(diǎn)使用節(jié)點(diǎn)排序后的序號(hào)整理為相應(yīng)的定位格-節(jié)點(diǎn)關(guān)系。當(dāng)計(jì)算x、y、z維度的周期邊界條件時(shí)，令該周期維度的定位格數(shù)Nx、Ny、Nz等于1 即可。

3) 建立主節(jié)點(diǎn)和定位格的關(guān)聯(lián)信息數(shù)組。該方法中需要針對(duì)主節(jié)點(diǎn)構(gòu)造2 類(lèi)數(shù)組：長(zhǎng)度為定位格數(shù)nb的定位格數(shù)組與長(zhǎng)度為主節(jié)點(diǎn)個(gè)數(shù)nm的節(jié)點(diǎn)數(shù)組。定位格數(shù)組記錄了當(dāng)前定位格第一個(gè)節(jié)點(diǎn)的序號(hào)，主節(jié)點(diǎn)數(shù)組記錄了除第一個(gè)節(jié)點(diǎn)之外的所有節(jié)點(diǎn)序號(hào)。

4) 按照如下流程將定位格中的節(jié)點(diǎn)依次插入定位格數(shù)組和節(jié)點(diǎn)數(shù)組：

a) 將定位格中第一個(gè)節(jié)點(diǎn)序號(hào)放入定位格數(shù)組的相應(yīng)位置。

b) 將該定位格中的第二個(gè)節(jié)點(diǎn)序號(hào)放入第一個(gè)節(jié)點(diǎn)序號(hào)為序號(hào)的節(jié)點(diǎn)數(shù)組中。

c) 對(duì)第二個(gè)節(jié)點(diǎn)之后的節(jié)點(diǎn)循環(huán)，將后一個(gè)節(jié)點(diǎn)序號(hào)放入以前一個(gè)節(jié)點(diǎn)序號(hào)為序號(hào)的節(jié)點(diǎn)數(shù)組中。

d) 最后將沒(méi)有對(duì)應(yīng)節(jié)點(diǎn)序號(hào)放入的定位格數(shù)組與節(jié)點(diǎn)數(shù)組的值設(shè)置為 ?1。

圖3 中展示了以2 號(hào)定位格的數(shù)據(jù)為例的節(jié)點(diǎn)數(shù)組與定位格數(shù)組的構(gòu)建順序。至此，主節(jié)點(diǎn)定位格的兩類(lèi)數(shù)組已經(jīng)構(gòu)建完成。

圖3 主節(jié)點(diǎn)和定位格的關(guān)聯(lián)信息數(shù)組Fig. 3 Association information array of master node and sub-patches

5) 將從節(jié)點(diǎn)與主節(jié)點(diǎn)進(jìn)行匹配。先使用式(25)計(jì)算從節(jié)點(diǎn)與定位格的匹配關(guān)系，再順序?qū)?jié)點(diǎn)數(shù)組和定位格數(shù)組中的節(jié)點(diǎn)按照插入順序進(jìn)行遍歷，直至找到相應(yīng)匹配節(jié)點(diǎn)或遍歷至節(jié)點(diǎn)序號(hào)為?1時(shí)無(wú)匹配節(jié)點(diǎn)停止匹配。

3.2 節(jié)點(diǎn)匹配算法討論

空間節(jié)點(diǎn)匹配問(wèn)題是多體動(dòng)力學(xué)、計(jì)算機(jī)視覺(jué)、邊界元方法等問(wèn)題的核心問(wèn)題。其中比較典型的為：快速多極邊界元問(wèn)題的自適應(yīng)樹(shù)結(jié)構(gòu)[20]，和本文引入的定位格節(jié)點(diǎn)匹配。這兩類(lèi)匹配算法的核心思想是簡(jiǎn)單且相同的，即通過(guò)增加中間層來(lái)降低節(jié)點(diǎn)匹配的計(jì)算量。同時(shí)，二者依托數(shù)據(jù)結(jié)構(gòu)的不同造成它們搜索效率的不同。

自適應(yīng)樹(shù)結(jié)構(gòu)根據(jù)其處理的二維/三維空間問(wèn)題，分別對(duì)應(yīng)四叉樹(shù)/八叉樹(shù)的數(shù)據(jù)格式，平均搜索速度為O(Nlog4N)/O(Nlog8N)，最劣搜索速度為O(N2)。定位格是一種哈希表數(shù)據(jù)格式，平均搜索速度為O(N) ，最劣搜索速度為O(N2)。在理想情況下，定位格方法的搜索速度更高。

對(duì)于實(shí)際的碰撞或邊界元仿真等數(shù)值問(wèn)題而言，節(jié)點(diǎn)在空間中必然分布不均。而樹(shù)狀數(shù)據(jù)結(jié)構(gòu)可以適應(yīng)空間節(jié)點(diǎn)分布不均的問(wèn)題，定位格方法則沒(méi)有該類(lèi)優(yōu)點(diǎn)。在本文所討論的有限元網(wǎng)格聲子晶體界面匹配問(wèn)題中，在單胞界面上的網(wǎng)格節(jié)點(diǎn)(在保證良好的網(wǎng)格質(zhì)量的前提下)在通常情況下是均布在單胞界面上的。這一現(xiàn)狀使得本文在進(jìn)行算法選擇時(shí)無(wú)需考慮節(jié)點(diǎn)在空間中分布不均的問(wèn)題，消除了定位格方法本身的劣勢(shì)并可以充分發(fā)揮其優(yōu)勢(shì)，因此針對(duì)該問(wèn)題定位格方法更為優(yōu)越。同時(shí)，其占用空間大小固定，避免了反復(fù)增加數(shù)組空間的時(shí)間；數(shù)組在內(nèi)存中儲(chǔ)存連續(xù)，也使得使用C++語(yǔ)言中的指針遍歷數(shù)組效率高于內(nèi)存中不連續(xù)的列表結(jié)構(gòu)。

4 聲子晶體能帶分析軟件模塊實(shí)現(xiàn)

本節(jié)使用前文描述的節(jié)點(diǎn)匹配方法及廣義特征值求解算法，設(shè)計(jì)實(shí)現(xiàn)聲子晶體的能帶分析軟件模塊org.sipesc.fems.phononic。聲子晶體的能帶結(jié)構(gòu)分析可歸納為3 個(gè)有別于傳統(tǒng)有限元分析的步驟，即：

a)匹配主從節(jié)點(diǎn)，處理約束關(guān)系；

b)施加Bloch 邊界條件；

c)在該邊界條件下進(jìn)行Hermitian 矩陣的廣義特征值分析。

SiPESC 科學(xué)計(jì)算平臺(tái)已經(jīng)擁有完善的線性有限元分析功能，只需在SiPESC 平臺(tái)上實(shí)現(xiàn)節(jié)點(diǎn)匹配和Hermitian 特征值求解，即可實(shí)現(xiàn)聲子晶體的能帶分析功能。

聲子晶體能帶分析的有限元方法具體流程如圖4 所示，共可劃分為11 個(gè)計(jì)算步驟。圖4 虛線線框內(nèi)為實(shí)現(xiàn)能帶分析功能新增加的分析模塊，其余計(jì)算步驟為SiPESC 科學(xué)計(jì)算平臺(tái)已經(jīng)具備的有限元分析功能，可分別使用MDofUpdate 自由度更新類(lèi)與MMain 聲子晶體能帶結(jié)構(gòu)分析主控流程類(lèi)進(jìn)行計(jì)算控制，其余的計(jì)算流程都可以上兩計(jì)算類(lèi)使用接口調(diào)用實(shí)現(xiàn)。二者繼承自SiPESC 平臺(tái)的線性流程基類(lèi)MTaskManager，使用平臺(tái)統(tǒng)一的工廠模式在有限元計(jì)算流程中調(diào)用。本節(jié)所實(shí)現(xiàn)的聲子晶體能帶分析軟件模塊org.sipesc.fems.phononic 整體UML 類(lèi)關(guān)系如圖5 如所示。

圖4 SiPESC 平臺(tái)能帶結(jié)構(gòu)分析流程Fig. 4 band structure analysis process of SiPESC platform

圖5 中org.sipesc.core 軟件包為SiPESC 科學(xué)計(jì)算平臺(tái)的基礎(chǔ)軟件包，包含平臺(tái)所需的各類(lèi)基礎(chǔ)功能、擴(kuò)展所需繼承的基類(lèi)、代數(shù)方程組求解器與其他矩陣計(jì)算類(lèi)。org.sipesc.core.engdbs 為平臺(tái)的數(shù)據(jù)庫(kù)軟件包，負(fù)責(zé)數(shù)據(jù)庫(kù)相關(guān)的IO 操作和數(shù)據(jù)庫(kù)相關(guān)類(lèi)型的創(chuàng)建與修改。org.sipesc.fems.elements 為平臺(tái)的單元軟件包，包含各類(lèi)實(shí)體單元、結(jié)構(gòu)單元、和用于自由度控制的單元控制類(lèi)。

圖5 SiPESC 科學(xué)計(jì)算平臺(tái)聲子晶體軟件模塊的UML 類(lèi)圖Fig. 5 UML class diagram of phononic crystal plugin in SiPESC platform

4.1 約束自由度處理模塊

MDofUpdate 類(lèi)繼承于SiPESC 平臺(tái)的線性計(jì)算流程類(lèi)MTaskManager，擁有初始化方法initialize和計(jì)算流程方法Start。該類(lèi)所實(shí)現(xiàn)的功能是將單胞的邊界節(jié)點(diǎn)劃分為界面主節(jié)點(diǎn)與界面從節(jié)點(diǎn)，對(duì)主從節(jié)點(diǎn)根據(jù)對(duì)應(yīng)的周期維度進(jìn)行匹配，得到一一對(duì)應(yīng)的主從節(jié)點(diǎn)關(guān)系。最后將節(jié)點(diǎn)自由度劃分為主從自由度，用以后續(xù)計(jì)算的自由度排序處理。

約束自由度處理的主要步驟為主從節(jié)點(diǎn)匹配，針對(duì)該問(wèn)題抽象出邊界節(jié)點(diǎn)匹配類(lèi)MBoundary。該類(lèi)通過(guò)readBoundary 方法識(shí)別主從邊界節(jié)點(diǎn)，用以后續(xù)的算法匹配計(jì)算。聲子晶體的周期性可分為1 維、2 維、3 維 這3 種，其中1 維、2 維周期邊界條件都可通過(guò)3 維周期邊界條件退化得到。針對(duì)具有三維周期性的聲子晶體問(wèn)題，可將主從邊界節(jié)點(diǎn)類(lèi)型分為3 類(lèi)：點(diǎn)主從約束、線主從約束、面主從約束。如圖6(a)、圖6(b)、圖6(c)所示，圖中α 為頂點(diǎn)主自由度，a、b、c為邊主自由度，A、C、E為面主自由度，其余為從自由度。圖6(d)、圖6(e)、圖6(f)為約束關(guān)系所整理的約束關(guān)系樹(shù)。首先討論圖6(a)面主從約束的處理，以最為常見(jiàn)的簡(jiǎn)立方體晶格為例，針對(duì)三維周期問(wèn)題該晶格存在3 主面、3 從面，在不考慮頂點(diǎn)和棱邊的情況下，主面節(jié)點(diǎn)與從面節(jié)點(diǎn)一一對(duì)應(yīng)。在該主從節(jié)點(diǎn)匹配類(lèi)中通過(guò)節(jié)點(diǎn)匹配方法match2D，調(diào)用MFaceMatch 類(lèi)進(jìn)行面節(jié)點(diǎn)的匹配。MFace-Match 類(lèi)為定位格策略節(jié)點(diǎn)匹配類(lèi)，該類(lèi)實(shí)現(xiàn)了第3 節(jié)中所敘述的定位格節(jié)點(diǎn)匹配策略，進(jìn)行快速的主從節(jié)點(diǎn)匹配。匹配后面主從節(jié)點(diǎn)可直接套用Bloch 邊界條件公式進(jìn)行計(jì)算。

線約束與點(diǎn)約束會(huì)存在某一節(jié)點(diǎn)同時(shí)為主面節(jié)點(diǎn)與從面節(jié)點(diǎn)的情況，屬于多層的主從約束匹配問(wèn)題。線約束的幾何關(guān)系如圖6(b)，可將其整理為如圖6(e)的主從約束樹(shù)。簡(jiǎn)立方體晶格棱線d、g上 的節(jié)點(diǎn)分別在x、z周 期維度方向上是a棱線上節(jié)點(diǎn)的從節(jié)點(diǎn)，同時(shí)在x、z周期維度方向上為i棱線的主節(jié)點(diǎn)。經(jīng)過(guò)簡(jiǎn)單的遞推，可將主從關(guān)系整理為單層主從約束。而后調(diào)用match1D 方法針對(duì)邊主從節(jié)點(diǎn)，進(jìn)行匹配，邊主從節(jié)點(diǎn)的匹配方式僅需要匹配非周期自由度的坐標(biāo)即可。

圖6 三類(lèi)約束幾何關(guān)系及其關(guān)系樹(shù)Fig. 6 Three kinds of constrained geometric relations and their relation trees

點(diǎn)約束同線約束在處理方法上是一致的，僅僅將2 層主從約束變?yōu)? 層主從約束，且各頂點(diǎn)無(wú)需進(jìn)行節(jié)點(diǎn)匹配。

4.2 能帶分析模塊

MMain 流程類(lèi)為聲子晶體能帶分析計(jì)算主控流程類(lèi)，其負(fù)責(zé)了從組裝無(wú)約束條件下的聲子晶體單胞總剛度矩陣至計(jì)算得到能帶色散曲線的全部流程具體流程如圖7 所示。

圖7 能帶結(jié)構(gòu)分析算法流程Fig. 7 Band structure analysis algorithm

有限元能帶結(jié)構(gòu)分析的計(jì)算主體為結(jié)構(gòu)的總體剛度矩陣K，總體質(zhì)量矩陣M。針對(duì)有限元總體矩陣對(duì)象，抽象出有限元總體矩陣類(lèi)MStiffMat。該類(lèi)通過(guò)formGlobalStiffMatrix 方法實(shí)現(xiàn)有限元的總體剛度矩陣組集，獲得有限元的總體剛度、質(zhì)量矩陣。繼而通過(guò)transformBoundary 方法實(shí)現(xiàn)了直接消去法的約束變換矩陣N的計(jì)算。使用addConstraint 方法實(shí)現(xiàn)直接消去法的約束變換矩陣乘法公式(8)。最后使用outputCon 方法將施加約束變換后的總體剛度矩陣和總體質(zhì)量矩陣輸出至MMain流程類(lèi)中。

針對(duì)Hermitian 矩陣廣義特征值求解問(wèn)題抽象出子空間迭代計(jì)算類(lèi)SubSpaceIter。該類(lèi)實(shí)現(xiàn)了第二節(jié)中的Hermitian 子空間迭代廣義特征值算法。通過(guò)該類(lèi)求解得到當(dāng)前波矢約束條件下的廣義特征值。

最后，針對(duì)結(jié)果輸出操作抽象出結(jié)果輸出對(duì)象exportorManager，該類(lèi)使用exportUnvModal 方法調(diào)用私有節(jié)點(diǎn)數(shù)據(jù)輸出方法exportUnvNodePart，私有單元數(shù)據(jù)輸出方法exportUnvElePart，私有模態(tài)數(shù)據(jù)輸出方法exportUnvModePart，將計(jì)算結(jié)果輸出為SiPESC 平臺(tái)的后處理文件類(lèi)型*.unv。

4.3 能帶分析算法討論

能帶結(jié)構(gòu)分析可以歸結(jié)為一個(gè)不斷對(duì)矩陣做攝動(dòng)的特征值重分析問(wèn)題，這類(lèi)攝動(dòng)問(wèn)題有一類(lèi)簡(jiǎn)單的計(jì)算技巧，即使用攝動(dòng)前的收斂解作為攝動(dòng)后的初值。很多研究人員都已經(jīng)將該技巧使用在其所做的數(shù)值研究中。這一技巧可以使特征值迭代時(shí)的試探向量較隨機(jī)值更接近于終值，可以極大地加快能帶計(jì)算時(shí)的收斂效率。而COMSOL中的能帶計(jì)算流程(參考COMSOL 官網(wǎng)的聲子晶體算例)是使用參數(shù)遍歷方法，在不同的波矢點(diǎn)上反復(fù)進(jìn)行初始向量隨機(jī)的LANCZOS 方法，并沒(méi)有繼承上一步的收斂解。這是SiPESC 與COMSOL效率差異的第一點(diǎn)原因。

SiPESC 平臺(tái)與其他商用有限元軟件也有顯著區(qū)別。這些結(jié)構(gòu)有限元軟件無(wú)法施加復(fù)數(shù)的邊界約束條件，只能使用實(shí)邊界條件施加約束。將剛度矩陣K劃分為實(shí)部和虛部：

進(jìn)行求解。這一算法將矩陣維度增廣為2 倍，其求解的同階特征對(duì)也會(huì)變?yōu)? 倍。例如，使用SiPESC 平臺(tái)求解10 階模態(tài)可獲得10 階特征對(duì)，使用上文所述辦法只能求解5 階特征對(duì)，每階特征對(duì)出現(xiàn)2 次。由于增加求解的特征對(duì)數(shù)量時(shí)，特征值算法的計(jì)算量并非線性增長(zhǎng)，因此二者同時(shí)計(jì)算10 階特征值時(shí)SiPESC 計(jì)算效率遠(yuǎn)高于使用同類(lèi)商業(yè)軟件。

5 數(shù)值算例

5.1 計(jì)算環(huán)境說(shuō)明

為驗(yàn)證所開(kāi)發(fā)數(shù)值仿真軟件模塊的計(jì)算能力，在本節(jié)中分別針對(duì)實(shí)體單元組成的聲子晶體模型，結(jié)構(gòu)單元組成的聲子晶體模型及大規(guī)模聲子晶體模型進(jìn)行了數(shù)值算例驗(yàn)證。選取文獻(xiàn)中應(yīng)用較多的多物理場(chǎng)數(shù)值仿真軟件COMSOL Multiphysics 5.5 軟件作為標(biāo)準(zhǔn)對(duì)照軟件。

受現(xiàn)實(shí)因素限制，無(wú)法使用同一計(jì)算環(huán)境對(duì)兩軟件進(jìn)行比較，現(xiàn)將它們部署的硬件環(huán)境列舉于表1。顯然，COMSOL 所部署的硬件環(huán)境的CPU 主頻(3.20 GHz)高于SiPESC 平臺(tái)所部屬硬件環(huán)境的CPU 主頻(2.00 GHz)，COMSOL 所部署的硬件環(huán)境的DDR4 內(nèi)存全面優(yōu)于SiPESC 平臺(tái)的DDR3 內(nèi)存。同時(shí)，通過(guò)將SiPESC 所部屬的工作站并行核數(shù)限制在6 核(通過(guò)OpenMP 設(shè)置)，使兩種計(jì)算環(huán)境最大并行線程數(shù)一致。因此，理論上，COMSOL 平臺(tái)在算例測(cè)試上計(jì)算效率更佳。為證明本文所實(shí)現(xiàn)軟件模塊的高效性，在此忽略SiPESC 平臺(tái)計(jì)算在環(huán)境上的劣勢(shì)，將二者視為同環(huán)境進(jìn)行對(duì)比。

表1 材料數(shù)據(jù)Table 1 Material data

為了驗(yàn)證本文所實(shí)現(xiàn)的軟件效率，選取了3 例彈性波聲子晶體算例進(jìn)行測(cè)試，因此在COMSOL計(jì)算中選用了固體力學(xué)模塊，施加Floquet 周期性條件進(jìn)行模型建立。為保證計(jì)算網(wǎng)格的一致性使用COMSOL 軟件進(jìn)行建模，MSC.Patran 前處理軟件進(jìn)行網(wǎng)格顯示。在折疊結(jié)構(gòu)算例中，額外使用了固體力學(xué)中殼體單元的物理場(chǎng)模塊，以進(jìn)行結(jié)構(gòu)單元模擬。在第5 節(jié)中所使用的計(jì)算環(huán)境都統(tǒng)一遵從5.1 節(jié)的說(shuō)明，不再額外在具體算例描述中敘述。

5.2 經(jīng)典局域共振聲子晶體算例

本節(jié)使用劉正猷于2000 年提出的經(jīng)典局域共振聲子晶體構(gòu)型[21]作為測(cè)試模型，選取COMSOL Multiphysics 5.5 軟件作為比較對(duì)象，以驗(yàn)證本文程序的計(jì)算速度與計(jì)算精度。該模型為3 組元聲子晶體模型，使用硅橡膠包裹鉛球體埋置入環(huán)氧樹(shù)脂基體，結(jié)構(gòu)如圖8 所示。具體材料參數(shù)如表2所示，結(jié)構(gòu)幾何尺寸見(jiàn)表3。

圖8 局域共振聲子晶體模型Fig. 8 Locally resonant sonic materials

表2 材料數(shù)據(jù)Table 2 Material data

表3 幾何尺寸Table 3 Geometric dimension

算例使用COMSOL 軟件默認(rèn)單元：10 節(jié)點(diǎn)二階四面體單元，計(jì)算網(wǎng)格兩軟件完全一致。該聲子晶體單胞有限元模型網(wǎng)格具有134 898 自由度，將不可約布里淵區(qū)邊界離散為120 個(gè)波矢點(diǎn)。共計(jì)算20 階模態(tài)，SiPESC 計(jì)算時(shí)間為4 h :31 min，COMSOL 計(jì)算時(shí)間為11 h : 39 min，SiPESC平臺(tái)的聲子晶體能帶分析程序在相同的計(jì)算環(huán)境下效率顯著高于COMSOL。

分析能帶結(jié)構(gòu)及其振動(dòng)模式(圖9)，顯示該結(jié)構(gòu)在遍歷不可約布里淵區(qū)中低階呈現(xiàn)出3 類(lèi)共振模式。其中L1、L2、L3共振模式為鉛球芯體在橡膠包覆層下圍繞三個(gè)坐標(biāo)軸的旋轉(zhuǎn)共振模式；L7則體現(xiàn)為橡膠包覆層的旋轉(zhuǎn)共振模式。這兩類(lèi)共振模式都呈現(xiàn)為旋轉(zhuǎn)振動(dòng)模式，難以與低頻長(zhǎng)行波產(chǎn)生耦合作用。L4、L5、L6與L′4、L′5、L6′都呈現(xiàn)出鉛球芯體與樹(shù)脂基體的反向振動(dòng)，與低頻長(zhǎng)行波有較長(zhǎng)的耦合作用，導(dǎo)致了較寬的第一帶隙的產(chǎn)生。其中 Γ 點(diǎn)處L4、L5、L6共振模式能帶穿透了L1、L2、L3共振模式能帶，由于繪圖工具順序繪圖的邏輯無(wú)法在圖上呈現(xiàn)，特此注明。

圖9 局域共振聲子晶體能帶結(jié)構(gòu)及振型Fig. 9 Band structure and mode shapes of local resonant sonic crystals

將二者能帶分析結(jié)果進(jìn)行對(duì)比，SiPESC 與COMSOL 能帶分析結(jié)果精度在低階一致，高階剪切振型有較小誤差。其原因大致為兩類(lèi)軟件的單元對(duì)剪切項(xiàng)的處理不同。本文使用單元為防止剪切自鎖現(xiàn)象，使用減縮積分進(jìn)行積分計(jì)算。

模態(tài)振型進(jìn)行對(duì)比，圖9 中振型圖中展示了兩組振型結(jié)果。在每一組中左圖為COMSOL 軟件計(jì)算結(jié)果，右圖中為SiPESC 軟件計(jì)算結(jié)果。通過(guò)振型對(duì)比兩軟件計(jì)算振型一致。SiPESC 后處理顯示時(shí)將二階單元處理為線性單元顯示，在計(jì)算中仍為二階單元計(jì)算，特此注明。

該算例證明，本文所實(shí)現(xiàn)軟件較商業(yè)軟件具有一致的計(jì)算精度和較高的計(jì)算效率，基本實(shí)現(xiàn)了高效的聲子晶體能帶結(jié)構(gòu)仿真功能。

5.3 大規(guī)模固/氣聲子晶體算例

為測(cè)試本能帶結(jié)構(gòu)分析軟件的大規(guī)模計(jì)算性能，使用2 組元60 萬(wàn)自由度的氣/固聲子晶體模型能帶結(jié)構(gòu)[22]，固體材料為樹(shù)脂，材料的物理性能見(jiàn)表4，幾何結(jié)構(gòu)如圖10 所示，幾何尺寸見(jiàn)表5。共計(jì)算20 階模態(tài)。

表4 材料數(shù)據(jù)Table 4 Material data

圖10 氣/固聲子晶體Fig. 10 Gas / solid phononic crystal

表5 幾何尺寸Table 5 Geometric dimension

本算例使用4 節(jié)點(diǎn)線性四面體單元，有限元模型具有659 952 自由度，將不可約布里淵區(qū)離散為100 個(gè)波矢點(diǎn)，遍歷計(jì)算能帶結(jié)構(gòu)，計(jì)算環(huán)境同算例1。計(jì)算時(shí)間為9 小時(shí)47 分鐘，能帶結(jié)構(gòu)如圖11 所示，與參考文獻(xiàn)[21]一致，該結(jié)構(gòu)存在低頻完全帶隙和高頻完全帶隙。其中低頻第一帶隙與經(jīng)典局域共振振型一致，芯體和框架呈現(xiàn)反向位移呈現(xiàn)局域共振特性，高頻帶隙帶隙中心頻率位于ct/2a附近，符合布拉格散射特性。圖中分別呈現(xiàn)了第一、第二帶隙兩側(cè)的單胞振型。左側(cè)一組振型為第一帶隙兩側(cè)頻率的振型：內(nèi)部埋置的樹(shù)脂球振動(dòng)和外框基體振動(dòng)。右側(cè)線框中為第二帶隙兩側(cè)頻率的振型：分別是橫波引起的外框的扭轉(zhuǎn)振動(dòng)，及縱波引起的漲縮振動(dòng)。

圖11 氣/固聲子晶體能帶結(jié)構(gòu)及振型Fig. 11 Band structure and mode shapes of gas/solid phononic crystal

該算例證明，本文所實(shí)現(xiàn)軟件具有大規(guī)模聲子晶體能能帶結(jié)構(gòu)分析計(jì)算能力，適用于大規(guī)模聲子晶體計(jì)算。

5.4 折疊結(jié)構(gòu)聲子晶體算例

為測(cè)試軟件的結(jié)構(gòu)單元計(jì)算能力，本節(jié)使用具有二維周期性的折疊結(jié)構(gòu)聲子晶體模型計(jì)算能帶結(jié)構(gòu)[23]。在文獻(xiàn)調(diào)研過(guò)程中，有文章[24]指出對(duì)該模型的單胞劃分會(huì)影響能帶分析結(jié)果。本文基于COMSOL 5.5 版本進(jìn)行了數(shù)值仿真，分別使用2 種單胞劃分方式(圖12 中劃分方式1 和劃分方式2)、2 種網(wǎng)格密度(1.5 k 自由度及20 k 自由度)組合進(jìn)行4 類(lèi)能帶結(jié)構(gòu)分析，使用COMSOL軟件計(jì)算其10 階能帶結(jié)構(gòu)，具體材料數(shù)據(jù)見(jiàn)表6，幾何尺寸見(jiàn)表7。SiPESC 平臺(tái)則使用劃分方式1 k～1.5 k 自由度網(wǎng)格進(jìn)行計(jì)算與COMSOL 軟件進(jìn)行對(duì)比。

表6 材料數(shù)據(jù)Table 6 Material data

表7 幾何尺寸Table 7 Geometric dimension

圖12 折疊結(jié)構(gòu)聲子晶體模型Fig. 12 Folded structure phononic crystal model

對(duì)比計(jì)算得到能帶結(jié)構(gòu)如圖13 所示，圖中實(shí)線為SiPESC 計(jì)算結(jié)果，虛線為COMSOL 計(jì)算結(jié)果。SiPESC計(jì)算使用80 s，COMSOL 使用4 小時(shí)30 秒以上時(shí)間不等。從圖中觀察得知，5 類(lèi)分析的能帶結(jié)果大致相符。為探究結(jié)果的精確度，此處取子繪圖13(a)與子繪圖13(b)討論。由于符合標(biāo)準(zhǔn)的有限元單元一定遵循分片實(shí)驗(yàn)，隨著網(wǎng)格加密其結(jié)果將逐漸逐漸收斂至解析解。在子繪圖13(a)中，劃分方式1、方式2 在網(wǎng)格加密后都趨向于SiPESC 平臺(tái)1.5 k 自由度的結(jié)果，因此可以得出結(jié)論，在針對(duì)該問(wèn)題的能帶結(jié)構(gòu)仿真中，SiPESC 平臺(tái)1.5 k 自由度時(shí)的板單元精度高與COMSOL 軟件20 k 自由度時(shí)的板單元能帶仿真精度。在子繪圖13(b)中只有SiPESC 平臺(tái)的結(jié)果中L1與L2兩類(lèi)共振模式的能帶簡(jiǎn)并于X點(diǎn)。觀察L1與L2兩類(lèi)共振模式可以得出二者是一組反對(duì)稱的共振模式，在周期結(jié)構(gòu)中二者的共振模式完全相同。因此L1與L2兩類(lèi)共振模式在X點(diǎn)的頻率也必然是相同的，兩能帶必須交于一點(diǎn)，因此證明了SiPESC 平臺(tái)結(jié)構(gòu)單元能帶分析精度高于COMSOL軟件。

圖13 折疊結(jié)構(gòu)聲子晶體能帶結(jié)構(gòu)Fig. 13 Band structure of folded phononic crystals

顯然，該折疊結(jié)構(gòu)不存在完全帶隙，RM、RY方向能帶結(jié)構(gòu)相似、能帶集中，符合各向異性周期結(jié)構(gòu)的特性。該結(jié)構(gòu)具有復(fù)雜的極化特性，已有諸多研究學(xué)者進(jìn)行過(guò)討論。

本節(jié)算例1 與算例3 與COMSOL 進(jìn)行了求解效率對(duì)比，其結(jié)果為，COMSOL 軟件的求解時(shí)間是SiPESC 平臺(tái)計(jì)算時(shí)間的2 倍以上，驗(yàn)證了SiPESC 平臺(tái)能帶分析軟件具有高效的仿真效率。求解算例1、3 時(shí)的內(nèi)存負(fù)載，二者大致相當(dāng)都在10 GB以下，在求解大規(guī)模模型算例2 時(shí)SiPESC 共使用內(nèi)存30 GB。證明SiPESC 平臺(tái)在內(nèi)存充足的情況下有能力求解更大的單胞結(jié)構(gòu)。精度二者大致相當(dāng)，在計(jì)算結(jié)構(gòu)單元時(shí)SiPESC 有明顯優(yōu)勢(shì)。

6 結(jié)論

本文基于SiPESC 平臺(tái)，實(shí)現(xiàn)了聲子晶體高效能帶結(jié)構(gòu)分析軟件，分別實(shí)現(xiàn)了用以廣義特征值分析的Hermitian 矩陣的子空間迭代法，與用以進(jìn)行主從節(jié)點(diǎn)匹配的定位格匹配方案。

(1)針對(duì)Hermitian 矩陣的廣義特征值計(jì)算問(wèn)題，發(fā)展了Hermitian 矩陣的子空間迭代法，并討論了對(duì)稱矩陣與Hermitian 矩陣在使用子空間迭代法時(shí)的算法差異。

(2)針對(duì)能帶計(jì)算中的節(jié)點(diǎn)匹配問(wèn)題，將用于碰撞檢測(cè)的定位格匹配方案，應(yīng)用于聲子晶體能帶結(jié)構(gòu)分析的主從約束匹配。

(3)基于SiPESC 工程與科學(xué)計(jì)算平臺(tái)有限元分析模塊，編寫(xiě)了用于聲子晶體能帶結(jié)構(gòu)仿真的軟件模塊。并介紹了該軟件模塊的架構(gòu)設(shè)計(jì)和多層主從約束處理等設(shè)計(jì)細(xì)節(jié)。

(4)通過(guò)將軟件模塊與經(jīng)典多物理場(chǎng)仿真軟件COMSOL Multiphysics 進(jìn)行對(duì)比，驗(yàn)證了該軟件模塊所實(shí)現(xiàn)算法具有可靠的數(shù)值精度和可觀的計(jì)算效率。在此基礎(chǔ)上計(jì)算了66 萬(wàn)自由度的氣固聲子晶體模型、100 個(gè)離散波矢點(diǎn)的20 階模態(tài)自振分析，進(jìn)一步驗(yàn)證了本文所實(shí)現(xiàn)軟件具有高效的大規(guī)模分析能力。