李紅豫， 滕　軍， 李祚華

(哈爾濱工業(yè)大學(xué) 深圳研究生院，深圳　518055)

?

鋼筋混凝土框架結(jié)構(gòu)非線性靜、動力分析的高效計算平臺HSNAS(GPU)
——Ⅰ程序開發(fā)

李紅豫， 滕軍， 李祚華

(哈爾濱工業(yè)大學(xué) 深圳研究生院，深圳518055)

摘要：基于傳統(tǒng)串行計算平臺的有限元分析面臨精度不足、耗時巨大的問題成為目前高層鋼筋混凝土框架結(jié)構(gòu)非線性動力時程分析面臨的瓶頸，利用GPU強(qiáng)大的并行計算能力，開發(fā)了一種結(jié)構(gòu)非線性有限元靜力、動力分析的高精度和高效率分析平臺HSNAS(GPU)。針對靜力問題提出了適用于GPU計算，且能有效解決結(jié)構(gòu)負(fù)剛度問題的位移增量迭代算法，開發(fā)了相應(yīng)的GPU線性方程組并行求解器；針對動力問題，開發(fā)了GPU基于Newmark時間積分算法的動力分析軟件平臺，結(jié)合纖維模型單元技術(shù)，引入扭轉(zhuǎn)、剪切變形以及材料非線性。算例模型表明，HSNAS(GPU)平臺在滿足精度條件下能有效地提高結(jié)構(gòu)非線性靜、動力分析的計算效率。

關(guān)鍵詞：鋼筋混凝土結(jié)構(gòu)；非線性有限元；模擬平臺；GPU；纖維模型

高層建筑結(jié)構(gòu)的復(fù)雜性以及精細(xì)化分析模型龐大的自由度數(shù)量，造成數(shù)值計算的規(guī)模龐大，對軟件分析平臺的計算速度要求更高。然而在傳統(tǒng)CPU(Central Processor Unit，中央處理器)平臺上完成計算任務(wù)，耗時多成為目前結(jié)構(gòu)有限元分析面臨的瓶頸[1]。如何實現(xiàn)性能高價格低的有限元分析平臺，成為迫切需要解決的問題。隨著計算機(jī)軟硬件的進(jìn)步，高性能計算正在向普及化發(fā)展。近年來計算機(jī)圖形處理器GPU(Graphics Processor Unit)已超過摩爾定律的速度而高速發(fā)展，很大程度地提高了與計算機(jī)圖形相關(guān)應(yīng)用領(lǐng)域的快速發(fā)展。由于其硬件構(gòu)造特殊，GPU的浮點運(yùn)算、并行計算能力提供數(shù)十倍乃至于上百倍于CPU的性能，可以很好地解決大規(guī)模的科學(xué)計算問題[2]。因此，利用GPU強(qiáng)大的并行計算能力，開發(fā)基于GPU的鋼筋混凝土框架結(jié)構(gòu)非線性靜力、動力計算的高精度高效率分析平臺，突破目前復(fù)雜大型結(jié)構(gòu)高度非線性數(shù)值模擬計算時間成本大的瓶頸，具有重要實用價值。

結(jié)構(gòu)非線性分析模型和分析方法是制約結(jié)構(gòu)靜力、動力分析應(yīng)用的重要因素。為了更精確地反映構(gòu)件的非線性行為，提出了一種精細(xì)的分析模型——纖維模型[3]，該模型直接從材料本構(gòu)關(guān)系出發(fā)獲得截面乃至整個構(gòu)件單元的非線性性能。然而工程中用到的纖維模型是基于平截面假定的，這顯然不適用于需要考慮構(gòu)件剪切或扭轉(zhuǎn)影響的情形，如何彌補(bǔ)傳統(tǒng)纖維模型單元剛度忽略剪切和扭轉(zhuǎn)的不足，成為研究精細(xì)化分析模型需要解決的問題之一[4]。在求解非線性問題的計算方法上，仍然存在求解精度和效率無法兼顧、負(fù)剛度等問題[5]。

