余丁浩， 李 鋼

(大連理工大學(xué) 建設(shè)工程學(xué)部海岸和近海工程國家重點(diǎn)實(shí)驗(yàn)室， 遼寧 大連 116024)

隨著工程結(jié)構(gòu)逐漸向大型化和復(fù)雜化方向發(fā)展，基于傳統(tǒng)變剛度法的結(jié)構(gòu)非線性地震反應(yīng)計(jì)算平臺(tái)所需巨大計(jì)算耗時(shí)已越發(fā)難以滿足工程建設(shè)需求，盡管目前計(jì)算機(jī)技術(shù)已獲得了長足發(fā)展，但大型結(jié)構(gòu)抗震分析所面臨的效率瓶頸依然是制約結(jié)構(gòu)工程發(fā)展的重要因素，如何在保證高精度計(jì)算模擬的同時(shí)有效提升結(jié)構(gòu)分析效率始終是工程界重點(diǎn)關(guān)注問題。

利用計(jì)算機(jī)的多線程計(jì)算能力或搭建集群化計(jì)算環(huán)境通過并行方式進(jìn)行結(jié)構(gòu)非線性地震反應(yīng)分析是提升計(jì)算效率的重要途徑，從編程角度來看可將當(dāng)前結(jié)構(gòu)分析領(lǐng)域常用的并行計(jì)算模式分為兩類：共享內(nèi)存式和分布內(nèi)存式。共享內(nèi)存式是指多個(gè)計(jì)算線程共用一個(gè)數(shù)據(jù)存儲(chǔ)空間，主要用于單機(jī)計(jì)算，如基于OpenMP的多核CPU(central processing unit)并行計(jì)算[1-2]，基于GPU的并行計(jì)算[3-11]等，具有維護(hù)和編程相對(duì)簡便的優(yōu)點(diǎn)；分布內(nèi)存式主要利用計(jì)算機(jī)集群實(shí)現(xiàn)高性能計(jì)算，不同計(jì)算設(shè)備之間通過必要的數(shù)據(jù)通訊進(jìn)行協(xié)調(diào)，此類并行計(jì)算通常采用消息傳遞的編程模式(message passing interface，MPI)實(shí)現(xiàn)，雖然其理論上可以實(shí)現(xiàn)超大規(guī)模問題的計(jì)算分析，但編程方式復(fù)雜、設(shè)備維護(hù)代價(jià)高，且不同計(jì)算設(shè)備之間無法避免的數(shù)據(jù)通訊導(dǎo)致計(jì)算效率存在瓶頸。在實(shí)際應(yīng)用中，上述兩種并行模式可以混合使用，以進(jìn)一步降低分析耗時(shí)，例如何強(qiáng)等[12]通過對(duì)整個(gè)求解域進(jìn)行剖分，將不同子域分配到不同計(jì)算設(shè)備上，并在單機(jī)上利用OpenMP進(jìn)行線程級(jí)別并行加速，提出了基于MPI+OpenMP混合編程模式的并行計(jì)算方法。此外，為充分發(fā)揮不同并行計(jì)算模式的效能，眾多學(xué)者有針對(duì)性的提出了多種便于并行實(shí)現(xiàn)的結(jié)構(gòu)分析方法，如單元接單元方法[13-16]、區(qū)域分解法[17-21]等，其中單元接單元方法不需總裝整體剛度矩陣和荷載向量，運(yùn)算主要集中在單元層面，該方法通常與預(yù)處理共軛梯度法相結(jié)合，多用于共享內(nèi)存式的并行計(jì)算；區(qū)域分解法采用分而治之的思想，將結(jié)構(gòu)劃分為不同的計(jì)算區(qū)域，其關(guān)鍵在于建立能夠?qū)崿F(xiàn)負(fù)載平衡的區(qū)域劃分方式和界面方程求解方法，多用于分布式計(jì)算集群。

