胡鄭州，吳明兒

（同濟(jì)大學(xué) 建筑工程系，上海200092）

根據(jù)式（3），得：

桿系結(jié)構(gòu)在極限荷載、強(qiáng)震等作用下，表現(xiàn)高度非線性，若按照線性理論分析，很難滿足工程設(shè)計(jì)要求.Zeris［1］和Spacone［2］提出的纖維梁單元，并成功解決鋼筋混凝土柱的材料非線性問題；文獻(xiàn)［3］基于柔度法并考慮了雙向彎曲與變軸力作用下對空間桿系結(jié)構(gòu)進(jìn)行了非彈性分析，其分析結(jié)果與試驗(yàn)結(jié)果相當(dāng)吻合；文獻(xiàn)［4］采用纖維梁單元分析了鋼－混凝土組合結(jié)構(gòu)在地震作用下的受力機(jī)理和破壞規(guī)律；文獻(xiàn)［5］采用纖維梁單元在火災(zāi)作用下模擬整體結(jié)構(gòu)連續(xù)性倒塌過程、機(jī)理及破壞規(guī)律；文獻(xiàn)［6］采用纖維梁單元分析了在考慮方鋼管局部屈曲情況下薄壁方鋼管混凝土梁柱結(jié)構(gòu)在單調(diào)和循環(huán)荷載下的性能；文獻(xiàn)［7－8］采用考慮幾何和材料非線性的纖維梁單元對橋墩在軸力和彎矩作用下進(jìn)行了二階分析，亦得出了滿意結(jié)果.可見纖維梁單元是分析桿系結(jié)構(gòu)高度非線性的最好單元之一.

本文首先根據(jù)連續(xù)性介質(zhì)力學(xué)和虛位移原理推導(dǎo)了通用的考慮大位移增量非線性有限元增量更新拉格朗日（Updated Lagrangian，UL）列式，同時(shí)修正大位移增量矩陣.據(jù)此理論推到了三維纖維梁大位移增量有限元UL列式.以上理論編制了分析程序，通過對幾個(gè)算例分析，證明該單元的精確性、通用性.

1 增量UL非線性列式

基于連續(xù)性介質(zhì)力學(xué)和虛位移原理［9］，參考現(xiàn)實(shí)時(shí)刻t構(gòu)形坐標(biāo)系定義得：

式中：kL表示小位移剛度矩陣；kg表示為幾何剛度矩陣；kU表示為大位移增量剛度矩陣；Δue表示節(jié)點(diǎn)增量位移；N表示單元局部坐標(biāo)插值函數(shù)矩陣；t＋Δttq表示單元節(jié)點(diǎn)參考t時(shí)刻構(gòu)形在t＋Δt時(shí)刻施加的荷載向量；t＋Δttρ表示單元參考t時(shí)刻構(gòu)形在t＋Δt時(shí)刻的密度；t＋Δttf表示單元參考t時(shí)刻構(gòu)形在t＋Δt時(shí)刻的體荷載向量；t＋Δt tt表示單元參考t時(shí)刻構(gòu)形在t＋Δt時(shí)刻的面荷載向量；ttσ表示t時(shí)刻的應(yīng)力向量；tdv表示t時(shí)刻單元積分體積；tda表示t時(shí)刻單元積分面積；tV、tA分別表示為t時(shí)刻單元的體積、面積.在基于上一時(shí)刻構(gòu)形時(shí)下面公式推導(dǎo)過程省略左上角t＋Δt和左下角t，其中kL、kg、kU表達(dá)式為

式（2）—（4）中：BL表示線性應(yīng)變—位移矩陣；BN表示非線性應(yīng)變—位移矩陣；DT表示材料本構(gòu)關(guān)系矩陣；G表示對形函數(shù)偏導(dǎo)矩陣；P表示應(yīng)力矩陣.式（4）仍然是不對稱矩陣，進(jìn)行下面簡化，

于是可使kU矩陣對稱化，從而使整個(gè)剛度矩陣是對稱的，這對單元切線剛度矩陣元素在計(jì)算機(jī)內(nèi)存中存儲(chǔ)是有利的.Bathe等［9］在推導(dǎo)UL 時(shí)將kU項(xiàng)忽略掉，本文保留該項(xiàng)對高度非線性的影響，利用在增量迭代步上一步迭代步的節(jié)點(diǎn)增量位移收斂解去求解大位移增量矩陣.

2 三維纖維梁單元?jiǎng)偠染仃囃茖?dǎo)

