范宣華,于晨陽，王柯穎，方 葉，肖世富

(中國工程物理研究院 總體工程研究所， 四川 綿陽 621900)

大規(guī)模單點基礎(chǔ)激勵隨機(jī)振動分析及并行計算

范宣華,于晨陽，王柯穎，方 葉，肖世富

(中國工程物理研究院 總體工程研究所， 四川 綿陽 621900)

基礎(chǔ)激勵作用下的隨機(jī)振動分析為結(jié)構(gòu)動力學(xué)有限元計算中的一種主要分析類型，在實際結(jié)構(gòu)分析中得到了廣泛應(yīng)用。基于模態(tài)疊加法，對基礎(chǔ)運(yùn)動激勵下的隨機(jī)振動核心算法理論進(jìn)行了系統(tǒng)推導(dǎo)?；贘AUMIN框架并行數(shù)據(jù)結(jié)構(gòu)，在PANDA平臺下完成了大規(guī)模算法設(shè)計和并行程序研發(fā)。結(jié)合數(shù)值算例對隨機(jī)振動響應(yīng)分析模塊進(jìn)行了對比驗證和大規(guī)模并行可擴(kuò)展性測試。研究結(jié)果表明：所研發(fā)的隨機(jī)振動分析軟件與商業(yè)軟件的計算結(jié)果完全吻合，在并行求解能力上可達(dá)到上億自由度、上萬核，遠(yuǎn)超出商業(yè)有限元軟件的計算水平。

單點基礎(chǔ)激勵；并行計算；隨機(jī)振動；PANDA平臺；并行可擴(kuò)展性；模態(tài)疊加法

隨機(jī)振動是自然界和工程結(jié)構(gòu)分析中常見的一類振動。諸如大型建筑物因地震、車輛因路面不平等引起的振動都屬于隨機(jī)振動分析的范疇。隨機(jī)振動不同于其他確定性振動，通常沒有固定周期，無法用確定性函數(shù)進(jìn)行描述，但具有一定的統(tǒng)計規(guī)律，一般可以用功率譜密度等統(tǒng)計量進(jìn)行描述[1]。

大型復(fù)雜裝備或建筑物經(jīng)常會受到來自約束部位的一致基礎(chǔ)激勵作用，有限元隨機(jī)振動分析多以平穩(wěn)隨機(jī)振動為前提，主要采用功率譜密度對結(jié)構(gòu)的響應(yīng)進(jìn)行描述[2]。在目前的商業(yè)有限元軟件中，串行分析求解能力基本維持在百萬至數(shù)百萬自由度量級，即便是國內(nèi)部分開放的商業(yè)有限元并行軟件，受技術(shù)禁運(yùn)(開放核數(shù)最多在數(shù)百個以內(nèi)，而且價格異常昂貴)和多年來串行編程機(jī)制下形成的內(nèi)核影響，僅計算規(guī)模略有增加，但并行性能往往難以提升，極大制約了復(fù)雜結(jié)構(gòu)的精細(xì)數(shù)值模擬水平[3]。

并行計算是解決這類大規(guī)模有限元問題的一個重要途徑。作為近20年來隨著計算機(jī)硬件快速發(fā)展而形成的一個熱門研究方向，并行計算對提升計算分析的規(guī)模和精度、縮短數(shù)值模擬時間具有重要意義。本文以單點基礎(chǔ)激勵隨機(jī)振動分析為研究對象，對基于自主并行軟件平臺PANDA[4]的大規(guī)模有限元并行計算研究進(jìn)行介紹，以期為大規(guī)模精細(xì)有限元分析提供借鑒。

1 基礎(chǔ)激勵隨機(jī)振動理論

隨機(jī)振動分析的輸入一般為作用于結(jié)構(gòu)約束部位的載荷功率譜密度曲線，載荷可以是位移、速度、加速度等基礎(chǔ)激勵。本文從最為常用的基礎(chǔ)加速度激勵入手，對隨機(jī)振動的核心理論進(jìn)行推導(dǎo)。

