并行計(jì)算技術(shù)在結(jié)構(gòu)顯式動(dòng)力分析中的應(yīng)用

侯曉武 李志山 喬保娟 劉春明 楊志勇

(1.中國建筑科學(xué)研究院，北京 100013;2.廣州大學(xué)廣東省地震工程與應(yīng)用技術(shù)重點(diǎn)實(shí)驗(yàn)室，廣州 510415)

引言

近年來，隨著越來越多復(fù)雜超限高層建筑結(jié)構(gòu)的出現(xiàn)，一般的彈性分析設(shè)計(jì)方法已經(jīng)很難滿足要求，我國是地震多發(fā)國家，結(jié)構(gòu)在大震作用下的彈塑性分析顯得越來越重要。它是對(duì)傳統(tǒng)分析設(shè)計(jì)方法的一個(gè)很好補(bǔ)充，通過在計(jì)算模型中引入塑性鉸模型、纖維模型和非線性殼單元，模擬結(jié)構(gòu)受荷后變形、開裂直至破壞的整個(gè)過程，能夠準(zhǔn)確預(yù)測(cè)結(jié)構(gòu)在罕遇地震作用下的彈塑性響應(yīng)，以此作為大震作用下結(jié)構(gòu)抗震性能評(píng)價(jià)的依據(jù)，并為超限復(fù)雜建筑結(jié)構(gòu)抗震設(shè)計(jì)提供指導(dǎo)，發(fā)現(xiàn)結(jié)構(gòu)的薄弱環(huán)節(jié)，有針對(duì)性的改善抗震設(shè)計(jì)。

彈塑性分析方法分為靜力彈塑性分析(Pushover)方法和動(dòng)力彈塑性分析方法，靜力彈塑性分析將地震力等效為作用在結(jié)構(gòu)上按一定規(guī)律分布的水平靜力荷載，不斷增大該荷載直至結(jié)構(gòu)倒塌。這種方法原理簡單，實(shí)現(xiàn)起來相對(duì)容易。但也存在如下幾個(gè)主要問題:1)水平荷載分布與實(shí)際地震荷載分布有出入;2)不能真實(shí)反映構(gòu)件在地震作用下卸載時(shí)的剛度退化以及內(nèi)力重分布等特性;3)主要適用于以第一振型為主，結(jié)構(gòu)周期較短的結(jié)構(gòu)。動(dòng)力彈塑性分析方法將地震波直接作用于建筑結(jié)構(gòu)，通過數(shù)值積分算法求解動(dòng)力學(xué)方程式，以得到結(jié)構(gòu)在地震作用下的響應(yīng)。該方法相對(duì)簡化較少，得到的結(jié)果更加準(zhǔn)確可靠，因而越來越多地應(yīng)用于復(fù)雜結(jié)構(gòu)的性能化設(shè)計(jì)。

對(duì)于復(fù)雜高層結(jié)構(gòu)，采用纖維束以及非線性分層殼等精細(xì)化模型，其單元自由度數(shù)一般可以達(dá)到幾百萬個(gè)甚至千萬個(gè)。對(duì)其進(jìn)行彈塑性時(shí)程分析意味著超大規(guī)模的計(jì)算量以及數(shù)據(jù)存儲(chǔ)，這對(duì)于軟件的計(jì)算能力提出了相當(dāng)高的要求。以往的并行計(jì)算多采用CPU 多核并行或計(jì)算機(jī)集群并行，對(duì)于軟硬件要求比較高，PKPM SAUSAGE 軟件采用目前技術(shù)領(lǐng)先的CPU+GPU 異構(gòu)并行計(jì)算技術(shù)，通過一系列高效獨(dú)到的數(shù)據(jù)訪存策略，使得動(dòng)力彈塑性分析高效快速，以往采用大型通用有限元程序并且需要相當(dāng)熟練的使用技巧才能完成的顆粒仿真分析工作，可以用一臺(tái)配備很低成本的計(jì)算顯卡為主要硬件實(shí)現(xiàn)，同時(shí)方便的得到結(jié)果，為結(jié)構(gòu)工程師進(jìn)行結(jié)構(gòu)概念設(shè)計(jì)提供了方便實(shí)用的工具。

