一種非線性顯式分層殼單元及其GPU并行計算實現(xiàn)

曹勝濤， 李志山， 楊志勇

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

鋼筋混凝土剪力墻是一種在高層建筑中被廣泛應(yīng)用的抗側(cè)力構(gòu)件[1]。在外荷載作用下剪力墻的損壞類別和程度是判斷抗震性能的重要指標(biāo)[2]。

在高層結(jié)構(gòu)中剪力墻空間形狀多樣且邊界條件復(fù)雜，采用非線性宏觀單元分析時，其參數(shù)往往不易確定。基于應(yīng)力-應(yīng)變層次的混凝土非線性本構(gòu)模型參數(shù)可由材料試驗獲得，與構(gòu)件形式無關(guān)；將其與分層殼元結(jié)合可較全面地反映剪力墻面外和面內(nèi)的耦合特性。Hu等[3]基于彈塑性力學(xué)和彌散開裂概念得到一種混凝土非線性本構(gòu)模型，與分層殼元結(jié)合對在單調(diào)荷載作用下鋼筋混凝土結(jié)構(gòu)進(jìn)行了分析。Polak等[4]給出了一種基于修正斜壓場理論[5]的鋼筋混凝土非線性分層殼元，并通過與試驗對比分析了單元的準(zhǔn)確性。Lu等[6]在非線性分層殼元中采用了一種損傷力學(xué)和彌散開裂概念結(jié)合的本構(gòu)模型模擬混凝土力學(xué)特性；并在OpenSees[7]完成開發(fā)，對一個框架-核心筒超高層結(jié)構(gòu)進(jìn)行了分析；通過與商業(yè)程序MSC.Marc的計算結(jié)果進(jìn)行對比，驗證了單元的準(zhǔn)確性。

將塑性力學(xué)理論[8]和損傷力學(xué)[9]結(jié)合得到的混凝土彈塑損傷模型可反映混凝土的損傷程度和殘余變形。根據(jù)損傷變量的數(shù)量，混凝土損傷本構(gòu)模型可以分為單標(biāo)量模型[10-11]和雙標(biāo)量模型[12-13]。Lee等[14]利用Lubliner等[15]給出的三維強度準(zhǔn)則將混凝土單軸拉伸和單軸壓縮的力學(xué)特性三維化得到了一種雙標(biāo)量的混凝土三維彈塑性損傷本構(gòu)模型(下文簡稱：塑性損傷模型)。塑性損傷模型可較好地反映混凝土多軸動力荷載作用下的強度軟化、剛度退化、殘余變形、單邊效應(yīng)、裂面效應(yīng)、靜水壓力效應(yīng)、中主應(yīng)力效應(yīng)等力學(xué)特性。商業(yè)通用有限元分析程序ABAQUS在顯式分析模塊Explicit中開發(fā)實現(xiàn)了塑性損傷本構(gòu)模型。ABAQUS將平面應(yīng)力條件下的塑性損傷模型與顯式分層殼單元結(jié)合可較好地反映剪力墻在拉、壓、彎、剪、扭及各種荷載共同作用下的損傷和破壞模式，克服了大規(guī)模非線性隱式分析計算量大、收斂困難的技術(shù)瓶頸，因此在工程實踐中被廣泛地應(yīng)用[16-20]。由于塑性損傷模型是由單軸方向力學(xué)特性三維化得到，所以只提供了損傷因子d+和d-(本文稱之為廣義受拉損傷因子和廣義受壓損傷因子)。而眾所周知，混凝土構(gòu)件損壞類型一般包括拉彎型、壓彎型、壓剪型等類型；只通過d+和d-不能準(zhǔn)確判斷剪力墻損傷和破壞模式。同時ABAQUS Explicit 目前只支持CPU并行計算，不易實現(xiàn)細(xì)粒度并行計算，無法充分提高動力方程顯式積分的計算速度，進(jìn)行大規(guī)模計算時效率較低。

