大規(guī)模多點基礎(chǔ)激勵隨機(jī)振動并行計算研究

寧民亮， 范宣華*， 王柯穎， 陳 璞

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

1 引 言

在對航空航天、大型軍事和民用設(shè)施等復(fù)雜裝備進(jìn)行可靠性評估時，有限元模型的自由度動輒數(shù)千萬甚至數(shù)億，需要使用具備超大規(guī)模計算能力的軟件來分析整個系統(tǒng)的結(jié)構(gòu)動力學(xué)特性。當(dāng)前市面上廣泛使用的有限元軟件如ANSYS和ABAQUS等，由于底層框架的限制，其自由度計算規(guī)模僅能達(dá)到百萬量級，且相關(guān)技術(shù)受到嚴(yán)格禁運，難以滿足對復(fù)雜裝備進(jìn)行高效數(shù)值模擬的迫切需求[1,2]。

并行計算是解決超大型裝備動力學(xué)分析等挑戰(zhàn)問題的關(guān)鍵技術(shù)之一[3,4]。2004年，美國發(fā)布白皮書將高性能計算作為提高和維系其制造業(yè)市場競爭力的王牌[5]。近年來并行技術(shù)不斷得到發(fā)展，國外已研發(fā)出多款具備高性能并行計算能力的軟件。如美國圣地亞實驗室基于SIERRA框架研發(fā)的結(jié)構(gòu)動力學(xué)模塊SIERRA/SD[6]，具備模態(tài)和振動分析常用類型的超大規(guī)模計算能力，最大自由度計算規(guī)模已經(jīng)可以上億，但相關(guān)軟件代碼對國內(nèi)完全禁運；西班牙瓦倫西亞大學(xué)結(jié)合PETSc數(shù)據(jù)框架研發(fā)的SLEPc軟件[7]，集成了大量的特征值并行計算功能，可實現(xiàn)不同行業(yè)的模態(tài)分析大規(guī)模并行計算能力。該軟件重點針對結(jié)構(gòu)動力學(xué)的模態(tài)分析功能，不具備其他結(jié)構(gòu)動力學(xué)分析功能。

面對技術(shù)封鎖，國內(nèi)近些年也在積極推動自主力學(xué)分析軟件的研發(fā)工作，國家多個部委也設(shè)置了軟件自主化項目予以傾斜資助。國內(nèi)相關(guān)研究以緊跟國際前沿、集成開發(fā)創(chuàng)新為主，但研發(fā)進(jìn)度和計算能力與發(fā)達(dá)國家還存在較大差距。如張林波等[8]基于PHG平臺開發(fā)的并行自適應(yīng)結(jié)構(gòu)分析模塊PHG-Solid，具備數(shù)億自由度以上的大規(guī)模自適應(yīng)網(wǎng)格劃分能力以及靜力學(xué)線性方程組求解能力，在結(jié)構(gòu)動力學(xué)方面支持模態(tài)分析并行計算，但尚不具備其他結(jié)構(gòu)動力學(xué)分析能力；徐良寅等[9]研發(fā)的SiPESC軟件，采用平臺+插件的運行架構(gòu)，方便用戶在此框架上進(jìn)行二次開發(fā)，其中SiPESC.FEMS模塊具備模態(tài)分析、頻響分析和隨機(jī)振動分析等功能，在多個領(lǐng)域結(jié)構(gòu)分析方面發(fā)揮了重要作用，目前SiPESC的程序運行主要以串行模式為主，最大分析能力仍受到一定限制。

中國工程物理研究院從2007年開始進(jìn)行結(jié)構(gòu)力學(xué)大規(guī)模并行計算研究，基于院內(nèi)JAUMIN框架[10]自主研發(fā)形成了并行非結(jié)構(gòu)網(wǎng)格有限元計算平臺PANDA[11]。并在該計算平臺上完成了模態(tài)分析、諧響應(yīng)分析和響應(yīng)譜分析等多個動力學(xué)分析模塊的開發(fā)和應(yīng)用研究[12,13]。本文在已有研究基礎(chǔ)上，針對多點基礎(chǔ)激勵隨機(jī)振動分析問題，開展了理論推導(dǎo)、算法設(shè)計以及模塊研發(fā)等相關(guān)工作。最終測試結(jié)果表明，研發(fā)的多點隨機(jī)振動分析模塊擁有商業(yè)軟件不具備的超大規(guī)模運算能力。