1 動(dòng)力方程求解方法

動(dòng)力彈塑性分析的本質(zhì)是求解如公式1 所示的動(dòng)力學(xué)方程式:

式中，［M］-質(zhì)量矩陣;

［C］-阻尼矩陣;

［K］-剛度矩陣;

求解動(dòng)力學(xué)方程式可以采用振型疊加法和直接積分法。振型疊加法是把多自由度體系結(jié)構(gòu)的整體振動(dòng)分解為與振型個(gè)數(shù)相對(duì)應(yīng)的單自由度體系，求得各個(gè)單自由度體系的動(dòng)力響應(yīng)后，再進(jìn)行疊加得到結(jié)構(gòu)整體響應(yīng)。直接積分法是用有限差分格式通過位移來近似表示速度和加速度，用數(shù)值積分的方法在時(shí)域上逐步對(duì)方程進(jìn)行積分，以得到結(jié)構(gòu)響應(yīng)。對(duì)于動(dòng)力彈塑性分析，由于材料本身的非線性特性和幾何非線性，疊加原理不再適用，因而不能采用振型疊加法，僅能采用直接積分法。

根據(jù)差分方法的不同，直接積分法可分為顯式與隱式兩類。顯式差分法可以根據(jù)前一時(shí)刻的平衡條件直接求解下一時(shí)刻的各參數(shù)解，而隱式差分法則必須對(duì)方程進(jìn)行迭代求解。顯式差分法有中心差分法等，隱式差分法有Newmark 法及HHT(Hilber-Hughes-Taylor)法等。采用隱式算法求解非線性動(dòng)力學(xué)問題時(shí)，對(duì)于每一個(gè)時(shí)間步都需要進(jìn)行大量的迭代計(jì)算，對(duì)于大規(guī)模結(jié)構(gòu)非線性分析不容易收斂，同時(shí)耗時(shí)較長。因而在求解此類問題時(shí)，顯式算法具有較大的優(yōu)越性。

目前市場(chǎng)上可以用于動(dòng)力彈塑性分析的軟件中，采用隱式算法的軟件占據(jù)絕大多數(shù)，如Perform 3D，SAP2000，midas Building，midas Gen，CANNY，OpenSees，EPDA 等。采用顯式算法的包括ABAQUS，LSDYNA 以及PKPM-SAUSAGE。

采用顯式算法求解動(dòng)力學(xué)問題時(shí)，對(duì)于時(shí)間步長有要求，屬于條件穩(wěn)定1。即只有當(dāng)時(shí)間步長小于臨界穩(wěn)定時(shí)間步長時(shí)，顯式算法才穩(wěn)定。一般情況下，高層結(jié)構(gòu)的臨界穩(wěn)定時(shí)間步長可達(dá)到10-5s～10-6s，這與隱式算法通常采用的時(shí)間步長(一般為0.01s～0.02s)相差1 000 甚至10 000 倍。時(shí)間步長的減小意味著計(jì)算量的增加以及模型數(shù)據(jù)量的增大，因而必須尋求與其相適應(yīng)的計(jì)算技術(shù)。

2 精細(xì)化彈塑性模型

SAUSAGE 軟件中對(duì)于梁、柱、樓板和剪力墻均采用了精細(xì)化模型，梁柱采用纖維模型，剪力墻和樓板采用分層殼模型，通過模型精細(xì)模擬和盡可能少的人為簡化，以保證彈塑性分析結(jié)果的準(zhǔn)確和客觀。

2.1 線單元模擬

構(gòu)建結(jié)構(gòu)動(dòng)力彈塑性分析模型時(shí)，梁、柱、支撐等線單元一般通過塑性鉸或纖維模型進(jìn)行模擬。塑性鉸模型是指對(duì)試驗(yàn)確定的力與變形關(guān)系進(jìn)行理想化，得到線單元的滯回曲線，如采用雙折線模型模擬鋼構(gòu)件，采用三折線模型模擬混凝土構(gòu)件以考慮混凝土的開裂，或采用四折線模型模擬混凝土構(gòu)件達(dá)到極限承載力后承載能力的降低。塑性鉸模型概念比較清晰，結(jié)果比較直觀，因而目前應(yīng)用比較廣泛。