分別采用SAUSAGE和ABAQUS完成了上海地區(qū)某框架核心筒體系的超高層結(jié)構(gòu)罕遇地震動力非線性時程分析；通過與ABAQUS計算結(jié)果對比，驗證了本文殼元在大規(guī)模分析中的精度和效率；根據(jù)本文殼元得到混凝土損傷和鋼筋、鋼材的塑性變形，對核心筒在動力荷載作用下的損壞演化規(guī)律和破壞模式進(jìn)行了分析。

1 基于塑性損傷模型的損傷分類

平面應(yīng)力條件下塑性損傷模型在有效應(yīng)力空間中屈服函數(shù)和塑性勢函數(shù)如式(1)～式(3)所示

(1)

(2)

(3)

根據(jù)損傷力學(xué)的基本假定，真實應(yīng)力σij為

(4)

(5)

本文通過損傷發(fā)展時主應(yīng)力所在象限的差異對損傷進(jìn)行進(jìn)一步分類。

如圖1所示，屈服函數(shù)F在主應(yīng)力空間中，可根據(jù)受力狀態(tài)劃分為拉-拉、壓-拉(拉-壓)、壓-壓三個區(qū)域，分別對應(yīng)拉彎、剪切、壓彎的受力狀態(tài)。

(6)

(7)

(8)

(9)

圖1 塑性損傷模型的屈服面Fig.1 Yield surface of elastoplastic damage model

圖2 區(qū)域I加載時的屈服面Fig.2 Yield surface of loading in zone Ⅰ

圖3 區(qū)域III加載時的屈服面Fig.3 Yield surface of loading in zone Ⅲ

圖4 區(qū)域II加載時的屈服面Fig.4 Yield surface of loading in zone Ⅱ

塑性損傷本構(gòu)模型應(yīng)力更新算法采用Lee等[21]給出的應(yīng)力回映算法。混凝土模型參數(shù)按《混凝土結(jié)構(gòu)設(shè)計規(guī)范》(GB 50010—2010)附錄C.2[22](簡稱：《混規(guī)》)的標(biāo)準(zhǔn)值取值，塑性應(yīng)變和損傷取值參考任曉丹等[23]的做法。

2 基于塑性損傷模型的非線性顯式分層殼單元

本文殼元為4節(jié)點24自由度四邊形單元。如圖5所示，每節(jié)點有6個自由度(括號內(nèi)自由度排序)，其中包括：①面內(nèi)膜元的兩個平動自由度；②面外板元的一個出平面平動自由度和兩個轉(zhuǎn)動自由度；③一個面內(nèi)旋轉(zhuǎn)自由度。

圖5 殼元自由度Fig.5 The degrees of freedom of the shell element

如圖6所示，三維空間全局坐標(biāo)系(下標(biāo)為g)和局部坐標(biāo)系均為直角坐標(biāo)系，遵循右手螺旋法則。局部坐標(biāo)系位于殼元平面和截面的形心處，局部坐標(biāo)系z軸與單元所在面垂直，且z軸正向與zg軸正向夾角不大于90°；局部坐標(biāo)系x軸平行于單元某條邊，且x軸正向與xg軸正向夾角不大于90°；局部坐標(biāo)系y軸可由z軸和x軸通過右手螺旋法則確定。T為局部坐標(biāo)系與

(a)三維空間

(b)截面平面圖6 分層殼元Fig.6 The layered shell element

整體坐標(biāo)系的轉(zhuǎn)換矩陣。如圖 6(b)，局部坐標(biāo)系中zi為殼元截面第i層形心到截面形心的距離，hi為第i層厚度，h為截面總厚度。

鋼筋層采用縱橫纖維模擬，其本構(gòu)模型為單軸雙線性隨動強化模型。鋼筋層與混凝土層變形協(xié)調(diào)，不考慮兩者間的相對滑移。

本文殼元形函數(shù)如式(10)所示，為保證計算精度和效率，對常應(yīng)變場采用單點積分，對線性應(yīng)變場采用物理沙漏控制。本文構(gòu)造的殼元分為6層，每層均設(shè)置積分點，其中：①面內(nèi)部分參考了Belytschko等[24-25]構(gòu)造的高精度膜元；②面外部分參考了Hughes等[26]構(gòu)造的選擇積分的Mindlin板元；③面內(nèi)旋轉(zhuǎn)剛度則參考了Zienkiewicz等[27]的構(gòu)造方法。

