禚 一，李忠獻(xiàn)，李 寧

(1. 天津大學(xué)建筑工程學(xué)院，天津 300072；2. 天津大學(xué)濱海土木工程結(jié)構(gòu)與安全教育部重點(diǎn)實(shí)驗(yàn)室，天津 300072)

近年來，隨著地震作用下橋梁結(jié)構(gòu)倒塌過程模擬研究的深入，鋼筋混凝土橋墩等重要構(gòu)件的精細(xì)化模擬問題已越來越多地受到人們的重視.目前，構(gòu)件精細(xì)化模擬分析模型主要有三維實(shí)體有限元模型和離散桿系單元模型.三維實(shí)體有限元模型雖然可以相當(dāng)精細(xì)地模擬構(gòu)件的一些重要非線性特征，但較高的計(jì)算成本及計(jì)算收斂問題在很大程度上限制了這種模型的發(fā)展，使之難以用于復(fù)雜橋梁結(jié)構(gòu)的模擬.相比而言，離散桿系單元模型既可從宏觀上模擬構(gòu)件的性能又能深入分析構(gòu)件的局部非線性特性，且模型簡單、無需耗費(fèi)大量機(jī)時(shí)，因而受到大多數(shù)研究人員的青睞[1].根據(jù)單元塑性鉸分布方式和截面滯回特性模擬方法的不同，離散桿系單元分析模型又分為集中塑性模型[2-3]、分布塑性模型[4]和梁柱纖維單元模型[5-6].由于纖維模型直接從截面內(nèi)纖維的本構(gòu)關(guān)系出發(fā)得到單元乃至整個(gè)構(gòu)件的非線性性能，可有效地模擬構(gòu)件的剛度退化、強(qiáng)度退化等損傷效應(yīng)及軸力和(單向、雙向)彎矩的多維耦合效應(yīng)等復(fù)雜非線性行為，因而已成為結(jié)構(gòu)精細(xì)化模擬的必要手段，并得到廣泛應(yīng)用．

現(xiàn)階段，采用纖維梁柱單元進(jìn)行數(shù)值模擬，大多利用國外已有軟件，如 OpenSEES[7]、DRAIN[8]等.大型通用有限元軟件ABAQUS具備強(qiáng)大的非線性求解能力以及友好的前后處理界面，已為大多數(shù)研究人員所采用，但國內(nèi)外基于 ABAQUS的纖維梁柱單元的實(shí)用開發(fā)卻不多見.因此，課題組基于ABAQUS開發(fā)了一種實(shí)用的鋼筋混凝土纖維梁柱單元模擬分析平臺(tái) FENAP，編制相應(yīng)的材料庫，開發(fā)了混凝土和鋼材的本構(gòu)模型，并利用 FENAP平臺(tái)模擬了鋼筋混凝土梁、墩柱構(gòu)件的滯回性能，取得了較好的模擬結(jié)果，驗(yàn)證了FENAP平臺(tái)進(jìn)行非線性靜力分析的有效性.

筆者在以往研究成果的基礎(chǔ)上，進(jìn)一步研究了纖維梁柱單元與通用有限元軟件ABAQUS相結(jié)合進(jìn)行非線性動(dòng)力分析的基本原理，開發(fā)了 FENAP平臺(tái)的非線性動(dòng)力分析模塊.同時(shí)，利用 FENAP平臺(tái)分別模擬了橋墩單向和雙向地震動(dòng)作用下的動(dòng)力特性，并將計(jì)算結(jié)果與OpenSEES及試驗(yàn)結(jié)果進(jìn)行對(duì)比分析，以驗(yàn)證FENAP平臺(tái)進(jìn)行非線性動(dòng)力分析的有效性.

圖1 纖維梁柱單元端力、位移及截面力、變形的定義Fig.1 Definition of terminal force， displacement，cross-section force and distoration of fiber beam-column element

1 FENAP平臺(tái)非線性動(dòng)力分析的基本原理

在 FENAP平臺(tái)中，采用了基于剛度法的三維兩節(jié)點(diǎn)纖維梁柱單元，如圖1所示.

根據(jù)經(jīng)典的 Euler-Bernoulli梁柱單元理論，單元的位移場分布 ()xu~ 可用式(1)所示的單元插值函數(shù)進(jìn)行描述.

