結構動力彈塑性與倒塌分析(I)*——滯回曲線改進、ABAQUS子程序開發(fā)與驗證

柳國環(huán)，練繼建，國 巍

(1.天津大學建筑工程學院，天津300072;2.天津大學水利工程仿真與安全國家重點實驗室，天津300072;3.中南大學土木工程學院，湖南長沙410075)

0 前言

當今，隨著經(jīng)濟的快速發(fā)展和社會的不斷進步，大跨和超高層等超限結構日益增多。根據(jù)當前的規(guī)范要求，超限結構的地震反應分析需要作為一項專門的工作來進行，如結構在大震作用下的動力彈塑性表現(xiàn)行為(葉獻國等，2003;宣鋼等，2003;杜修力等，2007)，然后對計算結果是否符合相關規(guī)范和(或)規(guī)程的相關具體要求進一步審查。一般情況下，對于這種特殊結構的計算通常需要用到非線性計算能力相對較高級的大型分析軟件，例如:ABAQUS、MSC.MARC和LSDyna。柳國環(huán)等(2014a)針對ABAQUS大型程序，做了滯回規(guī)則改進、相關程序開發(fā)以及試驗結果對比，并通過帶有大跨—輸電實際工程計算進一步檢驗了相關程序的可行性、計算精度與現(xiàn)實可靠性。

就ABAQUS大型軟件而言，核心程序涉及到隱式和顯式兩種算法，隱式算法需要迭代計算以滿足相應的容差，而顯式算法無需迭代，但一般需要將步長設置為比較小(與結構體系中單元尺寸有關)以得到可靠可信的計算結果。在ABAQUS中，可以編制子程序并與主程序對接來實現(xiàn)結構的鋼筋混凝土動力彈塑性分析，相應隱式和顯式算法的接口子程序頭文件分別為UMAT和VUMAT，采用UMAT和VUMAT中的已有變量參數(shù)直接與程序中需要額外定義變量共同實現(xiàn)相應本構子程序的代碼。為此，本文首先根據(jù)已有資料詳細回顧了鋼材和混凝土材料的單軸本構模型，并改進了混凝土卸載后再加載曲線的規(guī)則，旨在使滯回曲線避免理論上存在加載曲線與卸載曲線在拐點前相交的可能性，并使得再加載曲線形狀更接近于試驗曲線形狀(先凹后凸的曲線狀)。然后，給出了ABAQUS中UMAT與VUMAT之間的區(qū)別與聯(lián)系，并按照UMAT與VUMAT的規(guī)則開發(fā)了直接能夠與ABAQUS主程序保持數(shù)據(jù)相互調(diào)用的顯式和隱式算法子程序ABA-Subroutine。最后，將編制的子程序分別用于鋼筋束、素混凝土、鋼筋混凝土、型鋼混凝土以及鋼管混凝土構件的彈塑性分析中，并將計算結果與試驗結果相對比，進一步考察和驗證程序的可行性和準確性。本文關于ABAQUS子程序開發(fā)工作具有一定基礎性，與柳國環(huán)等(2014b)的研究成果共同作為系列性工作，可為后續(xù)工作中計算結果的可靠性提供可行可信的前提條件。

1 材料本構關系

1.1 常用鋼材滯回規(guī)則的3種情形

本節(jié)重點詳述常用鋼筋滯回曲線3種情形的加載卸載規(guī)則。骨架曲線與滯回曲線如圖1所示(GB 50010－2010)。對于如圖1b所示的滯回曲線的滯回規(guī)則如式(1)所示，相關參數(shù)意義如圖1b中標注。

圖1 鋼材本構的骨架(a)和滯回(b)曲線示意圖Fig.1 Envelop(a)and hysteristic(b)curves of steel constitutive

加卸載規(guī)則分3種情形說明:

(1)對于所有卸載:均按照彈性模量斜率計算。

(2)對于穿越平衡位置后的反向加載:εa和σa表示卸載后再加載路徑起點的應變和應力，一般取為滯回曲線與橫坐標的交點，這時σa=0。

εb和σb表示卸載后再加載路徑終點的應變和應力(初始值取為初始屈服點對應量值)，如果再加載方向的鋼筋在加載歷史進程中尚未發(fā)生過屈服，則εb和σb取鋼筋初始屈服點的應變εb和應力fy，如果再加載方向的鋼筋曾經(jīng)有過屈服，則取該加載方向鋼筋歷史進程中的最大應變和應力。

(3)對于在平衡位置一側的卸載后再加載:

式中，σ和ε表示滯回曲線函數(shù)的因變量(應力)和自變量(應變)。k表示等效硬化斜率，一般取為峰值強度和屈服強度兩點所確定直線的斜率tanβ(圖1a)。上述處理方法旨在通過確定式(1)～(3)中的參數(shù)來最終確定此函數(shù)的表達式，以描述接近于實際鋼材的反復加載下的物理行為。