纖維模型是將構(gòu)件沿縱向分割成若干段，以每一段中間某一截面的特性代表該段的特性。然后將該特征截面分割為若干個(gè)纖維(按材料又可區(qū)分為混凝土纖維、鋼筋纖維和鋼材纖維等)，纖維與纖維之間滿足平截面假定，通過定義鋼筋、混凝土和鋼材的非線性本構(gòu)關(guān)系來模擬構(gòu)件的非線性特性。纖維模型相對(duì)簡化較少，各種材料的本構(gòu)關(guān)系一般可以通過規(guī)范或者試驗(yàn)得到，因而相對(duì)而言適用性更廣泛。

SAUSAGE 軟件采用纖維模型模擬，對(duì)于常用的矩形柱和矩形梁單元，程序默認(rèn)劃分為6 ×6 和2 ×6 個(gè)纖維單元2，如圖1 所示。

圖1 矩形柱和矩形梁截面纖維劃分

2.2 剪力墻和樓板單元模擬

剪力墻單元主要通過纖維束或分層殼模型進(jìn)行模擬。采用纖維束模型模擬剪力墻時(shí)，與梁單元類似，沿剪力墻平面方向劃分橫向和豎向纖維，并通過積分得到水平方向和豎向的剛度。分層殼模型是沿剪力墻厚度方向進(jìn)行分層，包括混凝土層和鋼筋層，對(duì)于含鋼板的混凝土剪力墻，還包括鋼板層，各層在厚度方向滿足平截面假定。SAUSAGE 軟件中，采用分層殼模型模擬剪力墻，混凝土剪力墻中包含6 個(gè)混凝土層，鋼筋網(wǎng)層數(shù)可由用戶自行設(shè)置，程序默認(rèn)設(shè)置為兩層。

地震荷載作用下，樓板具有傳遞水平荷載、協(xié)調(diào)豎向構(gòu)件變形的作用。大震作用時(shí)應(yīng)該保證樓板的承載能力，以保證結(jié)構(gòu)的完整性。目前大多數(shù)軟件中，一般忽略樓板的彈塑性變化，通過剛性樓板假定或彈性板模型模擬樓板進(jìn)行彈塑性分析，這不利于掌握地震作用下樓板的實(shí)際破壞狀態(tài)，為結(jié)構(gòu)設(shè)計(jì)留下了隱患。SAUSAGE 軟件中對(duì)于樓板采用和剪力墻相同的分層殼模型進(jìn)行模擬，混凝土沿板厚方向不再進(jìn)行分層，鋼筋網(wǎng)層數(shù)默認(rèn)設(shè)置為2。

如圖2 所示，某框架核心筒結(jié)構(gòu)，樓板構(gòu)件數(shù)為4 611，剪力墻構(gòu)件數(shù)為4 784，樓板與剪力墻構(gòu)件數(shù)量基本相當(dāng)，對(duì)樓板和剪力墻采用相同的網(wǎng)格劃分尺寸以后，考慮樓板彈塑性與否對(duì)于計(jì)算量有將近一倍左右的差距。

3 并行計(jì)算

如第1 節(jié)和第2 節(jié)所述，由于SAUSAGE 軟件采取了精細(xì)化模型，考慮了樓板的彈塑性特性，增加了單元數(shù)量，從而增加模型的自由度數(shù)量。對(duì)于大規(guī)模高層建筑結(jié)構(gòu)，自由度數(shù)可達(dá)到幾百萬個(gè)。同時(shí)，SAUSAGE 軟件中采用顯式積分求解動(dòng)力學(xué)方程式，需要保證計(jì)算時(shí)間步長小于臨界穩(wěn)定時(shí)間步長，所以實(shí)際的步長一般控制在10-5s 級(jí)。假設(shè)單條地震波計(jì)算時(shí)長為30s，則動(dòng)力分析需要計(jì)算3百萬步。龐大的自由度數(shù)量以及百萬級(jí)的計(jì)算步驟數(shù)決定了軟件必須采取相應(yīng)的計(jì)算方法，以提高彈塑性分析的計(jì)算效率。

