基于SiPESC 平臺的聲子晶體能帶結構分析算法研究及軟件實現

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

(大連理工大學工業裝備結構分析國家重點實驗室，遼寧，大連 116024)

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

聲子晶體的有限元方法能帶結構分析通常借助于諸如COMSOL、ABAQUS 等成熟的商用有限元軟件[16]。此類軟件沒有直接的能帶分析模塊，研究人員通過設定多個模態分析步，分別調用商業軟件中的廣義特征值計算模塊以完成聲子晶體仿真功能。COMSOL 軟件支持Hermitian 矩陣的廣義特征值求解，只能定義網格的疏密性，無法進行獨立單元的修改操作，因而無法使用該軟件進行模型拓撲優化計算。ABAQUS 等結構有限元軟件可對有限元網格單獨修改，卻不支持Hermitian矩陣的廣義特征值求解。

本文軟件基于工程與科學計算平臺SiPESC 的有限元分析模塊SiPESC.FEMS[17]開發。SiPESC平臺為“微核心+插件”的設計模式，所有有限元分析功能被封裝為名為插件的軟件模塊(在后文中通稱軟件模塊)。本文依托于SiPESC.FEMS 軟件模塊中的模型文件導入功能、多類型單元庫及后處理功能，實現了有限元方法的聲子晶體能帶結構分析功能，完成了聲子晶體能帶結構分析軟件。該功能可同樣用于各類周期隔振結構的仿真計算。通過將該軟件與COMSOL 對比，驗證了其實體單元能帶結構計算能力、結構單元能帶結構計算能力、大規模結構能帶結構計算能力。

1 有限元能帶結構分析的基本理論

聲子晶體能帶結構體現了彈性波在聲子晶體介質中的傳播能力。通過將波矢遍歷不可約布里淵區的邊界，計算得到對應波長的波矢在聲子晶體中的傳播頻率，將其繪制為波矢頻率關系曲線即為聲子晶體的能帶結構，又稱色散關系曲線。

波在理想無限周期聲子晶體中的傳播可以看作是一種平穩狀態，Bloch 證明了當簡諧波在介質中傳播時介質位移符合簡諧波位移場u(r,t)。

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

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

顯然，用于模態分析的結構線性有限元分析中的剛度矩陣K與質量矩陣M為實對稱矩陣。Bloch 邊界條件中顯含復系數，將原實廣義特征值問題變為復數矩陣廣義特征值問題。同時，由于約束的施加方法為酉合同變換，可通過Hermitian矩陣得定義公式證明約束后的矩陣K? 、M?為Hermitian矩陣，符合Hermitian 矩陣的定義式，即：

因此，聲子晶體能帶結構分析，經過上述公式推導，可整理為Hermitian 矩陣的廣義特征值問題。能帶結構分析往往只關心低階的前幾十階特征對，而非全部特征對，需通過稀疏Hermitian 矩陣得廣義特征值求解算法對其求解。

2 Hermitian 矩陣子空間迭代算法

Hermitian 矩陣是一類復共軛對稱矩陣。結構有限元分析中常用的實廣義特征值算法無法直接用于復Hermitian 矩陣的廣義特征值問題的求解。本節在實數子空間迭代法的基礎上發展了用于Hermitian 矩陣子空間迭代算法。

子空間迭代法是一類典型的特征值迭代求解方法，由BATHE[18]為求解結構動力學問題引入力學領域。子空間迭代法代表了求解大型稀疏矩陣特征值的經典思路，即通過多向量的逆迭代方法構造模態變換基底，使用正交的變換基底構建一個隨著子空間迭代步數增加逐步逼近原問題特征值的特征子空間。

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

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

算法的計算細節描述如下：

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

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

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

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

并且，由于子空間變換矩陣Qn、特征向量Rn塊內正交性，一組正交向量經過正交變換所產生向量是必然正交的。因此，在Hermitian 矩陣的廣義特征值求解中，同樣無需做塊內正交化。將Hermitian矩陣的共軛轉置退化為轉置，則可以推導得到用于實數子空間迭代法的歸一化公式。

3 定位格節點匹配方案

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

3.1 定位格方案描述

圖2 定位格節點關系Fig. 2 sub-patches-node relationship

[]int為取整操作。將主節點使用節點排序后的序號整理為相應的定位格-節點關系。當計算x、y、z維度的周期邊界條件時，令該周期維度的定位格數Nx、Ny、Nz等于1 即可。

3) 建立主節點和定位格的關聯信息數組。該方法中需要針對主節點構造2 類數組：長度為定位格數nb的定位格數組與長度為主節點個數nm的節點數組。定位格數組記錄了當前定位格第一個節點的序號，主節點數組記錄了除第一個節點之外的所有節點序號。