2 相關(guān)理論和求解優(yōu)化

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

多點基礎(chǔ)激勵結(jié)構(gòu)動力學(xué)方程[14]為

(1)

式中M，C和K為結(jié)構(gòu)自由節(jié)點的質(zhì)量陣、阻尼陣和剛度陣，MC，CC和KC為支承和結(jié)構(gòu)耦聯(lián)部分的質(zhì)量陣、阻尼陣和剛度陣，MR，CR和KR為支承約束節(jié)點的質(zhì)量陣、阻尼陣和剛度陣，u為結(jié)構(gòu)的絕對位移，可以分解成準(zhǔn)靜態(tài)位移us和動力相對位移ud兩部分，uR為支承的絕對位移，R為支撐反力。

結(jié)合文獻(xiàn)[13]相關(guān)推導(dǎo)過程，式(1)最終可表示為準(zhǔn)靜態(tài)響應(yīng)方程和動態(tài)響應(yīng)方程兩部分。

(2)

式中uR,l為第l個基礎(chǔ)激勵位移，al為位移影響系數(shù)向量，表示在第l處約束沿激勵方向施加單位位移引起的結(jié)構(gòu)自由節(jié)點各自由度的位移向量，N為基礎(chǔ)激勵的個數(shù)。

對大規(guī)模裝備進(jìn)行動力學(xué)并行計算時，一般采用模態(tài)疊加法進(jìn)行分析，以避免直接矩陣分解帶來災(zāi)難性存儲和大量并行通信。用模態(tài)分析獲得的m階模態(tài)特征對{(ω1,φ1)，(ω2,φ2),…,(ωm，φm)}進(jìn)行解耦，可得式(2)動態(tài)部分對應(yīng)的模態(tài)解耦方程為

(3)

(4)

式中Sp q(ω)為第p處和第q處基礎(chǔ)加速度激勵對應(yīng)的功率譜密度，Sd,k(ω)，Ss,k(ω)和St,k(ω)分別對應(yīng)第k個自由度的動態(tài)位移自功率譜密度、靜態(tài)位移自功率譜密度和總位移功率譜密度，Sd s,k(ω)為動態(tài)和準(zhǔn)靜態(tài)位移功率譜密度的交叉項，Hj(w)為式(3)對應(yīng)的頻響函數(shù)，φj,k為第j階模態(tài)對應(yīng)的第k個元素值，ap,k為第p處單位基礎(chǔ)位移引起的第k個自由度的準(zhǔn)靜態(tài)位移值。式(4)中Hj(w)和Sd s,k(ω)的具體表達(dá)式為

(5)

式(4)建立了位移功率譜密度表達(dá)式，相應(yīng)的速度和加速度自功率譜密度分別為位移自功率譜密度的ω2和ω4倍。獲得相應(yīng)關(guān)注點的響應(yīng)功率譜密度曲線后，對曲線對應(yīng)頻段進(jìn)行積分開方后可獲得對應(yīng)關(guān)注點的等效均方根響應(yīng)。

2.2 隨機(jī)振動應(yīng)力分量和Mises等效應(yīng)力求解

隨機(jī)振動的Mises等效應(yīng)力求解相對于靜動力學(xué)的確定性問題求解更為困難，這是因為載荷輸入給出的是統(tǒng)計特性，而且響應(yīng)本身也是統(tǒng)計性數(shù)據(jù)[15]。此外，Mises應(yīng)力是各應(yīng)力分量的非線性函數(shù)，前面計算位移、速度、加速度以及應(yīng)力分量響應(yīng)的理論方法無法直接應(yīng)用于Mises應(yīng)力的計算[16]。有限元分析中，節(jié)點的應(yīng)力分量可以寫成矩陣形式σ={σx，σy，σz，τx y，τy z，τx z}T，而Mises應(yīng)力σM與應(yīng)力分量σ的關(guān)系滿足：

(6)

式中矩陣B的表達(dá)式為

(7)