2.1 小位移剛度矩陣

三維纖維梁單元的位移函數(shù)采用Euler－Bernoulli梁的位移函數(shù)，局部坐標(biāo)如圖1所示，截面上任意一點(diǎn)的軸向應(yīng)變?yōu)?/p>

式中：ε11為截面上任意一點(diǎn)的軸向應(yīng)變；ε0為截面形心軸向應(yīng)變；y為纖維在局部坐標(biāo)y－z上y軸坐標(biāo)；z為纖維在局部坐標(biāo)y－z上z軸坐標(biāo)；v″為局部坐標(biāo)系沿y方向位移兩階偏導(dǎo)；w″為局部坐標(biāo)系沿z方向位移兩階偏導(dǎo).

圖1 三維纖維梁單元Fig.1 Three dimensional fiber beam element

利用虛位移原理便可推導(dǎo)截面切線剛度矩陣ks［1－2］：

式中：Etan為纖維彈性模量；A為截面面積，dA截面積分面積變量，其中L表達(dá)式如下：

得出ks表達(dá)式如下：

式中：n為纖維截面中纖維總數(shù)；Ei，tan表示第i纖維彈性模量；Ai表示第i個(gè)纖維面積；yi，zi定義見圖1所示.

根據(jù)式（2）得到單元的小位移剛度矩陣：

2.2 幾何剛度矩陣和大位移增量剛度矩陣

根據(jù)式（3），得：

其中，P，G 表達(dá)式如下：

其中：

其中，σ11表示截面形心出軸應(yīng)力，x表示單元纖維截面形心到單元端部距離.

根據(jù)式（5）大位移增量剛度矩陣：

其中：由式（9）—（11）得到單元的局部坐標(biāo)下參考t時(shí)刻構(gòu)形在t＋Δt時(shí)刻構(gòu)形的切線剛度矩陣，并通過單元坐標(biāo)轉(zhuǎn)換矩陣，將局部坐標(biāo)系的切線剛度矩陣轉(zhuǎn)化為整體坐標(biāo)系下，最后得到U.L.列式的大位移增量方程如下：

3 混凝土和鋼筋材料本構(gòu)

3.1 混凝土本構(gòu)關(guān)系

本文混凝土本構(gòu)模型采用修正的Kent－Park模型［12］，該模型通過修改素混凝土（無約束混凝土）受壓骨架曲線應(yīng)變軟化段斜率來考慮箍筋對混凝土的約束影響，其骨架曲線如圖2所示，σc，εc為是混凝土的應(yīng)力和應(yīng)變；ε0為混凝土強(qiáng)度峰值應(yīng)變；K為考慮箍筋對混凝土約束作用所引起混凝土強(qiáng)度增大系數(shù)；εC20為約束混凝土強(qiáng)度等級(jí)為C20對應(yīng)的極限應(yīng)變；ε′C20為素混凝土強(qiáng)度等級(jí)為C20 對應(yīng)的極限應(yīng)變；fc為混凝土軸心抗壓強(qiáng)度，MPa.分段函數(shù)加以描述為

式中：εu為混凝土強(qiáng)度極限應(yīng)變；fcu為混凝土極限強(qiáng)度，MPa；ε0＝0.002K，K＝1＋ρsvfyv／fc，ρsv為從箍筋外邊緣計(jì)算核心混凝土的體積配箍率；fyv為箍筋屈服強(qiáng)度.其中：Z為混凝土應(yīng)變軟化段斜率；b′為從箍筋外邊緣計(jì)算的核心混凝土寬度；sv為從箍筋中心算起的箍筋間距.

3.2 鋼筋本構(gòu)關(guān)系

本文采用的鋼筋本構(gòu)模型最初是Menegotto等［13］提出的，后經(jīng)Filippou等［12］修正已考慮等向強(qiáng)化影響的本構(gòu)模型（圖3）.Minegotto等所建議的模型如下，

圖2 約束和素混凝土本構(gòu)關(guān)系Fig.2 Constitutive relationship of confined concrete and unconfined

式中，R0，a1，a2均為參數(shù).

Filippou等［12］建議將線性的屈服漸近線進(jìn)行應(yīng)力平移，平移大小取決塑性應(yīng)變的最大值，表達(dá)式如下：

式中：σst為鋼筋拉應(yīng)力；fy為鋼筋屈服應(yīng)力；εmax為鋼筋最大應(yīng)變；εmin為鋼筋最小應(yīng)變；εy為鋼筋屈服應(yīng)變.