雙線性形函數(shù)為

(10)

式中:ξ和η為母元坐標(biāo)；ξi和ηi為節(jié)點i在參考系中的取值(i=1,2,3,4)，如圖7所示。

圖7 等參元Fig.7 The isoparametric element

形函數(shù)對局部坐標(biāo)系x,y偏導(dǎo)為

(11)

式中:J為等參元雅克比矩陣。(注：此處下標(biāo)中“，”表示偏導(dǎo)，其它處均為標(biāo)號)

局部坐標(biāo)系中殼元節(jié)點j對應(yīng)的面內(nèi)變形子矩陣塊、面外彎曲變形子矩陣塊和面外剪切變形子矩陣塊，如式(12)～式(14)所示。

(12)

(13)

(14)

殼元面內(nèi)、面外變形子矩陣與局部坐標(biāo)系中整體變形矩陣的轉(zhuǎn)換矩陣如式(15)～式(18)所示。

(15)

(16)

(17)

(18)

單元面內(nèi)旋轉(zhuǎn)剛度子矩陣與局部坐標(biāo)系中整體剛度矩陣的轉(zhuǎn)換矩陣如式(19)和式(20)所示。

Tmd,j=[0 0 0 0 0 1]

(19)

(20)

本文殼元面內(nèi)變形矩陣、面外彎曲變形矩陣和面外剪切變形矩陣如式(21)～式(23)所示。

Bm=[Bm,1Bm,2Bm,3Bm,4]Tm

(21)

Bsb=[Bsb,1Bsb,2Bsb,3Bsb,4]Ts

(22)

Bss=[Bss,1Bss,2Bss,3Bss,4]Ts

(23)

將ξ=0和η=0代入式(21)和式(22)，得到單點積分處單元面內(nèi)和面外彎曲常應(yīng)變場矩陣B0m和B0sb。單元面外符合平截面假定，由面內(nèi)和面外的耦合可得到殼元第i層積分點處應(yīng)變增量

Δεi=(ziB0sb+B0m)TΔde

(24)

式中:Δde單元節(jié)點位移增量。

將Δεi代入塑性損傷模型中完成應(yīng)力更新，得到第i層積分點處應(yīng)力增量Δσi、混凝土剛度退化系數(shù)di及相關(guān)歷史量。

對Δσi進(jìn)行積分得到此部分對殼元面內(nèi)節(jié)點力增量的貢獻(xiàn)

(25)

式中:Ae為單元面積。(注：本文矩陣上標(biāo)“t”表示矩陣轉(zhuǎn)置)

(26)

通過截面積分得到單點積分對單元面外彎曲節(jié)點力增量的貢獻(xiàn)

(27)

(28)

本文殼元不考慮面外剪切非線性。為避免面外剪切閉鎖，殼元面外剪切部分采用單點積分，對單元節(jié)點力增量貢獻(xiàn)

(29)

(30)

式中:Dss,G和ks分別為殼元截面剪切剛度矩陣、材料剪切模量和剪切校正因子(矩形一般取5/6)。

為避免計算奇異并保持單點積分，本文殼元面內(nèi)旋轉(zhuǎn)自由度參考了Zienkiewicz等[27]的構(gòu)造方法。面內(nèi)旋轉(zhuǎn)自由度與其它自由度變形解耦，但本文考慮了d對面內(nèi)旋轉(zhuǎn)剛度的折減，此部分對單元節(jié)點力增量的貢獻(xiàn)

(31)

(32)

式中:Dmd為面內(nèi)旋轉(zhuǎn)剛度矩陣；E0,i為第i層混凝土初始彈性模量；kd為面內(nèi)旋轉(zhuǎn)剛度系數(shù)(本文取0.027)。

綜上所述，全局坐標(biāo)系中本文殼元單元節(jié)點力增量

(33)

3 本文殼元的開發(fā)實現(xiàn)及驗證