1.2 反復受壓卸載規(guī)則回顧與再加載規(guī)則改進

混凝土材料骨架曲線在低周反復荷載作用下的表現(xiàn)比較穩(wěn)定，且與在單調(diào)荷載作用下的表現(xiàn)有很好的一致性。但是，由于混凝土材料自身的特性與組成相對復雜等原因，在低周反復荷載作用下的滯回曲線的表現(xiàn)與鋼材相比更難以被準確量化。例如，由于加載過程中混凝土中局部發(fā)生裂紋等損傷現(xiàn)象，而且損傷會逐漸積累，會直接導致混凝土隨后表現(xiàn)出來的剛度逐漸減小，即剛度退化。目前描述混凝土滯回曲線的規(guī)則頗為多見，限于一次加載卸載過程試驗不能夠考慮到多次反復荷載作用下混凝土的表現(xiàn)行為，然而這種反復加載卸載作用下表現(xiàn)行為更符合隨機荷載作用下(例如:地震)實際情況。圖2b為混凝土反復荷載作用下滯回曲線示意圖，圖中σ與ε分別表示應力和應變，下標un，s和un，e表示卸載始點和終點，下標l，m表示加載路徑中的拐點，i和j表示第i次由骨架曲線卸載后第j次加載或卸載次數(shù)。

圖2 混凝土單軸本構骨架(a)和滯回(b)曲線示意圖Fig.2 Envelop(a)and hysteristic(b)curves of concrete umiaxial constitutive

(1)往復受壓混凝土的卸載規(guī)則

為能夠體現(xiàn)殘余應變隨加卸載次數(shù)增加而逐漸累計增大以及卸載剛度的逐漸退化(GB 50010－2010;聶建國，陶慕軒，2011)，按照如下計算:

式中，εc0為混凝土的峰值壓應變，根據(jù)(GB 50010－2010)取值。

本文未依據(jù)混凝土結構設計規(guī)范(GB 50010－2010)(GB 50010－2010)中所述直線卸載規(guī)則，而按照式(6)確定的更新割線剛度Eir后直接計算當前的應力。根據(jù)聶建國和陶慕軒(2011)通過卸載初始點應力、初始點應變、殘余應變以及當前步的應變ε共同來描述卸載曲線，如式(7)所示。

(2)往復受壓再加載規(guī)則改進

與卸載路徑不同，考慮到卸載后再加載路徑曲線存在拐點(再加載過程中曲線先表現(xiàn)為凹形曲線而后形狀發(fā)生變化的點)。本文同樣分段考慮加載規(guī)則，但為了更好描述再加載曲線路徑，這里定義為:第一段曲線凹(非二次曲線)，第二段曲線凸(非直線)，并給出拐點位置的應變和應力的規(guī)則如下:

式(12)中，κ=0.9i定義為卸載后再加載段拐點之后剛度退化參數(shù)，與由骨架曲線上卸載的次數(shù)i有關，該值隨i增大而減小。

需強調(diào)的是，編制程序時需要調(diào)用ABAQUS主程序中參數(shù)，需注意與上述表達式中物理參數(shù)符號統(tǒng)一。

2 子程序開發(fā)與實現(xiàn)

ABAQUS主程序具有用于二次開發(fā)材料接口的子程序，對應于隱式與顯式算法的頭文件為UMAT和VUMAT。用戶可利用相應接口編制子程序，并按照本構關系滯回規(guī)則定義變量，通過調(diào)用主程序材料參數(shù)并與主程序狀態(tài)變量實現(xiàn)無縫對接。關于子程序與ABAQUS主程序鏈接的基本框架如圖3所示。