圖3 鋼筋本構(gòu)關(guān)系Fig.3 Constitutive relationship of steel

4 大轉(zhuǎn)動(dòng)計(jì)算方法

本文采用Yang［14］更新單元幾何坐標(biāo)的大轉(zhuǎn)動(dòng)的研究成果，將每步增量位移分解為剛體位移ur和自然變形un，

t時(shí)刻構(gòu)形的軸tx，ty，tz將轉(zhuǎn)化為t＋Δt時(shí)刻的構(gòu)形t＋Δtx，t＋Δty，t＋Δtz，如圖4所示，自然變形un參考t＋Δt時(shí)刻的構(gòu)形的表達(dá)式可以寫成如下，

式中，Δu可以表示為

式中：t＋ΔtL表示t＋Δt時(shí)刻構(gòu)形的單元長度；tL表示在t時(shí)刻構(gòu)形的單元長度.單元節(jié)點(diǎn)i轉(zhuǎn)動(dòng)位移表示為

式中：φ和ni表示如下：

其中，αi、βi、γi表示單元i端截面形心坐標(biāo)系（αi，βi，γi）向量.

圖4 單元?jiǎng)傮w位移和自然變形Fig.4 Rigid body displacements and natural deformations

5 算例分析與討論

5.1 鋼筋混凝土橋墩柱push－over分析

鋼筋混凝土橋墩柱，墩高2 265mm，矩形截面，沿橫橋向?qū)?00mm，順橋向?qū)?00mm，將橋墩支承的上部結(jié)構(gòu)的重量等效為墩頂?shù)募辛，大小為393.2kN.本文采用一個(gè)纖維梁單元模擬此橋墩柱，每個(gè)單元截面由32塊核心混凝土纖維、24塊保護(hù)層混凝土纖維和18根鋼筋纖維組成，如圖5 所示.保護(hù)層和核心混凝土纖維的本構(gòu)模型采用本文的3.1節(jié)的模型，材料特性參數(shù)［15］包括混凝土峰值壓應(yīng)力fc、峰值壓應(yīng)變?chǔ)?、極限壓應(yīng)力fcu、極限壓應(yīng)變?chǔ)與u，其值見表1.鋼筋纖維采用3.2節(jié)的本構(gòu)模型，初始彈性模量Es0＝200GPa，屈服強(qiáng)度fy＝357 MPa，屈服后剛度硬化系數(shù)b＝0.01，式（16）—（17）參數(shù)取值為，R0＝20，a1＝18.5，a3＝0.01，a4＝7.

表1 計(jì)算結(jié)果比較混凝土材料特性Tab.1 Material properties of concrete

對圖5的鋼筋混凝土橋墩結(jié)構(gòu)，采用本文建議的三維纖維梁單元并對截面的核心混凝土分成32、40、和64 塊三種情況進(jìn)行位移控制的增量pushover分析，得到圖6所示的該橋墩的墩底反力—位移關(guān)系曲線.

由圖6所示，當(dāng)位移增量增加到0.438m 時(shí)，由于單元出現(xiàn)大量的復(fù)特征值而首先出現(xiàn)不收斂狀況，程序自動(dòng)停止求解；核心混凝土分別劃分32、40和64塊計(jì)算結(jié)果近似，通過試算知一般鋼筋混凝土梁柱核心混凝土截面劃分到30～40塊基本可以滿足工程需求；當(dāng)頂點(diǎn)水平位移到0.012m，保護(hù)層混凝土已經(jīng)退出工作狀態(tài)，引起橋墩構(gòu)件的剛度部分退化，但考慮箍筋對核心混凝土約束作用而提高構(gòu)件的承載力，故在0.012～0.438m 范圍內(nèi)構(gòu)件內(nèi)混凝土、鋼筋出現(xiàn)應(yīng)力強(qiáng)化.鋼筋混凝土橋墩柱在push－over分析中，墩頂此時(shí)發(fā)生較大的位移和轉(zhuǎn)角，當(dāng)達(dá)到0.438m 時(shí)結(jié)構(gòu)承載力達(dá)到極限而發(fā)生推覆破壞.特別是在橋墩位移在0.012～0.438m 區(qū)間，鋼筋混凝土橋墩出現(xiàn)強(qiáng)非線性特性，需要考慮大位移、大轉(zhuǎn)動(dòng)幾何非線性.

5.2 William 肘式框架分析