隨機(jī)振動算法理論的核心就是推導(dǎo)基礎(chǔ)激勵自功率譜密度與結(jié)構(gòu)關(guān)注點自功率譜密度響應(yīng)之間的傳遞關(guān)系。對于基礎(chǔ)運(yùn)動情形，以結(jié)構(gòu)和地基的相對位移xr作為變量，可建立如下多自由度系統(tǒng)的運(yùn)動方程：

(1)

式(1)中xr為各節(jié)點自由度對應(yīng)的相對位移向量，與絕對位移x以及基礎(chǔ)運(yùn)動位移u之間滿足如下關(guān)系：

x=xr+du

(2)

對于運(yùn)動方程(1)，采用模態(tài)疊加法[5]進(jìn)行方程解耦。首先進(jìn)行模態(tài)分析，與式(1)對應(yīng)的廣義特征值方程為

Mφ=ω2Kφ

(3)

i=1,2,…,m

(4)

對于式(4)，在頻域內(nèi)的穩(wěn)態(tài)解可以寫為

(5)

其中Hi(ω)為第i階模態(tài)對應(yīng)的頻響函數(shù)，即

(6)

對于平穩(wěn)隨機(jī)振動，根據(jù)功率譜密度定義，相應(yīng)的模態(tài)功率譜密度可以表示為：

(7)

式(7)給出了基礎(chǔ)運(yùn)動加速度激勵自譜密度和模態(tài)功率譜密度之間的遞推關(guān)系。根據(jù)相對位移和模態(tài)位移的關(guān)系，可以得到有限元模型中第k個自由度相對位移功率譜密度Sxrkxrk和基礎(chǔ)激勵加速功率譜密度之間的關(guān)系：

φikφjkSηiηj=

(8)

其中φik為第i階振型在第k個自由度處的取值，上標(biāo)*表示共軛。結(jié)合虛擬激勵法[6]思想，可對式(8)做進(jìn)一步簡化，得到

(9)

再根據(jù)相對位移和絕對位移的關(guān)系表達(dá)式(2)，得到第k個自由度的絕對位移自功率譜密度：

(10)

式中dk為方向向量d在第k個自由度處的取值。式(9)(10)分別給出了第k個自由度相對位移和絕對位移自功率譜密度與基礎(chǔ)加速度激勵自功率譜密度之間的遞推表達(dá)式。根據(jù)求解得到的位移自譜密度，可以直接得到速度和加速度自功率譜密度，三者之間滿足如下轉(zhuǎn)換關(guān)系：

(11)

2 算法設(shè)計

根據(jù)本文的討論和理論推導(dǎo)，可以建立如下單點基礎(chǔ)加速度激勵作用下的隨機(jī)振動分析算法：

算法1 基礎(chǔ)加速度激勵單點隨機(jī)振動分析算法輸入:結(jié)構(gòu)有限元模型,模態(tài)分析階數(shù)m,模態(tài)阻尼比ξi,加速度激勵功率譜密度曲線S¨u¨u(ω),激勵方向;輸出:指定節(jié)點自由度的位移、速度或加速度功率譜密度。1.根據(jù)結(jié)構(gòu)有限元模型,在JAUMIN框架下進(jìn)行離散,并行生成分布式質(zhì)量矩陣M和剛度矩陣K;2.根據(jù)質(zhì)量矩陣和剛度矩陣構(gòu)造廣義特征值問題,進(jìn)行模態(tài)分析,獲取的前m階特征對(ωi,?i);3.根據(jù)基礎(chǔ)運(yùn)動作用方向在所有非約束自由度上構(gòu)造單位方向向量d;4.結(jié)合整體質(zhì)量矩陣、振型和方向向量求解基礎(chǔ)激勵對應(yīng)的各階模態(tài)參與系數(shù);5.以指定計算的節(jié)點自由度作為第一重循環(huán),以加速度輸入自譜曲線離散后的頻率作為第二重循環(huán): a.獲取每個計算自由度在各階模態(tài)振型中的取值和對應(yīng)方向向量中的值; b.計算模態(tài)坐標(biāo)系下的各階頻響函數(shù)Hi(ω)(式(6)); c.按照式(9)和式(10)進(jìn)行模態(tài)疊加,分別計算相對位移和絕對位移功率譜密度曲線; d.根據(jù)需要按照式(11)分別計算速度和加速度的自譜密度曲線。6.功率譜密度曲線計算結(jié)果輸出。