對于第k個節(jié)點，結(jié)合文獻(xiàn)[15，16]處理方法，可以給出與式(4)對應(yīng)的動態(tài)和準(zhǔn)靜態(tài)Mises應(yīng)力自功率譜密度為

(8)

2.3 求解過程的優(yōu)化

將N個基礎(chǔ)載荷的功率譜密度以譜矩陣形式給出

(9)

式中每個元素代表一條譜曲線，對角線元素對應(yīng)各基礎(chǔ)激勵的自譜，非對角線元素對應(yīng)不同基礎(chǔ)激勵之間的互譜，不相關(guān)載荷之間的互譜為0。在前述核心理論公式的計算過程中，不考慮零元素，只針對非零元素進(jìn)行循環(huán)，可減少計算量。

在對全場均方根響應(yīng)進(jìn)行求解時，按照均方根響應(yīng)的定義，需求解完全部自由度對應(yīng)的自功率譜密度曲線后，再對自功率譜密度曲線進(jìn)行積分求解。以動態(tài)位移自功率譜密度為例，其均方根響應(yīng)表達(dá)式為

(10)

對于n自由度系統(tǒng)，p個頻率離散點，全場均方根響應(yīng)對應(yīng)的浮點運算量為(n×m×m×p×N×N)，如此計算需要的浮點運算量非常巨大，甚至超出模態(tài)分析計算時間的數(shù)倍之多。優(yōu)化過程中發(fā)現(xiàn)，均方根積分求解的模態(tài)振型與積分頻率無關(guān)，因而可以將積分項進(jìn)行分離變換。以動態(tài)位移響應(yīng)為例，將式(10)第k個自由度的動態(tài)位移均方根積分公式變換為

(11)

式中Qi j實際是模態(tài)坐標(biāo)空間中的一個m×m階協(xié)方差矩陣對應(yīng)的第i行第j列的元素值。經(jīng)過處理變換后得到的式(11)，分別以二重求和作為基本運算單位取代式(10)的四重求和運算，只需要(n×m×m+p×N×N)次浮點運算即可。此外需要存儲一個模態(tài)空間下m×m階的模態(tài)協(xié)方差矩陣，這對于大規(guī)模計算來說占用的存儲量幾乎可以忽略。以上優(yōu)化過程可以理解為,將整體積分項中與頻率相關(guān)積分變量項和空間變量進(jìn)行分離，在較小的模態(tài)空間內(nèi)獲得積分結(jié)果，并作為模態(tài)貢獻(xiàn)在物理空間內(nèi)進(jìn)行疊加，從而將原來的多個乘法浮點運算轉(zhuǎn)換成兩個較小乘法浮點運算量之和，大大降低了計算量。

對于上述優(yōu)化，舉例說明如下，假設(shè)結(jié)構(gòu)自由度數(shù)n為10000，頻率離散點p為1000，多點激勵個數(shù)N為10個，模態(tài)參與階數(shù)m為100階，則優(yōu)化前的基本運算單位對應(yīng)的浮點計算量為1013，優(yōu)化后浮點計算量僅為108+105，減少的單位浮點計算量可以呈現(xiàn)接近5個量級的減少，即原來需要100000秒完成的計算，優(yōu)化后只需要1秒左右便可完成。

3 隨機(jī)振動的算法設(shè)計和并行實現(xiàn)

3.1 算法設(shè)計

根據(jù)第2節(jié)推導(dǎo)的理論公式，設(shè)計多點基礎(chǔ)激勵隨機(jī)振動分析的算法如下。

(1) 建模和矩陣離散。簡化實際工程結(jié)構(gòu)，建立有限元模型并劃分網(wǎng)格。在JAUMIN框架下進(jìn)行區(qū)域分解和必要的網(wǎng)格自適應(yīng)加密；在PANDA平臺下建立相應(yīng)的質(zhì)量陣和剛度陣。

(2) 獲取模態(tài)信息。建立廣義特征值問題并開展模態(tài)分析，獲得各階模態(tài)頻率和正交歸一化的模態(tài)振型，該部分工作已在PANDA中實現(xiàn)[11]。

(3) 求解位移影響系數(shù)向量al。借助PANDA靜力學(xué)并行分析模塊，分別求解N個不同約束沿激勵方向施加單位位移引起的結(jié)構(gòu)響應(yīng)，獲得一系列位移向量al。