3.1 并行計(jì)算簡介

并行計(jì)算(Parallel Computing)是指同時(shí)使用多種計(jì)算資源解決計(jì)算問題的過程，是提高計(jì)算機(jī)系統(tǒng)計(jì)算速度和處理能力的一種有效手段。它的基本思想是用多個(gè)處理器來協(xié)同求解同一問題，即將被求解的問題分解成若干個(gè)部分，各部分均由一個(gè)獨(dú)立的處理器來并行計(jì)算。

傳統(tǒng)的并行計(jì)算以多核CPU 并行或多機(jī)并行為特征，其實(shí)質(zhì)是若干個(gè)串行任務(wù)的并行執(zhí)行。通用有限元軟件ABAQUS 和ANSYS 提供多核CPU 并行計(jì)算功能，這種并行方式對(duì)軟硬件配置有較高要求，需采用大中型服務(wù)器或多機(jī)集群，實(shí)施起來相對(duì)復(fù)雜，成本也比較高。GPU(Graphics Processing Unit，顯卡的處理器)具有多個(gè)高內(nèi)存帶寬驅(qū)動(dòng)的計(jì)算內(nèi)核，具有浮點(diǎn)計(jì)算能力強(qiáng)、帶寬高、性價(jià)比高、能耗低等優(yōu)點(diǎn)，以前主要用于圖形處理，目前已經(jīng)廣泛應(yīng)用于通用計(jì)算領(lǐng)域。以Nvidia 公司的GTX TITAN X 為例，其處理器數(shù)為3 072 個(gè)，在進(jìn)行大規(guī)模并行計(jì)算時(shí)，其效率可以達(dá)到CPU 的數(shù)百倍。

圖2 某框架核心筒結(jié)構(gòu)

由于CPU 在分支處理及隨機(jī)內(nèi)存讀取方面有優(yōu)勢(shì)，適合于處理串聯(lián)工作。而GPU 在處理大量有浮點(diǎn)運(yùn)算的并行計(jì)算時(shí)有優(yōu)勢(shì)。為了充分發(fā)揮CPU 和GPU 各自在計(jì)算方面的優(yōu)勢(shì)，使用CPU 進(jìn)行程序中的串行部分，而對(duì)于程序中的并行部分，通過GPU 進(jìn)行加速，這就是“CPU+GPU 異構(gòu)計(jì)算”的核心思想。

3.2 SAUSAGE 軟件并行計(jì)算實(shí)現(xiàn)及效率

為了充分利用GPU 的并行計(jì)算性能，SAUSAGE 軟件構(gòu)造了具有針對(duì)性的數(shù)據(jù)結(jié)構(gòu)和分析算法，對(duì)于隱式計(jì)算和顯式計(jì)算，均實(shí)現(xiàn)了并行求解。對(duì)于模態(tài)分析、最大頻率分析以及靜力分析等隱式分析，并行計(jì)算單元?jiǎng)偠染仃嚥⒔M集總體剛度矩陣，而后通過自行開發(fā)的基于GPU 并行計(jì)算的求解器實(shí)現(xiàn)代數(shù)方程組的求解。SAUSAGE 動(dòng)力計(jì)算采用的顯式積分求解方法，對(duì)于單元坐標(biāo)轉(zhuǎn)換矩陣、單元應(yīng)變以及應(yīng)力等，均可以實(shí)現(xiàn)單元級(jí)別的并行計(jì)算。對(duì)于每個(gè)時(shí)間步，求解初始加速度時(shí)，僅需對(duì)質(zhì)量矩陣進(jìn)行求逆，質(zhì)量采用集中質(zhì)量矩陣后，其為對(duì)角單元矩陣，求逆更加高效快速。阻尼考慮方法選擇瑞利阻尼時(shí)，［C］=α［M］+β［K］，α、β 分別為質(zhì)量阻尼系數(shù)和剛度阻尼系數(shù)。如果考慮剛度阻尼，一般會(huì)使顯式積分時(shí)間步長減小1～2 個(gè)數(shù)量級(jí)，整體計(jì)算時(shí)間會(huì)增加幾十倍，因而一般不予考慮。忽略剛度阻尼以后，阻尼矩陣僅與質(zhì)量矩陣有關(guān)，由于質(zhì)量矩陣為對(duì)角陣，因而阻尼矩陣也為對(duì)角陣，方程組可以很方便的進(jìn)行解耦求解。阻尼考慮方法選擇振型阻尼時(shí)，由于振型對(duì)于質(zhì)量矩陣和剛度矩陣均具有正交性，可利用這種特點(diǎn)對(duì)阻尼矩陣進(jìn)行變換，也可以實(shí)現(xiàn)方程組的解耦。對(duì)于解耦后的方程便于實(shí)現(xiàn)并行計(jì)算，同時(shí)避免了剛度矩陣組裝以及矩陣求逆等工作，提升了工作效率。