在以上算法設(shè)計中，基礎(chǔ)激勵輸入可以是位移或速度的功率譜密度曲線，三者輸入之間同樣滿足式(11)的換算關(guān)系。此外，總體質(zhì)量矩陣和剛度矩陣的生成、模態(tài)分析過程、響應(yīng)求解等多個環(huán)節(jié)均可借助JAUMIN框架和PANDA平臺進(jìn)行并行求解。

3 并行實現(xiàn)

3.1 JAUNMIN框架和PANDA平臺

以上單點基礎(chǔ)激勵隨機(jī)振動的并行實現(xiàn)主要借助中國工程物理研究院自主研發(fā)的JAUMIN并行計算框架[7]和PANDA平臺實現(xiàn)。JAUMIN是根據(jù)超大機(jī)群硬件結(jié)構(gòu)特點面向大規(guī)模非結(jié)構(gòu)網(wǎng)格計算而研發(fā)的并行計算框架，提供基本的底層并行數(shù)據(jù)結(jié)構(gòu)和矩陣向量操作等數(shù)學(xué)運(yùn)算，負(fù)責(zé)不同數(shù)據(jù)塊之間的并行通信和數(shù)據(jù)管理等操作，同時為各類應(yīng)用軟件提供各類接口。有關(guān)JAUMIN的詳細(xì)介紹參見文獻(xiàn)[7]。

PANDA是筆者團(tuán)隊基于JAUMIN框架并行數(shù)據(jù)結(jié)構(gòu)研發(fā)的結(jié)構(gòu)力學(xué)有限元并行分析軟件平臺。PANDA平臺基本架構(gòu)以及與JAUMIN框架之間的關(guān)系如圖1所示。

圖1 PANDA平臺的基本結(jié)構(gòu)

目前整個PANDA平臺包括靜力學(xué)、模態(tài)和振動、沖擊動力學(xué)、多物理場耦合分析等多個有限元分析模塊，主要采用 C++ 和MPI編寫語言，包含代碼10萬余行。對于復(fù)雜工程結(jié)構(gòu)，JAUMIN框架結(jié)合前后處理軟件完成結(jié)構(gòu)的有限元建模和網(wǎng)格區(qū)域分解，PANDA平臺根據(jù)結(jié)構(gòu)材料屬性和外載荷條件等完成矩陣并行組裝和求解。單點基礎(chǔ)激勵隨機(jī)振動分析只是模態(tài)和振動分析軟件中的一個基本分析類型，此外還包含模態(tài)分析、地震響應(yīng)譜分析、諧響應(yīng)分析以及多點隨機(jī)振動分析等多個動力學(xué)分析類型。

3.2 并行實現(xiàn)概述

首先，根據(jù)工程結(jié)構(gòu)特點建立結(jié)構(gòu)有限元模型。在基于JAUMIN框架的并行計算中，多采用自主前處理軟件Supermesh進(jìn)行建模，此外也支持商業(yè)有限元軟件建模方式。對于建立的有限元模型，JAUMIN框架采用圖剖分功能，將有限元網(wǎng)格進(jìn)行區(qū)域分解，分成多個子區(qū)域。在區(qū)域分解時綜合應(yīng)用負(fù)載平衡技術(shù)，將有限元模型網(wǎng)格信息均勻分配到各個CPU計算節(jié)點，在各個CPU內(nèi)，PANDA平臺將結(jié)合模型本構(gòu)和材料、邊界等物理參數(shù)，并行生成分布式質(zhì)量矩陣和剛度矩陣。

其次，利用生成的質(zhì)量矩陣和剛度矩陣，PANDA平臺調(diào)用模態(tài)分析模塊進(jìn)行模態(tài)分析并行計算，獲取模態(tài)疊加所需的固有頻率和模態(tài)振型等。模態(tài)分析在整個動力學(xué)分析過程中是最耗時間和資源的環(huán)節(jié)，其計算能力也基本決定了后續(xù)振動分析的計算能力。由于已經(jīng)在PANDA平臺下進(jìn)行了相應(yīng)并行求解的集成實現(xiàn)，故目前支持Krylov-Schur算法[8]和Jacobi-Davidson算法[9]等開展大規(guī)模模態(tài)分析并行求解[4,10]。