式中：U為單元節(jié)點(diǎn)位移矢量；s()xN 為插值函數(shù)矩陣；軸向位移 ()u x采用拉格朗日線性插值；橫向位移()v x、()w x分別采用三次埃米特多項(xiàng)式插值.在平截面假定條件下，單元截面變形 ()xd 與節(jié)點(diǎn)位移U的變形協(xié)調(diào)關(guān)系可表示為

式中：iε為第i根纖維的應(yīng)變，n表示纖維總數(shù)；iy、iz表示第i根纖維的中心位置坐標(biāo).L為線性幾何變換矩陣，即

根據(jù)應(yīng)變所在位置處纖維的材料本構(gòu)關(guān)系便可計(jì)算得到相應(yīng)纖維的材料切線模量iE及應(yīng)力值寫成矩陣形式為

式中diag()表示對(duì)角陣．對(duì)截面內(nèi)所有纖維進(jìn)行積分，可得到截面剛度矩陣s()xk 和抗力矢量 ()xD ，分別為

將式(2)代入式(9)，兩端消去虛加位移δU后，得到單元抗力矢量為

將式(2)和式(8)代入式(10)，得到單元?jiǎng)偠染仃嚍?/p>

FENAP平臺(tái)是基于 ABAQUS/standard求解器開發(fā)的，在求解器中動(dòng)力時(shí)程分析時(shí)，采用了Hilbert-Hughes-Taylor(HHT)積分算法進(jìn)行迭代計(jì)算[9].HHT算法通過建立結(jié)構(gòu)整體等效剛度矩陣和等效荷載矢量來計(jì)算結(jié)構(gòu)的節(jié)點(diǎn)位移，最終得到結(jié)構(gòu)的響應(yīng).因此，在單元層次上，需要進(jìn)一步建立單元的等效剛度矩陣eqK 和等效抗力矢量epP.根據(jù) HHT算法對(duì)等效剛度矩陣eqK 和等效抗力矢量eqP的定義，可得到

式中：M為單元質(zhì)量矩陣，這里選用集中質(zhì)量矩陣形式；C為阻尼矩陣，選用 Rayleigh阻尼形式；表示單元內(nèi)力貢獻(xiàn)；α、β、γ均為分析參數(shù)，1/2 γα=?，為了加快收斂，選用

2 FENAP平臺(tái)材料庫本構(gòu)模型

FENAP平臺(tái)的材料庫中包含的單軸本構(gòu)模型包括適用于混凝土材料的 Mohd-Yassin[10]提出的混凝土損傷本構(gòu)模型和適用于鋼筋材料的理想彈塑性本構(gòu)模型、雙線性等向強(qiáng)化本構(gòu)模型和Filippou等[11]修正的 Menegotto-Pinto本構(gòu)模型[12]等.其中，修正的Menegotto-Pinto本構(gòu)模型為課題組近期開發(fā)的，可考慮鋼筋反復(fù)加載過程中的Bauschinger效應(yīng)和等向強(qiáng)化效應(yīng)，并且在計(jì)算過程中，采用應(yīng)變的顯函數(shù)表達(dá)形式求解應(yīng)力和切線模量，因而具有較高的計(jì)算精度，與鋼筋反復(fù)加載試驗(yàn)結(jié)果吻合較好，其循環(huán)加載滯回曲線如圖2所示.

圖2 Menegotto-Pinto鋼筋本構(gòu)模型循環(huán)加載滯回曲線Fig.2 Cyclic hysteretic curves of Menegotto-Pinto steel Fig.2 constitutive model

3 模擬分析

為驗(yàn)證應(yīng)用FENAP平臺(tái)進(jìn)行非線性動(dòng)力分析的有效性，將模擬分析鋼筋混凝土橋墩在單向地震激勵(lì)下的數(shù)值試驗(yàn)和鋼筋混凝土橋墩在雙向地震激勵(lì)下的振動(dòng)臺(tái)試驗(yàn).

3.1 鋼筋混凝土橋墩在單向地震激勵(lì)下的數(shù)值試驗(yàn)

