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

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

0 前言

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

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

1 材料本構(gòu)關(guān)系

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

(2)往復(fù)受壓再加載規(guī)則改進(jìn)

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

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

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

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

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

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

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

圖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ā)程序的可靠性，首先需要分別對鋼筋與混凝土代碼的準(zhǔn)確性進(jìn)行驗證。

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

圖4 鋼筋束加載制度與試驗結(jié)果對比(a)對應(yīng)于試驗應(yīng)變加載的等幅值位移加載曲線;(b)數(shù)值計算與試驗結(jié)果對比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

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

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

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

4 試驗驗證——鋼管/型鋼/鋼混構(gòu)件

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

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

圖6 模型示意圖與加載歷程(a)圓鋼管;(b)型鋼;(c)鋼筋;(d)圓鋼管混凝土構(gòu)件加載曲線;(e)型鋼混凝土構(gòu)件加載曲線圖;(f)鋼筋混凝土構(gòu)件加載曲線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轉(zhuǎn)化后的ABAQUS模型(b)Fig.7 Comparison between SAP model(a)and ABAQUS model obtained by using SAP2ABAQUS(b)對比

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

圖8 圓鋼管混凝土構(gòu)件計算結(jié)果(a)顯式與隱式結(jié)果對比;(b)顯式、隱式與試驗結(jié)果對比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 型鋼混凝土構(gòu)件計算結(jié)果(a)顯式與隱式結(jié)果對比;(b)顯式、隱式與試驗結(jié)果三者對比Fig.9 Calculation results for SRC column(a)comparison between explicit and implicit results;(b)comparison among explicit，inplicit and test results

圖10 鋼筋混凝土構(gòu)件分析結(jié)果對比(a)顯式與隱式結(jié)果對比;(b)顯式、隱式與試驗結(jié)果三者對比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 結(jié)語

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

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

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

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

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

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

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

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

李敏.2012材料的率相關(guān)性對鋼筋混凝土結(jié)構(gòu)動力性能的影響［D］.大連:大連理工大學(xué).

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

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

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

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

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

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

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

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

GB 50010－2010，混凝土結(jié)構(gòu)設(shè)計規(guī)范［S］.