為此，本文基于GPU開發(fā)了一種鋼筋混凝土框架結(jié)構(gòu)非線性有限元分析高效計算平臺HSNAS(GPU) (High-Rise Structure Nonlinear Analysis Software based on GPU, 軟件著作權(quán)編號：2014SR086353)。針對靜力問題提出了有效地解決結(jié)構(gòu)負(fù)剛度問題且適用于GPU計算的位移增量迭代算法，在此基礎(chǔ)上開發(fā)了GPU的線性方程組并行求解器；針對動力問題，開發(fā)了基于GPU的Newmark時間積分算法的結(jié)構(gòu)動力彈塑性時程分析軟件，該分析軟件結(jié)合纖維模型單元技術(shù)，可以考慮扭轉(zhuǎn)、剪切變形和材料非線性的影響，采用了多個適用于GPU計算的并行計算策略，實現(xiàn)了結(jié)構(gòu)非線性分析的加速，為工程應(yīng)用中梁柱結(jié)構(gòu)非線性有限元分析提供了一種高精度、高效率的計算分析平臺。

1GPU并行計算平臺設(shè)計

1.1GPU平臺的任務(wù)劃分

在實際應(yīng)用程序中，大量的執(zhí)行時間都耗費(fèi)在循環(huán)迭代上?；贙nuth的統(tǒng)計，一個應(yīng)用程序中規(guī)模不到4%的循環(huán)迭代耗費(fèi)了約50%的應(yīng)用程序執(zhí)行時間。因此發(fā)掘應(yīng)用程序的并行性是提高應(yīng)用程序或應(yīng)用算法執(zhí)行性能的關(guān)鍵。

對一棟6層2×2跨的鋼筋混凝土框架進(jìn)行靜力和動力作用下的模擬，統(tǒng)計各計算步驟的耗時，結(jié)果表明求解有限元方程組占整個問題求解時間的70%以上，動力問題相對靜力問題占的比重更大。因此本文在HSNAS(GPU)平臺上對算法流程在CPU和GPU進(jìn)行任務(wù)劃分(如圖1所示)。將靜力和動力計算程序交付給GPU執(zhí)行，而初始化計算及前后處理則由CPU端完成。

圖1　HSNAS(GPU)平臺程序任務(wù)劃分Fig.1 HSNAS(GPU) task partitioning

1.2GPU平臺的核心架構(gòu)

基于GPU架構(gòu)的計算本質(zhì)上是一種并行尋優(yōu)的計算過程。在鋼筋混凝土結(jié)構(gòu)非線性問題的求解過程中，GPU計算模式需要遵循以下幾條原則對整體有限元程序架構(gòu)進(jìn)行重新設(shè)計：① 充分考慮開發(fā)程序的并行性和異構(gòu)性，以達(dá)到最優(yōu)化的執(zhí)行配置模式；② 能夠識別計算任務(wù)中各任務(wù)的并行性并調(diào)度任務(wù)的執(zhí)行；③ 能夠使具有不同計算類型的計算資源能互相進(jìn)行數(shù)據(jù)通信，協(xié)調(diào)運(yùn)行；④ 最終目的是讓整個應(yīng)用程序的執(zhí)行效率最高，計算耗時最短。

由此可見，基于GPU分析平臺開發(fā)的關(guān)鍵是合理分配CPU與GPU的計算任務(wù)，因此本文設(shè)計和開發(fā)的分析平臺讓CPU執(zhí)行數(shù)據(jù)讀寫和邏輯控制等順序型任務(wù)，而計算密集的大規(guī)模數(shù)據(jù)運(yùn)算則讓GPU并行執(zhí)行。

在HSNAS(GPU)平臺上實現(xiàn)GPU并行計算的步驟如下：首先對需要并行計算的數(shù)據(jù)進(jìn)行CPU端的預(yù)處理，由此轉(zhuǎn)換成能被GPU端識別的數(shù)據(jù)文件；接著在GPU顯存中開辟數(shù)據(jù)空間，將處理好的數(shù)據(jù)傳輸?shù)紾PU的顯存中；然后讀取計算配置文件，調(diào)用編寫的內(nèi)核函數(shù)(kernel)，開啟GPU細(xì)粒度的并行執(zhí)行任務(wù)；待GPU并行計算完成，再通過數(shù)據(jù)傳輸，將運(yùn)算結(jié)果從GPU顯存中拷貝回CPU內(nèi)存里，接著繼續(xù)進(jìn)行余下任務(wù)的操作。