動力方程求解格式采用王進(jìn)廷和杜修力[28]給出的一種適用于一般阻尼體系且精度為兩階的顯式差分格式。本文上述本構(gòu)模型、殼元和動力方程求解格式在計算過程中數(shù)據(jù)相關(guān)性低，可利用顯卡計算內(nèi)核多的特點實現(xiàn)細(xì)粒度并行計算，提高分析效率。SAUSAGE利用NVIDIA提供的CUDA實現(xiàn)了結(jié)構(gòu)非線性顯式動力分析的CPU+GPU異構(gòu)并行計算，步驟為：①將相關(guān)數(shù)據(jù)由CPU內(nèi)存拷貝到GPU顯存中；②單元內(nèi)力更新和動力方程求解分別在單元和自由度層次實現(xiàn)了細(xì)粒度的GPU并行計算，如圖 8所示；③將計算結(jié)果由GPU顯存拷貝到CPU內(nèi)存，進(jìn)而寫入硬盤存儲。需要說明的是，目前SAUSAGE只支持單顯卡并行計算，因此其計算能力會受到單顯卡性能的限制。

圖8 非線性顯式動力計算流程圖Fig.8 Diagram of nonlinear explicit dynamic analysis

為測試開發(fā)內(nèi)容的正確性，在SAUSAGE和ABAQUS(采用S4R單元)中建立了測試benchmark算例，如圖 9所示，單個四邊形殼元邊長1 m×1 m、厚度0.1 m。底部固定，對單元頂部節(jié)點進(jìn)行先壓剪后拉剪的位移加載，加載時程如圖 10所示。

圖9 測試模型Fig.9 Test model

圖10 測試模型的加載路徑Fig.10 Load path of test model

SAUSAGE和ABAQUS計算得到單元y向應(yīng)變-y向應(yīng)力關(guān)系和單元x向基底剪力時程，如圖11所示，兩者結(jié)果一致，驗證了本文開發(fā)內(nèi)容的正確性。

4 單片剪力墻損傷與破壞分析

Lefas等[29]對鋼筋混凝土剪力墻進(jìn)行了單調(diào)加載試驗，利用本文殼元對試驗中編號為SW-16和SW-23的剪力墻進(jìn)行擬靜力非線性有限元分析并與試驗結(jié)果進(jìn)行了對比。

圖11 x向基底剪力時程Fig.11 Base shear history in x direction

SW-16和SW-23軸壓比0.2，厚度分別70 mm和65 mm，幾何尺寸和有限元模型如圖12和圖13所示，配筋詳見文獻(xiàn)[29]。SW-16和SW-23混凝土塑性損傷本構(gòu)參數(shù)分別取《混規(guī)》中C80和C75的標(biāo)準(zhǔn)值；鋼筋單軸隨動強化模型參數(shù)取《混規(guī)》中HRB500的標(biāo)準(zhǔn)值。

如圖14所示，SAUSAGE計算得到頂點水平位移和水平反力的關(guān)系與試驗結(jié)果基本一致。相對于SW-23而言，SW-16剛度大而延性小。SW-16和SW-23從損傷到破壞的過程，如表1所示，分別屬于典型的剪切破壞和彎曲破壞[30]。

表1 剪力墻的損傷和破壞Tab.1 Damage and failure of the shear walls

圖12 剪力墻SW-16Fig.12 Shear wall SW-16

圖13 剪力墻SW-23Fig.13 Shear wall SW-23

圖14 頂點水平位移和水平反力Fig.14 The horizontal displacement and force

試驗過程中的裂縫類型和塑性損傷模型損傷類型的對應(yīng)關(guān)系，如表2所示。通過對比試驗，SAUSAGE利用本文殼元較準(zhǔn)確地反映了剪力墻剪切破壞模式和彎曲破壞模式中混凝土各類損傷、鋼筋屈服發(fā)生順序和特點，總結(jié)如表3所示，驗證了本文殼元及損傷分類的合理性。

表2 剪力墻裂縫類型和損傷類型Tab.2 The crack type and damage types of shear walls