如圖3a所示，對于采用UMAT的隱式算法需要更新應力、應變和相關的狀態(tài)變量，還要更新本構的雅克比矩陣(應力增量對應變增量的導數(shù)(dΔσ/dΔε)，然后通過當前更新的雅克比矩陣、應變矩陣以及幾何屬性可得到單元剛度矩陣的增量，可得到當前步局部坐標系下的單元剛度矩陣，進而結合單元的局部與整體坐標系的關系形成結構的整體剛度矩陣，隱式運算需要迭代以滿足相應的容差。

如圖3b所示，對于采用VUMAT的顯式算法子程序內(nèi)部直接更新的物理量是應力和相關狀態(tài)變量，因此不需要本構的雅克比矩陣。顯式算法不需要迭代，但需要事先通過ABAQUS設置的計算步長應足夠小以得到可靠可信的計算結果。這里需要強調(diào)，適為降低顯式算法的VUMAT子程序可能的數(shù)值誤差，相應變量宜被定義為雙精度，同時使用real*8和include vaba_param.inc定義，若只采用include vaba_param.inc語句不能確保變量被承認為雙精度，而與real*8共同使用具有雙保險意義。此外Job中顯式精度(ABAQUS/Explicit Precision)選為 double。對于 UMAT或 VUMAT，考慮應變率因素對本構動態(tài)強度提高影響，均需要通過應變增量DSTRAN與DTIME計算應變率絕對值DSTRAN/DTIME，然后根據(jù)本構動態(tài)提高系數(shù)公式進行計算。

圖3 ABAQUS主程序與子程序UMAT(a)與VUMAT(b)鏈接基本流程圖Fig.3 The basic flowchart of main program ABAQUS link to subroutine program UMAT(a)and VUAMT(b)

3 試驗驗證——鋼筋和素混凝土

為檢驗上述開發(fā)程序的可靠性，首先需要分別對鋼筋與混凝土代碼的準確性進行驗證。

驗證鋼筋計算結果:如圖4a所示的鋼筋位移加載制度根據(jù)應變—位移加載曲線以及試驗樣本標距給出(李敏，2012;李敏，李宏男，2010)，用于ABAQUS位移方式加載。分析時，ABAQUS中的隱式算法為采用靜力分析，顯式算法采用加大阻尼的動力計算方法(顯式算法只有動力分析，下同)以得到靜力計算的效果，與李敏(2012)給出的進入屈服階段后三周的滯回環(huán)結果相對應，相應給出了后3周的數(shù)值結果，兩種算法的結果與試驗結果對比曲線如圖4b所示。本文中鋼筋與其他算例相比尺寸較小，故顯式算法最大步長較取為1.0E－7。

圖4 鋼筋束加載制度與試驗結果對比(a)對應于試驗應變加載的等幅值位移加載曲線;(b)數(shù)值計算與試驗結果對比Fig.4 Loading system of steel bar HRB400 and test results comparison(a)amplitude displacement loading curve corresponding to the test strain loading;(b)comparison betweeen numerical calculation and test results

驗證素混凝土代碼計算結果:以編號BII－6試件為例，為100 mm×100 mm×300 mm棱柱體(過鎮(zhèn)海，張秀琴，1981)。試驗應變增量控制為1.0E－3(每次位移加載的增量為3.0E－1 mm)，每次卸載至應力為零結束后再反向加載，在第9次(應變?yōu)?.0E－3)開始卸載至應力為零時結束。加卸載曲線依據(jù)1.2節(jié)所述規(guī)則。為考察混凝土加卸載規(guī)則以及代碼的計算精度，應力—應變曲線(骨架曲線)按照過鎮(zhèn)海和張秀琴(1981)取為

式中，峰值應力σ1取值為30.664 MPa，峰值應變ε1取值為1.486E－3;a和α為參數(shù)，據(jù)過鎮(zhèn)海和張秀琴(1981)表2取值分別為1.2和1.64，b=3－2a，c=a－2。

兩種算法的結果如圖5b所示，兩者結果相近，與圖5a所示的試驗結果對比具有較好的一致性。

4 試驗驗證——鋼管/型鋼/鋼混構件