分別利用 FENAP平臺(tái)和 OpenSEES軟件對(duì)一RC橋墩在單向加載下的地震響應(yīng)進(jìn)行了數(shù)值模擬.該橋墩高 2.265,m，矩形截面 400,mm×800,mm.在FENAP建模過程中，將橋墩劃分為10個(gè)纖維梁柱單元，共計(jì)11個(gè)節(jié)點(diǎn)，每個(gè)單元采用兩個(gè)Gauss積分截面，截面內(nèi)由 24根保護(hù)層混凝土纖維、32根核心混凝土纖維以及24根鋼筋纖維組成.其中，混凝土纖維均采用 FENAP材料庫中的混凝土損傷本構(gòu)模型，其材料參數(shù)包括混凝土的初始彈性模量c0E 、峰值壓應(yīng)力cf′、相應(yīng)的壓應(yīng)變0ε、極限壓應(yīng)變uε、峰值拉應(yīng)力和抗拉剛化模量tsE .箍筋對(duì)核心混凝土的約束作用通過提高核心混凝土的材料本構(gòu)參數(shù)來實(shí)現(xiàn).表 1給出了保護(hù)層和核心混凝土的材料特性參數(shù).鋼筋纖維采用 FENAP材料庫中的雙線性等向強(qiáng)化本構(gòu)模型，其材料特性參數(shù)包括初始彈性模量 Es、屈服強(qiáng)度fy和屈服后剛度系數(shù)b，如表 2所示.圖 3給出了單元?jiǎng)澐旨敖孛胬w維離散化方式.OpenSEES建模過程中，單元及截面纖維的劃分方式與 FENAP一致，每個(gè)單元設(shè)為5個(gè)積分截面，并選用Concrete01模型模擬混凝土纖維的本構(gòu)，Steel01模型模擬縱筋纖維的本構(gòu)，相應(yīng)材料的強(qiáng)度和應(yīng)變參數(shù)也與FENAP相同.

表1 混凝土材料特性Tab.1 Material properties of concrete

表2 鋼筋材料特性Tab.2 Material properties of steel

圖3 單元?jiǎng)澐旨敖孛胬w維離散化Fig.3 Element division and cross-section fiber discretization

分析過程中，選用 El-centro波南北向在墩底進(jìn)行單向加載，為了得到結(jié)構(gòu)的非線性響應(yīng)，將地震動(dòng)峰值調(diào)整為0.312g，圖4、圖5分別給出了FENAP平臺(tái)和 OpenSEES計(jì)算得到的墩頂位移時(shí)程和墩底截面彎矩-曲率關(guān)系對(duì)比曲線，比較結(jié)果可看出，F(xiàn)ENAP平臺(tái)得到的位移時(shí)程響應(yīng)的相位變化和幅值大小與OpenSEES的結(jié)果十分吻合．FENAP計(jì)算得到的墩頂位移最大值為 34.1,mm，發(fā)生在 6.268,s處，而OpenSEES計(jì)算得到的墩頂位移最大值為33.47,mm，發(fā)生在6.26,s處．另外，即使進(jìn)入非線性狀態(tài)后，截面彎矩-曲率關(guān)系曲線的形狀和包圍面積的耗能，二者也給出了十分接近的分析結(jié)果．進(jìn)一步研究發(fā)現(xiàn)，對(duì)于該算例OpenSEES計(jì)算用時(shí)為908,s，而FENAP僅用時(shí) 746,s就得到了收斂的計(jì)算結(jié)果，這更加驗(yàn)證了FENAP平臺(tái)在保證求解精度的同時(shí)還具備較高的計(jì)算效率.

圖4 墩頂相對(duì)位移時(shí)程對(duì)比曲線Fig.4 Comparison of relative displacement time history at top of bridge pier

圖5 墩底截面彎矩-曲率對(duì)比曲線Fig.7 Comparison of moment-curvature curves at bottom Fig.7 of bridge pier

3.2 橋墩在雙向地震激勵(lì)下的振動(dòng)臺(tái)試驗(yàn)