剪力墻剪跨比、軸壓比、配筋等因素均會影響剪力墻破壞模式。一定條件下，剪力墻可能出現(xiàn)介于彎曲破壞和剪切破壞之間的彎剪破壞，可結(jié)合表3，根據(jù)試驗或數(shù)值分析結(jié)果進(jìn)行判斷。

5 框架核心筒結(jié)構(gòu)動力損傷和破壞分析

上海地區(qū)某框架核心筒結(jié)構(gòu)，如圖15所示，共49層，總高228 m，配筋取SATWE計算結(jié)果，配筋后模型總重1.52×105t；抗震設(shè)防烈度為7度，場地類別為Ⅳ，場地分組為第一組。有限元模型中梁單元數(shù)53 298、殼元數(shù)140 640、自由度數(shù)890 490。在SAUSAGE計算模型的基礎(chǔ)上，通過接口程序SAUSAGE-SSG2ABA得到了單元和自由度數(shù)量相同的ABAQUS計算模型。通過與ABAQUS對比，驗證SAUSAGE大規(guī)模非線性動力求解的正確性和計算效率。

表3 剪力墻的剪切和彎曲破壞Tab.3 Shearing and bending failure of shear walls

圖15 某框架核心筒結(jié)構(gòu)Fig.15 A frame-core tube structure

如表4所示，兩軟件模態(tài)分析結(jié)果中前10階振型特征周期對比差異小于3%，驗證了兩軟件計算模型的一致性。

表4 SAUSAGE和ABAQUS模態(tài)分析結(jié)果對比Tab.4 Comparison between modal analysisresults of ABAQUS and SAUSAGE

采用的地震動加速度時程記錄持時35 s，如圖 16所示。x向為主震方向，幅值2.0 m/s2；y向為次震方向，幅值1.7 m/s2；z向為豎向，時程與x相同，但幅值為1.3 m/s2。

圖16 地震動加速度時程Fig.16 Seismic acceleration history

5.1 罕遇地震非線性計算結(jié)果

使用同一臺計算機(CPU為i7-4770，4核8線程，內(nèi)存32 G；GPU為NVIDIA GTX970，流處理器單元數(shù)1 664，顯存4 G)，分別采用SAUSAGE和ABAQUS對此結(jié)構(gòu)進(jìn)行罕遇地震非線性時程分析。SAUSAGE采用GPU并行計算，耗時132 min；ABAQUS采用8核CPU并行計算，耗時751 min。如圖17和圖18所示，兩個軟件計算得到結(jié)構(gòu)最大層間位移和最大層間剪力基本一致，其中最大層間位移角小于1/100，滿足《建筑抗震設(shè)計規(guī)范》(GB 50011—2010)[31]中“大震不倒”抗震性能的要求。

如圖19所示，兩個軟件得到廣義受拉損傷d+和廣義受壓損傷d-的分布和數(shù)值基本一致?；炷翉V義受拉損傷范圍較大，主要出現(xiàn)在連梁、核心筒底部、核心筒收進(jìn)處、高階振型影響較大的上部等位置；廣義受壓損傷范圍較小，主要出現(xiàn)在連梁和核心筒收進(jìn)位置。在保證計算規(guī)模、計算結(jié)果和ABAQUS一致的前提下，SAUSAGE計算速度約為ABAQUS計算速度的5.89倍；驗證了SAUSAGE大規(guī)模非線性動力求解精度和計算效率。

圖17 最大層間位移角Fig.17 Maximun story drift

圖18 最大層間剪力Fig.18 Maximun story shear

如圖19所示，ABAQUS塑性損傷模型只能給出廣義受拉損傷因子和廣義受壓損傷因子，無法直接反映剪力墻的受剪破壞模式。

圖19 核心筒損傷(左：ABAQUS，右：SAUSAGE)Fig.19 Damage of the core tube(Left：ABAQUS，Right：SAUSAGE)

5.2 核心筒損傷演化分析

SAUSAGE計算得到核心筒受拉損傷、受壓損傷和受剪損傷在地震動作用0 s,10 s,20 s和30 s的損傷云圖，如表5所示。