圖6a～c依次給出了圓鋼管混凝土、型鋼混凝土以及鋼筋混凝土構件的示意圖，并給出了相應的邊界條件、幾何參數(shù)與材料參數(shù)(李敏，2012;李敏，李宏男，2010;過鎮(zhèn)海，張秀琴，1981;劉文淵等，2012，周緒紅等(2009)，加載制度如圖6d～f所示。

首先采用SAP2000軟件分別建立上述3種構件的有限元模型，然后采用柳國環(huán)等(2014a)中SAP2ABAQUS接口程序分別ABAQUS的inp文件，文件格式包括相應隱式和顯式格式。進而應用ABAQUS讀入inp文件調(diào)整模型，轉換前后的有限元模型與相應說明如圖7所示。分別采用隱式和顯式算法分析圖中算例，并將模擬結果分別與文獻中試驗結果相比較，比較結果見圖8。應該說明，圖7～8和9中模型試驗數(shù)據(jù)由提取數(shù)據(jù)的Getdata軟件獲取得到，顯式和隱式方法計算最大步長設置為1.0E－6和0.02。

圖6 模型示意圖與加載歷程(a)圓鋼管;(b)型鋼;(c)鋼筋;(d)圓鋼管混凝土構件加載曲線;(e)型鋼混凝土構件加載曲線圖;(f)鋼筋混凝土構件加載曲線Fig.6 Model diagrammatic sketch and its corresponding loading history(a)circular steel tube;(b)section reinforced;(c)reinforced;(d)loading curve of concrete filled steel tube column;(e)loading curve of SRC column;(f)loading cauve of RC beam

圖7 SAP模型(a)與采用SAP2ABAQUS轉化后的ABAQUS模型(b)Fig.7 Comparison between SAP model(a)and ABAQUS model obtained by using SAP2ABAQUS(b)對比

從圖8a、9a和10a中可以看出:采用顯式與隱式方法的計算結果比較一致，說明在合理選擇步長情況下，不同算法對計算結果的影響較小。從圖8b、9b和10b中，可以看出:數(shù)值計算與試驗獲取數(shù)據(jù)結果接近程度令人滿意。需要說明的是，考慮到李敏和李宏男(2010)研究中型鋼混凝土給出了位移加載曲線周期數(shù)(6周)和重復循環(huán)次數(shù)(每2周采用同一等幅加載)，因此在提取其試驗數(shù)據(jù)時也只提取至位移加載的第6周(所顯示的第3個滯回環(huán))加載結束。

圖8 圓鋼管混凝土構件計算結果(a)顯式與隱式結果對比;(b)顯式、隱式與試驗結果對比Fig.8 Calculation results for concrete filled circular steel-tube column(a)comparison between explicit and implicit results;(b)comparison among explicit，inplicit and test results

圖9 型鋼混凝土構件計算結果(a)顯式與隱式結果對比;(b)顯式、隱式與試驗結果三者對比Fig.9 Calculation results for SRC column(a)comparison between explicit and implicit results;(b)comparison among explicit，inplicit and test results

圖10 鋼筋混凝土構件分析結果對比(a)顯式與隱式結果對比;(b)顯式、隱式與試驗結果三者對比Fig.10 Comparison of analysis results for RC beam(a)comparison between explicit and implicit results;(b)comparison among explicit，inplicit and test results

5 結語

本文作為系列研究的第I部分工作，主要涉及本構與子程序開發(fā)等方面內(nèi)容，簡要總結如下:

(1)改進了混凝土卸載再加載曲線規(guī)則并給出具體表達式，定義了再加載段拐點之后剛度退化參數(shù)。

(2)給出并重點論述了UMAT和VUMAT與ABAQUS鏈接的框架示意圖，并給出相關注記。

(3)開發(fā)并實現(xiàn)了適用于隱式算法UMAT的子程序以及適用于顯式算法的雙精度VUMAT的子程序。

(4)驗證了本文改進的混凝土再加載規(guī)則與所開發(fā)程序的可行性和計算精度。首先考察了兩種算法(隱式和顯式法)對單一材料(鋼筋和素混凝土)構件的計算結果，進而考察了對混合材料(鋼筋混凝土、鋼管混凝土和型鋼混凝土)的計算結果，并與試驗結果相對比。

本文的理論基礎和本構改進、程序開發(fā)以及試驗驗證，可為實際工程計算提供可靠性保證。

杜修力，賈鵬，趙均.2007.鋼筋混凝土核心筒不同軸壓比作用下抗震性能試驗研究［J］.哈爾濱工業(yè)大學學報，39(2):567－572.

過鎮(zhèn)海，張秀琴.1981.混凝土在反復荷載作用下的應力—應變?nèi)€［J］.冶金建筑，9:14 －17.

李敏.2012材料的率相關性對鋼筋混凝土結構動力性能的影響［D］.大連:大連理工大學.

李敏，李宏男.2010.建筑鋼筋動態(tài)試驗及本構模型［J］.土木工程學報，43(4):70－75.

劉文淵，冷杰，段文峰.2012.往復荷載下圓鋼管混凝土柱的數(shù)值模擬［J］.吉林建筑工程學院學報，29(1):1－4.

柳國環(huán)，練繼建，國巍.2014a.結構動力彈塑性和倒塌分析(II):SAP2ABAQUS接口技術、程序開發(fā)與驗證［J］.地震研究，37(1):132－140.

柳國環(huán)，練繼建，孫雪艷，等.2014b.結構動力彈塑性和倒塌分析(III):地震差動作用下輸電塔—線體系的彈塑性與倒塌分析［J］.地震研究，37(1):141 －150.

聶建國，陶慕軒.2011.采用纖維梁單元分析鋼—混凝土組合結構地震反應的原理［J］.建筑結構學報，32(10):1－10.

宣鋼，顧祥林，呂西林.2003.強震作用下混凝土框架結構倒塌過程的數(shù)值分析［J］.地震工程與工程振動，23(6):24－30.

葉獻國，徐勤，李康寧，等.2003.地震中受損鋼筋混凝土建筑彈塑性時程分析與振動臺試驗研究［J］.土木工程學報，36(12):20－25.

周緒紅，張小東，劉屆鵬.2009.鋼管約束鋼筋混凝土柱與型鋼混凝土柱滯回性能試驗研究［J］.建筑結構學報，30(supp1):121－128.

GB 50010－2010，混凝土結構設計規(guī)范［S］.