再次，根據(jù)模態(tài)分析計算結(jié)果，結(jié)合JAUMIN框架提供的矩陣向量等并行操作運(yùn)算，獲取各階模態(tài)參與系數(shù)。

最后，將模態(tài)參與系數(shù)、計算自由度對應(yīng)的各階振型值、模態(tài)固有頻率和模態(tài)阻尼比等作為響應(yīng)計算輸入，設(shè)計相應(yīng)的模態(tài)疊加C++類，完成從基礎(chǔ)激勵到結(jié)構(gòu)計算自由度直接頻響函數(shù)的求解，并按照式(9)(10)設(shè)計隨機(jī)振動分析C++類，實現(xiàn)相對位移和絕對位移自功率譜密度曲線的求解和計算結(jié)果輸出。

除矩陣組裝和模態(tài)分析階段的并行環(huán)節(jié)以外，在隨機(jī)振動分析實現(xiàn)中的并行環(huán)節(jié)主要有2個：① 求解模態(tài)參與系數(shù); ② 計算各自由度的自譜密度曲線。在求解模態(tài)參與系數(shù)時，主要借助模態(tài)分析的質(zhì)量矩陣以及模態(tài)振型在各個進(jìn)程的分布式數(shù)據(jù)進(jìn)行并行求解，而在計算各個節(jié)點自由度的自功率譜密度時，其計算過程是一個天然并行模式，每個CPU進(jìn)程只負(fù)責(zé)本進(jìn)程內(nèi)節(jié)點自由度的自功率譜密度計算，相鄰進(jìn)程之間不需要任何數(shù)據(jù)通信。

4 數(shù)值算例

為驗證研發(fā)模塊的正確性和并行可擴(kuò)展性，本節(jié)以光機(jī)靶球隨機(jī)振動分析作為算例開展并行計算研究。光機(jī)靶球的有限元模型如圖2所示。采用四面體單元進(jìn)行劃分，初始模型自由度數(shù)為340萬，模態(tài)分析提取前100階模態(tài)，采用Jacobi-Davidson算法求解，各階模態(tài)阻尼比取0.01，在支架4個腳底施加沿水平方向的基礎(chǔ)加速度激勵，激勵曲線為0～50 Hz范圍內(nèi)的白直譜(加速度功率譜密度為1(m/s2)2/Hz的直線譜)。計算靶球上方桿件最上端的水平方向的位移功率譜密度。

圖2 光機(jī)靶球有限元模型

對于340萬自由度初始模型，采用商業(yè)有限元軟件ANSYS進(jìn)行對比驗證，ANSYS和PANDA計算得到的前100階固有頻率保持小數(shù)點后3位有效數(shù)字一致。圖3給出了靶球頂端關(guān)注點的絕對位移功率譜密度對比情況，可以看出，ANSYS和PANDA計算得到的功率譜密度曲線幾乎完全重合，驗證了PANDA單點基礎(chǔ)激勵程序模塊的正確性。

為進(jìn)一步驗證PANDA單點激勵隨機(jī)振動分析的并行可擴(kuò)展性，在百萬億次大型機(jī)群上對靶球340萬自由度模型進(jìn)行了自適應(yīng)網(wǎng)格加密后的隨機(jī)振動計算，經(jīng)過1次網(wǎng)格自適應(yīng)加密后達(dá)到 1 700 萬自由度，經(jīng)過2次加密后達(dá)到1.3億自由度。對于如此上千萬乃至上億自由度的有限元模型，已超出國內(nèi)通用商業(yè)有限元軟件的計算能力，目前只能借助并行計算軟件完成。

圖3 絕對位移自功率譜密度曲線對比