如圖7所示肘式框架（P 為集中力），兩端嵌固，桿件為截面為矩形，截面的寬為19.13mm，高6.17 mm，彈性模量為71 018.5 MPa.William［16］對該結(jié)構(gòu)進(jìn)行試驗(yàn)研究和非線性分析，其分析結(jié)果與試驗(yàn)數(shù)據(jù)相當(dāng)吻合.William 在進(jìn)行分析時(shí)，考慮了構(gòu)件的大位移、軸力對彎曲剛度影響及彎曲引起構(gòu)件縮短等非線性耦合關(guān)系.本文在分析該結(jié)構(gòu)時(shí)，對每個(gè)構(gòu)件采用6個(gè)纖維單元，每個(gè)纖維單元截面沿高度方向劃分3塊纖維，沿寬度方向劃分2塊纖維進(jìn)行模擬分析，并與William 對該結(jié)構(gòu)的實(shí)驗(yàn)數(shù)據(jù)進(jìn)行對比分析，如圖8所示，表明本文計(jì)算結(jié)果在前、后屈曲階段與William 的實(shí)驗(yàn)結(jié)果基本一致.

5.3 Remseth空間框架穹頂分析

如圖9—10所示Remseth空間框架穹頂結(jié)構(gòu)，支座嵌固，桿件的截面為矩形，其截面尺寸為0.76 m×1.22 m.截面材料特性：彈性模量為20 690 MPa，剪切模量為8 621 MPa.本文計(jì)算分析結(jié)果與Shi等［17］對該結(jié)構(gòu)進(jìn)行彈塑性大變形、大轉(zhuǎn)動(dòng)分析結(jié)果和Park等［18］對該結(jié)構(gòu)大位移分析結(jié)果進(jìn)行了對比，如圖11所示.本文對該結(jié)構(gòu)每個(gè)構(gòu)件用一個(gè)纖維梁單元模擬，單元截面沿寬度劃分8塊纖維和沿高度劃分14塊纖維；Shi對該結(jié)構(gòu)分析時(shí)：采用能量法和單元任意分布的塑性鉸模型，并考慮彎曲—拉伸高度耦合非線性、大變形、大轉(zhuǎn)動(dòng)作用，用一個(gè)單元來模擬一個(gè)構(gòu)件；Park在分析該結(jié)構(gòu)時(shí)：構(gòu)件所有截面沿寬度采用8積分點(diǎn)和沿高度采用14積分點(diǎn)，每個(gè)構(gòu)件采用16個(gè)單元進(jìn)行該結(jié)構(gòu)的大位移、大轉(zhuǎn)動(dòng)前屈曲和空間框架大變形倒塌分析.圖11顯示，結(jié)果基本一致，表明本文建議的三維纖維梁單元可以精確分析空間框架考慮大位移、大轉(zhuǎn)動(dòng)的屈曲分析和倒塌分析.圖11中橫坐標(biāo)值是根據(jù)本文計(jì)算結(jié)果每3 000個(gè)子增量步計(jì)算結(jié)果來描述的，是根據(jù)區(qū)域來劃分?jǐn)?shù)值而非常規(guī)等間距.

6 結(jié)語

本文以三維連續(xù)體介質(zhì)力學(xué)虛功原理和虛位移原理為基礎(chǔ)，推導(dǎo)了大位移增量UL 列式，此增量列式中極大程度上保留大位移增量剛度矩陣kU.然后基于UL增量列式推導(dǎo)出小應(yīng)變、大位移、大轉(zhuǎn)動(dòng)三維纖維梁大位移增量非線性有限元UL 列式.該梁單元的切線剛度矩陣考慮了復(fù)合材料的非線性、大位移大轉(zhuǎn)動(dòng)高度幾何非線性.通過對幾個(gè)算例分析，得到非常有用的以下兩個(gè)結(jié)論：

（1）算例表明，本文建議的三維纖維梁單元是一種比較精確的單元，單元截面可以模擬真實(shí)的非線性復(fù)合材料本構(gòu)關(guān)系；能通過大荷載步增量、相對大的位移增量、每步增量步較少的子迭代步精確求解大位移、大轉(zhuǎn)動(dòng)問題.

（2）此外本文建議的三維纖維梁可以用于pushover分析；可精確分析結(jié)構(gòu)的前屈曲和后屈服狀態(tài)荷載與變形高度非線性特性；同時(shí)該單元也可用于結(jié)構(gòu)倒塌分析.