(4) 計算振型參與系數(shù)γl j。利用質(zhì)量陣、模態(tài)振型和位移向量al計算不同模態(tài)坐標(biāo)系下的振型參與系數(shù)γl j。

(5) 計算自功率譜密度和全場均方根響應(yīng)。以關(guān)注節(jié)點的每個自由度為一次循環(huán)，以給定的頻譜曲線離散頻率為二次循環(huán)，進(jìn)行每個離散頻率下的自功率譜密度和均方根響應(yīng)計算。

① 利用曲線離散頻率及各階模態(tài)阻尼比，建立解耦方程，算出相應(yīng)的頻響函數(shù)。

② 根據(jù)功率譜密度、頻響函數(shù)值和振型參與系數(shù)，按照前面的理論公式計算動態(tài)位移、準(zhǔn)靜態(tài)位移以及總位移的自功率譜密度。

③ 按照模態(tài)位移振型獲得的應(yīng)力振型計算各節(jié)點應(yīng)力分量的自功率譜密度。

④ 計算模態(tài)坐標(biāo)系下的協(xié)方差矩陣，并根據(jù)推導(dǎo)的理論公式分別計算各部分位移(速度和加速度)的等效均方根響應(yīng)以及不同節(jié)點位置的各應(yīng)力分量和Mises應(yīng)力的等效均方根響應(yīng)。

(6) 輸出關(guān)注節(jié)點的自功率譜密度曲線以及全場等效均方根值。

以上的算法設(shè)計中，可以求解結(jié)構(gòu)部分關(guān)注點或全部節(jié)點自由度的自功率譜密度和響應(yīng)均方根值。各自由度的求解根據(jù)并行區(qū)域分解的結(jié)果可以在各自進(jìn)程內(nèi)分別求解，屬于天然并行模式。其中，區(qū)域分解基于JAUMIN框架完成，通過在不同子區(qū)域設(shè)置映像區(qū)的方式實現(xiàn)相鄰子區(qū)域之間的數(shù)據(jù)交換和通信，具體實現(xiàn)方式參見文獻(xiàn)[10]。

3.2 并行實現(xiàn)

本文開展的大規(guī)模多點基礎(chǔ)激勵隨機(jī)振動并行實現(xiàn)主要基于JAUMIN框架和PANDA平臺，采用C++和MPI編程完成。JAUMIN底層框架具有優(yōu)異的數(shù)據(jù)并行處理和矩陣運算能力，在并行實現(xiàn)過程中，搭配前后處理軟件，負(fù)責(zé)完成數(shù)據(jù)塊之間的通訊和管理、建模以及區(qū)域分解等工作；PANDA平臺則負(fù)責(zé)構(gòu)建相應(yīng)的質(zhì)量陣和剛度陣并組裝，然后調(diào)用研發(fā)的隨機(jī)振動分析模塊開展求解工作。關(guān)于JAUMIN框架和PANDA平臺的具體介紹可參見文獻(xiàn)[10,13]，在此不做贅述。

本文進(jìn)行的隨機(jī)振動相關(guān)算法設(shè)計是基于模態(tài)疊加法開展的，其總體架構(gòu)設(shè)計如圖1所示。隨機(jī)振動分析結(jié)合模態(tài)參與系數(shù)、頻響函數(shù)以及譜曲線輸入等，通過模態(tài)疊加計算響應(yīng)，然后導(dǎo)向輸出流。其中各階振型參與系數(shù)根據(jù)求解的準(zhǔn)靜態(tài)位移向量，并結(jié)合各階模態(tài)振型來確定；譜線管理器管理著多個自譜和互譜，為模態(tài)疊加過程提供譜曲線輸入。

圖1 多點隨機(jī)振動分析架構(gòu)設(shè)計

在多點基礎(chǔ)激勵隨機(jī)振動并行求解模塊中，首先通過JAUMIN框架進(jìn)行區(qū)域分解，構(gòu)建分布式質(zhì)量矩陣和剛度矩陣，不同區(qū)域?qū)?yīng)的矩陣數(shù)據(jù)發(fā)送到不同的CPU進(jìn)程，PANDA對這些分布式矩陣數(shù)據(jù)調(diào)用模態(tài)分析并行計算功能獲得模態(tài)頻率和振型等相關(guān)信息，之后結(jié)合載荷曲線和模態(tài)信息構(gòu)建模態(tài)參與系數(shù)，并在不同子區(qū)域內(nèi)進(jìn)行結(jié)構(gòu)隨機(jī)振動響應(yīng)的計算。隨機(jī)振動模塊主要包含以下四個C++類。