結(jié)構(gòu)在地震作用下的非線性變形通常集中于局部區(qū)域，如部分梁柱端部、墻底等部位，已有研究表明，充分利用結(jié)構(gòu)的局部非線性特征有助于大幅提升非線性地震反應(yīng)分析效率，當(dāng)前已發(fā)展出多種高效局部非線性分析方法，如擬力法[22-26]、數(shù)值子結(jié)構(gòu)法[27-30]、多尺度類方法[31-32]、快速非線性分析法[33-34]等，這些方法多通過保持結(jié)構(gòu)剛度矩陣恒定或縮減問題規(guī)模的方式避免結(jié)構(gòu)大規(guī)模剛度矩陣實(shí)時(shí)更新及由此導(dǎo)致的計(jì)算效率降低。在現(xiàn)有的局部非線性分析方法中，基于Woodbury公式的結(jié)構(gòu)非線性分析方法[35-37](簡稱Woodbury方法)因具有較好的數(shù)學(xué)理論基礎(chǔ)和較廣泛的適用性，近年來得到了越來越多的關(guān)注，此類方法將任意時(shí)刻的結(jié)構(gòu)切向剛度改寫為初始彈性剛度的低秩攝動(dòng)形式，進(jìn)而引入線性代數(shù)領(lǐng)域中的Woodbury公式求解結(jié)構(gòu)切向位移反應(yīng)，能夠避免了大規(guī)模整體切向剛度的迭代更新，僅需對(duì)代表局部非線性行為的小規(guī)模矩陣進(jìn)行實(shí)時(shí)更新和分解，可顯著提升非線性問題計(jì)算效率。Deng等通過將結(jié)構(gòu)切向剛度中局部時(shí)變元素單獨(dú)取出，建立了以Woodbury公式為高效求解器的虛擬力分析法。近期Li等[38-40]以擬力法中的變形分解思想為出發(fā)點(diǎn)，通過對(duì)局部區(qū)域內(nèi)進(jìn)入塑性的材料進(jìn)行應(yīng)變分解，在單元層面構(gòu)造非線性應(yīng)變場函數(shù)，提出了隔離非線性有限元法基本理論，該方法采用Woodbury公式進(jìn)行控制方程求解，且具有適用性廣、通用性好等優(yōu)點(diǎn)，能夠?qū)崿F(xiàn)一般結(jié)構(gòu)的高效非線性地震反應(yīng)分析。隨后，眾多學(xué)者圍繞該方法的算法優(yōu)化、計(jì)算應(yīng)用等開展了系列研究[41-48]，進(jìn)一步豐富了其理論體系。盡管Woodbury方法能夠從算法理論層面解決傳統(tǒng)非線性分析方法由于剛度矩陣時(shí)變導(dǎo)致的效率瓶頸，但現(xiàn)有研究均采用串行計(jì)算模式，并未充分利用計(jì)算硬件的并行加速能力，因而其高效性優(yōu)勢并未得到充分發(fā)掘。

本文將OpenMP并行計(jì)算模式引入到Woodbury方法中，提出了一種基于Woodbury+OpenMP的高效結(jié)構(gòu)非線性地震反應(yīng)分析方法。首先采用隔離非線性法構(gòu)造結(jié)構(gòu)切向剛度方程的攝動(dòng)展開表達(dá)式及相應(yīng)Woodbury求解公式，隨后，分別針對(duì)Woodbury公式計(jì)算、單元狀態(tài)確定、局部非線性系數(shù)矩陣更新這3個(gè)非線性迭代求解過程中主要計(jì)算步驟建立了基于OpenMP的并行優(yōu)化策略，實(shí)現(xiàn)了該方法的全過程并行加速，最后，通過一個(gè)高層結(jié)構(gòu)算例驗(yàn)證了本文方法的高效性和準(zhǔn)確性。

1 Woodbury方法基本理論

本文采用隔離非線性法構(gòu)造用于結(jié)構(gòu)地震反應(yīng)分析的Woodbury求解公式。使用圖1所示結(jié)構(gòu)對(duì)該方法中控制方程的基本推演過程進(jìn)行說明，設(shè)在某個(gè)計(jì)算步下單元i進(jìn)入非線性，圖1給出了該單元中任一點(diǎn)材料的應(yīng)力應(yīng)變關(guān)系，其中σ和ε分別為當(dāng)前的材料應(yīng)力和應(yīng)變，根據(jù)材料初始彈性模量Ee可將應(yīng)變?chǔ)欧纸鉃閺椥院头蔷€性兩部分，其中彈性應(yīng)變?chǔ)拧涞扔诩僭O(shè)材料保持初始彈性狀態(tài)時(shí)達(dá)到應(yīng)力σ對(duì)應(yīng)的應(yīng)變(即ε′=σ/Ee)，非線性應(yīng)變?chǔ)拧宓扔诳倯?yīng)變?chǔ)排c彈性應(yīng)變?chǔ)拧涞牟?即ε″=ε-ε′)。進(jìn)一步利用插值方法在單元層面構(gòu)造材料非線性應(yīng)變分布場函數(shù)，并結(jié)合虛功原理，即可建立增量形式的單元隔離非線性控制方程