［1］ Zeris C A，Mahin S A.Analysis of reinforced concrete beamcolumns under uniaxial excitation［J］.Journal of Structural Engineering，ASCE，1988，114（4）：804.

［2］ Spacone E，F(xiàn)ilippou F C，Taucer F F.Fiber beam－column model for nonlinear analysis of R／C frames： part I.Formulation［J］.Earthquake Engineering and Structural Dynamics，1996，25：711.

［3］ 陳滔，黃宗明.基于有限單元柔度法的材料與幾何雙重非線性空間梁柱單元［J］.計(jì)算力學(xué)學(xué)報(bào)，2006，23（5）：524.CHEN Tao，HUANG Zongming.Material and geometrically nonlinear spatial beam－column element based on the finite element flexibility method ［J］. Chinese Journal of Computational Mechanics，2006，23（5）：524.

［4］ 聶建國，陶慕軒.采用纖維梁單元分析鋼－混凝土組合結(jié)構(gòu)地震反應(yīng)的原理［J］.建筑結(jié)構(gòu)學(xué)報(bào)，2011，32（10）：1.NIE Jianguo，TAO Muxuan.Theory of seismic response analysis of steel－concrete composite structures using fiber beam elements［J］.Journal of Building Structures，2011，32（10）：1.

［5］ 李易，陸新征，葉列平，等.混凝土框架結(jié)構(gòu)火災(zāi)連續(xù)倒塌數(shù)值分析模型［J］.工程力學(xué)，2012，29（4）：96.LI Yi，LU Xinzheng，YE Lieping，et al.Numerical models of fire induced progressive collapse analysis for reinforced concrete frame structures［J］.Engineering Mechanics，2012，29（4）：96.

［6］ Zubydan A H，ElSabbagh A I.Monotonic and cyclic behavior of concrete－filled steel－tube beam－columns considering local buckling effect［J］.Thin－Walled Structures，2011，49（4）：465.

［7］ Thai H T，Kim S E.Second－order inelastic analysis of cablestayed bridges［J］.Finite Elements in Analysis and Design，2012，53：48.

［8］ Zupan D，Saje M.The linearized three－dimensional beam theory of naturally curved and twisted beams：the strain vectors formulation ［J］. Computer Methods in Applied Mechanics and Engineering，2006，195（33／36）：4557.

［9］ Bathe K J，Wilson E L.Nonsap—a nonlinear structural analysis program［J］.Nuclear Engineering and Design，1974，29：266.

［10］ Riks E.An incremental approach to the solution of snapping and buckling problems［J］.International Journal of Solids and Structures，1979，15（7）：529.

［11］ Crisfield M A.An arc－length method including line searches and accelerations［J］.International Journal for Numerical Methods in Engineering，1983，19（9）：1269.

［12］ Taucer F F，Spacone E，F(xiàn)ilippou F C.A fiber beam－column element for seismic response analysis of reinforced concrete structures［R］.Berkeley：Earthquake Engineering Research Center of University of California，Berkeley，1991.

［13］ Menegotto M，Pinto P E，Method of analysis for cyclically loaded reinforced concrete Plane frames including changes in geometry and non－elastic behavior of elements under combined normal force and bending ［C］／／ Proceedings，IABSE Symposium on Resistance and Ultimate Deformability of Structures Acted on by Well－Defined Repeated Loads.Lisbon：IABSE，1973：15－20.

［14］ YANG Yeongbin.Theory and analysis of nonlinear framed structures［M］.［S.l.］：Prentice－Hall，Inc，1994.

［15］ 禚一，李忠獻(xiàn).基于顯式算法的纖維梁柱單元模型［J］.工程力學(xué)，2011，28（12）：39.ZHUO Yi，LI Zhongxian.Explicit algorithm－based fiber beam column element model［J］.Engineering Mechanics，2011，28（12）：39.

［16］ William F W.An approach to the non－linear behaviour of the members of a rigid jointed plane framework with finite deflections［J］.The Quarterly Journal of Mechanics ＆Applied Mathematics，1964，17（4）：451.

［17］ Shi G，Atluri S N.Elasto－plastic large deformation analysis of space－frames： a plastic－h(huán)inge and stress－based explicit derivation of tangent stiffnesses［J］.International Journal for Numerical Methods in Engineering，1988，26（3）：589.

［18］ Park M S，Lee B C.Geometrically non－linear and elastoplastic three dimensional shear flexible beam element of von－misestype hardending material［J］.International Journal for Numerical Methods in Engineering，1996，99（3）：383.