文獻(xiàn)[13]對(duì)一方形截面鋼筋混凝土橋墩進(jìn)行了雙向加載振動(dòng)臺(tái)試驗(yàn)，筆者采用 FENAP平臺(tái)對(duì)其進(jìn)行數(shù)值模擬分析.建模過程中，沿橋墩縱向劃分為10個(gè)纖維梁柱單元和2個(gè)剛體單元，并在墩頂配重的形心位置布置一個(gè)質(zhì)量點(diǎn)單元，用于模擬上部結(jié)構(gòu)的等效質(zhì)量，如圖 6所示.在纖維單元中，將截面劃分為 144根核心混凝土纖維、112根保護(hù)層混凝土纖維和 48根鋼筋纖維，混凝土纖維采用 FENAP材料庫中的混凝土損傷本構(gòu)模型，強(qiáng)度為 34.1,MPa；鋼筋纖維采用FENAP材料庫中的修正Menegotto-Pinto本構(gòu)模型，強(qiáng)度為384,MPa，其他材料參數(shù)按照文獻(xiàn)[13]選取.加載時(shí)，選用Kobe波的東西和南北向同時(shí)在墩底輸入，并按時(shí)間相似比0.5對(duì)原記錄進(jìn)行壓縮，地震波東西和南北向峰值加速度分別取為0.666g和0.642g.

圖7、圖8分別給出了FENAP平臺(tái)和文獻(xiàn)[13]中作者數(shù)值模擬結(jié)果及試驗(yàn)結(jié)果的對(duì)比曲線．與文獻(xiàn)[13]的模擬結(jié)果比較，F(xiàn)ENAP平臺(tái)模擬結(jié)果與試驗(yàn)曲線較為接近.在3,s之前，F(xiàn)ENAP平臺(tái)與試驗(yàn)曲線基本吻合.對(duì)于x和y方向的峰值響應(yīng)出現(xiàn)時(shí)間和大小兩者并無較大差別.因此，對(duì)于結(jié)構(gòu)在未進(jìn)入塑性狀態(tài)前，F(xiàn)ENAP平臺(tái)可以給出較為準(zhǔn)確的預(yù)測(cè)結(jié)果；經(jīng)過峰值點(diǎn)，構(gòu)件完全進(jìn)入到塑性狀態(tài)后，F(xiàn)ENAP平臺(tái)結(jié)果與試驗(yàn)結(jié)果偏差較大，試驗(yàn)得到的位移明顯大于分析結(jié)果，其主要原因在于進(jìn)入塑性階段后，鋼筋和混凝土間發(fā)生的黏結(jié)滑移會(huì)降低構(gòu)件的剛度，增大位移響應(yīng)，而 FENAP中所采用的纖維模型基本原理是在平截面假定基礎(chǔ)上建立的，并未考慮黏結(jié)滑移的影響，因而得到的位移響應(yīng)較試驗(yàn)結(jié)果偏小.另一方面，由雙向位移耦合曲線的對(duì)比結(jié)果可看出，與文獻(xiàn)[13]得到的雙向位移耦合曲線相比，F(xiàn)ENAP平臺(tái)得到的雙向加載下的位移發(fā)展方向和曲線的形狀與試驗(yàn)是較為吻合的.

圖6 橋墩單元?jiǎng)澐旨敖孛胬w維離散化Fig.6 Element division and cross-section fiber discretization of bridge pier

圖7 墩頂相對(duì)位移時(shí)程響應(yīng)曲線Fig.7 Time history of relative displacement at top of bridge pier

圖8 雙向位移耦合曲線Fig.8 Coupled curves of bilateral displacement

為進(jìn)一步比較在FENAP平臺(tái)中考慮雙向彎曲耦合作用對(duì)動(dòng)力結(jié)果的影響，將上述橋墩試件所采用的東西向和南北向地震波，分別在試件的兩個(gè)垂直方向單獨(dú)施加，以模擬單向地震動(dòng)加載情況，從而對(duì)橋墩在單向加載與雙向加載條件下的數(shù)值模擬結(jié)果進(jìn)行對(duì)比研究.計(jì)算得到的單向和雙向耦合的相對(duì)位移和絕對(duì)加速度時(shí)程對(duì)比曲線如圖9和圖10所示.

圖9 單向和雙向加載位移時(shí)程響應(yīng)對(duì)比Fig.9 Comparison of time histories of displacement between unilateral and bilateral excitation