4) 按照如下流程將定位格中的節點依次插入定位格數組和節點數組：

a) 將定位格中第一個節點序號放入定位格數組的相應位置。

b) 將該定位格中的第二個節點序號放入第一個節點序號為序號的節點數組中。

c) 對第二個節點之后的節點循環，將后一個節點序號放入以前一個節點序號為序號的節點數組中。

d) 最后將沒有對應節點序號放入的定位格數組與節點數組的值設置為 ?1。

圖3 中展示了以2 號定位格的數據為例的節點數組與定位格數組的構建順序。至此，主節點定位格的兩類數組已經構建完成。

圖3 主節點和定位格的關聯信息數組Fig. 3 Association information array of master node and sub-patches

5) 將從節點與主節點進行匹配。先使用式(25)計算從節點與定位格的匹配關系，再順序對節點數組和定位格數組中的節點按照插入順序進行遍歷，直至找到相應匹配節點或遍歷至節點序號為?1時無匹配節點停止匹配。

3.2 節點匹配算法討論

空間節點匹配問題是多體動力學、計算機視覺、邊界元方法等問題的核心問題。其中比較典型的為：快速多極邊界元問題的自適應樹結構[20]，和本文引入的定位格節點匹配。這兩類匹配算法的核心思想是簡單且相同的，即通過增加中間層來降低節點匹配的計算量。同時，二者依托數據結構的不同造成它們搜索效率的不同。

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

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

4 聲子晶體能帶分析軟件模塊實現

本節使用前文描述的節點匹配方法及廣義特征值求解算法，設計實現聲子晶體的能帶分析軟件模塊org.sipesc.fems.phononic。聲子晶體的能帶結構分析可歸納為3 個有別于傳統有限元分析的步驟，即：

a)匹配主從節點，處理約束關系；

b)施加Bloch 邊界條件；

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

SiPESC 科學計算平臺已經擁有完善的線性有限元分析功能，只需在SiPESC 平臺上實現節點匹配和Hermitian 特征值求解，即可實現聲子晶體的能帶分析功能。

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

圖4 SiPESC 平臺能帶結構分析流程Fig. 4 band structure analysis process of SiPESC platform

圖5 中org.sipesc.core 軟件包為SiPESC 科學計算平臺的基礎軟件包，包含平臺所需的各類基礎功能、擴展所需繼承的基類、代數方程組求解器與其他矩陣計算類。org.sipesc.core.engdbs 為平臺的數據庫軟件包，負責數據庫相關的IO 操作和數據庫相關類型的創建與修改。org.sipesc.fems.elements 為平臺的單元軟件包，包含各類實體單元、結構單元、和用于自由度控制的單元控制類。

圖5 SiPESC 科學計算平臺聲子晶體軟件模塊的UML 類圖Fig. 5 UML class diagram of phononic crystal plugin in SiPESC platform

4.1 約束自由度處理模塊

MDofUpdate 類繼承于SiPESC 平臺的線性計算流程類MTaskManager，擁有初始化方法initialize和計算流程方法Start。該類所實現的功能是將單胞的邊界節點劃分為界面主節點與界面從節點，對主從節點根據對應的周期維度進行匹配，得到一一對應的主從節點關系。最后將節點自由度劃分為主從自由度，用以后續計算的自由度排序處理。

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

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

圖6 三類約束幾何關系及其關系樹Fig. 6 Three kinds of constrained geometric relations and their relation trees

點約束同線約束在處理方法上是一致的，僅僅將2 層主從約束變為3 層主從約束，且各頂點無需進行節點匹配。

4.2 能帶分析模塊

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

圖7 能帶結構分析算法流程Fig. 7 Band structure analysis algorithm

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

針對Hermitian 矩陣廣義特征值求解問題抽象出子空間迭代計算類SubSpaceIter。該類實現了第二節中的Hermitian 子空間迭代廣義特征值算法。通過該類求解得到當前波矢約束條件下的廣義特征值。

最后，針對結果輸出操作抽象出結果輸出對象exportorManager，該類使用exportUnvModal 方法調用私有節點數據輸出方法exportUnvNodePart，私有單元數據輸出方法exportUnvElePart，私有模態數據輸出方法exportUnvModePart，將計算結果輸出為SiPESC 平臺的后處理文件類型*.unv。

4.3 能帶分析算法討論

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

SiPESC 平臺與其他商用有限元軟件也有顯著區別。這些結構有限元軟件無法施加復數的邊界約束條件，只能使用實邊界條件施加約束。將剛度矩陣K劃分為實部和虛部：