在地震作用的初始時刻(即在施工荷載作用完成時刻)，核心筒的收進(jìn)處和上部樓層連梁的兩端出現(xiàn)一定的受拉損傷。在地震作用下，受拉損傷在原樓層開始蔓延，同時也向下部樓層的筒體發(fā)展，主要出現(xiàn)在核心筒邊緣的墻肢、連梁端部、核心筒根部等受拉區(qū)域。地震作用至10 s后，連梁中部和核心筒收進(jìn)處開始出現(xiàn)受剪損傷。動力響應(yīng)過程中，除局部外，基本未出現(xiàn)明顯的受壓損傷。

5.3 典型樓層核心筒的損傷和破壞分析

如圖20所示，核心筒第1層的受拉損傷集中在筒體根部；受壓損傷范圍和數(shù)值均較小，主要分布主震向墻肢的底部邊緣處；受剪損傷分布在連梁處；筒體鋼筋未屈服。

表5 核心筒損傷演化Tab.5 Damage evolution of the core tube

圖20 底部樓層混凝土損傷和 鋼筋塑性應(yīng)變(SAUSAGE)Fig.20 Concrete damage and plastic strain of base story by SAUSAGE

核心筒20～22層位于第1次收進(jìn)位置，如圖21所示，為提高結(jié)構(gòu)抗震性能，20層的3片剪力墻設(shè)置為鋼板剪力墻。筒體的受拉損傷主要分布在連梁和框架的端部墻肢上；受壓損傷區(qū)域很小，集中在鋼板剪力墻的下邊緣；受剪損傷主要分布在連梁、鋼板剪力墻及周圍墻肢；鋼板剪力墻周圍墻肢的部分鋼筋出現(xiàn)屈服，但鋼板未屈服。

圖21 第20～22層混凝土損傷和鋼材(鋼筋)塑性應(yīng)變(SAUSAGE)Fig.21 Concrete damage and steel(rebar) plastic strain of 20～22 story by SAUSAGE

6 結(jié) 論

(1) 為更直接、準(zhǔn)確地反映混凝土損傷和破壞模式，本文根據(jù)混凝土塑性損傷模型中損傷發(fā)展時主應(yīng)力象限的差異，將損傷因子細(xì)化為拉彎型損傷因子、拉剪型損傷因子、壓彎型損傷因子和壓剪型損傷因子；結(jié)合工程實踐需求又將其分為受拉、受壓和受剪損傷。

(2) 基于可損傷分類的塑性損傷模型，將單點積分與本文給出的非線性沙漏力、面內(nèi)旋轉(zhuǎn)力的構(gòu)造方法結(jié)合，給出了含面內(nèi)旋轉(zhuǎn)自由度的4節(jié)點24自由度四邊形非線性顯式分層殼單元。本文殼元可較合理地反映剪力墻在面內(nèi)彎曲、面內(nèi)剪切和面外彎曲的耦合作用下的非線性發(fā)展。將本文殼元在SAUSAGE中完成基于CPU+GPU異構(gòu)并行計算的開發(fā)；通過與ABAQUS benchmark模型和單片剪力墻單調(diào)加載試驗對比，驗證了本文單元的正確性和合理性。

(3) 分別利用SAUSAGE和ABAQUS完成某框架核心筒結(jié)構(gòu)罕遇地震非線性動力時程分析，結(jié)果表明：SAUSAGE和ABAQUS的計算結(jié)果基本一致，SAUSAGE計算速度為ABAQUS計算速度的5.89倍；驗證了本文殼元大規(guī)模非線性動力求解的精度和計算效率。由SAUSAGE計算結(jié)果可知：①核心筒墻肢處于彎曲型受力狀態(tài)，未出現(xiàn)明顯損壞；②連梁作為耗能構(gòu)件，出現(xiàn)一定程度的剪切損壞，連梁端部出現(xiàn)受拉損傷；③核心筒收進(jìn)樓層的鋼板剪力墻及周圍墻肢的混凝土出現(xiàn)一定程度的剪切損壞，但鋼板未屈服。采用本文殼元可更直接地反映了剪力墻損壞模式，有助于了解和優(yōu)化結(jié)構(gòu)抗震性能。

致謝

感謝徐培福研究員對本文工作的指導(dǎo)。