千萬自由度量級有限元模態分析并行計算研究

第一作者范宣華男，博士，高級工程師，1981年生

千萬自由度量級有限元模態分析并行計算研究

范宣華1,肖世富1,陳璞2

(1. 中國工程物理研究院總體工程研究所，四川綿陽621999; 2.北京大學工學院力學與工程科學系，北京100871)

摘要：大規模有限元模態分析在一些重大裝備研制過程中有迫切需求,對于實現裝置系統級分析具有重要意義?；陔[式重啟動Arnoldi、Krylov-Schur和Jacobi-Davidson三種主流算法和PANDA并行計算框架，構建了大規模模態分析并行計算體系；將并行求解體系應用于某光機主體結構，實現了其上千萬自由度、數千核的模態分析并行計算；結合算例對三種主流算法的適應性和并行可擴展性進行了評估。研究結果表明，基于三種算法構建的并行求解體系均可在1小時內求解千萬自由度量級的大規模模態分析問題,并行可擴展性非常優異。

關鍵詞：有限元法; 大規模模態分析; 并行計算

基金項目：國家自然科學基金面上項目(11472256)；中國工程物理研究院發展基金項目(2012A0202008);中國工程物理研究院重點學科項目“計算固體力學”;中國工程物理研究院院長基金項目

收稿日期：2014-05-07修改稿收到日期：2014-08-27

中圖分類號:TB123; TP311

文獻標志碼：A

DOI:10.13465/j.cnki.jvs.2015.17.013

Abstract:In the development of large equipments, the demand of large-scale finite element modal analysis is very strong due to its significance in realizing the systemic analysis of the entire large structure. Based on three predominant algorithms, i.e., implicitly restarted Arnoldi method, Krylov-Schur method and Jacobi-Davidson method and PANDA framework, a parallel computing system for a large-scale modal analysis was established. As a typical application, it was used to get solutions to the main structure of a certain optical machinery, a parallel modal analysis with over ten-million-DOF was performed and thousands of CPU processors were applied. The adaptability and parallel scalability of the three algorithms with numerical examples. Results showed that the parallel solving system can perform the modal analysis with over ten million-DOF within one hour and its parallel scalability is excellent.

Parallel computing for finite element modal analysis with over ten-million-DOF

FANXuan-hua1,XIAOShi-Fu1,CHENPu2(1. Institute of Systems Engineering, CAEP, Mianyang 621900, China;2. Department of Mechanics and Engineering Science, College of Engineering, Peking University, Beijing 100871, China)

Key words:finite element method; large-scale modal analysis; parallel computation

在一些特殊大型裝備如巨型光機裝置、復雜武器系統以及大型飛行器等的研制過程中，往往既需要關注整體動力學特性，也需要關注局部，涉及大規模復雜動力學系統求解，這對傳統有限元計算方法和分析工具形成挑戰。一方面傳統的簡化建模會影響結構局部關鍵細節的預測能力；另一方面大型裝備結構復雜，往往包含大量的密集模態，需要精細建模才能得以反映。目前國內現有的商業有限元軟件模態分析能力維持在100萬自由度量級、數小時的計算水平，無法滿足大型裝備系統級高精度數值分析的需求。這樣以來，高效的大規模有限元計算就顯得特別重要。

在動力學有限元分析過程中，模態分析不僅是最耗時間的計算環節，而且也是其他動力分析的基礎，其求解精度和效率將直接影響到后續動力分析的可行性和可靠性[1-2]。隨著超級計算機的快速發展和大量先進算法的涌現，開展大規模模態分析并行計算的時機也逐漸成熟。本文主要從特征值算法、模態分析并行計算體系構建以及某光機裝置大型裝備應用等方面介紹我們所取得的一些進展。

1大規模特征值問題求解算法

模態分析在數學上可以歸結為大型稀疏矩陣特征值問題，對于該問題的研究大多基于子空間投影技術[3-4]。依據子空間序列產生的方式不同，子空間類方法大致可以分為兩類——固定維數子空間法和變化維數子空間法。總體說來變化維數比固定維數的方式更為高效，最具代表性的是Krylov類子空間方法和Davidson類子空間方法，以下分別進行簡述。

1.1Krylov類子空間方法

Krylov類子空間方法可以追溯到20世紀50年代提出的用于對稱特征值求解的Lanczos算法和用于非對稱特征值求解的Arnoldi算法[2, 5]，基本思想是通過在Krylov子空間上建立正交基并利用這組正交基將矩陣A縮減成Hessenberg矩陣形式。后來很多學者對基本的Arnoldi/Lanczos算法進行了一系列重啟動改進[6-7]，比較著名的是Sorensen提出的隱式重啟動Arnoldi/Lanczos算法[6]，其具體代碼實現于ARPACK軟件包[7]，已經得到廣泛應用。Stewart[8]在隱式重啟動Arnoldi/Lanczos算法基礎上進行了改進，將Arnoldi分解代之以更為寬泛的Krylov分解，提出了更容易代碼實現的Krylov-Schur算法。隱式重啟動Arnoldi/Lanczos算法和Krylov-Schur算法在數學上具有等價性，是目前Krylov類子空間算法中的兩個主流數值算法。

對模態分析而言，其數學本質是求解由質量矩陣M和剛度矩陣K組成的廣義特征值方程Kx=λMx，所關心的是低階特征對。而Krylov類算法多收斂于最大特征值，為此需要進行譜變換。常用的譜變換有Shift-Invert變換和Cayley變換[9]。以Shift-Invert變換為例，對應的變換形式為：

(1)

式中:σ為移位值。通過譜變換可以將原廣義特征值問題可以化為標準特征值問題Asx=μsx，其中As=(K-σM)-1M，μs=1/(λ-σ)。在擴展Krylov子空間時，只需要將矩陣A替換成平移-逆處理矩陣As即可，由此帶來的問題就是在每個迭代步中都需要求解一個線性方程組，即在增加Krylov子空間序列時，w=Avi變成了w= (K-σM)-1Mvi，要得到w，即相當于求解線性系統(K-σM)w=Mvi。該大型線性方程組的求解是整個特征值算法的關鍵，一般采用矩陣分解的直接法進行求解，矩陣分解將占據整個特征值問題求解的大部分時間和存儲空間。

1.2Davidson類子空間方法

Davidson類子空間方法最初源于由Davidson[10]提出的一種求解對角占優的對稱矩陣特征值問題算法，后來逐漸發展成求解一般矩陣特征值問題的廣義Davidson算法[1,3,11]。該類算法的一個重大進展是由Sleijpen和Vorst[12]共同提出的Jacobi-Davidson算法。

Jacobi-Davidson方法的核心在于計算“修正方程”，對于廣義特征值問題，其修正方程形式為：

(I-MuuT)(K-θM) (I-uuTM)t=-r

(2)

式中:θ為某階近似特征值(里茲值)，u為和θ相對應的近似特征向量(里茲向量)，r為留數向量，t為待求的正交補向量。和Krylov類子空間算法的一個重要不同在于，Jacobi-Davidson方法通過得到t的近似解來擴展子空間，而Krylov類子空間方法通過將算子A的冪形式施加到初始迭代向量上來擴展子空間。Jacobi-Davidson方法采用內-外層迭代方式求解特征系統，外層迭代擴展子空間，內層迭代求解“修正方程”。該算法的優勢在于計算過程中涉及的運算主要是矩陣-向量積和向量內積，原特征矩陣的稀疏特性得以保留，對內存需求低，并行通訊少；但對不同問題的適應性不如采用Krylov類子空間方法。

2大規模模態分析并行求解體系

圍繞隱式重啟動Arnoldi、Krylov-Schur和Jacobi-Davidson三種主流算法，我們基于團隊研發的PANDA框架[13-14]構建了一套模態分析并行求解體系。Panda框架旨在為大規模并行計算力學應用提供基礎支撐服務組件和集成開發環境。該框架基于層次化、模塊化及面向對象的設計思路，相應的層次結構見圖1。處于底層的基礎支撐層形成了Panda框架的核心數據結構，除數學解法器本身的并行機制外，其他并行通訊操作主要在這一部分定義和實現；中間一層為數值共性層，該部分主要提供具有通用性的一些數學計算服務，可根據問題特點和不同需要進行選擇使用；最頂層為應用個性層，包含一些單元庫、材料模型庫及其接口服務等?；谠摽蚣芸梢詾楦鞣N有限元計算提供并行分布式數據結構，一方面通過Metis區域分解并行生成有限元計算所需質量和剛度矩陣，另一方面通過集成或研發相應解法器實現大規模有限元并行計算。

圖1　Panda框架的基本結構 Fig.1 The basic configuration of PANDA frame

利用PANDA框架提供的并行數據結構以及接口等功能形成的模態分析求解體系描述見圖2。整個體系分成MSC.Patran前處理建模，PANDA框架下的并行計算以及計算結果的后處理三個部分。具體介紹如下：

圖2　模態分析并行計算體系 Fig.2 Parallel computing system for modal analysis

(1)前處理建模。模態分析的前處理借助MSC.Patran實現。在MSC.Patran中，根據不同單元類型、不同材料和邊界約束等對工程結構的有限元模型進行分組處理，然后將這些單元節點信息以MSC.Patran 的中性文件(neutral)形式進行輸出，為PANDA并行計算提供基本的有限元模型數據。目前采用MSC.Patran可以在普通工作站上建立千萬自由度以上的有限元模型。

(2)PANDA并行計算。該部分是整個模態分析體系的核心，一方面，利用PANDA框架的接口功能將MSC.Patran的網格模型轉換成PANDA本身支持的網格模型，通過將美國圣地亞實驗室開發的有限元數據庫Exodus II接入，實現了模型文本文件到exo格式二進制網格文件的高效轉換；另一方面，結合PANDA集成的圖剖分軟件Metis[15]對有限元網格進行區域分解，利用不同子區域生成質量和剛度矩陣的并行分布式數據結構。生成矩陣后，通過編制的解法器接口調用相應的解法器進行模態分析求解。目前我們在Panda框架下集成了PETSc[16]、SLEPc[17]、ARPACK[7]以及MUMPS[18]等多個數值軟件包來實現隱式重啟動Arnoldi、Krylov-Schur和Jacobi-Davidson三種主流算法的模態分析并行求解。其中PETSc和MUMPS主要負責矩陣向量的基本運算以及線性方程組的并行求解；SLEPc負責Krylov-Schur算法和Jacobi-Davdison算法的具體實現；APPACK負責隱式重啟動Arnoldi算法的具體實現。同時為保證大規模模態分析的收斂性以及計算效率，在各算法和軟件中添加了譜變換、內部大型線性方程組求解技術以及內外層優化等關鍵求解技術，算法詳細設計和實現過程可以參見文獻[14]，限于篇幅在此不過多闡述。

(3)模態分析結果的后處理。主要是模態頻率、振型和求解過程中的一些統計信息顯示。

3某光機裝備典型應用

基于搭建的模態分析并行求解體系，我們開展了大量的工程結構的模態分析計算研究，并對計算軟件的正確性與商業軟件進行了大量的對比驗證，相關工作可參見文獻[14]。在此主要給出我國某光機主體結構的模態分析應用。大型光機裝置的研制對解決人類未來能源問題及開拓前沿科學技術研究具有重要意義，該裝置定位精度和穩定性設計等方面的苛刻要求對大規模數值模擬提出了嚴峻挑戰，數值模擬必須充分反映裝置的結構細節響應，由此導致整體結構的有限元計算規模至少要達到數千萬自由度以上。

本文所建立的光機主體結構有限元模型見圖3，該模型包含了靶球裝置、混凝土樓層、鋼架等不同結構形式和材料類型，幾乎囊括了大部分工程結構類型，圖2所示的有限元模型基本能夠反映結構本身的復雜特性，此外考慮到脫密等因素，在建模時對裝置各部分進行了必要的等效處理，并對各部分材料屬性做了相應改變。整個有限元模型采用六面體單元進行劃分，約束底部地基部分三個方向自由度。整體結構被劃分成百萬自由度到千萬自由度三種規模，三種情形的相關統計信息見表1。所有規模均求解前50階模態頻率。所有算例的并行測試均在某百萬億次大型集群上進行，該集群包含了數千個計算節點，每個計算節點含兩個六核CPU，12個核共享48 GB內存。

圖3　某光機主體結構有限元模型 Fig.3 FEM of Shenguang Structure

測試規模描述自由度規模剛度矩陣非零元個數平均帶寬Case11,236,98172,678,84359Case210,018,296713,214,22571Case320,189,3071,497,573,99974

Case 1的小規模算例主要用來對三種算法的準確性進行驗證和比較。首先對三種算法的求解精度采用了留數誤差和相對誤差雙重誤差進行評判，這兩個誤差能夠較為充分地反映出各階頻率和振型的求解精度。相應的誤差定義如下所示：

εm=‖Kxm-λmMxm‖2(留數誤差)

(4)

(5)

三種算法求解前50階模態頻率所對應的留數誤差和模態誤差見圖4。

圖4　三種算法對應的模態誤差和留數誤差(Case 1) Fig.4 The modal errors and residue errors for the three methods (Case 1)

從圖4中可以看出隱式重啟動Arnoldi算法和Krylov-Schur算法各階頻率所對應的留數誤差在10-6～10-4之間，模態誤差在10-10～10-7之間；而Jacobi-Davidson算法的模態誤差和留數誤差相對較低，在10-4～10-2之間，這主要是受Jacobi-Davidson算法本身影響，該算法通過內-外層迭代近似求解，為保證內層迭代的收斂性并盡可能減少迭代步數，降低了迭代控制精度。

本文中所采用的三種算法前50階模態頻率計算結果幾乎完全一致，三種算法間的相對誤差在0.05%以內。此外，還將Case1算例采用相同的單元類型在ANSYS中采用塊Lanczos方法進行了對比計算，結果發現模態頻率誤差也在0.05%以內，各階振型保持一致。

對于Case2的1000萬自由度規模，首先在512個CPU核內對三種算法進行了測試，測試過程發現隱式重啟動Arnoldi算法(IRA)和Krylov-Schur算法(K-S)的內存占用大致為400GB，Jacobi-Davidson算法(J-D)的內存占用大致為100GB。三種算法的并行計算時間和效率如表2所示，相應的加速比曲線見圖5。

圖5　三種算法的加速比曲線(Case 2) Fig.5 The speedup curves for the three algorithms (Case 2)

圖6　J-D算法測試得到的加速比曲線 Fig.6 The speedup curves for the J-D algorithms

對于Case3的2 000萬自由度模型，隱式重啟動Arnoldi算法(IRA)和Krylov-Schur算法(K-S)因矩陣分解失敗而無法計算，僅測試了Jacobi-Davidson算法，連同Case2，我們將最大并行CPU核數擴展到了4 800個。表3給出了部分測試結果，相應的加速比曲線見圖6。

通過這些測試結果，可以看出以下幾點：

(1)對于Krylov-Schur算法和隱式重啟動Arnoldi算法(Case2)，在400個CPU核附近出現計算了拐點，即隨著核數增加，并行計算時間不再降低。這主要是因為兩種Krylov類子空間算法在每個子空間迭代過程中都采用直接法求解一個大型線性方程組，矩陣分解帶來了較大的內存占用和較多的并行通訊，導致其并行可擴展性相對較弱；而Jacobi-Davidson算法采用內-外層迭代的方式求解特征對，并行通訊相對較少，使得其并行可擴展性非常優異，在512個CPU核內的并行加速比解決理想加速比，且未出現計算拐點。

(2)對于Jacobi-Davidson算法，在Case2的1000萬自由度規模時，并行計算拐點出現在2500個CPU核附近；在Case3的2 000萬自由度模型時，并行計算拐點出現在3 500個CPU核附近，兩種規模的可擴展CPU核數均達到數千個。

(3)綜合比較三種算法在各種并行進程下的實際計算時間可以看出，Jacobi-Davidson算法在少量CPU核數的計算時間高出Krylov類兩種算法一個數量級以上，隨著CPU核數的增加，三種算法的計算時間均縮減在1h內。這主要是在少量CPU核數時，Krylov類子空間算法在線性方程組求解時以存儲空間換時間，能快速準確得到擴展子空間向量，體現出一定優勢；而Jacobi-Davidson算法的內層迭代使得其獲取子空間向量的精度較弱，迭代也相對緩慢，導致整個外層迭代收斂較慢，但隨著CPU核數的增加，其優異的并行性能使得計算時間迅速下降，最終和Krylov類子空間算法求解時間相當。

表2　Case 2模型在某大型集群上的計算時間與效率比較

表3　Case 2和Case 3兩種規模通過J-D算法測試得到的計算時間和并行效率

此外，為揭示該復雜裝置大規模計算的必要性，我們對Case1～Case3的模態頻率結果進行了比較分析，表4給出了三種規模前10階模態頻率的對比情況。眾所周知，對于工程結構的有限元分析，存在剛度矩陣的“硬化”效應，這會導致在較小自由度規模(有限元網格較稀疏)計算時得到的模態頻率偏高，隨著網格加密(矩陣規模增大)模態頻率會逐漸接近真實頻率。我們從Case 1到Case 3的計算結果可以明顯發現這一現象，其中其中Case 2的變化率相對于Case 1而言，Case 3的變化率相對于Case 2而言。從表4可以看出，從Case 1的123萬自由度規模到Case 2的1000萬自由度規模，各階模態頻率發生了比較明顯的變化，最大變化率超過10%；而從Case 2的1 000萬自由度規模到Case 3的2 000萬自由度規模，各階模態頻率的變化率相對低的多。這說明對于光機結構這樣的大型復雜結構，百萬自由度的計算規模根本無法滿足其真實模態頻率預測的精度要求，相應的計算規模至少要達到千萬自由度以上。

表4　光機結構三種不同規模前10階

4結論

本文針對大型復雜裝備的模態分析問題，通過對三種主流特征值算法進行研究和代碼集成開發，在PANDA并行計算框架基礎上構建了一套大規模模態分析并行求解體系，該體系覆蓋Krylov子空間和Davidson兩大類主流算法，并且均具備上千萬自由度規模的并行計算能力，基本可以滿足大型裝備模態分析的需求。

數值算例結果表明，Krylov類子空間算法與Jacobi-Davidson算法具有很好的互補性， Krylov子空間類算法采用直接矩陣分解求解線性方程組，在少量CPU核時會凸顯出計算時間上的優勢，較快地得到計算結果，但內存占用相對較多，并行可擴展性不如Jacobi-Davidson算法；而Jacobi-Davidson算法采用內-外層迭代方式求解，內存占用比較少，并行可擴展性非常優異，但內層迭代過程會使得少量CPU核時的計算時間遠高于Krylov類子空間算法。

參考文獻

[1]Saad Y. Numerical methods for large eigenvalue problems[M]. Halsted Press, Div. of John Wiley &Sons, Inc., New York:1992.

[2]尹家聰. 模態分析的自動多重子結構法與動力重分析研究[D].北京：北京大學，2012.

[3]Bai Z, Demmel J, Dongarra J, et al. Templates for the solution of algebraic eigenvalue problems: A practical guide[M]. SIAM, Philadelphia, 2000.

[4]Watkins D S. The matrix eigenvalue problem: GR and krylov subspace methods[M]. SIAM, Philadelphia, 2007.

[5]Arnoldi W E. The principle of minimized iterations in the solution of the matrix eigenvalue problem[J]. Quart. Appl. Math., 1951(9):17-29.

[6]Sorensen D C. Implicit application of polynomial filters in a k-step Arnoldi method[J]. SIAM J. Matrix Anal. Appl., 1992(13):357-385.

[7]Lehoucq R B, Sorensen D C, Yang C. ARPACK Users’ guide: solution of large-scale eigenvalue problems by implicitly restarted arnoldi methods[M]. SIAM, Philadelphia, PA,1998.

[8]Stewart G W. A Krylov-Schur algorithm for large eigenproblems[J]. SIAM J. Matrix Anal. Appl., 2001,23(3): 601-614.

[9]Meerbergen K, Spence A, Roose D. Shift-invert and cayley transforms for detection of rightmost eigenvalues of nonsymmetric matrices[J]. BIT Numerical Mathematics, 1994,3(34): 409-423.

[10]Davidson E R. The iterative calculation of a few of the lowest eigenvalues and corresponding eigenvectors of large real symmetric matrices[J]. J. Comp. Phy., 1975(17): 87-94.

[11]Crouzeix M, Philippe B, Sadkane M. The davidson method[J]. SIAM J. Sci. Comput., 1994(15):62-67.

[12]Sleijpen G L G, van der Vorst H A. A Jacobi-Davidson iteration method for linear eigenvalue problems[J]. SIAM J. Matrix Anal. Appl., 1996, 17(2):401-425.

[13]Shi G M, Wu R, Wang K Y. Discussion about the design for mesh data structure within the parallel framework[C]. 9th World Congress on Computational Mechanics and 4th Asian Pacific Congress on Computational Mechanics, Sydney, Australia, 2010.

[14]范宣華. 基于Panda框架的大規模有限元模態分析并行計算及應用[D].北京：北京大學, 2013.

[15]Karypis G， Kumar V. Metis-a software package for partitioning unstructured graphs, partitioning meshes, and computing fill-reducing orderings of sparse matrices[R]. Technical Report, University of Minnesota, Minneapolis, 1998.

[16]Balay S, Buschelman K, Eijkhout V,et al. PETSc users manual[R]. Technical Report ANL-95/11-Revision 3.0.0, Argonne National Laboratory, 2008.

[17]Hernandez V, Roman J, Tomas A, et al. SLEPc: a scalable and flexible toolkit for the solution of eigenvalue problems[J]. ACM Trans. Math. Softw., 2005, 31(3):351-362.

[18]Amestoy P R, Duff I S, L’Excellent J Y, et al. MUMPS: A general purpose distributed memory sparse solver[J]. Lecture Notes in Computer Science, 2001(1947):121-131.