進行求解。這一算法將矩陣維度增廣為2 倍，其求解的同階特征對也會變為2 倍。例如，使用SiPESC 平臺求解10 階模態可獲得10 階特征對，使用上文所述辦法只能求解5 階特征對，每階特征對出現2 次。由于增加求解的特征對數量時，特征值算法的計算量并非線性增長，因此二者同時計算10 階特征值時SiPESC 計算效率遠高于使用同類商業軟件。

5 數值算例

5.1 計算環境說明

為驗證所開發數值仿真軟件模塊的計算能力，在本節中分別針對實體單元組成的聲子晶體模型，結構單元組成的聲子晶體模型及大規模聲子晶體模型進行了數值算例驗證。選取文獻中應用較多的多物理場數值仿真軟件COMSOL Multiphysics 5.5 軟件作為標準對照軟件。

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

表1 材料數據Table 1 Material data

為了驗證本文所實現的軟件效率，選取了3 例彈性波聲子晶體算例進行測試，因此在COMSOL計算中選用了固體力學模塊，施加Floquet 周期性條件進行模型建立。為保證計算網格的一致性使用COMSOL 軟件進行建模，MSC.Patran 前處理軟件進行網格顯示。在折疊結構算例中，額外使用了固體力學中殼體單元的物理場模塊，以進行結構單元模擬。在第5 節中所使用的計算環境都統一遵從5.1 節的說明，不再額外在具體算例描述中敘述。

5.2 經典局域共振聲子晶體算例

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

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

表2 材料數據Table 2 Material data

表3 幾何尺寸Table 3 Geometric dimension

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

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

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

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

模態振型進行對比，圖9 中振型圖中展示了兩組振型結果。在每一組中左圖為COMSOL 軟件計算結果，右圖中為SiPESC 軟件計算結果。通過振型對比兩軟件計算振型一致。SiPESC 后處理顯示時將二階單元處理為線性單元顯示，在計算中仍為二階單元計算，特此注明。

該算例證明，本文所實現軟件較商業軟件具有一致的計算精度和較高的計算效率，基本實現了高效的聲子晶體能帶結構仿真功能。

5.3 大規模固/氣聲子晶體算例

為測試本能帶結構分析軟件的大規模計算性能，使用2 組元60 萬自由度的氣/固聲子晶體模型能帶結構[22]，固體材料為樹脂，材料的物理性能見表4，幾何結構如圖10 所示，幾何尺寸見表5。共計算20 階模態。

表4 材料數據Table 4 Material data

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

表5 幾何尺寸Table 5 Geometric dimension

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

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

該算例證明，本文所實現軟件具有大規模聲子晶體能能帶結構分析計算能力，適用于大規模聲子晶體計算。

5.4 折疊結構聲子晶體算例

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

表6 材料數據Table 6 Material data

表7 幾何尺寸Table 7 Geometric dimension

圖12 折疊結構聲子晶體模型Fig. 12 Folded structure phononic crystal model

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

圖13 折疊結構聲子晶體能帶結構Fig. 13 Band structure of folded phononic crystals

顯然，該折疊結構不存在完全帶隙，RM、RY方向能帶結構相似、能帶集中，符合各向異性周期結構的特性。該結構具有復雜的極化特性，已有諸多研究學者進行過討論。

本節算例1 與算例3 與COMSOL 進行了求解效率對比，其結果為，COMSOL 軟件的求解時間是SiPESC 平臺計算時間的2 倍以上，驗證了SiPESC 平臺能帶分析軟件具有高效的仿真效率。求解算例1、3 時的內存負載，二者大致相當都在10 GB以下，在求解大規模模型算例2 時SiPESC 共使用內存30 GB。證明SiPESC 平臺在內存充足的情況下有能力求解更大的單胞結構。精度二者大致相當，在計算結構單元時SiPESC 有明顯優勢。

6 結論

本文基于SiPESC 平臺，實現了聲子晶體高效能帶結構分析軟件，分別實現了用以廣義特征值分析的Hermitian 矩陣的子空間迭代法，與用以進行主從節點匹配的定位格匹配方案。

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

(2)針對能帶計算中的節點匹配問題，將用于碰撞檢測的定位格匹配方案，應用于聲子晶體能帶結構分析的主從約束匹配。

(3)基于SiPESC 工程與科學計算平臺有限元分析模塊，編寫了用于聲子晶體能帶結構仿真的軟件模塊。并介紹了該軟件模塊的架構設計和多層主從約束處理等設計細節。

(4)通過將軟件模塊與經典多物理場仿真軟件COMSOL Multiphysics 進行對比，驗證了該軟件模塊所實現算法具有可靠的數值精度和可觀的計算效率。在此基礎上計算了66 萬自由度的氣固聲子晶體模型、100 個離散波矢點的20 階模態自振分析，進一步驗證了本文所實現軟件具有高效的大規模分析能力。