與單純采用CPU多核并行架構(gòu)[6]的并行模式相比，本文開發(fā)的HSNAS(GPU)平臺架構(gòu)相當(dāng)于引入了更多數(shù)量的線程(thread)進(jìn)行數(shù)據(jù)運(yùn)算，GPU也根據(jù)相應(yīng)的優(yōu)化準(zhǔn)則進(jìn)行并行化計算，因而其計算效率能夠得到明顯地提高。

1.3CUDA并行編程模式

目前主流的GPU編程模式有CUDA、OpenGL、Direct Compute等。CUDA[7]是一種并行編程模型和軟件開發(fā)環(huán)境，它提供了細(xì)粒度的數(shù)據(jù)并行化和線程并行化計算能力。本文開發(fā)的HSNAS(GPU)分析平臺采用CUDA編程模式，且具備以下三個基本能力：① 任務(wù)管理：在CPU和GPU之間進(jìn)行任務(wù)劃分，并調(diào)度任務(wù)的執(zhí)行；② 并行性描述：描述CPU和GPU之間的任務(wù)級并行以及CPU和GPU內(nèi)部細(xì)粒度的線程級并行；③ 數(shù)據(jù)管理：能夠管理CPU和GPU的程序數(shù)據(jù)，并在兩者之間進(jìn)行數(shù)據(jù)通信。

基于CUDA架構(gòu)，對線程(thread)、線程塊(block)以及內(nèi)核函數(shù)(kernel)進(jìn)行程序設(shè)計相當(dāng)于對GPU的硬件單元進(jìn)行控制設(shè)計。

2基于GPU的纖維模型非線性有限元理論

2.1虛加剛性彈簧的位移增量迭代算法

結(jié)構(gòu)非線性靜力分析的過程實際上是對非線性方程組求解的過程。已有的非線性方程組求解方法以力為增量的迭代法居多，其中Newton-Raphson法在求解結(jié)構(gòu)力-位移曲線上升段時，適應(yīng)性比較強(qiáng)。但是當(dāng)臨近極限強(qiáng)度以及軟化下降段時，由于剛度矩陣失去正定性，采用以力為增量的迭代法導(dǎo)致計算不收斂，此稱之為負(fù)剛度問題。目前已有一些負(fù)剛度問題的求解方法[8-9]，但這些方法也受到許多條件限制。本文提出虛加剛性彈簧的位移增量迭代算法，是一種有效解決負(fù)剛度問題的新方法。

非線性靜力問題可歸結(jié)為解代數(shù)方程組，一般采用增量格式：

[K]{Δδ}={ΔF}

(1)

式中:[K]為總剛度陣；{ΔF}節(jié)點荷載增量；{Δδ}為節(jié)點位移增量。

在原結(jié)構(gòu)體系施加變載的各自由度方向上虛設(shè)一個剛度較大的彈簧，并且假設(shè)彈簧的剛度在受力過程中保持不變。原結(jié)構(gòu)體系虛設(shè)彈簧后的合成體系如圖2所示。

圖2　原結(jié)構(gòu)虛設(shè)彈簧示意圖Fig.2 Structure add virtual spring

虛加彈簧作用之后，式(1)改寫為,

(2)

(3)

式中:[K]s為虛加彈簧的剛度矩陣，表達(dá)式如下,

[K]s=diag(ξ1ks1,ξ2ks2,…,ξnksn)

(4)

式中:ksi表示在第i個自由度上虛加的彈簧剛度，一般取為剛性彈簧，剛度較大，保證在整個求解過程中，彈簧的內(nèi)力與位移關(guān)系為線性。ξi為彈簧剛度的取值系數(shù)，即：

(5)

此外，虛加彈簧后的結(jié)構(gòu)節(jié)點荷載向量也要考慮由于節(jié)點位移引起的彈簧虛加荷載增量：

(6)

因此式(2)展開寫成：

[K]{Δδ}+[K]s{Δδ}={ΔF}+{ΔF}s

(7)

