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

曹勝濤， 李志山， 楊志勇

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

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

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

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

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

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

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

(1)

(2)

(3)

根據損傷力學的基本假定，真實應力σij為

(4)

(5)

本文通過損傷發展時主應力所在象限的差異對損傷進行進一步分類。

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

(6)

(7)

(8)

(9)

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

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

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

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

塑性損傷本構模型應力更新算法采用Lee等[21]給出的應力回映算法。混凝土模型參數按《混凝土結構設計規范》(GB 50010—2010)附錄C.2[22](簡稱：《混規》)的標準值取值，塑性應變和損傷取值參考任曉丹等[23]的做法。

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

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

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

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

(a)三維空間

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

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

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

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

雙線性形函數為

(10)

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

圖7 等參元Fig.7 The isoparametric element

形函數對局部坐標系x,y偏導為

(11)

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

局部坐標系中殼元節點j對應的面內變形子矩陣塊、面外彎曲變形子矩陣塊和面外剪切變形子矩陣塊，如式(12)～式(14)所示。

(12)

(13)

(14)

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

(15)

(16)

(17)

(18)

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

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

(19)

(20)

本文殼元面內變形矩陣、面外彎曲變形矩陣和面外剪切變形矩陣如式(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)，得到單點積分處單元面內和面外彎曲常應變場矩陣B0m和B0sb。單元面外符合平截面假定，由面內和面外的耦合可得到殼元第i層積分點處應變增量

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

(24)

式中:Δde單元節點位移增量。

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

對Δσi進行積分得到此部分對殼元面內節點力增量的貢獻

(25)

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

(26)

通過截面積分得到單點積分對單元面外彎曲節點力增量的貢獻

(27)

(28)

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

(29)

(30)

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

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

(31)

(32)

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

綜上所述，全局坐標系中本文殼元單元節點力增量

(33)

3 本文殼元的開發實現及驗證

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

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

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

圖9 測試模型Fig.9 Test model

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

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

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

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

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

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

如圖14所示，SAUSAGE計算得到頂點水平位移和水平反力的關系與試驗結果基本一致。相對于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

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

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

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

5 框架核心筒結構動力損傷和破壞分析

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

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

圖15 某框架核心筒結構Fig.15 A frame-core tube structure

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

表4 SAUSAGE和ABAQUS模態分析結果對比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 罕遇地震非線性計算結果

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

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

圖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所示。

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

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

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

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

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

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

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

6 結 論

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

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

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

致謝

感謝徐培福研究員對本文工作的指導。