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

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

1 引 言

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

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

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

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

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

2.1 多點基礎激勵隨機振動理論

多點基礎激勵結構動力學方程[14]為

(1)

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

結合文獻[13]相關推導過程，式(1)最終可表示為準靜態(tài)響應方程和動態(tài)響應方程兩部分。

(2)

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

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

(3)

(4)

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

(5)

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

2.2 隨機振動應力分量和Mises等效應力求解

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

(6)

式中矩陣B的表達式為

(7)

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

(8)

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

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

(9)

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

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

(10)

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

(11)

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

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

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

3.1 算法設計

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

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

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

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

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

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

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

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

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

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

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

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

3.2 并行實現(xiàn)

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

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

圖1 多點隨機振動分析架構設計

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

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

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

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

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

在整個分析流程中，網(wǎng)格區(qū)域的劃分、質量陣和剛度陣的構建以及模態(tài)分析和響應求解是主要的并行階段。其中模態(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所示。模型全部以六面體單元劃分，采用單一材料，對應的彈性模量為210 GPa，密度為7800 kg/m3，泊松比為0.3。將模型兩端面固定，作為施加加速度激勵的基礎，沿端面法向施加兩個大小不同的加速度激勵，頻率范圍為1 Hz～2000 Hz。激勵曲線采用白直譜，右端幅值為1，左端幅值為2。兩端激勵互譜實部和虛部幅值分別為0.7和1.2，取前20階模態(tài)進行模態(tài)疊加。選取模型中部五角星標記點作為關注點。

圖2 扇形體有限元模型

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

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

圖4 關注點絕對位移自功率譜密度對比

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

圖5 全場位移等效均方根響應對比

圖6 Mises等效應力均方根響應對比

4.2 大型光機測試算例

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

圖7 某大型光機有限元模型

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

圖9 總位移自譜曲線對比

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

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

表1 光機裝置并行可擴展性測試結果

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

5 結 論

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

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