從式(7)左邊項看出，彈簧剛度矩陣[K]s在整個求解過程中保持不變，結(jié)構(gòu)的剛度矩陣[K]隨著結(jié)構(gòu)非線性受力變形而發(fā)生改變。當(dāng)結(jié)構(gòu)達(dá)到繼續(xù)承載力時，[K]將從正定轉(zhuǎn)為負(fù)定，有限元剛度方程的系數(shù)矩陣出現(xiàn)奇異，也就是前文提到的負(fù)剛度問題。為了克服這一問題，修正式(7)，在方程兩邊同時加上[K]0{Δδ}，其中[K]0為結(jié)構(gòu)的初始剛度矩陣，這時式(7)變?yōu)椋?/p>

([K]0+[K]s){Δδ}={ΔF}+{ΔF}s+{ΔR}

(8)

式中{ΔR}為不平衡力向量,

{ΔR}=([K]0-[K]){Δδ}

(9)

由式(8)可見，通過改寫的方程左端系數(shù)矩陣在分析過程中始終保持正定。

式(8)的迭代格式,

(10)

式中:i為迭代次數(shù)；j為荷載增量步數(shù)。

如果結(jié)構(gòu)只受位移控制，上式可簡化為,

(11)

迭代的收斂準(zhǔn)則采用位移增量為判斷標(biāo)準(zhǔn),

(12)

式中:tol為迭代收斂的容差，根據(jù)計算要求精度來確定。

上述修正后的算法對系數(shù)矩陣只需要在初始時進(jìn)行一次分解，可保證節(jié)省存儲空間的同時，顯著提高計算效率。

2.2Newmark動力時程積分算法

非線性動力問題平衡方程的增量格式為：

(13)

上式求解方法常采用逐步積分法，本文采用Newmark-β法，時間微段的速度增量和位移增量為：

(14)

{Δu}=

(15)

將式(14)和式(15)代入式(13)：

[K*]{Δu}={ΔF*}

(16)

式中:[K*]為等效剛度，{ΔF*}為等效荷載增量。

有關(guān)Newmark-β法的更多詳細(xì)介紹，可參考文獻(xiàn)[10]，在此本文不作贅述。

2.3纖維模型單元技術(shù)

2.3.1基于Timoshenko梁的纖維單元模型

基于Timoshenko梁的纖維模型基本假定：① 滿足平截面假定，變形前垂直于梁中心線的截面變形后仍保持平面，但不再與梁中心線保持垂直，而是發(fā)生翹曲；② 一個單元沿長度方向劃分若干個積分點，各纖維處于單軸應(yīng)力狀態(tài)；③ 忽略纖維的黏結(jié)滑移影響；④ 剪力和剪應(yīng)變沿截面均勻分布，在截面層次上考慮剪切，與其他變形分量不耦合。

截面的剪切應(yīng)變?yōu)椋?/p>

γ=dω/dx-φ

(17)

梁截面的曲率為：

χ=-dφ/dx

(18)

因此，考慮剪切變形影響的梁單元任一截面的變形為：

ds(x)=

(19)

根據(jù)Timoshenko梁垂直于中性面的截面變形后仍保持為平面的假定，纖維的應(yīng)變?yōu)椋?/p>

εf(x,yif,zif)=b(x,yif,zif)ds(x)

(20)

式中:b(x,yif,zif)為考慮剪切效應(yīng)的幾何變換矩陣。

(21)

式中:zif和yif為纖維的幾何中心坐標(biāo)，如圖3所示。

圖3　纖維模型單元Fig.3 Fiber element

纖維的應(yīng)力表示為：

σf(x,yif,zif)=Df(x,yif,zif)εf(x,yif,zif)

(22)

式中:Df(x,yif,zif)為纖維的材料剛度。

(23)

式中:Eif和Gif分別為纖維的切線模量和剪切模量。

對截面內(nèi)所有纖維積分，獲得截面剛度矩陣和截面抗力矩陣：

(24)

(25)

在單元層次上，單元剛度矩陣和單元抗力向量：

(26)

(27)

式中：B(x)為應(yīng)變矩陣。

2.3.2鋼筋混凝土纖維本構(gòu)模型

纖維模型的計算精度和效率還取決于材料本構(gòu)模型的選取。如下所述，本文選取的材料本構(gòu)模型均是較為精細(xì)卻又不失較高計算效率的本構(gòu)模型。

(1) 鋼筋本構(gòu)模型