(1)

圖1 隔離非線性法基本原理

(2)

在任一迭代計(jì)算步中，將式(2)中的非線性自由度進(jìn)行凝聚可得攝動(dòng)式剛度方程式(3)。

(3)

式中:n為結(jié)構(gòu)位移自由度數(shù)；m為當(dāng)前計(jì)算步下結(jié)構(gòu)的非線性自由度數(shù)，代表了結(jié)構(gòu)非線性區(qū)域規(guī)模的大小。式(3)等號(hào)左邊括號(hào)中表達(dá)式的計(jì)算結(jié)果與結(jié)構(gòu)當(dāng)前的切向剛度Kt等效，該表達(dá)式代表了切向剛度的低秩攝動(dòng)。當(dāng)結(jié)構(gòu)局部非線性區(qū)域的規(guī)模較小時(shí)，非線性自由度數(shù)目m將遠(yuǎn)低于結(jié)構(gòu)的位移自由度數(shù)數(shù)目n(即m?n)，可利用式(4)對(duì)式(3)進(jìn)行高效求解。

(4)

(5)

當(dāng)結(jié)構(gòu)規(guī)模較大時(shí)，雖然式(4)中整體剛度Ke的規(guī)模較大，但其為常數(shù)矩陣，可僅在分析前合成并進(jìn)行一次分解，可見，Woodbury公式的使用可以避免傳統(tǒng)非線性分析方法中最為耗時(shí)的大規(guī)模整體剛度矩陣反復(fù)更新和分解運(yùn)算，取而代之的是僅對(duì)小規(guī)模Schur補(bǔ)矩陣進(jìn)行相應(yīng)運(yùn)算，實(shí)現(xiàn)了計(jì)算降維，從而大幅降低非線性分析所需計(jì)算消耗。對(duì)于地震反應(yīng)分析問題，需結(jié)合運(yùn)動(dòng)方程并引入Newmark方法構(gòu)造動(dòng)力控制方程，其對(duì)應(yīng)的Woodbury求解公式具有與式(4)完全一致的矩陣表達(dá)形式和計(jì)算特點(diǎn)。

2 Woodbury+OpenMP并行分析方法

2.1 OpenMP并行機(jī)制

OpenMP是適用于共享內(nèi)存多處理器并行計(jì)算的程序設(shè)計(jì)標(biāo)準(zhǔn)，具有靈活度高、易于編程實(shí)現(xiàn)等優(yōu)點(diǎn)。OpenMP采用Fork-Join的并行模式(Fork即創(chuàng)建或喚醒派生線程，Join即多線程的匯合)，如圖2所示，其基本思路為：程序開始時(shí)首先由一個(gè)主線程執(zhí)行，當(dāng)遇到某個(gè)并行任務(wù)時(shí)，程序?qū)⑴缮龆鄠€(gè)線程，此時(shí)主線程和派生線程將同步進(jìn)行程序計(jì)算，當(dāng)該并行任務(wù)結(jié)束時(shí)，派生線程退出工作，主線程繼續(xù)執(zhí)行后續(xù)程序計(jì)算任務(wù)，遇到下一個(gè)并行任務(wù)后派生線程將被再次喚醒并參與計(jì)算，整個(gè)計(jì)算程序中可以包含多個(gè)并行任務(wù)。

圖2 OpenMP并行機(jī)制

2.2 Woodbury方法的并行實(shí)現(xiàn)

結(jié)構(gòu)進(jìn)入非線性狀態(tài)后，需通過迭代完成每個(gè)增量計(jì)算步的求解?；贜ewton-Raphson迭代格式的Woodbury方法非線性分析流程，如圖3所示。由圖3可知，在開始進(jìn)行迭代求解之前需首先對(duì)結(jié)構(gòu)初始彈性剛度矩陣Ke進(jìn)行LDLT分解，將分解結(jié)果進(jìn)行存儲(chǔ)，以便于在迭代求解過程中隨時(shí)調(diào)用。在每次迭代時(shí)，主要包含3個(gè)計(jì)算步驟，分別是Woodbury公式計(jì)算、單元狀態(tài)確定、非線性相關(guān)系數(shù)矩陣更新，見圖3。本節(jié)通過分析上述3個(gè)步驟的計(jì)算特點(diǎn)，分別設(shè)計(jì)了基于OpenMP的并行加速方法，各部分具體并行實(shí)現(xiàn)方式詳述如下。