對于靶球模型加密一次后得到的1 700萬自由度規(guī)模，CPU核數(shù)從64個一直測試到 8 192個，得到不同并行CPU核數(shù)上的并行計算時間和并行效率，如表1所示?？梢钥闯觯簩τ谠撚嬎阋?guī)模，在8 192核內(nèi)，隨著核數(shù)的增加，計算時間持續(xù)下降，并未出現(xiàn)計算拐點。8 192核相對于64核的并行效率為24%，說明PANDA平臺下的單點隨機(jī)振動程序具有優(yōu)異的并行可擴(kuò)展性。

表1 靶球1 700萬自由度模型隨機(jī)振動分析計算時間和并行效率

對于加密2次得到的1.3億自由度規(guī)模，分別進(jìn)行1 024核、2 048核、4 096核、8 192核以及10 000核的并行可擴(kuò)展性測試，均成功算出。相應(yīng)計算時間情況如表2所示?？梢钥闯觯涸?8 192核內(nèi)隨著核數(shù)的增加，計算時間近乎線性減少，8 192核較1 024核的并行效率提高達(dá)84%，表現(xiàn)出非常優(yōu)異的并行可擴(kuò)展性；超過8 192核后，由于大型機(jī)群系統(tǒng)自身網(wǎng)絡(luò)分組配置原因?qū)е驴缬虿⑿型ㄐ牛⑿行视兴陆?，出現(xiàn)計算拐點，但萬核級相對于千核級并行效率仍高達(dá)60%以上。為了更直觀地表示分析模塊的并行可擴(kuò)展性，繪制1 700萬自由度和1.3億自由度的并行加速比曲線，如圖4所示。

表2 靶球1.3億自由度模型隨機(jī)振動分析計算時間和并行效率

圖4 兩種規(guī)模的計算加速比曲線

從表1～ 2以及圖4可以看出：

1) 計算規(guī)模越大，PANDA隨機(jī)振動分析功能模塊的并行可擴(kuò)展性越優(yōu)異。這主要是因為隨著模型的增大，每個節(jié)點計算區(qū)域內(nèi)的計算自由度數(shù)較邊界通信部分的自由度數(shù)大大增加，通信帶來的計算量相對變小。同理，對于同一計算規(guī)模，隨CPU核數(shù)的增加，各區(qū)域帶來的并行通信量增大，并行效率也逐漸下降。

2) 對于光機(jī)靶球之類復(fù)雜的算例，千萬自由度計算時間從64核的26 h減少到8 192核的1 h左右。上億自由度規(guī)模從1 024核的40余h減少到上萬核的6 h左右。這相對于串行計算而言是幾乎無法實現(xiàn)的，并行計算則不但大幅提升了計算規(guī)模，還大大縮短了數(shù)值模擬的計算時間，凸顯了并行計算的優(yōu)越性；

3) 本文的數(shù)值算例表明隨機(jī)振動分析模塊具備“上億自由度、上萬核”的并行可擴(kuò)展能力。

5 結(jié)束語

本文基于模態(tài)疊加法，利用中國工程物理研究院自主研發(fā)的JAUMIN框架和PANDA平臺，對單點基礎(chǔ)激勵隨機(jī)振動分析進(jìn)行了算法設(shè)計和并行實現(xiàn)，結(jié)合光機(jī)靶球數(shù)值算例驗證了隨機(jī)振動分析模塊的正確性和并行可擴(kuò)展性，取得了商業(yè)有限元分析軟件無法達(dá)到的大規(guī)模并行計算能力。

本文研發(fā)的單點基礎(chǔ)激勵隨機(jī)振動分析模塊僅是PANDA平臺下模態(tài)和振動分析軟件中的一個基本分析模塊，旨在說明基于JAUMIN框架研發(fā)的PANDA平臺具有超強(qiáng)的并行可擴(kuò)展性。在單點激勵基礎(chǔ)上，筆者所在團(tuán)隊近期已將隨機(jī)振動分析從單點分析擴(kuò)展到了多點激勵分析，并將分析規(guī)模進(jìn)一步提升到了10億自由度以上，限于篇幅，在此不做深入介紹。

[1] CHRISTIAN L..Random Vibration:Mechanical Vibration and Shock Analysis[C]//3rdedition.Published in Great Britain and the United States by ISTE Ltd and John Wiley & Sons,2014.