本文所采用的修正Menegotto-Pinto鋼筋本構(gòu)模型[11]，如圖4所示，能夠考慮鋼筋屈服、等向強(qiáng)化、包興格效應(yīng)等特性，該模型不僅分析效率高并且與試驗數(shù)據(jù)吻合較好。

圖4　修正的Menegotto-Pinto鋼筋本構(gòu)Fig.4 Menegotto-Pinto steel model

(2) 混凝土本構(gòu)模型

本文采用如圖5(a)所示的混凝土受壓滯回模型綜合了Kent-Park模型[12]、Blakeley兩折線模型[13]和焦點模型的優(yōu)點。一方面，通過修改骨架曲線中的強(qiáng)度增大系數(shù)和應(yīng)變軟化段斜率，可以模擬箍筋對混凝土的約束作用；另一方面，采用的兩折線卸載路徑和焦點模型，可以模擬混凝土在滯回過程中的剛度退化和滯回耗能性能?；炷潦芾瓬啬Ｐ筒捎萌鐖D5(b)所示的Yassin模型[14]。該模型較精細(xì)，可以考慮受拉剛化效應(yīng)，模型參數(shù)也較少，能兼顧計算精度和計算效率。

圖5　混凝土滯回本構(gòu)Fig.5 Hysteretic concrete stress-strain relation

3基于GPU的非線性分析并行化實現(xiàn)

3.1計算數(shù)據(jù)與線程之間的映射關(guān)系

從有限元離散化角度出發(fā)，一個精確的有限元結(jié)構(gòu)模型包含數(shù)十萬個單元和節(jié)點，因此整個模型將有大量的全局自由度。將計算數(shù)據(jù)和線程一一對應(yīng)是一種可行且高效的GPU并行策略。并行策略可分為以下三種計算模式：① 線程與單元對應(yīng)；② 線程與節(jié)點對應(yīng)；③ 線程與自由度對應(yīng)。前兩種模式的計算數(shù)據(jù)對象可以分割并分發(fā)到每個處理器的并行計算平臺上。然而，GPU的并行模型是一個細(xì)粒度模型，用相鄰的線程去映射相鄰的數(shù)據(jù)是最有效的。因此本文將每一個全局自由度視為一個獨立的計算單位，它的變量數(shù)據(jù)更新是獨立的，即在全局坐標(biāo)下的矩陣/向量(剛度矩陣、力、位移等)的元素是獨立的。密集的算術(shù)運(yùn)算使這些數(shù)據(jù)特別適合GPU線程級的并行實現(xiàn)。

計算數(shù)據(jù)和線程存在如圖6所示的一一映射關(guān)系。圖中示例為一個線程塊(block)設(shè)置了128個線程(thread)，根據(jù)數(shù)據(jù)規(guī)模([n]+128-1)/128配置線程塊數(shù)量，保證線程塊數(shù)量為一個整數(shù)，因此所有線程總數(shù)大于或等于實際計算數(shù)據(jù)總量。此外，GPU的并行計算能力還依賴于GPU上固定集成的硬件資源，因此對硬件資源合理分配，才能保證計算程序效率的最大化。

圖6　計算數(shù)據(jù)與線程之間一一映射關(guān)系Fig.6 Mapping relationship between calculated data and threads

3.2基于GPU的線性方程組并行求解器

PCG(Preconditioned Conjugate Gradients, 預(yù)處理共軛梯度法)法是目前求解大型對稱正定稀疏線性方程組最為有效的算法之一[15]。該算法存在內(nèi)在并行性，易于實現(xiàn)并行化，且要求的存儲空間和計算量相對較小。圖7給出了基于GPU并行PCG求解器的流程圖。流程圖中空心箭頭指向的表示在GPU中執(zhí)行的內(nèi)核函數(shù)。

圖7　基于GPU的線性方程組并行求解器流程Fig.7 PCG parallel solver based on GPU

從圖7看出，除了在迭代開始前和結(jié)束后，需要進(jìn)行CPU和GPU的數(shù)據(jù)傳輸，在迭代程序中，只需要在CPU中進(jìn)行少量計算以及進(jìn)行收斂控制外，其余計算基本都在GPU中執(zhí)行。該求解器的設(shè)計既可以充分發(fā)揮CPU的任務(wù)調(diào)度與邏輯判斷能力以及GPU強(qiáng)大的浮點運(yùn)算能力，又盡可能地減少CPU和GPU之間的數(shù)據(jù)通信，減少時間開銷，提高并行效率。