譜曲線管理類。負(fù)責(zé)處理多個不同基礎(chǔ)激勵的譜曲線，包括曲線的描述和離散等。曲線離散方式支持均勻離散和根據(jù)模態(tài)固有頻率及模態(tài)阻尼比的自適應(yīng)離散，這些離散結(jié)果在振動響應(yīng)求解時調(diào)用。

振型參與系數(shù)類。根據(jù)模態(tài)振型以及基礎(chǔ)載荷系數(shù)向量計算各階振型參與系數(shù)，這些振型參與系數(shù)以數(shù)組的形式存儲，在響應(yīng)求解時調(diào)用。

隨機(jī)振動響應(yīng)計算類。該類用于各關(guān)注點的PSD響應(yīng)曲線計算以及全場的等效均方根響應(yīng)計算，該類的輸入包括模態(tài)參與系數(shù)、頻響函數(shù)、各階模態(tài)振型和輸入載荷PSD曲線。

輸出類。該類用于關(guān)注點PSD曲線和全場等效均方根響應(yīng)的輸出。相關(guān)計算數(shù)據(jù)按照J(rèn)AUMIN框架后處理軟件規(guī)定的數(shù)據(jù)格式進(jìn)行輸出，便于曲線繪制和全場云紋圖顯示。

在整個分析流程中，網(wǎng)格區(qū)域的劃分、質(zhì)量陣和剛度陣的構(gòu)建以及模態(tài)分析和響應(yīng)求解是主要的并行階段。其中模態(tài)分析主要依靠PANDA平臺的模態(tài)分析功能實現(xiàn)求解，目前支持Krylov-Schur[17]與Jacobi-Davidson[18]兩種主流并行求解算法。

4 數(shù)值算例

4.1 Benchmark算例

采用不規(guī)則扇體作為Benchmark算例。扇面尺寸左邊寬為1.75 m，右邊低端寬為1.25 m，體厚度為0.6 m，有限元模型如圖2所示。模型全部以六面體單元劃分，采用單一材料，對應(yīng)的彈性模量為210 GPa，密度為7800 kg/m3，泊松比為0.3。將模型兩端面固定，作為施加加速度激勵的基礎(chǔ)，沿端面法向施加兩個大小不同的加速度激勵，頻率范圍為1 Hz～2000 Hz。激勵曲線采用白直譜，右端幅值為1，左端幅值為2。兩端激勵互譜實部和虛部幅值分別為0.7和1.2，取前20階模態(tài)進(jìn)行模態(tài)疊加。選取模型中部五角星標(biāo)記點作為關(guān)注點。

圖2 扇形體有限元模型

圖3和圖4是利用PANDA和ANSYS算出的標(biāo)記節(jié)點的位移自功率譜密度曲線對比情況?？梢钥闯觯还苁莿討B(tài)位移還是絕對位移，兩個軟件算出來的自功率譜密度曲線都基本重合。關(guān)注點的速度響應(yīng)以及加速度響應(yīng)曲線和位移曲線僅差一個頻率平方的倍數(shù)，兩個軟件計算結(jié)果亦完全一致，在此不再列出。

圖3 關(guān)注點動態(tài)位移自功率譜密度對比

圖4 關(guān)注點絕對位移自功率譜密度對比

圖5和圖6是PANDA和ANSYS關(guān)于全場等效位移云圖和全場Mises應(yīng)力云圖的對比。可以看出，兩種軟件關(guān)于等效均方根響應(yīng)的計算結(jié)果在數(shù)值分布、最大值以及最大值出現(xiàn)的位置都非常吻合。這說明前面推導(dǎo)的均方根響應(yīng)公式以及對應(yīng)的算法程序是正確的。計算時間方面，在同一配置的臺式機(jī)上，本算例PANDA串行計算的時間為1.05 s，ANSYS的運行時間為1.92 s。這是因為，PANDA軟件基于Linux系統(tǒng)安裝，運行過程可以充分利用機(jī)器內(nèi)存完成計算；ANSYS等商業(yè)有限元軟件基于Windows系統(tǒng)，在軟件運行過程中需要從硬盤讀取數(shù)據(jù)和存放中間變量信息，內(nèi)存利用率較低。這使得PANDA計算時間遠(yuǎn)低于ANSYS。隨著計算規(guī)模的增大，PANDA的計算時間優(yōu)勢較ANSYS更為明顯。