如第1 節(jié)所述，顯式分析的時(shí)間步長應(yīng)小于臨界時(shí)間步長，而臨界時(shí)間步長與結(jié)構(gòu)的最高階圓頻率有關(guān)系，結(jié)構(gòu)的最高階圓頻率越大，則臨界時(shí)間步長越小。結(jié)構(gòu)的最高階圓頻率與單元尺寸有關(guān)，因而模型的網(wǎng)格劃分會(huì)影響臨界時(shí)間步長。一方面單元尺寸不能過大以保證足夠的精度，另一方面單元尺寸也不能過小導(dǎo)致臨界時(shí)間步長過小，這就對(duì)網(wǎng)格劃分提出了比較高的要求。SAUSAGE 開發(fā)了一套高效的網(wǎng)格自動(dòng)生成技術(shù)，使用鋪砌法生成以四邊形為主、三角形為過渡的較高質(zhì)量混合網(wǎng)格，以確保顯式分析時(shí)間步長在合理范圍內(nèi)。某模型樓板網(wǎng)格劃分如圖3 所示。

圖3 樓板網(wǎng)格劃分

通過數(shù)十個(gè)實(shí)際工程彈塑性分析對(duì)比，在同等計(jì)算規(guī)模下，采用CPU +GPU 異構(gòu)并行計(jì)算技術(shù)，SAUSAGE 的單機(jī)動(dòng)力彈塑性分析速度可以比ABAQUS 單機(jī)CPU 并行計(jì)算快3～6 倍。

4 結(jié)語

本文首先介紹了SAUSAGE 軟件求解動(dòng)力學(xué)方程式時(shí)所采用的顯式積分算法以及使用該算法時(shí)時(shí)間步長的控制。然后介紹了SAUSAGE 軟件所采用的梁、柱、樓板和剪力墻彈塑性模型，對(duì)梁柱采用纖維模型，對(duì)剪力墻和樓板采用分層殼模型進(jìn)行模擬以及精細(xì)化模型帶來自由度數(shù)量的增加。最后介紹了SAUSAGE 軟件采用CPU+GPU 異構(gòu)并行計(jì)算技術(shù)，顯著提升大規(guī)模建筑結(jié)構(gòu)動(dòng)力彈塑性分析的計(jì)算效率，適用于大規(guī)模高層建筑結(jié)構(gòu)的抗震性能分析。

［1］上?，F(xiàn)代建筑設(shè)計(jì)(集團(tuán))有限公司技術(shù)中心，動(dòng)力彈塑性時(shí)程分析在建筑結(jié)構(gòu)抗震設(shè)計(jì)中的應(yīng)用［M］.上海:上?？茖W(xué)技術(shù)出版社，2013.

［2］廣州建研數(shù)力建筑科技有限公司，PKPM-SAUSAGE 高性能彈塑性動(dòng)力時(shí)程分析軟件使用手冊(cè)［M］.2013.