圖3 Woodbury方法的迭代求解流程

2.2.1 非線性系數(shù)矩陣的并行更新策略

(6)

(7)

對(duì)于系數(shù)矩陣Kinf，其為對(duì)稱矩陣，將式(6)代入式(5)，可得該矩陣的計(jì)算表達(dá)式

(8)

式中，Aji為與單元j和i有關(guān)的子矩陣塊，代表單元j對(duì)單元i的影響程度，其計(jì)算表達(dá)式為

1≤i≤p

(9)

圖4 矩陣Kinf更新過程示意圖

假設(shè)在當(dāng)前迭代計(jì)算步中僅單元p為新進(jìn)入非線性的單元，同時(shí)考慮矩陣Kinf的對(duì)稱特性，則根據(jù)式(8)可知此時(shí)僅需計(jì)算其第p列數(shù)據(jù)(包含子矩陣塊A1p，A2p，…，App)并對(duì)計(jì)算結(jié)果進(jìn)行稀疏近似即可。根據(jù)式(9)可知各子矩陣的計(jì)算互不干擾，為此，本文采用如下步驟實(shí)現(xiàn)這一計(jì)算過程的并行化。

步驟1首先計(jì)算式(10)?？紤]到分析前已完成初始彈性剛度矩陣Ke的分解運(yùn)算，可直接調(diào)用矩陣分解結(jié)果通過回代完成式(10)計(jì)算，回代過程的計(jì)算消耗極小，且可實(shí)現(xiàn)基于OpenMP并行加速，本文方法利用Intel MKL中Pardiso求解器實(shí)現(xiàn)并行回代。

(10)

式中，Sp為計(jì)算Kinf第p列數(shù)據(jù)時(shí)所需中間變量。

步驟2計(jì)算對(duì)角子矩陣App

(11)

步驟3在程序中創(chuàng)建并行區(qū)并添加OpenMP指令對(duì)，將子矩陣A1p，A2p，…，A(p-1), p的求解及其稀疏判定過程分發(fā)到不同線程上進(jìn)行并行計(jì)算。稀疏判定是通過將每個(gè)計(jì)算出的子矩陣與式(11)計(jì)算結(jié)果(即App)進(jìn)行對(duì)比實(shí)現(xiàn)的，如圖5所示。圖5中：‖·‖為矩陣范數(shù)；τ為預(yù)定義的稀疏判定閾值。由圖5可知，若某個(gè)計(jì)算出的子矩陣Aip(1≤i≤p-1)的范數(shù)與相應(yīng)對(duì)角塊App范數(shù)之比小于閾值τ，則可將其近似置零。在程序?qū)嶋H執(zhí)行時(shí)，將這些近似置零的子矩陣直接刪除，以降低數(shù)據(jù)存儲(chǔ)需求和后續(xù)Woodbury公式的計(jì)算復(fù)雜度。

圖5 系數(shù)矩陣矩陣Kinf的并行更新過程