[2] CHOPRA A K.Dynamics of Structures-Theory and Applications to Earthquake Engineering[M]//4th edition.New Jersey:Prentice Hall,2012.

[3] FAN X,WANG K,XIAO S,et al.Some progress on parallel modal and vibration analysis using the JAUMIN framework[J].Mathematical Problems in Engineering,2015(2):1-8.

[4] 范宣華.基于Panda框架的大規(guī)模有限元模態(tài)分析并行計算及應(yīng)用[D].北京：北京大學(xué),2013.

[5] ITOH T.Damped vibration mode superposition method for dynamic response analysis[J].Earthquake Engineering & Structure dynamics,1973,2(1):47-57.

[6] 林家浩，張亞輝.隨機(jī)振動的虛擬激勵法[M].北京：科學(xué)出版社,2006.

[7] LIU Q K,ZHAO W B,CHENG J,et al.A programming framework for large scale numerical simulations on unstructured mesh[C]//IEEE 2nd International Conference on High Performance and Smart Computing.New York:[s.n.]， 2016:310-315.

[8] STEWART G W.A Krylov-Schur algorithm for large eigenproblems[J].SIAM Journal on Matrix Analysis and Applications,2001,23(3):601-614.

[9] SLEIJPEN G L G，VORST H A V.A Jacobi-Davidson iteration method for linear eigenvalue problems[J].SIAM Journal on Matrix Analysis and Applications,1996,17(2):401- 425.

[10] FAN X H,CHEN P,WU R,et al.Parallel computing study for the large-scale generalized eigenvalue problems in modal analysis[J].Science China Physics,Mechanics and Astronomy,2014,57(3):477- 489.

(責(zé)任編輯楊黎麗)

ParallelComputationofLarge-ScaleRandomVibrationAnalysisUnderSingle-PointMotion-BasedExcitation

FAN Xuanhua， YU Chenyang, WANG Keying, FANG Ye, XIAO Shifu

(Institute of Systems Engineering, China Academy of Engineering Physics, Mianyang 621900, China)

The random vibration analysis under the basic excitation is a kind of main analysis type in the finite element calculation of structural dynamics, and it has been widely used in practical structural analysis. Based on the modal superposition method, the random vibration algorithm theories under the basic motion excitation were systematically deduces. Based on the parallel data structure of JAUMIN framework, large-scale algorithm design and parallel program development were carried out under the PANDA platform. Numerical examples are given to verify the rightness as well as the large-scale parallel scalabilities of the random vibration response analysis. The results show that the random vibration analysis software developed by us is in good agreement with the calculated results of commercial software, and can reach hundreds of millions of degrees of freedom, tens of thousands of CPU cores in the parallel solution scalability, going far beyond the calculation level of commercial finite element software.

single-point base excitation; parallel computation; random vibration; PANDA platform; parallel scalability; mode superposition method

2017-06-21

國家自然科學(xué)基金面上資助項目(11472256)；科技部“高性能計算”重大專項課題(2016YFB0201005)；國防基礎(chǔ)科研計劃項目(C1520110002)；中國工程物理研究院院長基金、院發(fā)展基金和雙百人才基金資助項目(YZ2015011，2014B0202025，ZX04003)

范宣華(1981—),男,山東日照人,博士,研究員,主要從事結(jié)構(gòu)動力學(xué)并行計算研究,E-mail：fanxh@caep.cn。

范宣華,于晨陽，王柯穎，等.大規(guī)模單點基礎(chǔ)激勵隨機(jī)振動分析及并行計算[J].重慶理工大學(xué)學(xué)報(自然科學(xué))，2017(10):56-61,89.

formatFAN Xuanhua，YU Chenyang, WANG Keying, et al.Parallel Computation of Large-Scale Random Vibration Analysis Under Single-Point Motion-Based Excitation[J].Journal of Chongqing University of Technology(Natural Science)，2017(10):56-61,89.

10.3969/j.issn.1674-8425(z).2017.10.009

TB132

A

1674-8425(2017)10-0056-06