圖10 單向和雙向加載加速度時(shí)程響應(yīng)對(duì)比Fig.10 Comparison of time histories of acceleration between unilateral and bilateral excitation

結(jié)果顯示，雙向地震下得到的位移響應(yīng)結(jié)果大于單向(x方向加載位移峰值僅為雙向加載時(shí)的63.3%，y方向僅為雙向加載時(shí)的 82.2%)，而雙向地震下得到的加速度響應(yīng)小于單向(x方向單向加載加速度峰值為雙向加載時(shí)的1.354倍，y方向?yàn)殡p向加載時(shí)的1.216倍)，可見，F(xiàn)ENAP平臺(tái)可有效模擬雙向彎曲耦合作用對(duì)結(jié)構(gòu)動(dòng)力響應(yīng)的影響，反映其真實(shí)抗震性能.

4 結(jié) 語

基于 FENAP平臺(tái)所開發(fā)的非線性動(dòng)力分析模塊，充分發(fā)揮了 ABAQUS求解器的強(qiáng)大非線性動(dòng)力分析功能，可有效模擬鋼筋混凝土橋墩的復(fù)雜非線性動(dòng)力行為，可考慮雙向彎矩和軸力的耦合效應(yīng)，且具有較高計(jì)算效率和求解精度，從而為長大橋梁結(jié)構(gòu)地震災(zāi)變過程模擬提供了一種實(shí)用分析手段.

［1］ Taucer F F，Spacone E，F(xiàn)ilippou F C. A Fiber Beamcolumn Element for Seismic Response Analysis of Reinforced Concrete Structures[R].Berkeley：Earthquake Engineering Research Center，University of California，1991.

［2］ Powell G H，Chen P F S．3D beam-column element with generalized plastic hinges[J]. Journal of Engineering Mechanics，ASCE，1986，112(7)：627-641.

［3］ Sfakianakis M G，F(xiàn)ardis M N. Bounding surface model for cyclic biaxial bending of RC sections[J]. Journal of Engineering，Mechanics，ASCE，1991，117(12)：2748-2769.

［4］ Roufaiel M S L，Meyer C. Analytical modeling of hysteretic behavior of R/C frames[J]. Journal of Structural Engineering，ASCE，1987，113(3)：429-444.

［5］ Spacone E，F(xiàn)ilippou F C，Taucer F F. Fibre beam-column model for non-linear analysis of R/C frames(Part I)：Formulation[J]. Earthquake Engineering and Structural Dynamics，ASCE，1996，25(7)：711-725.

［6］ Spacone E，F(xiàn)ilippou F C，Taucer F F. Fibre beam-column model for non-linear analysis of R/C frames(Part Ⅱ)：Applications[J]. Earthquake Engineering and Structural Dynamics，1996，25(7)：727-742.

［7］ Mckenna F，F(xiàn)enves G L. The OpenSEES command language manual[EB/OL]. http：//opensees. berkeley. edu/wiki/index. php/Main_Page，2010-03-31.

［8］ Prakash V，Powell G H，Campbell S D. DRAIN-2DX：Static and Dynamic Analysis of Inelastic Plane Structures[R]. Berkeley：Earthquake Engineering Research Center，University of California，1993.

［9］ ABAQUS Inc. ABAQUS User Subroutines Reference Manual[M]. Rhode Island：ABAQUS Inc，2006.

［10］ Mohd-Yassin M H. Nonlinear Analysis of Prestressed Concrete Structures under Monotonic and Cyclic Loads[D]. Berkeley：School of Civil and Environmental Engineering，University of California，1994.

［11］ Filippou F C，Popov E P，Bertero V V. Effects of Bond Deterioration on Hysteretic Behavior of Reinforced Concrete Joints[R]. Bekerley：Earthquake Engineering Research Center，University of California，1983.

［12］ Menegotto M，Pinto P E，Slender R C. Compressed members in biaxial bending[J]. Journal of Structural Division，ASCE，1977，103(3)：587-605.

［13］ Nishida H，Unjoh S. Dynamic response characteristic of reinforced concrete column subjected to bilateral earthquake ground motions[C]// 13th World Conference on Earthquake Engineering．Vancouver，Canada，2004：576-587.