3.3基于GPU的動力時程分析程序

逐步積分法的每個時間步計算將消耗大量的計算時間。因此本文開發(fā)的動力分析程序?qū)⒄麄€時間步積分在GPU中完成，計算流程如圖8所示。每個時間步下等效剛度和等效荷載采用內(nèi)核函數(shù)計算，位移增量則采用并行PCG求解器加速求解，新的速度和加速度也采用相應(yīng)的內(nèi)核函數(shù)計算，GPU中的計算直至整個時間步循環(huán)結(jié)束為止。

圖8　基于GPU的動力時程計算程序流程Fig.8 Flowchart of dynamic analysis program based on GPU

4平臺測試

4.1平臺參數(shù)

計算平臺：CPU為Intel i5-2300，頻率為2.8 GHz，內(nèi)存為4.00 GB；GPU為NVIDIA GeForce GTX 460，336個GPU流處理器核心數(shù)目，計算能力2.1，流處理器頻率1.4 GHz，顯存為1.0 GB，顯存帶寬為115.2 GB/s。

軟件開發(fā)環(huán)境：① Microsoft Windows 7 64位操作系統(tǒng)；② Microsoft Visual Studio 2012開發(fā)環(huán)境；③ 全面支持CUDA Fortran的PGI Accelerator Visual Fortran編譯器；④ CUDA 5.0驅(qū)動(driver)。

4.2算例模型

為進(jìn)行初步驗證，本文算例模型保持為彈性。非線性驗證參考論文的第Ⅱ篇。

算例為一棟10層鋼筋混凝土框架結(jié)構(gòu)，結(jié)構(gòu)平面和立面如圖9所示?？蚣軐痈邽? m，跨度為6 m，框架柱截面450 mm×450 mm，混凝土采用C40，框架梁截面300 mm×500 mm，混凝土采用C30，樓面承受恒荷載3.0 kN/m2，活荷載2.0 kN/m2。

圖9　框架模型Fig.9 Frame model

4.3靜力計算

為了對比計算精度，采用本文開發(fā)的HSNAS(GPU)平臺計算得到的①軸和軸相交的角柱豎向變形與有限元軟件ABAQUS得到的計算結(jié)果對比，如表1所示。從表中看出HSNAS(GPU)平臺與ABAQUS計算結(jié)果完全吻合，兩者之間的誤差幾乎為0，初步驗證了本文開發(fā)HSNAS(GPU)平臺以及所提算法的正確性，可用來進(jìn)行下一步的分析。

采用不同網(wǎng)格劃分的計算模型，圖10給出計算耗時及加速比與自由度大小的關(guān)系。由圖中看出HSNAS(GPU)平臺計算耗時遠(yuǎn)遠(yuǎn)少于傳統(tǒng)CPU串行平臺，計算效率非常明顯，隨著自由度增大，獲得了30倍以上的加速比。

表1?、佥S和軸角柱豎向變形計算結(jié)果對比

Tab.1 Comparison of vertical deformation of column in axis ① and 　mm

表1　①軸和軸角柱豎向變形計算結(jié)果對比

樓層HSNAS(GPU)ABAQUS誤差(%)1-0.425-0.4250.0002-0.809-0.8090.0003-1.152-1.1520.0004-1.454-1.4540.0005-1.713-1.7130.0006-1.929-1.9290.0007-2.102-2.1020.0008-2.231-2.2310.0009-2.317-2.3170.00010-2.360-2.3600.000

圖10　計算耗時對比及HSNAS(GPU)平臺加速比Fig.10 Comparison of computation time and speed up

4.4動力計算

采用本文開發(fā)的HSNAS(GPU)平臺進(jìn)行動力時程分析，地震波采用El Centro(N-S)波，加速度峰值調(diào)整為35 cm/s2(7度多遇烈度)，結(jié)構(gòu)采用瑞利阻尼，前兩階阻尼比取0.05。HSNAS(GPU)平臺得到的頂層位移、速度、加速度時程曲線與ABAQUS計算結(jié)果對比如圖11所示。