當(dāng)某個(gè)迭代步中有多個(gè)新進(jìn)入非線性的單元時(shí)，對(duì)每個(gè)單元的相關(guān)子矩陣列均執(zhí)行上述過程，即可實(shí)現(xiàn)系數(shù)矩陣Kinf的更新。從論述可知，Woodbury公式中與結(jié)構(gòu)非線性相關(guān)的系數(shù)矩陣均可劃分為若干相互獨(dú)立的子矩陣，各子矩陣在計(jì)算時(shí)不存在數(shù)據(jù)交互，因此這些矩陣的更新過程具有較好的可并行特征。上述對(duì)于矩陣Kinf的稀疏化處理將導(dǎo)致每次迭代時(shí)Woodbury公式計(jì)算結(jié)果存在一定近似誤差，但該誤差可通過少量增加迭代次數(shù)或引入適當(dāng)?shù)恼`差修正策略予以消除，本文采用Yu等提出的基于縮減基的誤差修正方法消除這一近似誤差帶來的不利影響。

本文方法無需在分析前預(yù)先設(shè)定非線性區(qū)域，而是在每次迭代時(shí)根據(jù)單元狀態(tài)確定結(jié)果實(shí)時(shí)確定進(jìn)入非線性的單元位置，根據(jù)論述可知，非線性相關(guān)系數(shù)矩陣更新過程所需計(jì)算量與非線性單元的位置無關(guān)，僅與這些單元的數(shù)量有關(guān)，因此該部分計(jì)算效率將隨非線性區(qū)域規(guī)模的擴(kuò)大逐漸降低。

2.2.2 Woodbury公式的并行計(jì)算

式(4)僅為Woodbury公式的數(shù)學(xué)表達(dá)，為充分發(fā)揮其高效性優(yōu)勢，在實(shí)際程序?qū)崿F(xiàn)時(shí)應(yīng)避免執(zhí)行耗時(shí)的矩陣與矩陣乘法運(yùn)算。為此，將式(4)從右向左拆解為如表1所示6個(gè)計(jì)算步，每次迭代時(shí)依次執(zhí)行這6個(gè)步驟即可得到Woodbury公式的計(jì)算結(jié)果。表1中，向量ΔX1～ΔX5均為這一計(jì)算過程的中間變量。由表1可知，各計(jì)算步僅涉及簡單的矩陣與向量或向量與向量之間的運(yùn)算，每個(gè)計(jì)算步的計(jì)算量均極小?？紤]到表1各計(jì)算步中涉及的系數(shù)矩陣均具有較高稀疏性，本文采用壓縮存儲(chǔ)格式CSR3(compressed sparse row)對(duì)矩陣進(jìn)行存儲(chǔ)并基于此開展相關(guān)矩陣運(yùn)算，以避免零元素相關(guān)運(yùn)算對(duì)計(jì)算資源的占用，降低計(jì)算復(fù)雜度。

表1 Woodbury公式計(jì)算步驟

對(duì)于表1中的第1個(gè)和第5個(gè)計(jì)算步，可直接調(diào)用矩陣Ke的分解結(jié)果通過回代計(jì)算完成，本文基于Pardiso求解器使用并行回代方法進(jìn)行這兩個(gè)計(jì)算步的快速求解；對(duì)于表1中的第2個(gè)和第4個(gè)計(jì)算步，其均為稀疏矩陣與向量之間的乘法運(yùn)算，根據(jù)矩陣運(yùn)算基本原理可知，等號(hào)右邊向量ΔX2和ΔX4中各元素的求解過程均互不干擾，因而可基于OpenMP并行方法實(shí)現(xiàn)這兩個(gè)計(jì)算步的加速求解；表1中的第3個(gè)計(jì)算步要求對(duì)以Schur補(bǔ)矩陣為系數(shù)矩陣的方程組進(jìn)行求解，可采用基于OpenMP的并行分解和回代方法進(jìn)行高效計(jì)算，本文方法利用Pardiso求解器實(shí)現(xiàn)該計(jì)算步的并行求解，由于Schur補(bǔ)矩陣的規(guī)模通常較小，且可近似化為具有高度稀疏性的矩陣，該步的計(jì)算耗時(shí)通常極??；表1中的第6個(gè)計(jì)算步為兩個(gè)向量的加法運(yùn)算，其計(jì)算量極小，且亦可直接實(shí)現(xiàn)OpenMP的并行加速。綜上可知，通過將Woodbury公式拆解為6個(gè)連續(xù)執(zhí)行的計(jì)算步，并基于OpenMP方法對(duì)每個(gè)計(jì)算步的求解進(jìn)行加速，可實(shí)現(xiàn)該公式計(jì)算的全過程并行化?？紤]到基于矩陣稀疏存儲(chǔ)進(jìn)行相關(guān)運(yùn)算可大幅降低計(jì)算過程的時(shí)間和空間復(fù)雜度，本文方法可極大提高每次迭代時(shí)Woodbury公式的實(shí)際程序運(yùn)行效率。

2.2.3 單元狀態(tài)確定過程的并行實(shí)現(xiàn)

單元狀態(tài)確定的主要作用是在每次迭代時(shí)根據(jù)計(jì)算出的結(jié)點(diǎn)位移數(shù)據(jù)對(duì)各單元材料應(yīng)力、應(yīng)變、切向模量等參量進(jìn)行更新，并計(jì)算各單元的結(jié)點(diǎn)恢復(fù)力，以便于更新不平衡力。對(duì)于Woodbury方法，還需在單元狀態(tài)確定過程中判定各單元是否進(jìn)入非線性，以便于在后續(xù)迭代計(jì)算中對(duì)Woodbury求解公式中非線性相關(guān)矩陣進(jìn)行更新。由于單元狀態(tài)確定過程中各單元之間不存在數(shù)據(jù)交互，具有典型的并行化特點(diǎn)，因而可直接在該部分程序代碼的開始和結(jié)束部位添加OpenMP并行指令對(duì)，以實(shí)現(xiàn)計(jì)算過程的并行加速。

3 算例分析

參考文獻(xiàn)[49]建立如圖6所示的巨型鋼框架結(jié)構(gòu)分析模型，結(jié)構(gòu)共計(jì)48層，層高4 m，總高度為192 m，平面尺寸為40 m×40 m，每邊均為5跨，跨度均為8 m，圖6(b)和圖6(c)分別給出了該結(jié)構(gòu)的平立面布置。結(jié)構(gòu)1～24層柱采用Q345鋼，其他層均采用Q235鋼。使用基于隔離非線性法的纖維梁單元建立該結(jié)構(gòu)分析模型，采用考慮隨動(dòng)硬化的雙線性彈塑性本構(gòu)模型模擬材料非線性行為，其中屈服后剛度系數(shù)設(shè)為0.01。為充分驗(yàn)證本文建立的并行Woodbury方法對(duì)于大規(guī)模結(jié)構(gòu)非線性問題的高效性，對(duì)模型網(wǎng)格進(jìn)行適當(dāng)加密，本例中，將各構(gòu)件每隔1 m劃分一個(gè)單元，共計(jì)劃分40 608個(gè)單元，結(jié)構(gòu)的總位移自由度數(shù)超過21萬(共計(jì)215 712)。選取El-Centro地震記錄對(duì)結(jié)構(gòu)進(jìn)行非線性地震反應(yīng)分析，考慮施加平面雙向地震激勵(lì)，其中X方向地面運(yùn)動(dòng)的峰值加速度調(diào)幅至6.0 m/s2，Y方向地面運(yùn)動(dòng)的峰值加速度調(diào)幅到5.1 m/s2。計(jì)算步長設(shè)為0.01 s，截取前20 s地震記錄進(jìn)行分析，共計(jì)分析2 000個(gè)增量計(jì)算步。為驗(yàn)證本文方法的準(zhǔn)確性，同時(shí)使用有限元軟件ABAQUS建立該結(jié)構(gòu)的數(shù)值分析模型并選用相同地震記錄進(jìn)行非線性分析，ABAQUS軟件模型的網(wǎng)格密度與本文方法相同，單元類型為B32。

(a) 三維模型

本文Woodbury+OpenMP方法、串行Woodbury法和ABAQUS軟件分析得到的結(jié)構(gòu)頂點(diǎn)沿X和Y方向的位移時(shí)程曲線，如圖7所示。由圖7可知，并行和串行Woodbury法的分析結(jié)果相同，且兩者與ABAQUS軟件分析結(jié)果基本一致。本文方法和ABAQUS軟件算得的結(jié)構(gòu)X方向最大頂點(diǎn)位移分別為0.429 9 m和0.425 3 m，僅相差1.08%，Y方向的結(jié)構(gòu)最大頂點(diǎn)位移分別為0.366 7 m 和0.368 5 m，僅相差0.48%，這表明本文方法計(jì)算結(jié)果具有與ABAQUS軟件一致的精度水平。本文方法分析時(shí)非線性自由度數(shù)目的時(shí)程曲線，如圖8所示。由于非線性自由度數(shù)等于Woodbury求解公式中需實(shí)時(shí)更新和分解的Schur補(bǔ)矩陣階數(shù)，該值越小表明Woodbury公式的求解效率越高，由圖8可知，結(jié)構(gòu)非線性自由度數(shù)的峰值為8 379，僅為結(jié)構(gòu)位移自由度數(shù)目的3.88%，遠(yuǎn)小于傳統(tǒng)非線性分析方法中需實(shí)時(shí)更新和分解的整體剛度矩陣階數(shù)(即位移自由度數(shù))，這說明對(duì)于局部非線性問題，使用Woodbury公式能夠在非線性分析時(shí)實(shí)現(xiàn)大幅度的計(jì)算降維，這也是其計(jì)算高效性的主要來源。

(a) X方向

圖8 非線性自由度時(shí)程曲線

本文方法的所需計(jì)算量與結(jié)構(gòu)非線性區(qū)域規(guī)模有關(guān)，對(duì)于和本例類似的局部非線性問題，其高效性優(yōu)勢尤為顯著，雖然隨結(jié)構(gòu)非線性區(qū)域規(guī)模的提高其計(jì)算效率將逐漸降低，但通過結(jié)合適當(dāng)?shù)挠?jì)算優(yōu)化策略，如本文采用的Schur補(bǔ)矩陣稀疏優(yōu)化方法，本文方法的高效性在非線性區(qū)域規(guī)模較大時(shí)依然能夠得到保證。若結(jié)構(gòu)出現(xiàn)全局非線性或大部分單元進(jìn)入非線性的情況，本文方法的適用性將降低，然而，對(duì)于工程結(jié)構(gòu)的抗震分析問題，這一極端情況通常較少出現(xiàn)。

為驗(yàn)證本文提出的并行加速策略的有效性，本文Woodbury+OpenMP方法、串行Woodbury方法和ABAQUS軟件的總計(jì)算耗時(shí)，如表2所示。使用的計(jì)算平臺(tái)為裝配Intel Core i9-10920X處理器的個(gè)人PC機(jī)，處理器核數(shù)為12?？梢钥闯?，本文方法僅耗時(shí)967 s，為串行Woodbury分析方法計(jì)算耗時(shí)的18%，這表明對(duì)于該算例，本文并行方法的加速比為5.5。此外，ABAQUS軟件的計(jì)算耗時(shí)約為9.5 h，分別為串行和并行Woodbury方法計(jì)算耗時(shí)的6.3倍和34.6倍，并行與串行Woodbury方法各部分的計(jì)算耗時(shí)，如表3所示?？梢钥闯?，本文方法在單元狀態(tài)確定部分的并行效率最高，這主要是由于該部分計(jì)算時(shí)各單元完全獨(dú)立，不存在串行計(jì)算代碼，無需執(zhí)行額外的進(jìn)程創(chuàng)建與銷毀工作。由表3可知，系數(shù)矩陣更新和Woodbury公式計(jì)算的并行效率雖然相對(duì)較低，但此兩部分的計(jì)算耗時(shí)依然能夠?qū)崿F(xiàn)較大幅度的降低(分別降低了60.4%和49.8%)。由于Woodbury求解公式本身僅需極小的計(jì)算消耗，串行分析程序的主要計(jì)算耗時(shí)集中于單元狀態(tài)確定部分和矩陣更新，因而從這兩部分節(jié)約出的較高計(jì)算時(shí)間有助于大幅提升本文方法的整體并行計(jì)算效率。

表2 計(jì)算耗時(shí)對(duì)比

表3 Woodbury方法不同計(jì)算部分耗時(shí)統(tǒng)計(jì)

4 結(jié) 論

本文利用OpenMP并行編程模式，針對(duì)基于Woodbury公式的結(jié)構(gòu)非線性動(dòng)力反應(yīng)分析方法，通過對(duì)其迭代求解過程中的主要3個(gè)主要計(jì)算部分(非線性相關(guān)系數(shù)矩陣更新、Woodbury公式計(jì)算和單元狀態(tài)確定)進(jìn)行并行加速，提出了一種用于高效非線性地震反應(yīng)分析的Woodbury+OpenMP方法。本文具體結(jié)論如下：

(2) Woodbury公式在實(shí)際應(yīng)用時(shí)可以拆解為依次執(zhí)行的6個(gè)計(jì)算步，各計(jì)算步僅涉及矩陣與向量或向量與向量之間的運(yùn)算，由于各計(jì)算步均可實(shí)現(xiàn)基于OpenMP并行化，本文方法能夠大幅提升該公式的實(shí)際程序執(zhí)行效率。

(3) 本文方法中，Woodbury公式的使用避免了傳統(tǒng)非線性地震反應(yīng)分析方法所需的整體剛度方程反復(fù)更新求解，進(jìn)一步結(jié)合并行計(jì)算技術(shù)，本文方法在保證計(jì)算精度與傳統(tǒng)非線性分析方法相當(dāng)?shù)那疤嵯驴纱蠓嵘Y(jié)構(gòu)地震反應(yīng)分析效率。