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

(1.中國建筑科學研究院，北京 100013;2.廣州大學廣東省地震工程與應用技術重點實驗室，廣州 510415)

引言

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

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

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

1 動力方程求解方法

動力彈塑性分析的本質是求解如公式1 所示的動力學方程式:

式中，［M］-質量矩陣;

［C］-阻尼矩陣;

［K］-剛度矩陣;

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

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

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

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

2 精細化彈塑性模型

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

2.1 線單元模擬

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

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

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

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

2.2 剪力墻和樓板單元模擬

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

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

如圖2 所示，某框架核心筒結構，樓板構件數為4 611，剪力墻構件數為4 784，樓板與剪力墻構件數量基本相當，對樓板和剪力墻采用相同的網格劃分尺寸以后，考慮樓板彈塑性與否對于計算量有將近一倍左右的差距。

3 并行計算

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

3.1 并行計算簡介

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

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

圖2 某框架核心筒結構

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

3.2 SAUSAGE 軟件并行計算實現及效率

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

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

圖3 樓板網格劃分

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

4 結語

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

［1］上?，F代建筑設計(集團)有限公司技術中心，動力彈塑性時程分析在建筑結構抗震設計中的應用［M］.上海:上?？茖W技術出版社，2013.

［2］廣州建研數力建筑科技有限公司，PKPM-SAUSAGE 高性能彈塑性動力時程分析軟件使用手冊［M］.2013.