圖11　HSNAS(GPU)平臺與ABAQUS計算結(jié)果對比Fig.11 Comparison of ABAQUS and HSNAS(GPU) platform

從圖11中看出，HSNAS(GPU)平臺得到的動力響應(yīng)與ABAQUS有限元計算結(jié)果符合較好。同樣，為了測試HSNAS(GPU)平臺的動力問題計算效率，采用不同網(wǎng)格劃分的計算模型，統(tǒng)計計算耗時及加速效果如圖12所示。能得到與靜力計算類似的曲線結(jié)果，隨著自由度增大，獲得了25倍以上的加速比。

圖12　計算耗時對比及HSNAS(GPU)平臺加速比Fig.12 Comparison of computation time and speed up

5結(jié)論

(1) 利用GPU并行計算能力，開發(fā)了基于GPU的鋼筋混凝土框架結(jié)構(gòu)非線性靜力、動力計算的高精度高效率分析平臺HSNAS(GPU)。讓CPU處理數(shù)據(jù)讀寫和邏輯控制等順序型任務(wù)，GPU則處理計算密集的大規(guī)模數(shù)據(jù)操作，將結(jié)構(gòu)有限元分析各步驟進(jìn)行了“粗粒度”的任務(wù)劃分。從有限元離散化角度出發(fā)，將每一個全局自由度視為一個獨立的計算單位，建立了計算數(shù)據(jù)與GPU線程一一對應(yīng)映射關(guān)系，實現(xiàn)了對硬件資源的合理分配，保證了計算程序效率的最大化。

(2) 對非線性求解中常見的負(fù)剛度問題，提出一種新的適用于GPU計算的位移增量迭代算法。針對動力問題，提出了基于GPU的Newmark時間積分算法。結(jié)合纖維模型單元技術(shù)，引入剪切、扭轉(zhuǎn)變形以及材料非線性的影響，采用了多個適用于GPU并行的計算策略并且編寫了內(nèi)核函數(shù)，在HSNAS(GPU)平臺上開發(fā)完成了鋼筋混凝土框架結(jié)構(gòu)非線性有限元程序。該程序能夠識別出數(shù)據(jù)并行程度高且計算密集的部分，通過配置語句由CPU主程序調(diào)用GPU的線性方程組并行求解器，從而實現(xiàn)了結(jié)構(gòu)非線性靜、動力分析的加速。

(3) 在HSNAS(GPU)平臺上對鋼筋混凝土框架結(jié)構(gòu)算例模型進(jìn)行了彈性的靜力、動力分析進(jìn)行初步驗證。HSNAS(GPU)平臺獲得的計算結(jié)果與有限元軟件對比吻合較好，隨著計算模型的增大，分別獲得了30和25倍以上的提速效率。

(4) 鑒于GPU的巨大計算潛能，為進(jìn)一步提高計算平臺體系的執(zhí)行效率，本文對HSNAS(GPU)平臺進(jìn)行了相應(yīng)的程序優(yōu)化策略、更深入地非線性分析以及平臺精度和效率的驗證，這些研究內(nèi)容將在論文的第Ⅱ篇詳述。

參 考 文 獻(xiàn)

[1] 李云貴. 工程結(jié)構(gòu)設(shè)計中的高性能計算[J]. 建筑結(jié)構(gòu)學(xué)報, 2010, 31(6): 89-95.

LI Yun-gui. High-performance computing in structural design [J]. Journal of Building Structures, 2010, 31(6): 89-95.

[2] 李紅豫, 滕軍, 李祚華. 基于CPU-GPU異構(gòu)平臺的高層結(jié)構(gòu)地震響應(yīng)分析方法研究[J]. 振動與沖擊,2014,33(13)：86-91.

LI Hong-yu, TENG Jun, LI Zuo-hua. Analysis method for seismic response of high-rise structure based on CPU-GPU heterogeneous platform [J]. Journal of Vibration and Shock, 2014, 33(13): 86-91.

[3] Spacone E, Fillippou F C, Taucer F F. Fiber beam-column model for nonlinear analysis of RC frames: Part I. Formulation [J]. Earthquake Engineering & Structure Dynamics, 1996, 25: 711-725.

[4] Petrangeli M, Pinto P E, Ciampi V. Fiber element for cyclic bending and shear of RC structures. I: Theory [J]. Journal of Engineering Mechanics, 1999, 125(9): 994-1001.