圖5 全場位移等效均方根響應(yīng)對比

圖6 Mises等效應(yīng)力均方根響應(yīng)對比

4.2 大型光機(jī)測試算例

任取模型上某一節(jié)點，圖8和圖9是PANDA和ANSYS關(guān)于該點Y方向動態(tài)位移和絕對總位移自功率譜密度曲線的對比?？梢钥闯觯瑑蓚€軟件計算的響應(yīng)曲線基本一致，再一次驗證了研發(fā)模塊的正確性。圖10給出了光機(jī)裝置全場Z向動態(tài)位移等效均方根位移云圖的計算結(jié)果。二者的云圖亦保持一致，略微差別主要來自不同的頻率離散策略帶來的積分誤差。

圖7 某大型光機(jī)有限元模型

圖8 動態(tài)位移自譜曲線對比

圖9 總位移自譜曲線對比

圖10 光機(jī)裝置Z向動態(tài)位移云紋圖對比

為測試多點隨機(jī)振動分析模塊的并行可擴(kuò)展性，對光機(jī)結(jié)構(gòu)進(jìn)行網(wǎng)格自適應(yīng)加密3次后，得到了11.88億自由度計算規(guī)模，并在天河超大并行機(jī)群上開展了多達(dá)上萬核的并行擴(kuò)展測試。表1展示了使用不同CPU核數(shù)對應(yīng)的計算時間和并行效率?？梢钥闯觯瑢τ?1.88億自由度，萬核相對模型能夠計算的2048核并行效率高達(dá)80%以上，說明了研發(fā)模塊具備優(yōu)異的并行可擴(kuò)展能力。表2 給出了光機(jī)裝置在10240個CPU核下各個主要求解環(huán)節(jié)的時間統(tǒng)計情況，可以看出,模態(tài)分析占據(jù)了整個求解時間的99%以上，隨機(jī)振動響應(yīng)求解只占了0.2%的時間，充分說明了2.3節(jié)中算法和模塊求解優(yōu)化帶來的巨大計算效益。

表1 光機(jī)裝置并行可擴(kuò)展性測試結(jié)果

表2 光機(jī)裝置萬核計算各主要環(huán)節(jié)時間統(tǒng)計

5 結(jié) 論

本文從隨機(jī)振動的基本理論入手，系統(tǒng)推導(dǎo)了基于模態(tài)疊加法的多點基礎(chǔ)激勵隨機(jī)振動分析理論體系，對求解流程進(jìn)行了必要的優(yōu)化，在此基礎(chǔ)上完成了相應(yīng)的算法設(shè)計和模塊研發(fā)。將算例結(jié)果與商業(yè)有限元軟件進(jìn)行對比，驗證了研發(fā)模塊的正確性；對該模塊進(jìn)行了并行可擴(kuò)展性測試，結(jié)果顯示該模塊具備優(yōu)異的并行計算可擴(kuò)展能力。

結(jié)合開展的算例測試結(jié)果可以看出，開發(fā)的多點基礎(chǔ)激勵隨機(jī)振動分析模塊能夠?qū)崿F(xiàn)位移、速度、加速度以及應(yīng)力等典型動力學(xué)響應(yīng)量的功率譜密度計算和等效均方根計算，并且可以達(dá)到與商業(yè)有限元軟件一致的計算精度，該模塊功能對單點基礎(chǔ)激勵的情形同樣適用。對于商業(yè)軟件無法計算的上億自由度規(guī)模，PANDA軟件通過并行求解技術(shù)完全可以勝任，目前測試得到的最大計算規(guī)模可以達(dá)到10億自由度以上，其并行計算可擴(kuò)展能力比現(xiàn)有商業(yè)軟件更強。