[5] 殷有泉. 非線性有限元基礎(chǔ)[M]. 北京：北京大學(xué)出版社，2007.

[6] 丁峻宏, 宋雅麗, 王惠, 等. 大規(guī)模柴油機(jī)動力學(xué)抗沖擊并行仿真計算[J]. 振動與沖擊, 2014, 33(2)：163-167.

DING Jun-hong, SONG Ya-li, WANG Hui, et al. Parallel computing for large scale anti-shock dynamic simulation of diesel engine [J]. Journal of Vibration and Shock，2014,33(2): 163-167.

[7] Nvidai Corporation. CUDA C Programming Guide [EB/OL]. http://docs.nvidia.com/cuda/pdf/CUDA_C_Programming_Guide.pdf. July 2013.

[8] Crisfield M A. An arc-length method including line searches and accelerations [J]. International Journal of Numerical Methods in Engineering, 1983, 19(9): 1269-1289.

[9] 杜修力, 曹惠, 金瀏. 力-變位關(guān)系全過程模擬的有限元位移控制新方法[J]. 工程力學(xué), 2012, 29(1): 1-6.

DU Xiu-li, CAO Hui, JIN Liu. A new finite element displacement control method of the whole process simulation of force-displacement relation [J]. Engineering Mechanics, 2012, 29(1): 1-6.

[10] 克拉夫 R, 彭津 J. 結(jié)構(gòu)動力學(xué)[M].2版.王光遠(yuǎn), 等,譯. 北京：高等教育出版社，2006.

[11] Menegotto M, Pinto P E, Slender R C. Compressed members in biaxial bending [J]. Journal of Structural Division, ASCE, 1977, 103(3): 587-605.

[12] Kent D C, Park R. Flexural members with confined concrete [J]. Journal of the Structural Division, ASCE,1971,97(7): 1969-1990.

[13] Blakely R W G, Park R. Prestressed concrete sections with cyclic flexure [J]. Journal of the Structural Division, ASCE, 1973, 99(8): 1717-1742.

[14] Yassin M H M. Nonlinear analysis of prestressed concrete structures under monotonic and cycling loads [D]. Berkeley，CA:University of California, 1994.

[15] Benzi M, Tuma M. A comparative study of sparse approximate inverse preconditioners [J]. Applied Numerical Mathematics, 1999, 30(2): 305-340.

基金項目：國家自然科學(xué)基金重大國際(中美)合作研究項目(51261120374)；國家自然科學(xué)基金面上項目(51278155;51378007)

收稿日期：2015-04-07修改稿收到日期：2015-07-07

通信作者滕軍 男，博士，教授，博士生導(dǎo)師，1962年生

中圖分類號:TU375

文獻(xiàn)標(biāo)志碼：A

DOI:10.13465/j.cnki.jvs.2016.14.008

An efficient platform HSNAS(GPU) for nonlinear static and dynamic analysis of reinforced concrete frames—Ⅰ. Program development

LI Hong-yu, TENG Jun, LI Zuo-hua

(Shenzhen Graduate School, Harbin Institute of Technology, Shenzhen 518055, China)

Abstract:The traditional serial computation platform has some disadvantages such as low accuracy and dramatically high time consuming, which have hindered the development of structural nonlinear dynamics analysis. In order to achieve a higher computing accuracy and save calculation time in the process of nonlinear dynamic analysis of reinforced concrete (RC) frames, a simulation platform HSNAS(GPU) based on the graphics processing unit (GPU) was developed. For static analysis, a GPU-based incremental displacement algorithm was introduced to deal with negative stiffness problems, and a parallel preconditioned conjugate gradients (PCG) solver was developed. For dynamic analysis, the GPU-based Newmark-beta algorithm was presented. The fiber beam model was improved by considering the effects of shear and torsion. In addition, the constitutive models of steel and concrete were developed. The results of numerical examples illustrate that the developed platform HSNAS(GPU) could improve the efficiency of nonlinear static and dynamic analysis besides satisfying the accuracy requirement.

Key words:reinforced concrete structure; nonlinear finite element analysis; simulation platform; GPU; fiber model

第一作者 李紅豫 女，博士生，1985年生